changeset 320:c540e8281d5a

Don't crash with non-square matrices
author hauberg
date Tue, 08 Apr 2008 18:24:35 +0000
parents 7077805b4fef
children df97404d1fb4
files src/__bwdist.cc
diffstat 1 files changed, 82 insertions(+), 65 deletions(-) [+]
line wrap: on
line diff
--- a/src/__bwdist.cc
+++ b/src/__bwdist.cc
@@ -32,7 +32,8 @@
 inline T square(const T &x) { return x*x; };
 
 /* dt of 1d function using squared distance */
-static float *dt(float *f, int n) {
+static float* dt (float *f, int n)
+{
   float *d = new float[n];
   int *v = new int[n];
   float *z = new float[n+1];
@@ -40,24 +41,27 @@
   v[0] = 0;
   z[0] = -INF;
   z[1] = +INF;
-  for (int q = 1; q <= n-1; q++) {
-    float s  = ((f[q]+square(q))-(f[v[k]]+square(v[k])))/(2*q-2*v[k]);
-    while (s <= z[k]) {
-      k--;
-      s  = ((f[q]+square(q))-(f[v[k]]+square(v[k])))/(2*q-2*v[k]);
+  for (int q = 1; q <= n-1; q++)
+    {
+      float s  = ((f[q] + square (q)) - (f[v[k]] + square (v[k])))/(2*q-2*v[k]);
+      while (s <= z[k])
+        {
+          k--;
+          s  = ((f[q] + square (q)) - (f[v[k]] + square (v[k])))/(2*q-2*v[k]);
+        }
+      k++;
+      v[k] = q;
+      z[k] = s;
+      z[k+1] = +INF;
     }
-    k++;
-    v[k] = q;
-    z[k] = s;
-    z[k+1] = +INF;
-  }
 
   k = 0;
-  for (int q = 0; q <= n-1; q++) {
-    while (z[k+1] < q)
-      k++;
-    d[q] = square(q-v[k]) + f[v[k]];
-  }
+  for (int q = 0; q <= n-1; q++)
+    {
+      while (z[k+1] < q)
+        k++;
+      d[q] = square (q-v[k]) + f[v[k]];
+    }
 
   delete [] v;
   delete [] z;
@@ -65,57 +69,67 @@
 }
 
 /* dt of 2d function using squared distance */
-static void dt(NDArray &im) {
-  const int width = im.dim2();
-  const int height = im.dim1();
+static void dt (NDArray &im)
+{
+  const int width = im.dim1 ();
+  const int height = im.dim2 ();
   float *f = new float[std::max(width,height)];
 
   // transform along columns
-  for (int x = 0; x < width; x++) {
-    for (int y = 0; y < height; y++) {
-      f[y] = im(x, y);
+  for (int x = 0; x < width; x++)
+    {
+      for (int y = 0; y < height; y++)
+        {
+          f[y] = im (x, y);
+        }
+      float *d = dt (f, height);
+      for (int y = 0; y < height; y++)
+        {
+          im (x, y) = d[y];
+        }
+      delete [] d;
     }
-    float *d = dt(f, height);
-    for (int y = 0; y < height; y++) {
-      im(x, y) = d[y];
-    }
-    delete [] d;
-  }
 
-  // transform along rows
-  for (int y = 0; y < height; y++) {
-    for (int x = 0; x < width; x++) {
-      f[x] = im(x, y);
-    }
-    float *d = dt(f, width);
-    for (int x = 0; x < width; x++) {
-      im(x, y) = d[x];
-    }
-    delete [] d;
-  }
+    // transform along rows
+    for (int y = 0; y < height; y++)
+      {
+        for (int x = 0; x < width; x++)
+          {
+            f[x] = im (x, y);
+          }
+        float *d = dt (f, width);
+        for (int x = 0; x < width; x++)
+          {
+            im (x, y) = d[x];
+          }
+        delete [] d;
+      }
 
-  delete f;
+    delete [] f;
 }
 
 
 /* dt of binary image using squared distance */
-void dt(boolNDArray &im, NDArray &out) {
-  const int width = im.dim2();
-  const int height = im.dim1();
+void dt (const boolNDArray &im, NDArray &out)
+{
+  const int width = im.dim1 ();
+  const int height = im.dim2 ();
 
-  for (int y = 0; y < height; y++) {
-    for (int x = 0; x < width; x++) {
-      if (im(x, y))
-        out(x, y) = 0;
-      else
-        out(x, y) = INF;
+  for (int y = 0; y < height; y++)
+    {
+      for (int x = 0; x < width; x++)
+        {
+          if (im (x, y))
+            out (x, y) = 0;
+          else
+            out (x, y) = INF;
+        }
     }
-  }
 
-  dt(out);
+  dt (out);
 }
 
-DEFUN_DLD(__bwdist, args, nargout, "\
+DEFUN_DLD (__bwdist, args, nargout, "\
 -*- texinfo -*-\n\
 @deftypefn {Function File} __bwdist (@var{bw})\n\
 Computes the Euclidian Distance Transform for a binary image @var{bw}.\n\
@@ -126,22 +140,25 @@
 {
   octave_value_list retval;
   
-  boolNDArray bw = args(0).bool_array_value();
-  if (error_state) {
-    error("__bwdist: input must be a boolean matrix");
-    return retval;
-  }
+  const boolNDArray bw = args(0).bool_array_value ();
+  if (error_state)
+    {
+      error ("__bwdist: input must be a boolean matrix");
+      return retval;
+    }
   
-  const dim_vector dims = bw.dims();
-  if (dims.length() != 2) {
-    error("__bwdist: currently only binary images are supported");
-    return retval;
-  }
+  const dim_vector dims = bw.dims ();
+  if (dims.length () != 2)
+    {
+      error ("__bwdist: currently only binary images are supported");
+      return retval;
+    }
   
-  NDArray *out = new NDArray(dims);
+  NDArray out (dims);
   
-  dt(bw, *out);
+  dt (bw, out);
   
-  retval.append(*out);
+  retval.append (out);
+printf("146\n"); fflush(stdout);
   return retval;
 }