/branches/R2-1-x/octave-forge/main/image/edge.m
MATLAB | 240 lines | 220 code | 19 blank | 1 comment | 16 complexity | b9e8853110860bf790bb81d29fd04f23 MD5 | raw file
Possible License(s): GPL-2.0, BSD-3-Clause, LGPL-2.1, GPL-3.0, LGPL-3.0
- 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: edge.m 187 2002-03-17 02:38:52Z aadler $
- [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
- #
- #
- #