octave_packages/image-1.0.15/fspecial.m

   1 ## Copyright (C) 2005 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 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{filter} = } fspecial(@var{type}, @var{arg1}, @var{arg2})
  18 ## Create spatial filters for image processing.
  19 ##
  20 ## @var{type} determines the shape of the filter and can be
  21 ## @table @t
  22 ## @item "average"
  23 ## Rectangular averaging filter. The optional argument @var{arg1} controls the
  24 ## size of the filter. If @var{arg1} is an integer @var{N}, a @var{N} by @var{N}
  25 ## filter is created. If it is a two-vector with elements @var{N} and @var{M}, the
  26 ## resulting filter will be @var{N} by @var{M}. By default a 3 by 3 filter is
  27 ## created.
  28 ## @item "disk"
  29 ## Circular averaging filter. The optional argument @var{arg1} controls the
  30 ## radius of the filter. If @var{arg1} is an integer @var{N}, a 2 @var{N} + 1
  31 ## filter is created. By default a radius of 5 is used.
  32 ## @item "gaussian"
  33 ## Gaussian filter. The optional argument @var{arg1} controls the size of the
  34 ## filter. If @var{arg1} is an integer @var{N}, a @var{N} by @var{N}
  35 ## filter is created. If it is a two-vector with elements @var{N} and @var{M}, the
  36 ## resulting filter will be @var{N} by @var{M}. By default a 3 by 3 filter is
  37 ## created. The optional argument @var{arg2} sets spread of the filter. By default
  38 ## a spread of @math{0.5} is used.
  39 ## @item "log"
  40 ## Laplacian of Gaussian. The optional argument @var{arg1} controls the size of the
  41 ## filter. If @var{arg1} is an integer @var{N}, a @var{N} by @var{N}
  42 ## filter is created. If it is a two-vector with elements @var{N} and @var{M}, the
  43 ## resulting filter will be @var{N} by @var{M}. By default a 5 by 5 filter is
  44 ## created. The optional argument @var{arg2} sets spread of the filter. By default
  45 ## a spread of @math{0.5} is used.
  46 ## @item "laplacian"
  47 ## 3x3 approximation of the laplacian. The filter is approximated as
  48 ## @example
  49 ## (4/(@var{alpha}+1))*[@var{alpha}/4,     (1-@var{alpha})/4, @var{alpha}/4; ...
  50 ##                (1-@var{alpha})/4, -1,          (1-@var{alpha})/4;  ...
  51 ##                @var{alpha}/4,     (1-@var{alpha})/4, @var{alpha}/4];
  52 ## @end example
  53 ## where @var{alpha} is a number between 0 and 1. This number can be controlled
  54 ## via the optional input argument @var{arg1}. By default it is @math{0.2}.
  55 ## @item "unsharp"
  56 ## Sharpening filter. The following filter is returned
  57 ## @example
  58 ## (1/(@var{alpha}+1))*[-@var{alpha},   @var{alpha}-1, -@var{alpha}; ...
  59 ##                 @var{alpha}-1, @var{alpha}+5,  @var{alpha}-1; ...
  60 ##                -@var{alpha},   @var{alpha}-1, -@var{alpha}];
  61 ## @end example
  62 ## where @var{alpha} is a number between 0 and 1. This number can be controlled
  63 ## via the optional input argument @var{arg1}. By default it is @math{0.2}.
  64 ## @item "motion"
  65 ## Moion blur filter of width 1 pixel. The optional input argument @var{arg1}
  66 ## controls the length of the filter, which by default is 9. The argument @var{arg2}
  67 ## controls the angle of the filter, which by default is 0 degrees.
  68 ## @item "sobel"
  69 ## Horizontal Sobel edge filter. The following filter is returned
  70 ## @example
  71 ## [ 1,  2,  1;
  72 ##   0,  0,  0;
  73 ##  -1, -2, -1 ]
  74 ## @end example
  75 ## @item "prewitt"
  76 ## Horizontal Prewitt edge filter. The following filter is returned
  77 ## @example
  78 ## [ 1,  1,  1;
  79 ##   0,  0,  0;
  80 ##  -1, -1, -1 ]
  81 ## @end example
  82 ## @item "kirsch"
  83 ## Horizontal Kirsch edge filter. The following filter is returned
  84 ## @example
  85 ## [ 3,  3,  3;
  86 ##   3,  0,  3;
  87 ##  -5, -5, -5 ]
  88 ## @end example
  89 ## @end table
  90 ## @end deftypefn
  91
  92 ## Remarks by Søren Hauberg (jan. 2nd 2007)
  93 ## The motion filter and most of the documentation was taken from Peter Kovesi's
  94 ## GPL'ed implementation of fspecial from
  95 ## http://www.csse.uwa.edu.au/~pk/research/matlabfns/OctaveCode/fspecial.m
  96
  97 function f = fspecial (type, arg1, arg2)
  98   if (!ischar (type))
  99     error ("fspecial: first argument must be a string");
 100   endif
 101
 102   switch lower(type)
 103     case "average"
 104       ## Get filtersize
 105       if (nargin > 1 && isreal (arg1) && length (arg1 (:)) <= 2)
 106         fsize = arg1 (:);
 107       else
 108         fsize = 3;
 109       endif
 110       ## Create the filter
 111       f = ones (fsize);
 112       ## Normalize the filter to integral 1
 113       f = f / sum (f (:));
 114
 115     case "disk"
 116       ## Get the radius
 117       if (nargin > 1 && isreal (arg1) && length (arg1 (:)) == 1)
 118         radius = arg1;
 119       else
 120         radius = 5;
 121       endif
 122       ## Create the filter
 123       [x, y] = meshgrid (-radius:radius, -radius:radius);
 124       r = sqrt (x.^2 + y.^2);
 125       f = (r <= radius);
 126       ## Normalize the filter to integral 1
 127       f = f / sum (f (:));
 128
 129     case "gaussian"
 130       ## Get hsize
 131       if (nargin > 1 && isreal (arg1))
 132         if (length (arg1 (:)) == 1)
 133           hsize = [arg1, arg1];
 134         elseif (length (arg1 (:)) == 2)
 135           hsize = arg1;
 136         else
 137           error ("fspecial: second argument must be a scalar or a vector of two scalars");
 138         endif
 139       else
 140         hsize = [3, 3];
 141       endif
 142       ## Get sigma
 143       if (nargin > 2 && isreal (arg2) && length (arg2 (:)) == 1)
 144         sigma = arg2;
 145       else
 146         sigma = 0.5;
 147       endif
 148       h1 = hsize (1)-1; h2 = hsize (2)-1;
 149       [x, y] = meshgrid(0:h2, 0:h1);
 150       x = x-h2/2; y = y-h1/2;
 151       gauss = exp( -( x.^2 + y.^2 ) / (2*sigma^2) );
 152       f = gauss / sum (gauss (:));
 153
 154     case "laplacian"
 155       ## Get alpha
 156       if (nargin > 1 && isscalar (arg1))
 157         alpha = arg1;
 158         if (alpha < 0 || alpha > 1)
 159           error ("fspecial: second argument must be between 0 and 1");
 160         endif
 161       else
 162         alpha = 0.2;
 163       endif
 164       ## Compute filter
 165       f = (4/(alpha+1))*[alpha/4,     (1-alpha)/4, alpha/4; ...
 166                          (1-alpha)/4, -1,          (1-alpha)/4;  ...
 167                          alpha/4,     (1-alpha)/4, alpha/4];
 168     case "log"
 169       ## Get hsize
 170       if (nargin > 1 && isreal (arg1))
 171         if (length (arg1 (:)) == 1)
 172           hsize = [arg1, arg1];
 173         elseif (length (arg1 (:)) == 2)
 174           hsize = arg1;
 175         else
 176           error ("fspecial: second argument must be a scalar or a vector of two scalars");
 177         endif
 178       else
 179         hsize = [5, 5];
 180       endif
 181       ## Get sigma
 182       if (nargin > 2 && isreal (arg2) && length (arg2 (:)) == 1)
 183         sigma = arg2;
 184       else
 185         sigma = 0.5;
 186       endif
 187       ## Compute the filter
 188       h1 = hsize (1)-1; h2 = hsize (2)-1;
 189       [x, y] = meshgrid(0:h2, 0:h1);
 190       x = x-h2/2; y = y = y-h1/2;
 191       gauss = exp( -( x.^2 + y.^2 ) / (2*sigma^2) );
 192       f = ( (x.^2 + y.^2 - 2*sigma^2).*gauss )/( 2*pi*sigma^6*sum(gauss(:)) );
 193
 194     case "motion"
 195       ## Taken (with some changes) from Peter Kovesis implementation
 196       ## (http://www.csse.uwa.edu.au/~pk/research/matlabfns/OctaveCode/fspecial.m)
 197       ## FIXME: The implementation is not quite matlab compatible.
 198       if (nargin > 1 && isreal (arg1))
 199         len = arg1;
 200       else
 201         len = 9;
 202       endif
 203       if (mod (len, 2) == 1)
 204         sze = [len, len];
 205       else
 206         sze = [len+1, len+1];
 207       end
 208       if (nargin > 2 && isreal (arg2))
 209         angle = arg2;
 210       else
 211         angle = 0;
 212       endif
 213
 214       ## First generate a horizontal line across the middle
 215       f = zeros (sze);
 216       f (floor (len/2)+1, 1:len) = 1;
 217
 218       # Then rotate to specified angle
 219       f = imrotate (f, angle, "bilinear", "loose");
 220       f = f / sum (f (:));
 221
 222     case "prewitt"
 223       ## The filter
 224       f = [1, 1, 1; 0, 0, 0; -1, -1, -1];
 225
 226     case "sobel"
 227       ## The filter
 228       f = [1, 2, 1; 0, 0, 0; -1, -2, -1];
 229
 230     case "kirsch"
 231       ## The filter
 232       f = [3, 3, 3; 3, 0, 3; -5, -5, -5];
 233
 234     case "unsharp"
 235       ## Get alpha
 236       if (nargin > 1 && isscalar (arg1))
 237         alpha = arg1;
 238         if (alpha < 0 || alpha > 1)
 239           error ("fspecial: second argument must be between 0 and 1");
 240         endif
 241       else
 242         alpha = 0.2;
 243       endif
 244       ## Compute filter
 245       f = (1/(alpha+1))*[-alpha,   alpha-1, -alpha; ...
 246                           alpha-1, alpha+5,  alpha-1; ...
 247                          -alpha,   alpha-1, -alpha];
 248
 249     otherwise
 250       error ("fspecial: filter type '%s' is not supported", type);
 251   endswitch
 252 endfunction