octave_packages/image-1.0.15/edge.m

   1 ## Copyright (C) 2008 Søren Hauberg
   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 3 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, write to the Free Software
  15 ## Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  16 ##
  17 ## Note: The implementation and help text for the 'andy' edge detector came with
  18 ## the following notice:
  19 ##   Copyright (C) 1999 Andy Adler
  20 ##   This code has no warrany whatsoever.
  21 ##   Do what you like with this code as long as you
  22 ##   leave this copyright in place.
  23
  24 ## -*- texinfo -*-
  25 ## @deftypefn {Function File} {@var{bw} =} edge (@var{im}, @var{method})
  26 ## @deftypefnx{Function File} {@var{bw} =} edge (@var{im}, @var{method}, @var{arg1}, @var{arg2})
  27 ## @deftypefnx{Function File} {[@var{bw}, @var{thresh}] =} edge (...)
  28 ## Detect edges in the given image using various methods. The first input @var{im}
  29 ## is the gray scale image in which edges are to be detected. The second argument
  30 ## controls which method is used for detecting the edges. The rest of the input
  31 ## arguments depend on the selected method. The first output @var{bw} is a
  32 ## @code{logical} image containing the edges. Most methods also returns an automatically
  33 ## computed threshold as the second output.
  34 ##
  35 ## The @var{method} input argument can any of the following strings (the default
  36 ## value is "Sobel")
  37 ##
  38 ## @table @asis
  39 ## @item "Sobel"
  40 ## Finds the edges in @var{im} using the Sobel approximation to the
  41 ## derivatives. Edge points are defined as points where the length of
  42 ## the gradient exceeds a threshold and is larger than it's neighbours
  43 ## in either the horizontal or vertical direction. The threshold is passed to
  44 ## the method in the third input argument @var{arg1}. If one is not given, a
  45 ## threshold is automatically computed as 4*@math{M}, where @math{M} is the mean
  46 ## of the gradient of the entire image. The optional 4th input argument controls
  47 ## the direction in which the gradient is approximated. It can be either
  48 ## "horizontal", "vertical", or "both" (default).
  49 ##
  50 ## @item "Prewitt"
  51 ## Finds the edges in @var{im} using the Prewitt approximation to the
  52 ## derivatives. This method works just like "Sobel" except a different aproximation
  53 ## the gradient is used.
  54 ##
  55 ## @item "Roberts"
  56 ## Finds the edges in @var{im} using the Roberts approximation to the
  57 ## derivatives. Edge points are defined as points where the length of
  58 ## the gradient exceeds a threshold and is larger than it's neighbours
  59 ## in either the horizontal or vertical direction. The threshold is passed to
  60 ## the method in the third input argument @var{arg1}. If one is not given, a
  61 ## threshold is automatically computed as 6*@math{M}, where @math{M} is the mean
  62 ## of the gradient of the entire image. The optional 4th input argument can be
  63 ## either "thinning" (default) or "nothinning". If it is "thinning" a simple
  64 ## thinning procedure is applied to the edge image such that the edges are only
  65 ## one pixel wide. If @var{arg2} is "nothinning", this procedure is not applied.
  66 ##
  67 ## @item "Kirsch"
  68 ## Finds the edges in @var{im} using the Kirsch approximation to the
  69 ## derivatives. Edge points are defined as points where the length of
  70 ## the gradient exceeds a threshold and is larger than it's neighbours
  71 ## in either the horizontal or vertical direction. The threshold is passed to
  72 ## the method in the third input argument @var{arg1}. If one is not given, a
  73 ## threshold is automatically computed as @math{M}, where @math{M} is the mean
  74 ## of the gradient of the entire image. The optional 4th input argument controls
  75 ## the direction in which the gradient is approximated. It can be either
  76 ## "horizontal", "vertical", or "both" (default).
  77 ##
  78 ## @item "LoG"
  79 ## Finds edges in @var{im} by convolving with the Laplacian of Gaussian (LoG)
  80 ## filter, and finding zero crossings. Only zero crossings where the
  81 ## filter response is larger than an automatically computed threshold are retained.
  82 ## The threshold is passed to the method in the third input argument @var{arg1}.
  83 ## If one is not given, a threshold is automatically computed as 0.75*@math{M},
  84 ## where @math{M} is the mean of absolute value of LoG filter response. The
  85 ## optional 4th input argument sets the spread of the LoG filter. By default
  86 ## this value is 2.
  87 ##
  88 ## @item "Zerocross"
  89 ## Finds edges in the image @var{im} by convolving it with the user-supplied filter
  90 ## @var{arg2} and finding zero crossings larger than the threshold @var{arg1}. If
  91 ## @var{arg1} is [] a threshold is computed as the mean value of the absolute
  92 ## filter response.
  93 ##
  94 ## @item "Canny"
  95 ## Finds edges using the Canny edge detector. The optional third input argument
  96 ## @var{arg1} sets the thresholds used in the hysteresis thresholding. If
  97 ## @var{arg1} is a two dimensional vector it's first element is used as the lower
  98 ## threshold, while the second element is used as the high threshold. If, on the
  99 ## other hand, @var{arg1} is a single scalar it is used as the high threshold,
 100 ## while the lower threshold is 0.4*@var{arg1}. The optional 4th input argument
 101 ## @var{arg2} is the spread of the low-pass Gaussian filter that is used to smooth
 102 ## the input image prior to estimating gradients. By default this scale parameter
 103 ## is 2.
 104 ##
 105 ## @item "Lindeberg"
 106 ## Finds edges using in @var{im} using the differential geometric single-scale edge
 107 ## detector given by Tony Lindeberg. The optional third input argument @var{arg1}
 108 ## is the scale (spread of Gaussian filter) at which the edges are computed. By
 109 ## default this 2.
 110 ##
 111 ## @item "Andy"
 112 ## A.Adler's idea (c) 1999. Somewhat based on the canny method. The steps are
 113 ## @enumerate
 114 ## @item
 115 ## Do a Sobel edge detection and to generate an image at
 116 ## a high and low threshold.
 117 ## @item
 118 ## Edge extend all edges in the LT image by several pixels,
 119 ## in the vertical, horizontal, and 45 degree directions.
 120 ## Combine these into edge extended (EE) image.
 121 ## @item
 122 ## Dilate the EE image by 1 step.
 123 ## @item
 124 ## Select all EE features that are connected to features in
 125 ## the HT image.
 126 ## @end enumerate
 127 ##
 128 ## The parameters for the method is given in a vector:
 129 ## @table @asis
 130 ## @item params(1)==0 or 4 or 8
 131 ## Perform x connected dilatation (step 3).
 132 ## @item params(2)
 133 ## Dilatation coeficient (threshold) in step 3.
 134 ## @item params(3)
 135 ## Length of edge extention convolution (step 2).
 136 ## @item params(4)
 137 ## Coeficient of extention convolution in step 2.
 138 ## @end table
 139 ## defaults = [8, 1, 3, 3]
 140 ##
 141 ## @end table
 142 ##
 143 ## @seealso{fspecial, nonmax_supress}
 144 ## @end deftypefn
 145
 146 function [bw, out_threshold, g45_out, g135_out] = edge (im, method, varargin)
 147   ## Get the image
 148   if (nargin == 0)
 149     error("edge: not enough input arguments");
 150   endif
 151   if ( !isgray(im) )
 152     error("edge: first input must be a gray-scale image");
 153   endif
 154
 155   ## Get the method
 156   if (nargin == 1)
 157     method = "Sobel";
 158   endif
 159   if (!ischar(method))
 160     error("edge: second argument must be a string");
 161   endif
 162   method = lower(method);
 163
 164   ## Perform the actual edge detection
 165   switch (method)
 166     #####################################
 167     ## S O B E L
 168     #####################################
 169     case "sobel"
 170       ## Get the direction argument
 171       direction = get_direction(varargin{:});
 172       ## Create filters;
 173       h1 = fspecial("sobel"); # horizontal
 174       h3 = h1'; # vertical
 175       ## Compute edge strength
 176       switch(direction)
 177         case "horizontal"
 178           strength = abs( conv2(im, h1, "same") );
 179         case "vertical"
 180           strength = abs( conv2(im, h3, "same") );
 181         case "both"
 182           strength = sqrt( conv2(im, h1, "same").^2 + ...
 183                            conv2(im, h3, "same").^2 );
 184       endswitch
 185       ## Get threshold
 186       if (nargin > 2 && isscalar(varargin{1}))
 187         thresh = varargin{1};
 188       else
 189         thresh = 2*mean(strength(:));
 190       endif
 191       ## Perform thresholding and simple thinning
 192       strength(strength<=thresh) = 0;
 193       bw = simple_thinning(strength);
 194
 195     #####################################
 196     ## P R E W I T T
 197     #####################################
 198     case "prewitt"
 199       ## Get the direction argument
 200       direction = get_direction(varargin{:});
 201       ## Create filters;
 202       h1 = fspecial("prewitt"); # vertical
 203       h3 = h1'; # horizontal
 204       ## Compute edge strength
 205       switch(direction)
 206         case "vertical"
 207           strength = abs( conv2(im, h1, "same") );
 208         case "horizontal"
 209           strength = abs( conv2(im, h3, "same") );
 210         case "both"
 211           strength = sqrt( conv2(im, h1, "same").^2 + ...
 212                            conv2(im, h3, "same").^2 );
 213       endswitch
 214       ## Get threshold
 215       if (nargin > 2 && isscalar(varargin{1}))
 216         thresh = varargin{1};
 217       else
 218         thresh = 4*mean(strength(:));
 219       endif
 220       ## Perform thresholding and simple thinning
 221       strength(strength<=thresh) = 0;
 222       bw = simple_thinning(strength);
 223
 224     #####################################
 225     ## K I R S C H
 226     #####################################
 227     case "kirsch"
 228       ## Get the direction argument
 229       direction = get_direction(varargin{:});
 230       ## Create filters;
 231       h1 = fspecial("kirsch"); # vertical
 232       h3 = h1'; # horizontal
 233       ## Compute edge strength
 234       switch(direction)
 235         case "vertical"
 236           strength = abs( conv2(im, h1, "same") );
 237         case "horizontal"
 238           strength = abs( conv2(im, h3, "same") );
 239         case "both"
 240           strength = sqrt( conv2(im, h1, "same").^2 + ...
 241                            conv2(im, h3, "same").^2 );
 242       endswitch
 243       ## Get threshold
 244       if nargin > 2 && isscalar(varargin{1})
 245         thresh = varargin{1};
 246       else
 247         thresh = mean(strength(:));
 248       endif
 249       ## Perform thresholding and simple thinning
 250       strength(strength<=thresh) = 0;
 251       bw = simple_thinning(strength);
 252
 253     #####################################
 254     ## R O B E R T S
 255     #####################################
 256     case "roberts"
 257       ## Get the thinning argument (option)
 258       if (nargin == 4)
 259         option = varargin{2};
 260         if (!ischar(option))
 261           error("edge: 'option' must be a string");
 262         endif
 263         option = lower(option);
 264         if (!any(strcmp(option, {"thinning", "nothinning"})))
 265           error("edge: 'option' must be either 'thinning', or 'nothinning'");
 266         endif
 267       else
 268         option = "thinning";
 269       endif
 270       ## Create filters;
 271       h1 = [1 0; 0 -1];
 272       h2 = [0 1; -1 0];
 273       ## Compute edge strength
 274       g45  = conv2(im, h1, "same");
 275       g135 = conv2(im, h2, "same");
 276       strength = abs( g45 ) + abs( g135 );
 277       ## Get threshold
 278       if (nargin > 2 && isscalar(varargin{1}))
 279         thresh = varargin{1};
 280       else
 281         thresh = 6*mean(strength(:));
 282       endif
 283       ## Perform thresholding and simple thinning
 284       strength(strength<=thresh) = 0;
 285       if (strcmp(option, "thinning"))
 286         bw = simple_thinning(strength);
 287       else
 288         bw = (strength > 0);
 289       endif
 290       ## Check if g45 and g135 should be returned
 291       if (nargout == 4)
 292         g45_out  = g45;
 293         g135_out = g135;
 294       endif
 295
 296     #####################################
 297     ## L A P L A C I A N   O F   G A U S S I A N
 298     #####################################
 299     case "log"
 300       ## Get sigma
 301       if (nargin == 4 && isscalar(varargin{2}))
 302         sigma = varargin{2};
 303       else
 304         sigma = 2;
 305       endif
 306       ## Create the filter
 307       s = ceil(3*sigma);
 308       %[x y] = meshgrid(-s:s);
 309       %f = (x.^2 + y.^2 - sigma^2) .* exp(-(x.^2 + y.^2)/(2*sigma^2));
 310       %f = f/sum(f(:));
 311       f = fspecial("log", 2*s+1, sigma);
 312       ## Perform convolution with the filter f
 313       g = conv2(im, f, "same");
 314       ## Get threshold
 315       if (nargin > 2 && isscalar(varargin{1}))
 316         thresh = varargin{1};
 317       else
 318         thresh = 0.75*mean(abs(g(:)));
 319       endif
 320       ## Find zero crossings
 321       zc = zerocrossings(g);
 322       bw = (abs(g) >= thresh) & zc;
 323
 324     #####################################
 325     ## Z E R O   C R O S S I N G
 326     #####################################
 327     case "zerocross"
 328       ## Get the filter
 329       if (nargin == 4 && ismatrix(varargin{2}))
 330         f = varargin{2};
 331       else
 332         error("edge: a filter must be given as the fourth argument when 'zerocross' is used");
 333       endif
 334       ## Perform convolution with the filter f
 335       g = conv2(im, f, "same");
 336       ## Get threshold
 337       if (nargin > 2 && isscalar(varargin{1}))
 338         thresh = varargin{1};
 339       else
 340         thresh = mean(abs(g(:)));
 341       endif
 342       ## Find zero crossings
 343       zc = zerocrossings(g);
 344       bw = (abs(g) >= thresh) & zc;
 345
 346     #####################################
 347     ## C A N N Y
 348     #####################################
 349     case "canny"
 350       ## Get sigma
 351       if (nargin == 4 && isscalar(varargin{2}))
 352         sigma = varargin{2};
 353       else
 354         sigma = 2;
 355       endif
 356
 357       ## Change scale
 358       J = imsmooth(double(im), "Gaussian", sigma);
 359
 360       ## Canny enhancer
 361       p = [1 0 -1]/2;
 362       Jx = conv2(J, p,  "same");
 363       Jy = conv2(J, p', "same");
 364       Es = sqrt( Jx.^2 + Jy.^2 );
 365       Eo = pi - mod (atan2 (Jy, Jx) - pi, pi);
 366
 367       ## Get thresholds
 368       if (nargin > 2 && isscalar(varargin{1}))
 369         thresh = [0.4*varargin{1}, varargin{1}];
 370       elseif (nargin > 2 && ismatrix (varargin{1}) && length (varargin{1}(:)) == 2)
 371         thresh = varargin{1}(:);
 372       else
 373         tmp = mean(abs(Es(:)));
 374         thresh = [0.4*tmp, tmp];
 375       endif
 376       bw = nonmax_supress(Es, Eo, thresh(1), thresh(2));
 377
 378     #####################################
 379     ## L I N D E B E R G
 380     #####################################
 381     case "lindeberg"
 382       ## In case the user asks for more then 1 output argument
 383       ## we define thresh to be -1.
 384       thresh = -1;
 385       ## Get sigma
 386       if (nargin > 2 && isscalar(varargin{1}))
 387         sigma = varargin{1};
 388       else
 389         sigma = 2;
 390       endif
 391       ## Filters for computing the derivatives
 392       Px   = [-1 0 1; -1 0 1; -1 0 1];
 393       Py   = [1 1 1; 0 0 0; -1 -1 -1];
 394       Pxx  = conv2(Px,  Px, "full");
 395       Pyy  = conv2(Py,  Py, "full");
 396       Pxy  = conv2(Px,  Py, "full");
 397       Pxxx = conv2(Pxx, Px, "full");
 398       Pyyy = conv2(Pyy, Py, "full");
 399       Pxxy = conv2(Pxx, Py, "full");
 400       Pxyy = conv2(Pyy, Px, "full");
 401       ## Change scale
 402       L = imsmooth(double(im), "Gaussian", sigma);
 403       ## Compute derivatives
 404       Lx   = conv2(L, Px,   "same");
 405       Ly   = conv2(L, Py,   "same");
 406       Lxx  = conv2(L, Pxx,  "same");
 407       Lyy  = conv2(L, Pyy,  "same");
 408       Lxy  = conv2(L, Pxy,  "same");
 409       Lxxx = conv2(L, Pxxx, "same");
 410       Lyyy = conv2(L, Pyyy, "same");
 411       Lxxy = conv2(L, Pxxy, "same");
 412       Lxyy = conv2(L, Pxyy, "same");
 413       ## Compute directional derivatives
 414       Lvv  = Lx.^2.*Lxx + 2.*Lx.*Ly.*Lxy + Ly.^2.*Lyy;
 415       Lvvv = Lx.^3.*Lxxx + 3.*Lx.^2.*Ly.*Lxxy ...
 416            + 3.*Lx.*Ly.^2.*Lxyy + 3.*Ly.^3.*Lyyy;
 417       ## Perform edge detection
 418       bw = zerocrossings(Lvv) & Lvvv < 0;
 419
 420     #####################################
 421     ## A N D Y
 422     #####################################
 423     case "andy"
 424       [bw, out_threshold] = andy (im, method, varargin{:});
 425
 426     otherwise
 427       error("edge: unsupported edge detector: %s", method);
 428   endswitch
 429
 430   if (nargout > 1)
 431     out_threshold = thresh;
 432   endif
 433 endfunction
 434
 435 ## An auxilary function that parses the 'direction' argument from 'varargin'
 436 function direction = get_direction(varargin)
 437   if (nargin >= 2)
 438     direction = varargin{2};
 439     if (!ischar(direction))
 440       error("edge: direction must be a string");
 441     endif
 442     direction = lower(direction);
 443     if (!any(strcmp(direction, {"horizontal", "vertical", "both"})))
 444       error("edge :direction must be either 'horizontal', 'vertical', or 'both'");
 445     endif
 446   else
 447     direction = "both";
 448   endif
 449 endfunction
 450
 451 ## An auxilary function that performs a very simple thinning.
 452 ## Strength is an image containing the edge strength.
 453 ## bw contains a 1 in (r,c) if
 454 ##  1) strength(r,c) is greater than both neighbours in the
 455 ##     vertical direction, OR
 456 ##  2) strength(r,c) is greater than both neighbours in the
 457 ##     horizontal direction.
 458 ## Note the use of OR.
 459 function bw = simple_thinning(strength)
 460   [r c] = size(strength);
 461   x = ( strength > [ zeros(r,1) strength(:,1:end-1) ] & ...
 462         strength > [ strength(:,2:end) zeros(r,1) ] );
 463   y = ( strength > [ zeros(1,c); strength(1:end-1,:) ] & ...
 464         strength > [ strength(2:end,:); zeros(1,c) ] );
 465   bw = x | y;
 466 endfunction
 467
 468 ## Auxilary function. Finds the zero crossings of the
 469 ## 2-dimensional function f. (By Etienne Grossmann)
 470 function z = zerocrossings(f)
 471   z0 = f<0;                 ## Negative
 472   [R,C] = size(f);
 473   z = zeros(R,C);
 474   z(1:R-1,:) |= z0(2:R,:);  ## Grow
 475   z(2:R,:) |= z0(1:R-1,:);
 476   z(:,1:C-1) |= z0(:,2:C);
 477   z(:,2:C) |= z0(:,1:C-1);
 478
 479   z &= !z0;                  ## "Positive zero-crossings"?
 480 endfunction
 481
 482 ## The 'andy' edge detector that was present in older versions of 'edge'.
 483 ## The function body has simply been copied from the old implementation.
 484 ##   -- Soren Hauberg, march 11th, 2008
 485 function [imout, thresh] = andy(im, method, thresh, param2)
 486    [n,m]= size(im);
 487    xx= 2:m-1;
 488    yy= 2:n-1;
 489
 490    filt= [1 2 1;0 0 0; -1 -2 -1]/8;  tv= 2;
 491    imo= conv2(im, rot90(filt), 'same').^2 + conv2(im, filt, 'same').^2;
 492    if nargin<3 || thresh==[];
 493       thresh= sqrt( tv* mean(mean( imo(yy,xx) ))  );
 494    end
 495 #     sum( imo(:)>thresh ) / prod(size(imo))
 496    dilate= [1 1 1;1 1 1;1 1 1]; tt= 1; sz=3; dt=3;
 497    if nargin>=4
 498       # 0 or 4 or 8 connected dilation
 499       if length(param2) > 0
 500          if      param2(1)==4 ; dilate= [0 1 0;1 1 1;0 1 0];
 501          elseif  param2(1)==0 ; dilate= 1;
 502          end
 503       end
 504       # dilation threshold
 505       if length(param2) > 2; tt= param2(2); end
 506       # edge extention length
 507       if length(param2) > 2; sz= param2(3); end
 508       # edge extention threshold
 509       if length(param2) > 3; dt= param2(4); end
 510
 511    end
 512    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;
 513                       0 .5 .5 0 0;0 1 0 0 0;.5 .5 0 0 0;1 0 0 0 0];
 514    fobliq= fobliq( 5-sz:5+sz, 3-ceil(sz/2):3+ceil(sz/2) );
 515
 516    xpeak= imo(yy,xx-1) <= imo(yy,xx) & imo(yy,xx) > imo(yy,xx+1) ;
 517    ypeak= imo(yy-1,xx) <= imo(yy,xx) & imo(yy,xx) > imo(yy+1,xx) ;
 518
 519    imht= ( imo >= thresh^2 * 2); # high threshold image
 520    imht(yy,xx)= imht(yy,xx) & ( xpeak | ypeak );
 521    imht([1,n],:)=0; imht(:,[1,m])=0;
 522
 523 %  imlt= ( imo >= thresh^2 / 2); # low threshold image
 524    imlt= ( imo >= thresh^2 / 1); # low threshold image
 525    imlt(yy,xx)= imlt(yy,xx) & ( xpeak | ypeak );
 526    imlt([1,n],:)=0; imlt(:,[1,m])=0;
 527
 528 # now we edge extend the low thresh image in 4 directions
 529
 530    imee= ( conv2( imlt, ones(2*sz+1,1)    , 'same') > tt ) | ...
 531          ( conv2( imlt, ones(1,2*sz+1)    , 'same') > tt ) | ...
 532          ( conv2( imlt, eye(2*sz+1)       , 'same') > tt ) | ...
 533          ( conv2( imlt, rot90(eye(2*sz+1)), 'same') > tt ) | ...
 534          ( conv2( imlt, fobliq            , 'same') > tt ) | ...
 535          ( conv2( imlt, fobliq'           , 'same') > tt ) | ...
 536          ( conv2( imlt, rot90(fobliq)     , 'same') > tt ) | ...
 537          ( conv2( imlt, flipud(fobliq)    , 'same') > tt );
 538 #  imee(yy,xx)= conv2(imee(yy,xx),ones(3),'same') & ( xpeak | ypeak );
 539    imee= conv2(imee,dilate,'same') > dt; #
 540
 541 %  ff= find( imht==1 );
 542 %  imout = bwselect( imee, rem(ff-1, n)+1, ceil(ff/n), 8);
 543    imout = imee;
 544
 545 endfunction