Commit [d13dba] default Maximize Restore History

New version of bwboundaries (from A Kelly)

hauberg hauberg 2010-12-22

added inst/fchcode.m
changed INDEX
changed inst/bwboundaries.m
changed src/Makefile
copied src/__imboundary__.cc -> src/__boundary__.cc
inst/fchcode.m Diff Switch to side-by-side view
Loading...
INDEX Diff Switch to side-by-side view
Loading...
inst/bwboundaries.m Diff Switch to side-by-side view
Loading...
src/Makefile Diff Switch to side-by-side view
Loading...
src/__imboundary__.cc to src/__boundary__.cc
--- a/src/__imboundary__.cc
+++ b/src/__boundary__.cc
@@ -1,162 +1,181 @@
-// Copyright (C) 2010 Soren Hauberg
-//
-// This program is free software; you can redistribute it and/or
-// modify it under the terms of the GNU General Public License
-// as published by the Free Software Foundation; either version 3
-// of the License, or (at your option) any later version.
-//
-// This program is distributed in the hope that it will be useful, but
-// WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-// General Public License for more details.
-//
-// You should have received a copy of the GNU General Public License
-// along with this program; if not, see <http://www.gnu.org/licenses/>.
+/*
+Copyright (C) 2010 Andrew Kelly, IPS Radio & Space Services
 
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 3 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program; If not, see <http://www.gnu.org/licenses/>.
+*/
+
+/**
+ *  Oct-file to trace the boundary of an object in a binary image.
+ *
+ *      b = boundary(region, conn=8)
+ */
 #include <octave/oct.h>
 
-inline int iseven (const int a)
+using namespace std;
+
+DEFUN_DLD(__boundary__, args, nargout,
+"-*- texinfo -*-\n\
+@deftypefn {Function File} {@var{b} = } boundary(@var{region})\n\
+@deftypefnx {Function File} {@var{b} = } boundary(@var{region}, @var{conn})\n\
+Trace the boundary of an object in a binary image.\n\
+\n\
+@code{boundary} computes the exterior clockwise boundary of the single \
+@var{conn}-connected object represented by the non-zero pixels \
+of @var{region}. It uses an algorithm based on Moore-neighbour tracing.\n\
+\n\
+@var{conn} can be either 8 (the default) or 4.\n\
+\n\
+@var{b} is an N-by-2 matrix containing the row/column coordinates of points \
+on the boundary. The first boundary point is the first non-zero \
+pixel of @var{region}, as determined by @code{find}. The last boundary \
+point is the same as the first.\n\
+@seealso{boundaries, bwlabel, find}\n\
+@end deftypefn")
 {
-  if ((a % 2) == 0)
-    return 1;
-  else
-    return 0;
-}
+  octave_value_list retval;
+  
+  enum { ROW, COL };
 
