octave_packages/image-1.0.15/bwmorph.m

   1 ## Copyright (C) 2004 Josep Mones i Teixidor
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify
   4 ## it under the terms of the GNU General Public License as published by
   5 ## the Free Software Foundation; either version 2 of the License, or
   6 ## (at your option) any later version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful,
   9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 ## GNU General Public License for more details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License
  14 ## along with this program; If not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ## @deftypefn {Function File} {@var{BW2} = } bwmorph (@var{BW},@var{operation})
  18 ## @deftypefnx {Function File} {@var{BW2} = } bwmorph (@var{BW},@var{operation},@var{n})
  19 ## Perform a morphological operation on a binary image.
  20 ##
  21 ## BW2=bwmorph(BW,operation) performs a morphological operation
  22 ## specified by @var{operation} on binary image @var{BW}. All possible
  23 ## operations and their meaning are specified in a table below.
  24 ##
  25 ## BW2=bwmorph(BW,operation,n) performs a morphological operation
  26 ## @var{n} times. Keep in mind that it has no sense to apply some
  27 ## operations more than once, since some of them return the same result
  28 ## regardless how many iterations we request. Those return a warning if
  29 ## are called with n>1 and they compute the result for n=1.
  30 ##
  31 ## @var{n}>1 is actually used for the following operations: diag,
  32 ## dilate, erode, majority, shrink, skel, spur, thicken and thin.
  33 ##
  34 ## @table @code
  35 ## @item 'bothat'
  36 ## Performs a bottom hat operation, a closing operation (which is a
  37 ## dilation followed by an erosion) and finally substracts the original
  38 ## image.
  39 ##
  40 ## @item 'bridge'
  41 ## Performs a bridge operation. Sets a pixel to 1 if it has two nonzero
  42 ## neighbours which are not connected, so it "bridges" them. There are
  43 ## 119 3-by-3 patterns which trigger setting a pixel to 1.
  44 ##
  45 ## @item 'clean'
  46 ## Performs an isolated pixel remove operation. Sets a pixel to 0 if all
  47 ## of its eight-connected neighbours are 0.
  48 ##
  49 ## @item 'close'
  50 ## Performs closing operation, which is a dilation followed by erosion.
  51 ## It uses a ones(3) matrix as structuring element for both operations.
  52 ##
  53 ## @item 'diag'
  54 ## Performs a diagonal fill operation. Sets a pixel to 1 if that
  55 ## eliminates eight-connectivity of the background.
  56 ##
  57 ## @item 'dilate'
  58 ## Performs a dilation operation. It uses ones(3) as structuring element.
  59 ##
  60 ## @item 'erode'
  61 ## Performs an erosion operation. It uses ones(3) as structuring element.
  62 ##
  63 ## @item 'fill'
  64 ## Performs a interior fill operation. Sets a pixel to 1 if all
  65 ## four-connected pixels are 1.
  66 ##
  67 ## @item 'hbreak'
  68 ## Performs a H-break operation. Breaks (sets to 0) pixels that are
  69 ## H-connected.
  70 ##
  71 ## @item 'majority'
  72 ## Performs a majority black operation. Sets a pixel to 1 if five
  73 ## or more pixels in a 3-by-3 window are 1. If not it is set to 0.
  74 ##
  75 ## @item 'open'
  76 ## Performs an opening operation, which is an erosion followed by a
  77 ## dilation. It uses ones(3) as structuring element.
  78 ##
  79 ## @item 'remove'
  80 ## Performs a iterior pixel remove operation. Sets a pixel to 0 if
  81 ## all of its four-connected neighbours are 1.
  82 ##
  83 ## @item 'shrink'
  84 ## Performs a shrink operation. Sets pixels to 0 such that an object
  85 ## without holes erodes to a single pixel (set to 1) at or near its
  86 ## center of mass. An object with holes erodes to a connected ring lying
  87 ## midway between each hole and its nearest outer boundary. It preserves
  88 ## Euler number.
  89 ##
  90 ## @item 'skel'
  91 ## Performs a skeletonization operation. It calculates a "median axis
  92 ## skeleton" so that points of this skeleton are at the same distance of
  93 ## its nearby borders. It preserver Euler number. Please read
  94 ## compatibility notes for more info.
  95 ##
  96 ## It uses the same algorithm as skel-pratt but this could change for
  97 ## compatibility in the future.
  98 ##
  99 ## @item 'skel-lantuejol'
 100 ## Performs a skeletonization operation as described in Gonzalez & Woods
 101 ## "Digital Image Processing" pp 538-540. The text references Lantuejoul
 102 ## as authour of this algorithm.
 103 ##
 104 ## It has the beauty of being a clean and simple approach, but skeletons
 105 ## are thicker than they need to and, in addition, not guaranteed to be
 106 ## connected.
 107 ##
 108 ## This algorithm is iterative. It will be applied the minimum value of
 109 ## @var{n} times or number of iterations specified in algorithm
 110 ## description. It's most useful to run this algorithm with @code{n=Inf}.
 111 ##
 112 ## @item 'skel-pratt'
 113 ## Performs a skeletonization operation as described by William K. Pratt
 114 ## in "Digital Image Processing".
 115 ##
 116 ## @item 'spur'
 117 ## Performs a remove spur operation. It sets pixel to 0 if it has only
 118 ## one eight-connected pixel in its neighbourhood.
 119 ##
 120 ## @item 'thicken'
 121 ## Performs a thickening operation. This operation "thickens" objects
 122 ## avoiding their fusion. Its implemented as a thinning of the
 123 ## background. That is, thinning on negated image. Finally a diagonal
 124 ## fill operation is performed to avoid "eight-connecting" objects.
 125 ##
 126 ## @item 'thin'
 127 ## Performs a thinning operation. When n=Inf, thinning sets pixels to 0
 128 ## such that an object without holes is converted to a stroke
 129 ## equidistant from its nearest outer boundaries. If the object has
 130 ## holes it creates a ring midway between each hole and its near outer
 131 ## boundary. This differ from shrink in that shrink converts objects
 132 ## without holes to a single pixels and thin to a stroke. It preserves
 133 ## Euler number.
 134 ##
 135 ## @item 'tophat'
 136 ## Performs a top hat operation, a opening operation (which is an
 137 ## erosion followed by a dilation) and finally substracts the original
 138 ## image.
 139 ## @end table
 140 ##
 141 ## Some useful concepts to understant operators:
 142 ##
 143 ## Operations are defined on 3-by-3 blocks of data, where the pixel in
 144 ## the center of the block. Those pixels are numerated as follows:
 145 ##
 146 ## @multitable @columnfractions 0.05 0.05 0.05
 147 ## @item X3 @tab X2 @tab X1
 148 ## @item X4 @tab  X @tab X0
 149 ## @item X5 @tab X6 @tab X7
 150 ## @end multitable
 151 ##
 152 ## @strong{Neighbourhood definitions used in operation descriptions:}
 153 ## @table @code
 154 ## @item 'four-connected'
 155 ## It refers to pixels which are connected horizontally or vertically to
 156 ## X: X1, X3, X5 and X7.
 157 ## @item 'eight-connected'
 158 ## It refers to all pixels which are connected to X: X0, X1, X2, X3, X4,
 159 ## X5, X6 and X7.
 160 ## @end table
 161 ##
 162 ## @strong{Compatibility notes:}
 163 ## @table @code
 164 ## @item 'fill'
 165 ## Checking MATLAB behaviour is needed because its documentation doesn't
 166 ## make clear if it creates a black pixel if all eight-connected pixels
 167 ## are black or if four-connected suffice (as we do currently following
 168 ## Pratt's book).
 169 ## @item 'skel'
 170 ## Algorithm used here is described in Pratt's book. When applying it to
 171 ## the "circles" image in MATLAB documentation, results are not the
 172 ## same. Perhaps MATLAB uses Blum's algoritm (for further info please
 173 ## read comments in code).
 174 ## @item 'skel-pratt'
 175 ## This option is not available in MATLAB.
 176 ## @item 'skel-lantuejoul'
 177 ## This option is not available in MATLAB.
 178 ## @item 'thicken'
 179 ## This implementation also thickens image borders. This can easily be
 180 ## avoided i necessary. MATLAB documentation doesn't state how it behaves.
 181 ## @end table
 182 ##
 183 ## References:
 184 ## W. K. Pratt, "Digital Image Processing"
 185 ## Gonzalez and Woods, "Digital Image Processing"
 186 ##
 187 ## @seealso{dilate, erode, makelut, applylut}
 188 ## @end deftypefn
 189
 190
 191 ## TODO: As soon as Octave doesn't segfault when assigning values to a
 192 ## TODO: bool matrix, remove all conversions from lut to logical and
 193 ## TODO: just create it as a logical from the beginning.
 194
 195 ## TODO: n behaviour should be tested in all cases for compatibility.
 196
 197 ## Author:  Josep Mones i Teixidor <jmones@puntbarra.com>
 198
 199 function BW2 = bwmorph(BW, operation, n)
 200   if(nargin<2 || nargin>3)
 201     usage("BW2=bwmorph(BW, operation [,n])");
 202   endif
 203   if(nargin<3)
 204     n=1;
 205   endif
 206   if(n<0)
 207     error("bwmorph: n should be > 0");
 208   elseif(n==0) ## we'll just return the same matrix (check this!)
 209     BW2=BW;
 210   endif
 211
 212   ## post processing command
 213   postcmd="";
 214
 215   switch(operation)
 216     case('bothat')
 217       se=ones(3);
 218       BW2=erode(dilate(BW, se), se)-BW;
 219       if(n>1)
 220         ## TODO: check if ignoring n>1 is ok. Should I just ignore it
 221         ## TODO: without a warning?
 222         disp("WARNING: n>1 has no sense here. Using n=1. Please fill a bug if you think this behaviour is not correct");
 223       endif
 224       return;
 225
 226     case('bridge')
 227       ## see __bridge_lut_fun__ for rules
 228       ## lut=makelut("__bridge_lut_fun__",3);
 229       lut=logical([0;0;0;0;0;0;0;0;0;0;0;0;1;1;0;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 230                    0;1;0;0;0;1;0;0;1;1;0;0;1;1;0;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 231                    0;0;1;1;1;1;1;1;0;0;0;0;1;1;0;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 232                    1;1;1;1;1;1;1;1;1;1;0;0;1;1;0;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 233                    0;1;1;1;1;1;1;1;0;0;0;0;1;1;0;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 234                    0;1;0;0;0;1;0;0;0;0;0;0;0;0;0;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 235                    0;1;1;1;1;1;1;1;0;0;0;0;1;1;0;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 236                    0;1;0;0;0;1;0;0;0;0;0;0;0;0;0;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 237                    0;1;1;1;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 238                    0;1;0;0;0;1;0;0;1;1;0;0;1;1;0;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 239                    0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 240                    1;1;1;1;1;1;1;1;1;1;0;0;1;1;0;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 241                    0;1;1;1;1;1;1;1;0;0;0;0;1;1;0;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 242                    0;1;0;0;0;1;0;0;0;0;0;0;0;0;0;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 243                    0;1;1;1;1;1;1;1;0;0;0;0;1;1;0;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 244                    0;1;0;0;0;1;0;0;0;0;0;0;0;0;0;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1]);
 245       BW2=applylut(BW, lut);
 246       if(n>1)
 247         ## TODO: check if ignoring n>1 is ok.
 248         disp("WARNING: n>1 has no sense here. Using n=1. Please fill a bug if you think this behaviour is not correct");
 249       endif
 250       return;
 251
 252     case('clean')
 253       ## BW(j,k)=X&&(X0||X1||...||X7)
 254       ## lut=makelut(inline("x(2,2)&&any((x.*[1,1,1;1,0,1;1,1,1])(:))","x"),3);
 255       ## which is the same as...
 256       lut=repmat([zeros(16,1);ones(16,1)],16,1);  ## identity
 257       lut(17)=0;                                  ## isolated to 0
 258       ## I'd prefer to create lut directly as a logical, but assigning a
 259       ## value to a logical segfaults 2.1.57. We'll change it as soon as
 260       ## it works.
 261       BW2=applylut(BW, logical(lut));
 262       if(n>1)
 263         ## TODO: check if ignoring n>1 is ok.
 264         disp("WARNING: n>1 has no sense here. Using n=1. Please fill a bug if you think this behaviour is not correct");
 265       endif
 266       return;
 267
 268     case('close')
 269       se=ones(3);
 270       BW2=erode(dilate(BW, se), se);
 271       if(n>1)
 272         ## TODO: check if ignoring n>1 is ok.
 273         disp("WARNING: n>1 has no sense here. Using n=1. Please fill a bug if you think this behaviour is not correct");
 274       endif
 275       return;
 276
 277     case('diag')
 278       ## see __diagonal_fill_lut_fun__ for rules
 279       ## lut=makelut("__diagonal_fill_lut_fun__",3);
 280       lut=logical([0;0;0;0;0;0;0;0;0;0;1;0;0;0;1;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 281                    0;0;1;1;0;0;0;0;0;0;1;1;0;0;1;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 282                    0;0;0;0;0;0;0;0;0;0;1;0;0;0;1;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 283                    0;0;1;1;0;0;0;0;0;0;1;1;0;0;1;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 284                    0;0;0;0;0;0;0;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 285                    1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 286                    0;0;0;0;0;0;0;0;0;0;1;0;0;0;1;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 287                    1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 288                    0;0;0;0;0;0;0;0;0;0;1;0;0;0;1;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 289                    0;0;1;1;0;0;0;0;0;0;1;1;0;0;1;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 290                    0;0;0;0;0;0;0;0;0;0;1;0;0;0;1;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 291                    0;0;1;1;0;0;0;0;0;0;1;1;0;0;1;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 292                    0;0;0;0;0;0;0;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 293                    0;0;1;1;0;0;0;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 294                    0;0;0;0;0;0;0;0;0;0;1;0;0;0;1;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 295                    0;0;1;1;0;0;0;0;0;0;1;1;0;0;1;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1]);
 296       cmd="BW2=applylut(BW, lut);";
 297
 298
 299     case('dilate')
 300       cmd="BW2=dilate(BW, ones(3));";
 301
 302     case('erode')
 303       cmd="BW2=erode(BW, ones(3));";
 304
 305     case('fill')
 306       ## lut=makelut(inline("x(2,2)||(sum((x&[0,1,0;1,0,1;0,1,0])(:))==4)","x"),3);
 307       ## which is the same as...
 308       lut=repmat([zeros(16,1);ones(16,1)],16,1); ## identity
 309       ## 16 exceptions
 310       lut([171,172,175,176,235,236,239,240,427,428,431,432,491,492,495,496])=1;
 311       BW2=applylut(BW, logical(lut));
 312       if(n>1)
 313         ## TODO: check if ignoring n>1 is ok.
 314         disp("WARNING: n>1 has no sense here. Using n=1. Please fill a bug if you think this behaviour is not correct");
 315       endif
 316       return;
 317
 318
 319     case('hbreak')
 320       ## lut=makelut(inline("x(2,2)&&!(all(x==[1,1,1;0,1,0;1,1,1])||all(x==[1,0,1;1,1,1;1,0,1]))","x"),3);
 321       ## which is the same as
 322       lut=repmat([zeros(16,1);ones(16,1)],16,1); ## identity
 323       lut([382,472])=0;                          ## the 2 exceptions
 324       BW2=applylut(BW, logical(lut));
 325       if(n>1)
 326         ## TODO: check if ignoring n>1 is ok.
 327         disp("WARNING: n>1 has no sense here. Using n=1. Please fill a bug if you think this behaviour is not correct");
 328       endif
 329       return;
 330
 331     case('majority')
 332       ## lut=makelut(inline("sum((x&ones(3,3))(:))>=5"),3);
 333       lut=logical([0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;
 334                    0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;0;0;0;0;0;0;0;1;0;0;0;1;0;1;1;1;
 335                    0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;0;0;0;0;0;0;0;1;0;0;0;1;0;1;1;1;
 336                    0;0;0;0;0;0;0;1;0;0;0;1;0;1;1;1;0;0;0;1;0;1;1;1;0;1;1;1;1;1;1;1;
 337                    0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;0;0;0;0;0;0;0;1;0;0;0;1;0;1;1;1;
 338                    0;0;0;0;0;0;0;1;0;0;0;1;0;1;1;1;0;0;0;1;0;1;1;1;0;1;1;1;1;1;1;1;
 339                    0;0;0;0;0;0;0;1;0;0;0;1;0;1;1;1;0;0;0;1;0;1;1;1;0;1;1;1;1;1;1;1;
 340                    0;0;0;1;0;1;1;1;0;1;1;1;1;1;1;1;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 341                    0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;0;0;0;0;0;0;0;1;0;0;0;1;0;1;1;1;
 342                    0;0;0;0;0;0;0;1;0;0;0;1;0;1;1;1;0;0;0;1;0;1;1;1;0;1;1;1;1;1;1;1;
 343                    0;0;0;0;0;0;0;1;0;0;0;1;0;1;1;1;0;0;0;1;0;1;1;1;0;1;1;1;1;1;1;1;
 344                    0;0;0;1;0;1;1;1;0;1;1;1;1;1;1;1;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 345                    0;0;0;0;0;0;0;1;0;0;0;1;0;1;1;1;0;0;0;1;0;1;1;1;0;1;1;1;1;1;1;1;
 346                    0;0;0;1;0;1;1;1;0;1;1;1;1;1;1;1;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 347                    0;0;0;1;0;1;1;1;0;1;1;1;1;1;1;1;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 348                    0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1]);
 349       cmd="BW2=applylut(BW, lut);";
 350
 351     case('open')
 352       se=ones(3);
 353       BW2=dilate(erode(BW, se), se);
 354       if(n>1)
 355         ## TODO: check if ignoring n>1 is ok.
 356         disp("WARNING: n>1 has no sense here. Using n=1. Please fill a bug if you think this behaviour is not correct");
 357       endif
 358       return;
 359
 360     case('remove')
 361       ## lut=makelut(inline("x(2,2)&&!(sum((x&[0,1,0;1,1,1;0,1,0])(:))==5)","x"),3);
 362       lut=repmat([zeros(16,1);ones(16,1)],16,1); ## identity
 363       ## 16 qualifying patterns
 364       lut([187,188,191,192,251,252,255,256,443,444,447,448,507,508,511,512])=0;
 365       BW2=applylut(BW, logical(lut));
 366       if(n>1)
 367         ## TODO: check if ignoring n>1 is ok.
 368         disp("WARNING: n>1 has no sense here. Using n=1. Please fill a bug if you think this behaviour is not correct");
 369       endif
 370       return;
 371
 372     case('shrink')
 373       ## lut1=makelut("__conditional_mark_patterns_lut_fun__",3,"S");
 374       lut1=logical([0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;1;1;1;0;1;1;1;1;0;1;0;0;1;1;
 375                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;0;0;1;1;0;1;1;0;0;0;0;0;0;0;1;
 376                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;0;0;0;0;0;0;0;1;1;0;1;0;0;0;1;
 377                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;
 378                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;
 379                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;0;0;0;0;0;0;0;0;0;0;0;
 380                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;0;0;0;0;0;0;0;1;1;0;1;0;0;0;1;
 381                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;0;0;0;0;0;0;0;0;0;0;0;
 382                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;
 383                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;0;1;1;1;0;1;1;0;0;0;0;0;0;0;1;
 384                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;
 385                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;
 386                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;
 387                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;0;0;0;1;0;1;1;0;0;0;0;0;0;0;0;
 388                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;0;0;0;0;0;0;0;1;1;0;1;0;0;0;1;
 389                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;0;0;0;1;0;1;1;1;1;0;0;1;1;0;0]);
 390       ## lut2=makelut(inline("!m(2,2)||__unconditional_mark_patterns_lut_fun__(m,'S')","m"),3);
 391       lut2=logical([1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;0;1;0;1;0;0;0;1;0;0;0;0;0;1;0;
 392                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;0;0;1;0;0;0;0;0;0;1;1;0;0;1;0;
 393                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;0;0;0;0;1;1;1;0;0;0;0;1;1;0;1;
 394                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;1;1;1;1;1;1;1;0;1;1;1;0;1;0;1;
 395                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;0;0;0;0;1;0;1;0;0;1;1;1;1;1;1;
 396                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;1;1;1;1;1;1;1;1;1;0;0;1;1;0;1;
 397                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;0;0;0;1;1;0;1;0;0;1;0;1;1;0;1;
 398                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;1;1;1;0;1;0;1;1;1;0;1;0;1;0;1;
 399                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;0;1;0;1;0;1;0;1;1;1;1;1;1;1;
 400                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;1;1;0;0;1;0;1;0;0;1;0;1;1;1;1;
 401                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;1;1;1;1;1;1;1;0;1;1;1;1;1;1;1;
 402                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;1;1;1;1;1;1;1;0;1;1;1;1;1;1;1;
 403                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;1;0;0;0;1;0;1;0;0;1;0;1;1;1;1;
 404                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 405                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;1;1;1;1;1;1;1;0;1;1;1;1;1;1;1;
 406                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1]);
 407       cmd="BW2=BW&applylut(applylut(BW, lut1), lut2);";
 408
 409     case({'skel','skel-pratt'})
 410       ## WARNING: Result doesn't look as MATLAB's sample. It has been
 411       ## WARNING: coded following Pratt's guidelines for what he calls
 412       ## WARNING: is a "reasonably close approximation". I couldn't find
 413       ## WARNING: any bug.
 414       ## WARNING: Perhaps MATLAB uses Blum's algorithm (which Pratt
 415       ## WARNING: refers to) in: H. Blum, "A Transformation for
 416       ## WARNING: Extracting New Descriptors of Shape", Symposium Models
 417       ## WARNING: for Perception of Speech and Visual Form, W.
 418       ## WARNING: Whaten-Dunn, Ed. MIT Press, Cambridge, MA, 1967.
 419
 420       ## lut1=makelut("__conditional_mark_patterns_lut_fun__",3,"K");
 421       lut1=logical([0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;0;0;1;0;0;0;0;1;
 422                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;0;0;0;0;1;0;0;0;0;0;0;0;1;
 423                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;0;1;0;0;0;1;
 424                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;
 425                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;0;0;0;0;0;0;0;
 426                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;
 427                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;0;1;0;0;0;1;
 428                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;
 429                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;
 430                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;0;1;1;0;0;0;0;0;0;0;1;
 431                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;
 432                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;
 433                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;
 434                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;0;1;1;0;0;0;0;0;0;0;1;
 435                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;0;0;0;0;0;0;0;1;1;0;1;0;0;0;1;
 436                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;0;0;0;1;0;1;1;1;1;0;1;1;1;1;0]);
 437
 438       ## lut2=makelut(inline("!m(2,2)||__unconditional_mark_patterns_lut_fun__(m,'K')","m"),3);
 439       lut2=logical([1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;1;1;0;1;0;0;0;1;0;1;1;0;0;0;1;
 440                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;1;0;0;0;1;1;0;0;1;1;0;0;1;1;
 441                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;0;0;0;1;0;1;0;0;0;1;0;1;0;1;
 442                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;1;1;1;0;1;1;1;0;1;1;1;0;1;1;1;
 443                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;0;0;0;1;0;1;1;0;1;1;1;1;1;1;
 444                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;0;1;1;1;1;1;1;1;1;1;1;1;
 445                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;0;0;0;0;1;0;1;1;1;1;1;1;1;1;1;
 446                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;1;1;1;0;1;1;1;1;1;1;1;1;1;1;1;
 447                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;0;0;0;1;0;1;0;0;1;1;1;1;1;1;
 448                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;0;0;0;0;1;1;1;0;0;1;1;1;1;1;1;
 449                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;1;1;1;1;1;1;1;0;1;1;1;1;1;1;1;
 450                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;1;1;1;1;1;1;1;0;1;1;1;1;1;1;1;
 451                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;0;0;0;0;1;0;1;0;0;1;1;1;1;1;1;
 452                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 453                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 454                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1]);
 455       cmd="BW2=BW&applylut(applylut(BW, lut1), lut2);";
 456       postcmd="BW2=bwmorph(BW2,'bridge');";
 457
 458     case('skel-lantuejoul')
 459       ## init values
 460       se=ones(3,3);              ## structuring element used everywhere
 461       BW2=zeros(size(BW));       ## skeleton result
 462       eBW=BW;                    ## eBW will hold k-times eroded BW
 463       i=1;
 464       while i<=n
 465         if(!any(eBW))            ## if erosion result is 0-matrix then
 466           break;                 ## we are over
 467         endif
 468         BW2|=eBW-dilate(erode(eBW, se), se); ## eBW - opening operation on eBW
 469                                              ## contributes to skeleton
 470         eBW=erode(eBW,se);
 471         i++;
 472       endwhile
 473       return;                    ## no general loop in this case
 474
 475     case('spur')
 476       ## lut=makelut(inline("xor(x(2,2),(sum((x&[0,1,0;1,0,1;0,1,0])(:))==0)&&(sum((x&[1,0,1;0,0,0;1,0,1])(:))==1)&&x(2,2))","x"),3);
 477       ## which is the same as
 478       lut=repmat([zeros(16,1);ones(16,1)],16,1); ## identity
 479       lut([18,21,81,273])=0; ## 4 qualifying patterns
 480       lut=logical(lut);
 481       cmd="BW2=applylut(BW, lut);";
 482
 483     case('thicken')
 484       ## This implementation also "thickens" the border. To avoid this,
 485       ## a simple solution could be to add a border of 1 to the reversed
 486       ## image.
 487       BW2=bwmorph(!BW,'thin',n);
 488       BW2=bwmorph(BW2,'diag');
 489       return;
 490
 491     case('thin')
 492       ## lut1=makelut("__conditional_mark_patterns_lut_fun__",3,"T");
 493       lut1=logical([0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;0;0;1;1;0;0;1;1;
 494                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;1;0;0;1;1;0;0;0;0;0;0;0;1;
 495                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;0;1;0;0;0;1;
 496                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;
 497                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;0;0;0;0;0;0;0;
 498                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;0;0;0;1;0;0;0;0;0;0;0;0;0;0;0;
 499                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;1;0;1;0;0;0;1;
 500                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;0;0;0;0;0;0;0;0;0;0;0;
 501                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;
 502                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;1;1;0;1;1;0;0;0;0;0;0;0;1;
 503                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;
 504                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;
 505                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;
 506                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;0;0;0;1;0;1;1;0;0;0;0;0;0;0;0;
 507                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;0;0;0;0;0;0;0;1;1;0;1;0;0;0;1;
 508                     0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;1;0;0;0;1;0;1;1;1;1;0;0;1;1;0;0]);
 509       ## lut2=makelut(inline("!m(2,2)||__unconditional_mark_patterns_lut_fun__(m,'T')","m"),3);
 510       lut2=logical([1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;0;1;1;1;0;1;0;1;1;0;0;0;0;1;0;
 511                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;0;0;1;1;0;0;0;0;0;1;1;0;0;1;0;
 512                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;0;0;0;0;1;1;1;1;0;0;0;1;1;0;1;
 513                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;1;1;1;1;1;1;1;0;1;1;1;0;1;0;1;
 514                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;0;0;0;0;1;0;1;0;0;1;1;1;1;1;1;
 515                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;1;1;1;1;1;1;1;1;1;0;0;1;1;0;1;
 516                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;0;0;1;1;0;1;0;0;1;0;1;1;0;1;
 517                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;1;1;1;0;1;0;1;1;1;0;1;0;1;0;1;
 518                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;0;1;0;1;0;1;0;1;1;1;1;1;1;1;
 519                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;0;1;0;1;0;0;1;0;1;1;1;1;
 520                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;1;1;1;1;1;1;1;0;1;1;1;1;1;1;1;
 521                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;1;1;1;1;1;1;1;0;1;1;1;1;1;1;1;
 522                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;0;0;1;0;1;0;0;1;0;1;1;1;1;
 523                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;
 524                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0;1;1;1;1;1;1;1;0;1;1;1;1;1;1;1;
 525                     1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1]);
 526       cmd="BW2=BW&applylut(applylut(BW, lut1), lut2);";
 527
 528
 529     case('tophat')
 530       se=ones(3);
 531       BW2=BW-dilate(erode(BW, se), se);
 532       if(n>1)
 533         ## TODO: check if ignoring n>1 is ok.
 534         disp("WARNING: n>1 has no sense here. Using n=1. Please fill a bug if you think this behaviour is not correct");
 535       endif
 536       return;
 537
 538
 539     otherwise
 540       error("bwmorph: unknown operation type requested.");
 541   endswitch
 542
 543   ## we use this assignment because of the swap operation inside the
 544   ## while.
 545   BW2=BW;
 546
 547   ## if it doesn't change we don't need to process it further
 548   i=1;
 549   while(i<=n) ## for wouldn't work because n can be Inf
 550     [BW,BW2]=swap(BW,BW2);
 551     eval(cmd);
 552     if(all((BW2==BW)(:)))
 553       break
 554     endif
 555     i+=1;
 556   endwhile
 557
 558   ## process post processing commands if needed
 559   if (isempty (postcmd))
 560     eval(postcmd);
 561   endif
 562
 563 endfunction
 564
 565 function [b, a] = swap (a, b)
 566 endfunction
 567
 568 %!demo
 569 %! bwmorph(ones(11),'shrink', Inf)
 570 %! # Should return 0 matrix with 1 pixel set to 1 at (6,6)
 571
 572 ## TODO: code tests
 573
 574
 575 ## Test skel-lantuejoul using Gozalez&Woods example (fig 8.39)
 576 %!shared slBW, rslBW
 577 %! uint8(0); # fail for 2.1.57 or less instead of crashing later
 578 %! slBW=logical(zeros(12,7));
 579 %! slBW(2,2)=true;
 580 %! slBW(3:4,3:4)=true;
 581 %! rslBW=slBW;
 582 %! slBW(5:6,3:5)=true;
 583 %! slBW(7:11,2:6)=true;
 584 %! rslBW([6,7,9],4)=true;
 585
 586 %!assert(bwmorph(slBW,'skel-lantuejoul',1),[rslBW(1:5,:);logical(zeros(7,7))]);
 587 %!assert(bwmorph(slBW,'skel-lantuejoul',2),[rslBW(1:8,:);logical(zeros(4,7))]);
 588 %!assert(bwmorph(slBW,'skel-lantuejoul',3),rslBW);
 589 %!assert(bwmorph(slBW,'skel-lantuejoul',Inf),rslBW);
 590
 591
 592 %
 593 % $Log$
 594 % Revision 1.4  2007/03/23 16:14:36  adb014
 595 % Update the FSF address
 596 %
 597 % Revision 1.3  2007/01/04 23:41:47  hauberg
 598 % Minor changes in help text
 599 %
 600 % Revision 1.2  2006/10/15 08:04:55  hauberg
 601 % Fixed texinfo in bwmorph
 602 %
 603 % Revision 1.1  2006/08/20 12:59:32  hauberg
 604 % Changed the structure to match the package system
 605 %
 606 % Revision 1.6  2004/09/16 02:14:40  pkienzle
 607 % Use frivolous uint8() call to block tests for version 2.1.57 and earlier
 608 %
 609 % Revision 1.5  2004/09/15 20:36:57  jmones
 610 % logical(1) => true
 611 %
 612 % Revision 1.4  2004/09/15 20:00:00  jmones
 613 % Updated tests to match Gonzalez&Woods example
 614 %
 615 % Revision 1.3  2004/09/15 13:51:10  pkienzle
 616 % Use logical in tests; reduce number of shared variables in tests.
 617 %
 618 % Revision 1.2  2004/09/01 22:35:47  jmones
 619 % Added Lantuejoul skeletonizing algorithm from Gonzalez&Woods
 620 %
 621 % Revision 1.1  2004/08/15 19:47:04  jmones
 622 % bwmorph added: Perform a morphological operation on a binary image
 623 %
 624 %