# HG changeset patch
# User carandraug
# Date 1366168147 0
# Node ID 5d87258abc1a6c30a6b388121c7d0b763560a421
# Parent  2d1db11b14705567a686135db3aa3f0597fbf9b6
bwdist: matlab compatibility fixes

* distance matrix must be class single
* index class is uint32 or uint64 dependent on matrix size
* object pixels are all non-zero (not larger than zero)
* added tests
* OCTAVE_QUIT on the main loops so we can abort
* sqrt when using euclidean distance is done after getting the distances
* dropped some templates from last time (not really a performance increase)

diff -r 2d1db11b1470 -r 5d87258abc1a main/image/NEWS
--- a/main/image/NEWS	Tue Apr 16 15:32:11 2013 +0000
+++ b/main/image/NEWS	Wed Apr 17 03:09:07 2013 +0000
@@ -71,6 +71,10 @@
  ** Bug fixes on the concavity, intermodes, maxlikelihood, and minimum methods
     of graythresh.
 
+ ** The `bwdist' function will now consider any non zero value as object pixels,
+    the class of the distance matrix has changed to single, and indexes an
+    uint dependent on the matrix size.
+
 Summary of important user-visible changes for image 2.0.0:
 -------------------------------------------------------------------
 
diff -r 2d1db11b1470 -r 5d87258abc1a main/image/src/bwdist.cc
--- a/main/image/src/bwdist.cc	Tue Apr 16 15:32:11 2013 +0000
+++ b/main/image/src/bwdist.cc	Wed Apr 17 03:09:07 2013 +0000
@@ -14,11 +14,6 @@
 // You should have received a copy of the GNU General Public License along with
 // this program; if not, see <http://www.gnu.org/licenses/>.
 
-// TODO output Matrix must be single class
-//      class of IDX is uint32 for numel () <= 2^32 - 1
-//                      uint64 for numel () >= 2^32
-//      check for any improvement in using boolean matrix
-
 #include <octave/oct.h>
 
 /*
@@ -41,17 +36,12 @@
  Edited in 2009 to form the foundation for Octave BWDIST:
  added #define-configurable distance measure and function name
  
- Edited in 2013 for C++ and removed the #define stuff
+ Edited in 2013 for C++, removed the #define stuff, and other
+ fixes for matlab compatibility.
  */
 
-/*
-  Using template for optimization.
-  Because func is called a lot of times in edtfunc, we instantiate edtfunc
-  multiple times for each of the different distance measuring functions.
-*/
-
-template<typename func>
-void edtfunc (const Matrix &img,
+void edtfunc (float (*func)(short int, short int),
+              const Matrix &img,
               short *distx,
               short *disty)
 {
@@ -60,7 +50,7 @@
   const int numel = img.numel ();
 
   int x, y, i;
-  double olddist2, newdist2, newdistx, newdisty;
+  float olddist2, newdist2, newdistx, newdisty;
   bool changed;
 
   // Initialize index offsets for the current image width
@@ -76,7 +66,7 @@
   // Initialize the distance images to be all large values
   for (i = 0; i < numel; i++)
     {
-      if(img(i) <= 0.0)
+      if(img(i) == 0.0)
         {
           distx[i] = 32000; // Large but still representable in a short, and
           disty[i] = 32000; // 32000^2 + 32000^2 does not overflow an int
@@ -96,19 +86,18 @@
       // Scan rows, except first row
       for (y = 1; y < h; y++)
         {
-
           // move index to leftmost pixel of current row
           i = y*w;
 
           /* scan right, propagate distances from above & left */
 
           /* Leftmost pixel is special, has no left neighbors */
-          olddist2 = func::apply(distx[i], disty[i]);
+          olddist2 = (*func)(distx[i], disty[i]);
           if(olddist2 > 0) // If not already zero distance
             {
               newdistx = distx[i+offset_u];
               newdisty = disty[i+offset_u]+1;
-              newdist2 = func::apply(newdistx, newdisty);
+              newdist2 = (*func)(newdistx, newdisty);
               if(newdist2 < olddist2)
                 {
                   distx[i]=newdistx;
@@ -119,7 +108,7 @@
 
               newdistx = distx[i+offset_ur]-1;
               newdisty = disty[i+offset_ur]+1;
-              newdist2 = func::apply(newdistx, newdisty);
+              newdist2 = (*func)(newdistx, newdisty);
               if(newdist2 < olddist2)
                 {
                   distx[i]=newdistx;
@@ -132,12 +121,13 @@
           /* Middle pixels have all neighbors */
           for(x=1; x<w-1; x++, i++)
             {
-              olddist2 = func::apply(distx[i], disty[i]);
+              OCTAVE_QUIT;
+              olddist2 = (*func)(distx[i], disty[i]);
               if(olddist2 == 0) continue; // Already zero distance
 
               newdistx = distx[i+offset_l]+1;
               newdisty = disty[i+offset_l];
-              newdist2 = func::apply(newdistx, newdisty);
+              newdist2 = (*func)(newdistx, newdisty);
               if(newdist2 < olddist2)
                 {
                   distx[i]=newdistx;
@@ -148,7 +138,7 @@
 
               newdistx = distx[i+offset_lu]+1;
               newdisty = disty[i+offset_lu]+1;
-              newdist2 = func::apply(newdistx, newdisty);
+              newdist2 = (*func)(newdistx, newdisty);
               if(newdist2 < olddist2)
                 {
                   distx[i]=newdistx;
@@ -159,7 +149,7 @@
 
               newdistx = distx[i+offset_u];
               newdisty = disty[i+offset_u]+1;
-              newdist2 = func::apply(newdistx, newdisty);
+              newdist2 = (*func)(newdistx, newdisty);
               if(newdist2 < olddist2)
                 {
                   distx[i]=newdistx;
@@ -170,7 +160,7 @@
 
               newdistx = distx[i+offset_ur]-1;
               newdisty = disty[i+offset_ur]+1;
-              newdist2 = func::apply(newdistx, newdisty);
+              newdist2 = (*func)(newdistx, newdisty);
               if(newdist2 < olddist2)
                 {
                   distx[i]=newdistx;
@@ -180,12 +170,12 @@
             }
 
           /* Rightmost pixel of row is special, has no right neighbors */
-          olddist2 = func::apply(distx[i], disty[i]);
+          olddist2 = (*func)(distx[i], disty[i]);
           if(olddist2 > 0) // If not already zero distance
             {
               newdistx = distx[i+offset_l]+1;
               newdisty = disty[i+offset_l];
-              newdist2 = func::apply(newdistx, newdisty);
+              newdist2 = (*func)(newdistx, newdisty);
               if(newdist2 < olddist2)
                 {
                   distx[i]=newdistx;
@@ -196,7 +186,7 @@
 
               newdistx = distx[i+offset_lu]+1;
               newdisty = disty[i+offset_lu]+1;
-              newdist2 = func::apply(newdistx, newdisty);
+              newdist2 = (*func)(newdistx, newdisty);
               if(newdist2 < olddist2)
                 {
                   distx[i]=newdistx;
@@ -207,7 +197,7 @@
 
               newdistx = distx[i+offset_u];
               newdisty = disty[i+offset_u]+1;
-              newdist2 = func::apply(newdistx, newdisty);
+              newdist2 = (*func)(newdistx, newdisty);
               if(newdist2 < olddist2)
                 {
                   distx[i]=newdistx;
@@ -224,12 +214,13 @@
           /* scan left, propagate distance from right */
           for(x=w-2; x>=0; x--, i--)
             {
-              olddist2 = func::apply(distx[i], disty[i]);
+              OCTAVE_QUIT;
+              olddist2 = (*func)(distx[i], disty[i]);
               if(olddist2 == 0) continue; // Already zero distance
               
               newdistx = distx[i+offset_r]-1;
               newdisty = disty[i+offset_r];
-              newdist2 = func::apply(newdistx, newdisty);
+              newdist2 = (*func)(newdistx, newdisty);
               if(newdist2 < olddist2)
                 {
                   distx[i]=newdistx;
@@ -242,18 +233,19 @@
       /* Scan rows in reverse order, except last row */
       for(y=h-2; y>=0; y--)
         {
+          OCTAVE_QUIT;
           /* move index to rightmost pixel of current row */
           i = y*w + w-1;
 
           /* Scan left, propagate distances from below & right */
 
           /* Rightmost pixel is special, has no right neighbors */
-          olddist2 = func::apply(distx[i], disty[i]);
+          olddist2 = (*func)(distx[i], disty[i]);
           if(olddist2 > 0) // If not already zero distance
             {
               newdistx = distx[i+offset_d];
               newdisty = disty[i+offset_d]-1;
-              newdist2 = func::apply(newdistx, newdisty);
+              newdist2 = (*func)(newdistx, newdisty);
               if(newdist2 < olddist2)
                 {
                   distx[i]=newdistx;
@@ -264,7 +256,7 @@
 
               newdistx = distx[i+offset_dl]+1;
               newdisty = disty[i+offset_dl]-1;
-              newdist2 = func::apply(newdistx, newdisty);
+              newdist2 = (*func)(newdistx, newdisty);
               if(newdist2 < olddist2)
                 {
                   distx[i]=newdistx;
@@ -277,12 +269,13 @@
           /* Middle pixels have all neighbors */
           for(x=w-2; x>0; x--, i--)
             {
-              olddist2 = func::apply(distx[i], disty[i]);
+              OCTAVE_QUIT;
+              olddist2 = (*func)(distx[i], disty[i]);
               if(olddist2 == 0) continue; // Already zero distance
 
               newdistx = distx[i+offset_r]-1;
               newdisty = disty[i+offset_r];
-              newdist2 = func::apply(newdistx, newdisty);
+              newdist2 = (*func)(newdistx, newdisty);
               if(newdist2 < olddist2)
                 {
                   distx[i]=newdistx;
@@ -293,7 +286,7 @@
 
               newdistx = distx[i+offset_rd]-1;
               newdisty = disty[i+offset_rd]-1;
-              newdist2 = func::apply(newdistx, newdisty);
+              newdist2 = (*func)(newdistx, newdisty);
               if(newdist2 < olddist2)
                 {
                   distx[i]=newdistx;
@@ -304,7 +297,7 @@
 
               newdistx = distx[i+offset_d];
               newdisty = disty[i+offset_d]-1;
-              newdist2 = func::apply(newdistx, newdisty);
+              newdist2 = (*func)(newdistx, newdisty);
               if(newdist2 < olddist2)
                 {
                   distx[i]=newdistx;
@@ -315,7 +308,7 @@
 
               newdistx = distx[i+offset_dl]+1;
               newdisty = disty[i+offset_dl]-1;
-              newdist2 = func::apply(newdistx, newdisty);
+              newdist2 = (*func)(newdistx, newdisty);
               if(newdist2 < olddist2)
                 {
                   distx[i]=newdistx;
@@ -324,12 +317,12 @@
                 }
             }
           /* Leftmost pixel is special, has no left neighbors */
-          olddist2 = func::apply(distx[i], disty[i]);
+          olddist2 = (*func)(distx[i], disty[i]);
           if(olddist2 > 0) // If not already zero distance
             {
               newdistx = distx[i+offset_r]-1;
               newdisty = disty[i+offset_r];
-              newdist2 = func::apply(newdistx, newdisty);
+              newdist2 = (*func)(newdistx, newdisty);
               if(newdist2 < olddist2)
                 {
                   distx[i]=newdistx;
@@ -340,7 +333,7 @@
 
               newdistx = distx[i+offset_rd]-1;
               newdisty = disty[i+offset_rd]-1;
-              newdist2 = func::apply(newdistx, newdisty);
+              newdist2 = (*func)(newdistx, newdisty);
               if(newdist2 < olddist2)
                 {
                   distx[i]=newdistx;
@@ -351,7 +344,7 @@
 
               newdistx = distx[i+offset_d];
               newdisty = disty[i+offset_d]-1;
-              newdist2 = func::apply(newdistx, newdisty);
+              newdist2 = (*func)(newdistx, newdisty);
               if(newdist2 < olddist2)
                 {
                   distx[i]=newdistx;
@@ -366,13 +359,14 @@
           i = y*w + 1;
           for(x=1; x<w; x++, i++)
             {
+              OCTAVE_QUIT;
               /* scan right, propagate distance from left */
-              olddist2 = func::apply(distx[i], disty[i]);
+              olddist2 = (*func)(distx[i], disty[i]);
               if(olddist2 == 0) continue; // Already zero distance
 
               newdistx = distx[i+offset_l]+1;
               newdisty = disty[i+offset_l];
-              newdist2 = func::apply(newdistx, newdisty);
+              newdist2 = (*func)(newdistx, newdisty);
               if(newdist2 < olddist2)
                 {
                   distx[i]=newdistx;
@@ -388,41 +382,49 @@
 
 // The different functions used to calculate distance, as a
 // class so its typename can be used for edtfunc template
-struct euclidean
-{
-  static double apply(short x, short y)
-  { return sqrt ((int)(x)*(x) + (y)*(y)); }
-};
-struct chessboard
+static float euclidean (short x, short y)
 {
-  static double apply(short x, short y)
-  { return std::max (abs (y), abs (x)); }
-};
-struct cityblock
-{
-  static double apply(short x, short y)
-  { return abs (x) + abs (y); }
-};
-struct quasi_euclidean
+  // the actual euclidean distance, is the square root of this. But
+  // squaring does not change the order of the distances, so we can
+  // do it in the end and only in the values that matter
+  return ((int)(x)*(x) + (y)*(y));
+}
+
+static float chessboard (short x, short y)
+{ return std::max (abs (y), abs (x)); }
+
+static float cityblock (short x, short y)
+{ return abs (x) + abs (y); }
+
+static float quasi_euclidean (short x, short y)
 {
-  static double apply(short x, short y)
-    {
-      static const double sqrt2_1 = sqrt (2) - 1;
-      return abs(x)>abs(y) ? (abs(x) + sqrt2_1 * abs(y)) :
-                             (sqrt2_1 * abs(x) + abs(y)) ;
-    }
-};
+  static const float sqrt2_1 = sqrt (2) - 1;
+  return abs(x)>abs(y) ? (abs(x) + sqrt2_1 * abs(y)) :
+                         (sqrt2_1 * abs(x) + abs(y)) ;
+}
 
-template<typename func>
-void loop_distance (const Matrix &bw,
-                    Matrix &dist,
-                    short *xdist,
-                    short *ydist)
+FloatMatrix calc_distances (float (*func)(short, short),
+                            Matrix bw,
+                            short *xdist,
+                            short *ydist)
 {
-  edtfunc<func> (bw, xdist, ydist);
+  FloatMatrix dist (bw.dims ());
+  edtfunc (func, bw, xdist, ydist);
+  const int numel = dist.numel ();
+  for (int i = 0; i < numel; i++)
+    dist(i) = (*func)(xdist[i], ydist[i]);
+  return dist;
+}
+
+template <class T>
+T calc_index (Matrix bw, short *xdist, short *ydist)
+{
+  T idx (bw.dims ());
   const int numel = bw.numel ();
-  for (int i = 0; i < numel; i++)
-    dist(i) = func::apply (xdist[i], ydist[i]);
+  const int rows  = bw.rows ();
+  for(int i = 0; i < numel; i++)
+    idx (i) = i+1 - xdist[i] - ydist[i]*rows;
+  return idx;
 }
 
 DEFUN_DLD (bwdist, args, nargout,
@@ -433,12 +435,14 @@
 Compute distance transform in binary image.\n\
 \n\
 The image @var{bw} must be a binary matrix  For @sc{matlab} compatibility, no\n\
-check is performed, all non-zero values are considered true, or object pixels.\n\
+check is performed, all non-zero values are considered object pixels.\n\
 The return value @var{dist}, is the distance of each background pixel to the\n\
-closest object pixel.\n\
+closest object pixel in a matrix of class @code{single}.\n\
 \n\
 @var{idx} is the linear index for the closest object, used to calculate the\n\
-distance for each of the pixels.\n\
+distance for each of the pixels.  Its class is dependent on the number of\n\
+elements in @var{bw}, @code{uint64} if less than 2^32 elements, @code{uint32}\n\
+otherwise.\n\
 \n\
 The distance can be measured through different @var{method}s:\n\
 \n\
@@ -468,9 +472,6 @@
 
   // for matlab compatibility, we do not actually check if the values are all
   // 0 and 1, any non-zero value is considered true
-
-  // FIXME const boolMatrix bw = args (0).bool_matrix_value();
-
   const Matrix bw = args (0).matrix_value ();
   if (error_state)
     {
@@ -505,15 +506,22 @@
   OCTAVE_LOCAL_BUFFER (short, xdist, numel);
   OCTAVE_LOCAL_BUFFER (short, ydist, numel);
 
-  Matrix dist (bw.dims ()); // the output distance matrix
+  FloatMatrix dist;
   if (method == "euclidean")
-    loop_distance<euclidean> (bw, dist, xdist, ydist);
+    {
+      dist = calc_distances (euclidean, bw, xdist, ydist);
+      const Array<octave_idx_type> positions = (!bw).find ();
+      const int zpos = positions.numel();
+      for (int i = 0; i < zpos; i++) {
+        dist (positions(i)) = sqrt(dist(positions(i)));
+      }
+    }
   else if (method == "chessboard")
-    loop_distance<chessboard> (bw, dist, xdist, ydist);
+    dist = calc_distances (chessboard,      bw, xdist, ydist);
   else if (method == "cityblock")
-    loop_distance<cityblock> (bw, dist, xdist, ydist);
+    dist = calc_distances (cityblock,       bw, xdist, ydist);
   else if (method == "quasi-euclidean")
-    loop_distance<quasi_euclidean> (bw, dist, xdist, ydist);
+    dist = calc_distances (quasi_euclidean, bw, xdist, ydist);
   else
     error ("bwdist: unknown METHOD '%s'", method.c_str ());
 
@@ -522,12 +530,10 @@
   // Compute optional 'index to closest object pixel', only if requested
   if (nargout > 1)
     {
-      Matrix idx (bw.dims ());
-      const int rows = bw.rows ();
-      for(int i = 0; i < numel; i++)
-        idx (i) = i+1 - xdist[i] - ydist[i]*rows;
-
-      retval(1) = idx;
+      if (numel >= pow (2, 32))
+        retval(1) = calc_index<uint64NDArray> (bw, xdist, ydist);
+      else
+        retval(1) = calc_index<uint32NDArray> (bw, xdist, ydist);
     }
 
   return retval;
@@ -553,6 +559,7 @@
 %!         0.00000   0.00000   0.00000   0.00000   1.00000   1.00000   1.00000   0.00000
 %!         0.00000   0.00000   0.00000   1.00000   1.41421   1.00000   0.00000   1.00000
 %!         1.00000   1.00000   0.00000   1.00000   2.00000   1.00000   0.00000   0.00000];
+%! out = single (out);
 %!
 %!assert (bwdist (bw), out, 0.0001);  # default is euclidean
 %!assert (bwdist (bw, "euclidean"), out, 0.0001);
@@ -566,6 +573,7 @@
 %!         0   0   0   0   1   1   1   0
 %!         0   0   0   1   1   1   0   1
 %!         1   1   0   1   2   1   0   0];
+%! out = single (out);
 %!
 %!assert (bwdist (bw, "chessboard"), out);
 %!
@@ -577,6 +585,7 @@
 %!         0   0   0   0   1   1   1   0
 %!         0   0   0   1   2   1   0   1
 %!         1   1   0   1   2   1   0   0];
+%! out = single (out);
 %!
 %!assert (bwdist (bw, "cityblock"), out);
 %!
@@ -588,11 +597,14 @@
 %!         0.00000   0.00000   0.00000   0.00000   1.00000   1.00000   1.00000   0.00000
 %!         0.00000   0.00000   0.00000   1.00000   1.41421   1.00000   0.00000   1.00000
 %!         1.00000   1.00000   0.00000   1.00000   2.00000   1.00000   0.00000   0.00000];
+%! out = single (out);
 %!
 %!assert (bwdist (bw, "quasi-euclidean"), out, 0.0001);
 %!
 %! bw(logical (bw)) = 3; # there is no actual check if matrix is binary or 0 and 1
 %!assert (bwdist (bw, "quasi-euclidean"), out, 0.0001);
+%! bw(logical (bw)) = -2; # anything non-zero is considered object
+%!assert (bwdist (bw, "quasi-euclidean"), out, 0.0001);
 %!
 %!error bwdist (bw, "not a valid method");
 */