-Matrix trace_boundary (const boolMatrix &im, const int N, const octave_idx_type r,
-                       const octave_idx_type c)
-{
-  // Get size information
-  const octave_idx_type rows = im.rows ();
-  const octave_idx_type cols = im.columns ();
+  // check number of arguments
+  const int nargin = args.length ();
+  if (nargin > 2 || nargout != 1)
+    {
+      error ("__boundary__: wrong number of input arguments");  
+      return retval;
+    }
 
-  // Create list of points
-  typedef std::pair<int, int> point;
-  std::list <point> P;
-  
-  // Add first point
-  const point P0 (r, c); // first list element
-  point P1 (-1, -1); // second list element
-  point Pn1 (-1, -1); // second last list element
-  P.push_back (P0);
-  int len = 1;
-  
-  // Create a simple lookup table that translates 'dir' to a row and column offset
-  const int dir2row4 [] = { 0, -1,  0, +1};
-  const int dir2col4 [] = {+1,  0, -1,  0};
-  const int dir2row8 [] = { 0, -1, -1, -1,  0, +1, +1, +1};
-  const int dir2col8 [] = {+1, +1,  0, -1, -1, -1,  0, +1};
-  
-  // Start searching from there...
-  int dir = N - 1;
-  int curr_dir = dir; //(N == 4) ? (dir + 3) % N : (dir + 6 + iseven (dir)) % N;
-  int delta_r, delta_c;
-  octave_idx_type row = r, col = c;
-  while (true)
-//  for (int z = 0; z < 1000; z++)
+  // extract arguments
+  const boolMatrix unpadded = args (0).bool_matrix_value ();
+  const int conn = (nargin > 1)
+      ? (int) args (1).scalar_value ()
+      : 8;
+    
+  if (error_state)
     {
-      OCTAVE_QUIT;
-      
-      // Get next search direction
-      if (N == 4)
-        {
-          //curr_dir = (dir + 3) % N;
-          delta_r = dir2row4 [curr_dir];
-          delta_c = dir2col4 [curr_dir];
-        }
-      else
-        {
-          //curr_dir = (dir + 6 + iseven (dir)) % N;
-          delta_r = dir2row8 [curr_dir];
-          delta_c = dir2col8 [curr_dir];
-        }
-        
-      // Is a pixel available at the search direction
-      const octave_idx_type curr_r = row + delta_r;
-      const octave_idx_type curr_c = col + delta_c;
-/*      std::cerr << " curr_r = " << curr_r
-                << " curr_c = " << curr_c
-                << " curr_dir = " << curr_dir
-                << " row = " << row
-                << " colr = " << col
-                << std::endl;
-*/      
-      if (curr_r >= 0 && curr_r < rows && curr_c >= 0 && curr_c < cols && im (curr_r, curr_c))
-        {
-          // Update 'dir'
-          dir = curr_dir;
-          curr_dir = (N == 4) ? (dir + 3) % N : (dir + 6 + iseven (dir)) % N;
-          
-          // Add point to list
-          const point Pn (curr_r, curr_c);
-          P.push_back (Pn);
-          len ++;
-          
-          // Update 'row' and 'col'
-          row = curr_r;
-          col = curr_c;
-          
-          // Save the second element of P for the stop criteria
-          if (len == 2)
-            {
-              P1.first = curr_r;
-              P1.second = curr_c;
-            }
-            
-          // Should we stop?
-          if (Pn.first ==  P1.first && Pn.second == P1.second &&
-              Pn1.first == P0.first && Pn1.second == P0.second)
-            break;
-            
-          // Save current point for next time
-          Pn1 = Pn;
-        }
-      else
-        {
-          // Update search direction
-          curr_dir = (curr_dir+1) % N;
-        }
-    } // end while
+      error ("__boundary__: internal error");
+      return retval;
+    }
 
-  // Copy data to output matrix
-  Matrix out (len-1, 2);
-  std::list<point>::const_iterator iter = P.begin ();
-  for (int idx = 0; idx < len-2; iter++, idx++)
-    {
-      out (idx, 0) = iter->second + 1;
-      out (idx, 1) = iter->first + 1;
-    }
-  out (len-2, 0) = P0.second + 1;
-  out (len-2, 1) = P0.first + 1;
-  
-  return out;
-}
+  // pad to avoid boundary issues
+  int rows = unpadded.rows ();
+  int cols = unpadded.columns ();
+  boolMatrix region (rows + 2, cols + 2, false);
+  for (int r = 0; r < rows; r++)
+    for (int c = 0; c < cols; c++)
+      region.elem (r+1, c+1) = unpadded (r, c);
 
