X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?p=CreaPhase.git;a=blobdiff_plain;f=octave_packages%2Fimage-1.0.15%2Fedge.m;fp=octave_packages%2Fimage-1.0.15%2Fedge.m;h=9ba0213372fdd20ee3c4b4387769b214263ee71b;hp=0000000000000000000000000000000000000000;hb=f5f7a74bd8a4900f0b797da6783be80e11a68d86;hpb=1705066eceaaea976f010f669ce8e972f3734b05

diff --git a/octave_packages/image-1.0.15/edge.m b/octave_packages/image-1.0.15/edge.m
new file mode 100644
index 0000000..9ba0213
--- /dev/null
+++ b/octave_packages/image-1.0.15/edge.m
@@ -0,0 +1,545 @@
+## Copyright (C) 2008 Søren Hauberg
+## 
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or
+## (at your option) any later version.
+## 
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+## 
+## You should have received a copy of the GNU General Public License
+## along with this program; if not, write to the Free Software
+## Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+##
+## Note: The implementation and help text for the 'andy' edge detector came with
+## the following notice:
+##   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.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{bw} =} edge (@var{im}, @var{method})
+## @deftypefnx{Function File} {@var{bw} =} edge (@var{im}, @var{method}, @var{arg1}, @var{arg2})
+## @deftypefnx{Function File} {[@var{bw}, @var{thresh}] =} edge (...)
+## Detect edges in the given image using various methods. The first input @var{im}
+## is the gray scale image in which edges are to be detected. The second argument
+## controls which method is used for detecting the edges. The rest of the input
+## arguments depend on the selected method. The first output @var{bw} is a 
+## @code{logical} image containing the edges. Most methods also returns an automatically
+## computed threshold as the second output.
+##
+## The @var{method} input argument can any of the following strings (the default
+## value is "Sobel")
+##
+## @table @asis
+## @item "Sobel"
+## Finds the edges in @var{im} using the Sobel approximation to the
+## derivatives. Edge points are defined as points where the length of
+## the gradient exceeds a threshold and is larger than it's neighbours
+## in either the horizontal or vertical direction. The threshold is passed to
+## the method in the third input argument @var{arg1}. If one is not given, a
+## threshold is automatically computed as 4*@math{M}, where @math{M} is the mean
+## of the gradient of the entire image. The optional 4th input argument controls
+## the direction in which the gradient is approximated. It can be either
+## "horizontal", "vertical", or "both" (default).
+##
+## @item "Prewitt"
+## Finds the edges in @var{im} using the Prewitt approximation to the
+## derivatives. This method works just like "Sobel" except a different aproximation
+## the gradient is used.
+##
+## @item "Roberts"
+## Finds the edges in @var{im} using the Roberts approximation to the
+## derivatives. Edge points are defined as points where the length of
+## the gradient exceeds a threshold and is larger than it's neighbours
+## in either the horizontal or vertical direction. The threshold is passed to
+## the method in the third input argument @var{arg1}. If one is not given, a
+## threshold is automatically computed as 6*@math{M}, where @math{M} is the mean
+## of the gradient of the entire image. The optional 4th input argument can be
+## either "thinning" (default) or "nothinning". If it is "thinning" a simple
+## thinning procedure is applied to the edge image such that the edges are only
+## one pixel wide. If @var{arg2} is "nothinning", this procedure is not applied.
+##
+## @item "Kirsch"
+## Finds the edges in @var{im} using the Kirsch approximation to the
+## derivatives. Edge points are defined as points where the length of
+## the gradient exceeds a threshold and is larger than it's neighbours
+## in either the horizontal or vertical direction. The threshold is passed to
+## the method in the third input argument @var{arg1}. If one is not given, a
+## threshold is automatically computed as @math{M}, where @math{M} is the mean
+## of the gradient of the entire image. The optional 4th input argument controls
+## the direction in which the gradient is approximated. It can be either
+## "horizontal", "vertical", or "both" (default).
+##
+## @item "LoG"
+## Finds edges in @var{im} by convolving with the Laplacian of Gaussian (LoG)
+## filter, and finding zero crossings. Only zero crossings where the 
+## filter response is larger than an automatically computed threshold are retained.
+## The threshold is passed to the method in the third input argument @var{arg1}.
+## If one is not given, a threshold is automatically computed as 0.75*@math{M},
+## where @math{M} is the mean of absolute value of LoG filter response. The
+## optional 4th input argument sets the spread of the LoG filter. By default
+## this value is 2.
+##
+## @item "Zerocross"
+## Finds edges in the image @var{im} by convolving it with the user-supplied filter
+## @var{arg2} and finding zero crossings larger than the threshold @var{arg1}. If
+## @var{arg1} is [] a threshold is computed as the mean value of the absolute
+## filter response.
+##
+## @item "Canny"
+## Finds edges using the Canny edge detector. The optional third input argument
+## @var{arg1} sets the thresholds used in the hysteresis thresholding. If 
+## @var{arg1} is a two dimensional vector it's first element is used as the lower
+## threshold, while the second element is used as the high threshold. If, on the
+## other hand, @var{arg1} is a single scalar it is used as the high threshold,
+## while the lower threshold is 0.4*@var{arg1}. The optional 4th input argument
+## @var{arg2} is the spread of the low-pass Gaussian filter that is used to smooth
+## the input image prior to estimating gradients. By default this scale parameter
+## is 2.
+##
+## @item "Lindeberg"
+## Finds edges using in @var{im} using the differential geometric single-scale edge
+## detector given by Tony Lindeberg. The optional third input argument @var{arg1}
+## is the scale (spread of Gaussian filter) at which the edges are computed. By
+## default this 2.
+##
+## @item "Andy"
+## A.Adler's idea (c) 1999. Somewhat based on the canny method. The steps are
+## @enumerate
+## @item
+## Do a Sobel edge detection and to generate an image at
+## a high and low threshold.
+## @item
+## Edge extend all edges in the LT image by several pixels,
+## in the vertical, horizontal, and 45 degree directions.
+## Combine these into edge extended (EE) image.
+## @item
+## Dilate the EE image by 1 step.
+## @item
+## Select all EE features that are connected to features in
+## the HT image.
+## @end enumerate
+## 
+## The parameters for the method is given in a vector:
+## @table @asis
+## @item params(1)==0 or 4 or 8
+## Perform x connected dilatation (step 3).
+## @item params(2)
+## Dilatation coeficient (threshold) in step 3.
+## @item params(3)
+## Length of edge extention convolution (step 2).
+## @item params(4)
+## Coeficient of extention convolution in step 2.
+## @end table
+## defaults = [8, 1, 3, 3]
+##
+## @end table
+##
+## @seealso{fspecial, nonmax_supress}
+## @end deftypefn
+
+function [bw, out_threshold, g45_out, g135_out] = edge (im, method, varargin)
+  ## Get the image
+  if (nargin == 0)
+    error("edge: not enough input arguments");
+  endif
+  if ( !isgray(im) )
+    error("edge: first input must be a gray-scale image");
+  endif
+
+  ## Get the method
+  if (nargin == 1)
+    method = "Sobel";
+  endif
+  if (!ischar(method))
+    error("edge: second argument must be a string");
+  endif
+  method = lower(method);
+
+  ## Perform the actual edge detection
+  switch (method)
+    #####################################
+    ## S O B E L
+    #####################################
+    case "sobel"
+      ## Get the direction argument
+      direction = get_direction(varargin{:});
+      ## Create filters;
+      h1 = fspecial("sobel"); # horizontal
+      h3 = h1'; # vertical
+      ## Compute edge strength
+      switch(direction)
+        case "horizontal"
+          strength = abs( conv2(im, h1, "same") );
+        case "vertical"
+          strength = abs( conv2(im, h3, "same") );
+        case "both"
+          strength = sqrt( conv2(im, h1, "same").^2 + ...
+                           conv2(im, h3, "same").^2 );
+      endswitch
+      ## Get threshold
+      if (nargin > 2 && isscalar(varargin{1}))
+        thresh = varargin{1};
+      else
+        thresh = 2*mean(strength(:));
+      endif
+      ## Perform thresholding and simple thinning
+      strength(strength<=thresh) = 0;
+      bw = simple_thinning(strength);
+
+    #####################################
+    ## P R E W I T T
+    #####################################
+    case "prewitt"
+      ## Get the direction argument
+      direction = get_direction(varargin{:});
+      ## Create filters;
+      h1 = fspecial("prewitt"); # vertical
+      h3 = h1'; # horizontal
+      ## Compute edge strength
+      switch(direction)
+        case "vertical"
+          strength = abs( conv2(im, h1, "same") );
+        case "horizontal"
+          strength = abs( conv2(im, h3, "same") );
+        case "both"
+          strength = sqrt( conv2(im, h1, "same").^2 + ...
+                           conv2(im, h3, "same").^2 );
+      endswitch
+      ## Get threshold
+      if (nargin > 2 && isscalar(varargin{1}))
+        thresh = varargin{1};
+      else
+        thresh = 4*mean(strength(:));
+      endif
+      ## Perform thresholding and simple thinning
+      strength(strength<=thresh) = 0;
+      bw = simple_thinning(strength);
+    
+    #####################################
+    ## K I R S C H
+    #####################################
+    case "kirsch"
+      ## Get the direction argument
+      direction = get_direction(varargin{:});
+      ## Create filters;
+      h1 = fspecial("kirsch"); # vertical
+      h3 = h1'; # horizontal
+      ## Compute edge strength
+      switch(direction)
+        case "vertical"
+          strength = abs( conv2(im, h1, "same") );
+        case "horizontal"
+          strength = abs( conv2(im, h3, "same") );
+        case "both"
+          strength = sqrt( conv2(im, h1, "same").^2 + ...
+                           conv2(im, h3, "same").^2 );
+      endswitch
+      ## Get threshold
+      if nargin > 2 && isscalar(varargin{1})
+        thresh = varargin{1};
+      else
+        thresh = mean(strength(:));
+      endif
+      ## Perform thresholding and simple thinning
+      strength(strength<=thresh) = 0;
+      bw = simple_thinning(strength);
+
+    #####################################
+    ## R O B E R T S
+    #####################################
+    case "roberts"
+      ## Get the thinning argument (option)
+      if (nargin == 4)
+        option = varargin{2};
+        if (!ischar(option))
+          error("edge: 'option' must be a string");
+        endif
+        option = lower(option);
+        if (!any(strcmp(option, {"thinning", "nothinning"})))
+          error("edge: 'option' must be either 'thinning', or 'nothinning'");
+        endif
+      else
+        option = "thinning";
+      endif
+      ## Create filters;
+      h1 = [1 0; 0 -1]; 
+      h2 = [0 1; -1 0]; 
+      ## Compute edge strength
+      g45  = conv2(im, h1, "same");
+      g135 = conv2(im, h2, "same");
+      strength = abs( g45 ) + abs( g135 );
+      ## Get threshold
+      if (nargin > 2 && isscalar(varargin{1}))
+        thresh = varargin{1};
+      else
+        thresh = 6*mean(strength(:));
+      endif
+      ## Perform thresholding and simple thinning
+      strength(strength<=thresh) = 0;
+      if (strcmp(option, "thinning"))
+        bw = simple_thinning(strength);
+      else
+        bw = (strength > 0);
+      endif
+      ## Check if g45 and g135 should be returned
+      if (nargout == 4)
+        g45_out  = g45;
+        g135_out = g135;
+      endif
+    
+    #####################################
+    ## L A P L A C I A N   O F   G A U S S I A N
+    #####################################
+    case "log"
+      ## Get sigma
+      if (nargin == 4 && isscalar(varargin{2}))
+        sigma = varargin{2};
+      else
+        sigma = 2;
+      endif
+      ## Create the filter
+      s = ceil(3*sigma);
+      %[x y] = meshgrid(-s:s);
+      %f = (x.^2 + y.^2 - sigma^2) .* exp(-(x.^2 + y.^2)/(2*sigma^2));
+      %f = f/sum(f(:));
+      f = fspecial("log", 2*s+1, sigma);
+      ## Perform convolution with the filter f
+      g = conv2(im, f, "same");
+      ## Get threshold
+      if (nargin > 2 && isscalar(varargin{1}))
+        thresh = varargin{1};
+      else
+        thresh = 0.75*mean(abs(g(:)));
+      endif
+      ## Find zero crossings
+      zc = zerocrossings(g);
+      bw = (abs(g) >= thresh) & zc;
+    
+    #####################################
+    ## Z E R O   C R O S S I N G 
+    #####################################
+    case "zerocross"
+      ## Get the filter
+      if (nargin == 4 && ismatrix(varargin{2}))
+        f = varargin{2};
+      else
+        error("edge: a filter must be given as the fourth argument when 'zerocross' is used");
+      endif
+      ## Perform convolution with the filter f
+      g = conv2(im, f, "same");
+      ## Get threshold
+      if (nargin > 2 && isscalar(varargin{1}))
+        thresh = varargin{1};
+      else
+        thresh = mean(abs(g(:)));
+      endif
+      ## Find zero crossings
+      zc = zerocrossings(g);
+      bw = (abs(g) >= thresh) & zc;
+
+    #####################################
+    ## C A N N Y 
+    #####################################
+    case "canny"
+      ## Get sigma
+      if (nargin == 4 && isscalar(varargin{2}))
+        sigma = varargin{2};
+      else
+        sigma = 2;
+      endif
+
+      ## Change scale
+      J = imsmooth(double(im), "Gaussian", sigma);
+
+      ## Canny enhancer
+      p = [1 0 -1]/2;
+      Jx = conv2(J, p,  "same");
+      Jy = conv2(J, p', "same");
+      Es = sqrt( Jx.^2 + Jy.^2 );
+      Eo = pi - mod (atan2 (Jy, Jx) - pi, pi);
+
+      ## Get thresholds
+      if (nargin > 2 && isscalar(varargin{1}))
+        thresh = [0.4*varargin{1}, varargin{1}];
+      elseif (nargin > 2 && ismatrix (varargin{1}) && length (varargin{1}(:)) == 2)
+        thresh = varargin{1}(:);
+      else
+        tmp = mean(abs(Es(:)));
+        thresh = [0.4*tmp, tmp];
+      endif
+      bw = nonmax_supress(Es, Eo, thresh(1), thresh(2));
+
+    #####################################
+    ## L I N D E B E R G 
+    #####################################
+    case "lindeberg"
+      ## In case the user asks for more then 1 output argument
+      ## we define thresh to be -1.
+      thresh = -1;
+      ## Get sigma
+      if (nargin > 2 && isscalar(varargin{1}))
+        sigma = varargin{1};
+      else
+        sigma = 2;
+      endif
+      ## Filters for computing the derivatives
+      Px   = [-1 0 1; -1 0 1; -1 0 1];
+      Py   = [1 1 1; 0 0 0; -1 -1 -1];
+      Pxx  = conv2(Px,  Px, "full");
+      Pyy  = conv2(Py,  Py, "full");
+      Pxy  = conv2(Px,  Py, "full");
+      Pxxx = conv2(Pxx, Px, "full");
+      Pyyy = conv2(Pyy, Py, "full");
+      Pxxy = conv2(Pxx, Py, "full");
+      Pxyy = conv2(Pyy, Px, "full");
+      ## Change scale
+      L = imsmooth(double(im), "Gaussian", sigma);
+      ## Compute derivatives
+      Lx   = conv2(L, Px,   "same");
+      Ly   = conv2(L, Py,   "same");
+      Lxx  = conv2(L, Pxx,  "same");
+      Lyy  = conv2(L, Pyy,  "same");
+      Lxy  = conv2(L, Pxy,  "same");
+      Lxxx = conv2(L, Pxxx, "same");
+      Lyyy = conv2(L, Pyyy, "same");
+      Lxxy = conv2(L, Pxxy, "same");
+      Lxyy = conv2(L, Pxyy, "same");
+      ## Compute directional derivatives
+      Lvv  = Lx.^2.*Lxx + 2.*Lx.*Ly.*Lxy + Ly.^2.*Lyy;
+      Lvvv = Lx.^3.*Lxxx + 3.*Lx.^2.*Ly.*Lxxy ...
+           + 3.*Lx.*Ly.^2.*Lxyy + 3.*Ly.^3.*Lyyy;
+      ## Perform edge detection
+      bw = zerocrossings(Lvv) & Lvvv < 0;
+
+    #####################################
+    ## A N D Y
+    #####################################
+    case "andy"
+      [bw, out_threshold] = andy (im, method, varargin{:});
+    
+    otherwise
+      error("edge: unsupported edge detector: %s", method);
+  endswitch
+  
+  if (nargout > 1)
+    out_threshold = thresh;
+  endif
+endfunction
+
+## An auxilary function that parses the 'direction' argument from 'varargin'
+function direction = get_direction(varargin)
+  if (nargin >= 2)
+    direction = varargin{2};
+    if (!ischar(direction))
+      error("edge: direction must be a string");
+    endif
+    direction = lower(direction);
+    if (!any(strcmp(direction, {"horizontal", "vertical", "both"})))
+      error("edge :direction must be either 'horizontal', 'vertical', or 'both'");
+    endif
+  else
+    direction = "both";
+  endif
+endfunction
+
+## An auxilary function that performs a very simple thinning.
+## Strength is an image containing the edge strength.
+## bw contains a 1 in (r,c) if
+##  1) strength(r,c) is greater than both neighbours in the
+##     vertical direction, OR
+##  2) strength(r,c) is greater than both neighbours in the
+##     horizontal direction.
+## Note the use of OR.
+function bw = simple_thinning(strength)
+  [r c] = size(strength);
+  x = ( strength > [ zeros(r,1) strength(:,1:end-1) ] & ...
+        strength > [ strength(:,2:end) zeros(r,1) ] );
+  y = ( strength > [ zeros(1,c); strength(1:end-1,:) ] & ...
+        strength > [ strength(2:end,:); zeros(1,c) ] );
+  bw = x | y;
+endfunction
+
+## Auxilary function. Finds the zero crossings of the 
+## 2-dimensional function f. (By Etienne Grossmann)
+function z = zerocrossings(f)
+  z0 = f<0;                 ## Negative
+  [R,C] = size(f);
+  z = zeros(R,C);
+  z(1:R-1,:) |= z0(2:R,:);  ## Grow
+  z(2:R,:) |= z0(1:R-1,:);
+  z(:,1:C-1) |= z0(:,2:C);
+  z(:,2:C) |= z0(:,1:C-1);
+
+  z &= !z0;                  ## "Positive zero-crossings"?
+endfunction
+
+## The 'andy' edge detector that was present in older versions of 'edge'.
+## The function body has simply been copied from the old implementation.
+##   -- Soren Hauberg, march 11th, 2008
+function [imout, thresh] = andy(im, method, thresh, param2)
+   [n,m]= size(im);
+   xx= 2:m-1;
+   yy= 2:n-1;
+
+   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;
+
+endfunction