1 ## Copyright (C) 2008 Søren Hauberg
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.
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.
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
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.
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.
35 ## The @var{method} input argument can any of the following strings (the default
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).
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.
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.
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).
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
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
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
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
112 ## A.Adler's idea (c) 1999. Somewhat based on the canny method. The steps are
115 ## Do a Sobel edge detection and to generate an image at
116 ## a high and low threshold.
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.
122 ## Dilate the EE image by 1 step.
124 ## Select all EE features that are connected to features in
128 ## The parameters for the method is given in a vector:
130 ## @item params(1)==0 or 4 or 8
131 ## Perform x connected dilatation (step 3).
133 ## Dilatation coeficient (threshold) in step 3.
135 ## Length of edge extention convolution (step 2).
137 ## Coeficient of extention convolution in step 2.
139 ## defaults = [8, 1, 3, 3]
143 ## @seealso{fspecial, nonmax_supress}
146 function [bw, out_threshold, g45_out, g135_out] = edge (im, method, varargin)
149 error("edge: not enough input arguments");
152 error("edge: first input must be a gray-scale image");
160 error("edge: second argument must be a string");
162 method = lower(method);
164 ## Perform the actual edge detection
166 #####################################
168 #####################################
170 ## Get the direction argument
171 direction = get_direction(varargin{:});
173 h1 = fspecial("sobel"); # horizontal
175 ## Compute edge strength
178 strength = abs( conv2(im, h1, "same") );
180 strength = abs( conv2(im, h3, "same") );
182 strength = sqrt( conv2(im, h1, "same").^2 + ...
183 conv2(im, h3, "same").^2 );
186 if (nargin > 2 && isscalar(varargin{1}))
187 thresh = varargin{1};
189 thresh = 2*mean(strength(:));
191 ## Perform thresholding and simple thinning
192 strength(strength<=thresh) = 0;
193 bw = simple_thinning(strength);
195 #####################################
197 #####################################
199 ## Get the direction argument
200 direction = get_direction(varargin{:});
202 h1 = fspecial("prewitt"); # vertical
203 h3 = h1'; # horizontal
204 ## Compute edge strength
207 strength = abs( conv2(im, h1, "same") );
209 strength = abs( conv2(im, h3, "same") );
211 strength = sqrt( conv2(im, h1, "same").^2 + ...
212 conv2(im, h3, "same").^2 );
215 if (nargin > 2 && isscalar(varargin{1}))
216 thresh = varargin{1};
218 thresh = 4*mean(strength(:));
220 ## Perform thresholding and simple thinning
221 strength(strength<=thresh) = 0;
222 bw = simple_thinning(strength);
224 #####################################
226 #####################################
228 ## Get the direction argument
229 direction = get_direction(varargin{:});
231 h1 = fspecial("kirsch"); # vertical
232 h3 = h1'; # horizontal
233 ## Compute edge strength
236 strength = abs( conv2(im, h1, "same") );
238 strength = abs( conv2(im, h3, "same") );
240 strength = sqrt( conv2(im, h1, "same").^2 + ...
241 conv2(im, h3, "same").^2 );
244 if nargin > 2 && isscalar(varargin{1})
245 thresh = varargin{1};
247 thresh = mean(strength(:));
249 ## Perform thresholding and simple thinning
250 strength(strength<=thresh) = 0;
251 bw = simple_thinning(strength);
253 #####################################
255 #####################################
257 ## Get the thinning argument (option)
259 option = varargin{2};
261 error("edge: 'option' must be a string");
263 option = lower(option);
264 if (!any(strcmp(option, {"thinning", "nothinning"})))
265 error("edge: 'option' must be either 'thinning', or 'nothinning'");
273 ## Compute edge strength
274 g45 = conv2(im, h1, "same");
275 g135 = conv2(im, h2, "same");
276 strength = abs( g45 ) + abs( g135 );
278 if (nargin > 2 && isscalar(varargin{1}))
279 thresh = varargin{1};
281 thresh = 6*mean(strength(:));
283 ## Perform thresholding and simple thinning
284 strength(strength<=thresh) = 0;
285 if (strcmp(option, "thinning"))
286 bw = simple_thinning(strength);
290 ## Check if g45 and g135 should be returned
296 #####################################
297 ## L A P L A C I A N O F G A U S S I A N
298 #####################################
301 if (nargin == 4 && isscalar(varargin{2}))
308 %[x y] = meshgrid(-s:s);
309 %f = (x.^2 + y.^2 - sigma^2) .* exp(-(x.^2 + y.^2)/(2*sigma^2));
311 f = fspecial("log", 2*s+1, sigma);
312 ## Perform convolution with the filter f
313 g = conv2(im, f, "same");
315 if (nargin > 2 && isscalar(varargin{1}))
316 thresh = varargin{1};
318 thresh = 0.75*mean(abs(g(:)));
320 ## Find zero crossings
321 zc = zerocrossings(g);
322 bw = (abs(g) >= thresh) & zc;
324 #####################################
325 ## Z E R O C R O S S I N G
326 #####################################
329 if (nargin == 4 && ismatrix(varargin{2}))
332 error("edge: a filter must be given as the fourth argument when 'zerocross' is used");
334 ## Perform convolution with the filter f
335 g = conv2(im, f, "same");
337 if (nargin > 2 && isscalar(varargin{1}))
338 thresh = varargin{1};
340 thresh = mean(abs(g(:)));
342 ## Find zero crossings
343 zc = zerocrossings(g);
344 bw = (abs(g) >= thresh) & zc;
346 #####################################
348 #####################################
351 if (nargin == 4 && isscalar(varargin{2}))
358 J = imsmooth(double(im), "Gaussian", sigma);
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);
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}(:);
373 tmp = mean(abs(Es(:)));
374 thresh = [0.4*tmp, tmp];
376 bw = nonmax_supress(Es, Eo, thresh(1), thresh(2));
378 #####################################
380 #####################################
382 ## In case the user asks for more then 1 output argument
383 ## we define thresh to be -1.
386 if (nargin > 2 && isscalar(varargin{1}))
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");
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;
420 #####################################
422 #####################################
424 [bw, out_threshold] = andy (im, method, varargin{:});
427 error("edge: unsupported edge detector: %s", method);
431 out_threshold = thresh;
435 ## An auxilary function that parses the 'direction' argument from 'varargin'
436 function direction = get_direction(varargin)
438 direction = varargin{2};
439 if (!ischar(direction))
440 error("edge: direction must be a string");
442 direction = lower(direction);
443 if (!any(strcmp(direction, {"horizontal", "vertical", "both"})))
444 error("edge :direction must be either 'horizontal', 'vertical', or 'both'");
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) ] );
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
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);
479 z &= !z0; ## "Positive zero-crossings"?
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)
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) )) );
495 # sum( imo(:)>thresh ) / prod(size(imo))
496 dilate= [1 1 1;1 1 1;1 1 1]; tt= 1; sz=3; dt=3;
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;
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
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) );
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) ;
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;
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;
528 # now we edge extend the low thresh image in 4 directions
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; #
541 % ff= find( imht==1 );
542 % imout = bwselect( imee, rem(ff-1, n)+1, ceil(ff/n), 8);