-DEFUN_DLD(__imboundary__, args, , "\
--*- texinfo -*-\n\
-@deftypefn {Function File} __imboundary__ (@var{bw}, @var{N}, @var{r}, @var{c})\n\
-Undocumented internal function.\n\
-User interface is available in @code{bwboundaries}.\n\
-@end deftypefn\n\
-")
-{
-  // Handle input
-  octave_value_list retval;
-  if (args.length () != 4) {
-    error ("__imboundary__: not enough input arguments");
+  // the padded size
+  rows += 2;
+  cols += 2;
+
+  // find the (first two) true pixels, if any
+  std::vector <int> pixels;
+  for (int i = 0; pixels.size () < 2 && i < region.numel (); ++i)
+    if (region.elem (i))
+      pixels.push_back (i);
+
+  if (pixels.empty ())
     return retval;
-  }
-    
-  const boolMatrix im = args (0).bool_matrix_value ();
-  const int N = (int) args (1).scalar_value ();
-  const octave_idx_type r = (octave_idx_type) args (2).scalar_value () - 1;
-  const octave_idx_type c = (octave_idx_type) args (3).scalar_value () - 1;
-  if (error_state || (N != 4)) // && N != 8))
-    error ("__imboundary__: invalid input arguments");
+
+  // the starting boundary point
+  const int start = pixels [0];
+  std::vector <int> bound;
+  bound.push_back (start);
+
+  // is this the only point?
+  if (pixels.size () == 1)
+    bound.push_back (start);
+
+  // otherwise, find the boundary by tracing the Moore neighbourhood of its pixels
+  //
+  //      8-connected:    7 0 1      4-connected:      0
+  //                      6 . 2                      3 . 1
+  //                      5 4 3                        2
   else
     {
-      Matrix out = trace_boundary (im, N, r, c);
-      retval.append (out);
+      // relative row/column positions
+      static const int row8 [] = {-1, -1,  0,  1,  1,  1,  0, -1};
+      static const int col8 [] = { 0,  1,  1,  1,  0, -1, -1, -1};
+      static const int row4 [] = {-1,      0,      1,      0    };
+      static const int col4 [] = { 0,      1,      0,     -1    };
+      const int* mr = (conn == 4) ? row4 : row8;
+      const int* mc = (conn == 4) ? col4 : col8;
+
+      // next after backing-up
+      static const int back8 [] = {7, 7, 1, 1, 3, 3, 5, 5};
+      static const int back4 [] = {3, 0, 1, 2};
+      const int* mBack = (conn == 4) ? back4 : back8;
+
+      // relative indexes into the region for the Moore neighbourhood pixels
+      int mi [conn];
+      for (int i = 0; i < conn; ++i)
+        mi[i] = mr[i] + (rows * mc [i]);
+
+      // next neighbourhood pixel
+      static const int next8 [] = {1, 2, 3, 4, 5, 6, 7, 0};
+      static const int next4 [] = {1, 2, 3, 0};
+      const int* mNext = (conn == 4) ? next4 : next8;
+
+      // the final boundary point to be visited
+      int finish = 0;
+      for (int i = 0; i < conn; ++i)
+        if (region.elem(start + mi [i]))
+          finish = start + mi [i];
+
+      // look for the next boundary point, starting at the next neighbour
+      int bp = start;
+      int mCurrent = mNext [0];
+      bool done = false;
+      while (!done)
+        {
+          // next neighbour
+          int cp = bp + mi [mCurrent];
+
+          // if this pixel is false, try the next one
+          if (!region.elem (cp))
+            {
+              mCurrent = mNext [mCurrent];
+            }
+          // otherwise, we have another boundary point
+          else
+            {
+              bound.push_back (cp);
+
+              // either we're back at the start for the last time
+              if (bp == finish && cp == start)
+                {
+                  done = true;
+                }
+              // or we step back to where we came in from, and continue
+              else
+                {
+                  bp = cp;
+                  mCurrent = mBack [mCurrent];
+                }
+            }
+        }
     }
-  return retval;    
+
+  // convert boundary points to row/column coordinates
+  Matrix b (bound.size (), 2);
+  for (unsigned int i = 0; i < bound.size (); i++)
+    {
+      const int point = bound [i];
+      b (i, ROW) = point % rows;
+      b (i, COL) = point / rows;
+    }
+  
+  retval.append (b);
+  return retval;
 }