new file mode 100644
--- /dev/null
+++ b/bwfill.cc
@@ -0,0 +1,215 @@
+/*
+ * BWFILL: fill a bw image starting at points
+ * imo= block(im, xregs, yregs);
+ *
+ * Copyright (C) 1999 Andy Adler
+ * This code has no warrany whatsoever.
+ * Do what you like with this code as long as you
+ *     leave this copyright in place.
+ *
+ * $Id$
+ */
+
+#include <octave/oct.h>
+
+#define   UP     (-1)
+#define   DN     (+1)
+#define   RT     (+ioM)
+#define   LF     (-ioM)
+
+/*
+ * check if the point needs to be filled, if so
+ * fill it and change the appropriate variables
+ */
+void checkpoint( int pt,
+      unsigned char* imo,
+               int * ptstack,
+               int * npoints ) {
+// printf("filling %d np=%d fill=%d\n",pt,*npoints, *(imo+pt)==0 );
+   if( *(imo+pt) != 0 ) return;
+
+   *(imo+pt) = 2;
+   *(ptstack + (*npoints))= pt;
+   (*npoints)++;
+}
+
+DEFUN_DLD (bwfill, args, ,
+  "[...] = bwfill (...)
+   [BW2,IDX] = BWFILL(BW1,Y,X,N) performs a flood-fill on BW1
+       
+       (X(k), Y(k)) are rows and columns of seed points
+
+   [BW2,IDX] = BWFILL(BW1,'holes',N) fills interior holes in BW1
+
+       N = 4 or 8(default) for neighborhood connectedness
+
+       IDX is the indices of the filled pixels
+       ")
+{
+   octave_value_list retval;
+   octave_value tmp;
+   ColumnVector xseed, yseed ;
+   int nargin = args.length ();
+
+   if (nargin < 2 ) {
+      print_usage ("bwfill");
+      return retval;
+   }
+
+   Matrix       im=    args(0).matrix_value();
+   int          imM=   im.rows();
+   int          imN=   im.columns();
+
+   int          nb=    8;
+   int          npoints= 0;
+   bool         fillmode= false;
+   if (args(1).is_string() && args(1).string_value() == "holes" ) {
+      fillmode= true;
+
+      npoints= 2*( imM + imN - 4 ); // don't start fill from corners
+
+      xseed= ColumnVector( npoints );
+      yseed= ColumnVector( npoints );
+      int idx= 0;
+      for (int j=2; j<= imN-1; j++) {
+         xseed( idx )= j;   yseed( idx++)= 1;   
+         xseed( idx )= j;   yseed( idx++)= imM;   
+      }
+
+      for (int i=2; i<= imM-1; i++) {
+         yseed( idx )= i;   xseed( idx++)= 1;   
+         yseed( idx )= i;   xseed( idx++)= imN;   
+      }
+
+      if (nargin >= 4 ) 
+         nb= (int) args(2).double_value();
+   } // holes mode? 
+   else {
+      {
+         ColumnVector tmp( args(2).vector_value() );
+         xseed= tmp;
+      }
+      {
+         ColumnVector tmp( args(1).vector_value() );
+         yseed= tmp;
+      }
+      npoints= xseed.length();
+      if (nargin >= 4 ) 
+         nb= (int) args(3).double_value();
+   } // holes mode? 
+
+/*
+ * put a one pixel thick boundary around the image
+ *  so that we can be more efficient in the main loop
+ */
+   int           ioM=   imM+2;
+   unsigned char imo[ (imM+2) * (imN+2) ];
+
+   for (int i=0; i<imM; i++) 
+      for (int j=0; j<imN; j++)
+         imo[(i+1) + ioM*(j+1)]= ( im(i,j) > 0 ) ;
+
+   for (int i=0; i<ioM; i++) 
+      imo[i]= imo[i + ioM*(imN+1)] = 3;
+
+   for (int j=1; j<imN+1; j++)
+      imo[ioM*j]= imo[imM+1 + ioM*j] = 3;
+
+// This is obviously big enough for the point stack, but I'm
+// sure it can be smaller. 
+   int     ptstack[ ioM*imN ];
+
+   int seedidx= npoints; 
+   npoints= 0;
+   while ( (--seedidx) >= 0 ) {
+// no need to add 1 to convert indexing style because we're adding a boundary
+      int pt= (int) xseed( seedidx )*ioM + (int) yseed( seedidx ); 
+      checkpoint( pt , imo, ptstack, &npoints );
+   }
+
+   while ( npoints > 0 ) {
+      npoints--;
+      int pt= ptstack[ npoints ];
+      
+      checkpoint( pt + LF, imo, ptstack, &npoints );
+      checkpoint( pt + RT, imo, ptstack, &npoints );
+      checkpoint( pt + UP, imo, ptstack, &npoints );
+      checkpoint( pt + DN, imo, ptstack, &npoints );
+      
+      if (nb==8) {
+         checkpoint( pt + LF + UP, imo, ptstack, &npoints );
+         checkpoint( pt + RT + UP, imo, ptstack, &npoints );
+         checkpoint( pt + LF + DN, imo, ptstack, &npoints );
+         checkpoint( pt + RT + DN, imo, ptstack, &npoints );
+      }
+   } // while ( npoints
+
+   Matrix       imout( imM, imN );
+   ColumnVector idxout (imM*imN );
+   int idx=0;
+
+   int notvalidpt= 0;
+   int idxpoint=   2;
+   if ( fillmode ) {
+      notvalidpt= 2;
+      idxpoint=   0;
+   }
+
+   for (int i=0; i<imM; i++) 
+      for (int j=0; j<imN; j++) {
+         imout(i,j) =    (double) ( imo[(i+1) + ioM*(j+1)] != notvalidpt );
+         if ( imo[(i+1) + ioM*(j+1)] == idxpoint )
+            idxout(idx++) = (double) (i + j*imM + 1);
+      }
+
+   /*
+   Matrix imout( imM+2, imN+2 );
+   for (int i=0; i<imM+2; i++) 
+      for (int j=0; j<imN+2; j++)
+         imout(i,j) = (double) imo[i + ioM*j];
+    */
+
+   retval(0)= imout;
+// we need to do this to be able to return a proper empty vector
+   if (idx > 0) 
+      retval(1)= idxout.extract(0, idx-1);
+   else
+      retval(1)= ColumnVector ( 0 );
+   return retval;
+}
+
+
+/*
+ * $Log$
+ * Revision 1.1  2002/03/17 02:38:51  aadler
+ * fill and edge detection operators
+ *
+ * Revision 1.9  2000/06/16 20:22:47  aadler
+ * mods for 2.1/2.0 compat
+ *
+ * Revision 1.8  2000/06/13 17:27:24  aadler
+ * mods for 2.1.30
+ *
+ * Revision 1.7  1999/06/10 19:42:12  aadler
+ * minor verbose fix
+ *
+ * Revision 1.6  1999/06/08 16:30:30  aadler
+ * bug fix. reversed r,c input arguments
+ *
+ * Revision 1.5  1999/06/08 15:41:02  aadler
+ * now fills in holes
+ *
+ * Revision 1.4  1999/06/08 15:21:02  aadler
+ * fixed bug that so specified points are only used if they can fill
+ *
+ * Revision 1.3  1999/06/08 15:05:08  aadler
+ * now returns 1 and gives index output
+ *
+ * Revision 1.2  1999/06/04 21:58:57  aadler
+ * fixed 8 vs 4 neighborhood
+ *
+ * Revision 1.1  1999/06/04 21:43:20  aadler
+ * Initial revision
+ *
+ *
+ */

new file mode 100644
--- /dev/null
+++ b/bwselect.m
@@ -0,0 +1,37 @@
+function [imout, idx] = bwselect( im, cols, rows, connect )
+# BWSELECT: select connected regions in a binary image
+# [imout, idx] = bwselect( im, cols, rows, connect )
+#
+#   im          -> binary input image
+#   [cols,rows] -> vectors of starting points (x,y)
+#   connect     -> connectedness 4 or 8. default is 8
+#   imout       -> the image of all objects in image im that overlap
+#                  pixels in (cols,rows)
+#   idx         -> index of pixels in imout
+
+# Copyright (C) 1999 Andy Adler
+# This code has no warrany whatsoever.
+# Do what you like with this code as long as you
+#     leave this copyright in place.
+#
+# $Id$
+
+if nargin<4
+   connect= 8;
+end
+
+[jnk,idx]= bwfill( ~im, cols,rows, connect );
+
+imout= zeros( size(jnk) );
+imout( idx ) = 1;
+
+# 
+# $Log$
+# Revision 1.1  2002/03/17 02:38:52  aadler
+# fill and edge detection operators
+#
+# Revision 1.1  1999/06/08 17:06:01  aadler
+# Initial revision
+#
+#
+

new file mode 100644
--- /dev/null
+++ b/edge.m
@@ -0,0 +1,240 @@
+function [imout, thresh] = edge( im, method, thresh, param2 )
+# EDGE: find image edges
+# [imout, thresh] = edge( im, method, thresh, param2 )
+#
+# OUTPUT
+#  imout  -> output image
+#  thresh -> output thresholds
+#
+# INPUT
+#  im     -> input image (greyscale)
+#  thresh -> threshold value (value is estimated if not given)
+#  
+# The following methods are based on high pass filtering the image in
+#   two directions, calculating a combined edge weight from and then thresholding
+#
+# method = 'roberts'
+#     filt1= [1 0 ; 0 -1];   filt2= rot90( filt1 )
+#     combine= sqrt( filt1^2 + filt2^2 )  
+# method = 'sobel'
+#     filt1= [1 2 1;0 0 0;-1 -2 -1];      filt2= rot90( filt1 ) 
+#     combine= sqrt( filt1^2 + filt2^2 )  
+# method = 'prewitt'
+#     filt1= [1 1 1;0 0 0;-1 -1 -1];      filt2= rot90( filt1 ) 
+#     combine= sqrt( filt1^2 + filt2^2 )  
+# method = 'kirsh'
+#     filt1= [1 2 1;0 0 0;-1 -2 -1];  filt2 .. filt8 are 45 degree rotations of filt1
+#     combine= max( filt1 ... filt8 )
+#
+# methods based on filtering the image and finding zero crossings
+#
+# method = 'log' -> Laplacian of Gaussians 
+#      param2 is the standard deviation of the filter, default is 2
+# method = 'zerocross' -> generic zero-crossing filter
+#      param2 is the user supplied filter
+# 
+# method = 'andy' -> my idea
+#      A.Adler's idea (c) 1999. somewhat based on the canny method
+#      Step 1: Do a sobel edge detection and to generate an image at
+#               a high and low threshold
+#      Step 2: Edge extend all edges in the LT image by several pixels,
+#               in the vertical, horizontal, and 45degree directions.
+#               Combine these into edge extended (EE) image
+#      Step 3: Dilate the EE image by 1 step
+#      Step 4: Select all EE features that are connected to features in
+#               the HT image
+#                
+#      Parameters:
+#        param2(1)==0 or 4 or 8 -> perform x connected dilatation (step 3)
+#        param2(2)    dilatation coeficient (threshold) in step 3
+#        param2(3)    length of edge extention convolution (step 2)
+#        param2(4)    coeficient of extention convolution in step 2
+#        defaults = [8 1 3 3]
+
+# Copyright (C) 1999 Andy Adler
+# This code has no warrany whatsoever.
+# Do what you like with this code as long as you
+#     leave this copyright in place.
+#
+# $Id$
+
+[n,m]= size(im);
+xx= 2:m-1;
+yy= 2:n-1;
+
+if   strcmp(method,'roberts') || strcmp(method,'sobel') || ...
+     strcmp(method,'prewitt') 
+     
+
+   if strcmp(method,'roberts') 
+      filt= [1 0;0 -1]/4;               tv= 6;
+   elseif strcmp(method,'sobel') 
+      filt= [1 2 1;0 0 0; -1 -2 -1]/8;  tv= 2;
+   elseif strcmp(method,'prewitt') 
+      filt= [1 1 1;0 0 0; -1 -1 -1]/6;  tv= 4;
+   end
+
+   imo= conv2(im, rot90(filt), 'same').^2 + conv2(im, filt, 'same').^2;
+   
+# check to see if the user supplied a threshold
+# if not, calculate one in the same way as Matlab
+
+   if nargin<3
+      thresh= sqrt( tv* mean(mean( imo(yy,xx) ))  );
+   end
+
+# The filters are defined for sqrt(imo), but since we calculated imo, compare
+#  to thresh ^2
+
+   imout= ( imo >= thresh^2 );   
+
+# Thin the wide edges
+   xpeak= imo(yy,xx-1) <= imo(yy,xx) & imo(yy,xx) > imo(yy,xx+1) ;
+   ypeak= imo(yy-1,xx) <= imo(yy,xx) & imo(yy,xx) > imo(yy+1,xx) ;
+   imout(yy,xx)= imout(yy,xx) & ( xpeak | ypeak );
+
+elseif strcmp(method,'kirsch')   
+
+   filt1= [1 2 1;0 0 0;-1 -2 -1];   fim1= conv2(im,filt1,'same');
+   filt2= [2 1 0;1 0 -1;0 -1 -2];   fim2= conv2(im,filt2,'same');
+   filt3= [1 0 -1;2 0 -2;1 0 -1];   fim3= conv2(im,filt3,'same');
+   filt4= [0 1 2;-1 0 1;-2 -1 0];   fim4= conv2(im,filt4,'same');
+
+   imo= reshape(max([abs(fim1(:)) abs(fim2(:)) abs(fim3(:)) abs(fim4(:))]'),n,m);
+
+   if nargin<3
+      thresh=  2* mean(mean( imo(yy,xx) )) ;
+   end
+
+   imout=  imo >= thresh ;   
+
+# Thin the wide edges
+   xpeak= imo(yy,xx-1) <= imo(yy,xx) & imo(yy,xx) > imo(yy,xx+1) ;
+   ypeak= imo(yy-1,xx) <= imo(yy,xx) & imo(yy,xx) > imo(yy+1,xx) ;
+   imout(yy,xx)= imout(yy,xx) & ( xpeak | ypeak );
+
+elseif  strcmp(method,'log') || strcmp(method,'zerocross') 
+
+   if strcmp(method,'log') 
+      if nargin >= 4;    sd= param2;
+      else               sd= 2;
+      end
+
+      sz= ceil(sd*3);
+      [x,y]= meshgrid( -sz:sz, -sz:sz );
+      filt = exp( -( x.^2 + y.^2 )/2/sd^2 ) .* ...
+                   ( x.^2 + y.^2 - 2*sd^2 ) / 2 / pi / sd^6 ;
+   else
+      filt = param2;
+   end
+   filt = filt - mean(filt(:));
+
+   imo= conv2(im, filt, 'same');
+
+   if nargin<3 || isempty( thresh )
+      thresh=  0.75* mean(mean( abs(imo(yy,xx)) )) ;
+   end
+
+   zcross= imo > 0;
+   yd_zc=  diff( zcross );
+   xd_zc=  diff( zcross' )';
+   yd_io=  abs(diff( imo ) ) > thresh;
+   xd_io=  abs(diff( imo')') > thresh;
+
+# doing it this way puts the transition at the <=0 point
+   xl= zeros(1,m);  yl= zeros(n,1);
+   imout= [    ( yd_zc ==  1 ) & yd_io ; xl] | ...
+          [xl; ( yd_zc == -1 ) & yd_io     ] | ...
+          [    ( xd_zc ==  1 ) & xd_io , yl] | ... 
+          [yl, ( xd_zc == -1 ) & xd_io     ];
+
+elseif  strcmp(method,'canny')  
+    error("method canny not implemented");
+
+elseif  strcmp(method,'andy')  
+
+   filt= [1 2 1;0 0 0; -1 -2 -1]/8;  tv= 2;
+   imo= conv2(im, rot90(filt), 'same').^2 + conv2(im, filt, 'same').^2;
+   if nargin<3 || thresh==[];
+      thresh= sqrt( tv* mean(mean( imo(yy,xx) ))  );
+   end
+#     sum( imo(:)>thresh ) / prod(size(imo))
+   dilate= [1 1 1;1 1 1;1 1 1]; tt= 1; sz=3; dt=3;
+   if nargin>=4
+      # 0 or 4 or 8 connected dilation
+      if length(param2) > 0
+         if      param2(1)==4 ; dilate= [0 1 0;1 1 1;0 1 0];
+         elseif  param2(1)==0 ; dilate= 1;
+         end
+      end
+      # dilation threshold
+      if length(param2) > 2; tt= param2(2); end
+      # edge extention length
+      if length(param2) > 2; sz= param2(3); end
+      # edge extention threshold
+      if length(param2) > 3; dt= param2(4); end
+      
+   end
+   fobliq= [0 0 0 0 1;0 0 0 .5 .5;0 0 0 1 0;0 0 .5 .5 0;0 0 1 0 0; 
+                      0 .5 .5 0 0;0 1 0 0 0;.5 .5 0 0 0;1 0 0 0 0];
+   fobliq= fobliq( 5-sz:5+sz, 3-ceil(sz/2):3+ceil(sz/2) );
+
+   xpeak= imo(yy,xx-1) <= imo(yy,xx) & imo(yy,xx) > imo(yy,xx+1) ;
+   ypeak= imo(yy-1,xx) <= imo(yy,xx) & imo(yy,xx) > imo(yy+1,xx) ;
+
+   imht= ( imo >= thresh^2 * 2); # high threshold image   
+   imht(yy,xx)= imht(yy,xx) & ( xpeak | ypeak );
+   imht([1,n],:)=0; imht(:,[1,m])=0;
+
+%  imlt= ( imo >= thresh^2 / 2); # low threshold image   
+   imlt= ( imo >= thresh^2 / 1); # low threshold image   
+   imlt(yy,xx)= imlt(yy,xx) & ( xpeak | ypeak );
+   imlt([1,n],:)=0; imlt(:,[1,m])=0;
+
+# now we edge extend the low thresh image in 4 directions
+
+   imee= ( conv2( imlt, ones(2*sz+1,1)    , 'same') > tt ) | ...
+         ( conv2( imlt, ones(1,2*sz+1)    , 'same') > tt ) | ...
+         ( conv2( imlt, eye(2*sz+1)       , 'same') > tt ) | ...
+         ( conv2( imlt, rot90(eye(2*sz+1)), 'same') > tt ) | ...
+         ( conv2( imlt, fobliq            , 'same') > tt ) | ...
+         ( conv2( imlt, fobliq'           , 'same') > tt ) | ...
+         ( conv2( imlt, rot90(fobliq)     , 'same') > tt ) | ...
+         ( conv2( imlt, flipud(fobliq)    , 'same') > tt );
+#  imee(yy,xx)= conv2(imee(yy,xx),ones(3),'same') & ( xpeak | ypeak );
+   imee= conv2(imee,dilate,'same') > dt; #
+
+%  ff= find( imht==1 );
+%  imout = bwselect( imee, rem(ff-1, n)+1, ceil(ff/n), 8);  
+   imout = imee;
+
+
+else
+
+   error (['Method ' method ' is not recognized']);
+
+end
+
+
+
+
+# 
+# $Log$
+# Revision 1.1  2002/03/17 02:38:52  aadler
+# fill and edge detection operators
+#
+# Revision 1.4  2000/11/20 17:13:07  aadler
+# works?
+#
+# Revision 1.3  1999/06/09 17:29:36  aadler
+# implemented 'andy' mode edge detection
+#
+# Revision 1.2  1999/06/08 14:26:50  aadler
+# zero-cross and LoG filters work
+#
+# Revision 1.1  1999/06/07 21:01:38  aadler
+# Initial revision
+#
+#
+
+#