X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=octave_packages%2Fimage-1.0.15%2Ffspecial.m;fp=octave_packages%2Fimage-1.0.15%2Ffspecial.m;h=540bb59edbe50ed8215e9eb2f1f4999475032470;hb=c880e8788dfc484bf23ce13fa2787f2c6bca4863;hp=0000000000000000000000000000000000000000;hpb=1705066eceaaea976f010f669ce8e972f3734b05;p=CreaPhase.git diff --git a/octave_packages/image-1.0.15/fspecial.m b/octave_packages/image-1.0.15/fspecial.m new file mode 100644 index 0000000..540bb59 --- /dev/null +++ b/octave_packages/image-1.0.15/fspecial.m @@ -0,0 +1,252 @@ +## Copyright (C) 2005 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 2 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, see . + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{filter} = } fspecial(@var{type}, @var{arg1}, @var{arg2}) +## Create spatial filters for image processing. +## +## @var{type} determines the shape of the filter and can be +## @table @t +## @item "average" +## Rectangular averaging filter. The optional argument @var{arg1} controls the +## size of the filter. If @var{arg1} is an integer @var{N}, a @var{N} by @var{N} +## filter is created. If it is a two-vector with elements @var{N} and @var{M}, the +## resulting filter will be @var{N} by @var{M}. By default a 3 by 3 filter is +## created. +## @item "disk" +## Circular averaging filter. The optional argument @var{arg1} controls the +## radius of the filter. If @var{arg1} is an integer @var{N}, a 2 @var{N} + 1 +## filter is created. By default a radius of 5 is used. +## @item "gaussian" +## Gaussian filter. The optional argument @var{arg1} controls the size of the +## filter. If @var{arg1} is an integer @var{N}, a @var{N} by @var{N} +## filter is created. If it is a two-vector with elements @var{N} and @var{M}, the +## resulting filter will be @var{N} by @var{M}. By default a 3 by 3 filter is +## created. The optional argument @var{arg2} sets spread of the filter. By default +## a spread of @math{0.5} is used. +## @item "log" +## Laplacian of Gaussian. The optional argument @var{arg1} controls the size of the +## filter. If @var{arg1} is an integer @var{N}, a @var{N} by @var{N} +## filter is created. If it is a two-vector with elements @var{N} and @var{M}, the +## resulting filter will be @var{N} by @var{M}. By default a 5 by 5 filter is +## created. The optional argument @var{arg2} sets spread of the filter. By default +## a spread of @math{0.5} is used. +## @item "laplacian" +## 3x3 approximation of the laplacian. The filter is approximated as +## @example +## (4/(@var{alpha}+1))*[@var{alpha}/4, (1-@var{alpha})/4, @var{alpha}/4; ... +## (1-@var{alpha})/4, -1, (1-@var{alpha})/4; ... +## @var{alpha}/4, (1-@var{alpha})/4, @var{alpha}/4]; +## @end example +## where @var{alpha} is a number between 0 and 1. This number can be controlled +## via the optional input argument @var{arg1}. By default it is @math{0.2}. +## @item "unsharp" +## Sharpening filter. The following filter is returned +## @example +## (1/(@var{alpha}+1))*[-@var{alpha}, @var{alpha}-1, -@var{alpha}; ... +## @var{alpha}-1, @var{alpha}+5, @var{alpha}-1; ... +## -@var{alpha}, @var{alpha}-1, -@var{alpha}]; +## @end example +## where @var{alpha} is a number between 0 and 1. This number can be controlled +## via the optional input argument @var{arg1}. By default it is @math{0.2}. +## @item "motion" +## Moion blur filter of width 1 pixel. The optional input argument @var{arg1} +## controls the length of the filter, which by default is 9. The argument @var{arg2} +## controls the angle of the filter, which by default is 0 degrees. +## @item "sobel" +## Horizontal Sobel edge filter. The following filter is returned +## @example +## [ 1, 2, 1; +## 0, 0, 0; +## -1, -2, -1 ] +## @end example +## @item "prewitt" +## Horizontal Prewitt edge filter. The following filter is returned +## @example +## [ 1, 1, 1; +## 0, 0, 0; +## -1, -1, -1 ] +## @end example +## @item "kirsch" +## Horizontal Kirsch edge filter. The following filter is returned +## @example +## [ 3, 3, 3; +## 3, 0, 3; +## -5, -5, -5 ] +## @end example +## @end table +## @end deftypefn + +## Remarks by Søren Hauberg (jan. 2nd 2007) +## The motion filter and most of the documentation was taken from Peter Kovesi's +## GPL'ed implementation of fspecial from +## http://www.csse.uwa.edu.au/~pk/research/matlabfns/OctaveCode/fspecial.m + +function f = fspecial (type, arg1, arg2) + if (!ischar (type)) + error ("fspecial: first argument must be a string"); + endif + + switch lower(type) + case "average" + ## Get filtersize + if (nargin > 1 && isreal (arg1) && length (arg1 (:)) <= 2) + fsize = arg1 (:); + else + fsize = 3; + endif + ## Create the filter + f = ones (fsize); + ## Normalize the filter to integral 1 + f = f / sum (f (:)); + + case "disk" + ## Get the radius + if (nargin > 1 && isreal (arg1) && length (arg1 (:)) == 1) + radius = arg1; + else + radius = 5; + endif + ## Create the filter + [x, y] = meshgrid (-radius:radius, -radius:radius); + r = sqrt (x.^2 + y.^2); + f = (r <= radius); + ## Normalize the filter to integral 1 + f = f / sum (f (:)); + + case "gaussian" + ## Get hsize + if (nargin > 1 && isreal (arg1)) + if (length (arg1 (:)) == 1) + hsize = [arg1, arg1]; + elseif (length (arg1 (:)) == 2) + hsize = arg1; + else + error ("fspecial: second argument must be a scalar or a vector of two scalars"); + endif + else + hsize = [3, 3]; + endif + ## Get sigma + if (nargin > 2 && isreal (arg2) && length (arg2 (:)) == 1) + sigma = arg2; + else + sigma = 0.5; + endif + h1 = hsize (1)-1; h2 = hsize (2)-1; + [x, y] = meshgrid(0:h2, 0:h1); + x = x-h2/2; y = y-h1/2; + gauss = exp( -( x.^2 + y.^2 ) / (2*sigma^2) ); + f = gauss / sum (gauss (:)); + + case "laplacian" + ## Get alpha + if (nargin > 1 && isscalar (arg1)) + alpha = arg1; + if (alpha < 0 || alpha > 1) + error ("fspecial: second argument must be between 0 and 1"); + endif + else + alpha = 0.2; + endif + ## Compute filter + f = (4/(alpha+1))*[alpha/4, (1-alpha)/4, alpha/4; ... + (1-alpha)/4, -1, (1-alpha)/4; ... + alpha/4, (1-alpha)/4, alpha/4]; + case "log" + ## Get hsize + if (nargin > 1 && isreal (arg1)) + if (length (arg1 (:)) == 1) + hsize = [arg1, arg1]; + elseif (length (arg1 (:)) == 2) + hsize = arg1; + else + error ("fspecial: second argument must be a scalar or a vector of two scalars"); + endif + else + hsize = [5, 5]; + endif + ## Get sigma + if (nargin > 2 && isreal (arg2) && length (arg2 (:)) == 1) + sigma = arg2; + else + sigma = 0.5; + endif + ## Compute the filter + h1 = hsize (1)-1; h2 = hsize (2)-1; + [x, y] = meshgrid(0:h2, 0:h1); + x = x-h2/2; y = y = y-h1/2; + gauss = exp( -( x.^2 + y.^2 ) / (2*sigma^2) ); + f = ( (x.^2 + y.^2 - 2*sigma^2).*gauss )/( 2*pi*sigma^6*sum(gauss(:)) ); + + case "motion" + ## Taken (with some changes) from Peter Kovesis implementation + ## (http://www.csse.uwa.edu.au/~pk/research/matlabfns/OctaveCode/fspecial.m) + ## FIXME: The implementation is not quite matlab compatible. + if (nargin > 1 && isreal (arg1)) + len = arg1; + else + len = 9; + endif + if (mod (len, 2) == 1) + sze = [len, len]; + else + sze = [len+1, len+1]; + end + if (nargin > 2 && isreal (arg2)) + angle = arg2; + else + angle = 0; + endif + + ## First generate a horizontal line across the middle + f = zeros (sze); + f (floor (len/2)+1, 1:len) = 1; + + # Then rotate to specified angle + f = imrotate (f, angle, "bilinear", "loose"); + f = f / sum (f (:)); + + case "prewitt" + ## The filter + f = [1, 1, 1; 0, 0, 0; -1, -1, -1]; + + case "sobel" + ## The filter + f = [1, 2, 1; 0, 0, 0; -1, -2, -1]; + + case "kirsch" + ## The filter + f = [3, 3, 3; 3, 0, 3; -5, -5, -5]; + + case "unsharp" + ## Get alpha + if (nargin > 1 && isscalar (arg1)) + alpha = arg1; + if (alpha < 0 || alpha > 1) + error ("fspecial: second argument must be between 0 and 1"); + endif + else + alpha = 0.2; + endif + ## Compute filter + f = (1/(alpha+1))*[-alpha, alpha-1, -alpha; ... + alpha-1, alpha+5, alpha-1; ... + -alpha, alpha-1, -alpha]; + + otherwise + error ("fspecial: filter type '%s' is not supported", type); + endswitch +endfunction