1 ## Copyright (C) 2005 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 2 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, see <http://www.gnu.org/licenses/>.
17 ## @deftypefn {Function File} {@var{filter} = } fspecial(@var{type}, @var{arg1}, @var{arg2})
18 ## Create spatial filters for image processing.
20 ## @var{type} determines the shape of the filter and can be
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
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.
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.
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.
47 ## 3x3 approximation of the laplacian. The filter is approximated as
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];
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}.
56 ## Sharpening filter. The following filter is returned
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}];
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}.
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.
69 ## Horizontal Sobel edge filter. The following filter is returned
76 ## Horizontal Prewitt edge filter. The following filter is returned
83 ## Horizontal Kirsch edge filter. The following filter is returned
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
97 function f = fspecial (type, arg1, arg2)
99 error ("fspecial: first argument must be a string");
105 if (nargin > 1 && isreal (arg1) && length (arg1 (:)) <= 2)
112 ## Normalize the filter to integral 1
117 if (nargin > 1 && isreal (arg1) && length (arg1 (:)) == 1)
123 [x, y] = meshgrid (-radius:radius, -radius:radius);
124 r = sqrt (x.^2 + y.^2);
126 ## Normalize the filter to integral 1
131 if (nargin > 1 && isreal (arg1))
132 if (length (arg1 (:)) == 1)
133 hsize = [arg1, arg1];
134 elseif (length (arg1 (:)) == 2)
137 error ("fspecial: second argument must be a scalar or a vector of two scalars");
143 if (nargin > 2 && isreal (arg2) && length (arg2 (:)) == 1)
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 (:));
156 if (nargin > 1 && isscalar (arg1))
158 if (alpha < 0 || alpha > 1)
159 error ("fspecial: second argument must be between 0 and 1");
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];
170 if (nargin > 1 && isreal (arg1))
171 if (length (arg1 (:)) == 1)
172 hsize = [arg1, arg1];
173 elseif (length (arg1 (:)) == 2)
176 error ("fspecial: second argument must be a scalar or a vector of two scalars");
182 if (nargin > 2 && isreal (arg2) && length (arg2 (:)) == 1)
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(:)) );
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))
203 if (mod (len, 2) == 1)
206 sze = [len+1, len+1];
208 if (nargin > 2 && isreal (arg2))
214 ## First generate a horizontal line across the middle
216 f (floor (len/2)+1, 1:len) = 1;
218 # Then rotate to specified angle
219 f = imrotate (f, angle, "bilinear", "loose");
224 f = [1, 1, 1; 0, 0, 0; -1, -1, -1];
228 f = [1, 2, 1; 0, 0, 0; -1, -2, -1];
232 f = [3, 3, 3; 3, 0, 3; -5, -5, -5];
236 if (nargin > 1 && isscalar (arg1))
238 if (alpha < 0 || alpha > 1)
239 error ("fspecial: second argument must be between 0 and 1");
245 f = (1/(alpha+1))*[-alpha, alpha-1, -alpha; ...
246 alpha-1, alpha+5, alpha-1; ...
247 -alpha, alpha-1, -alpha];
250 error ("fspecial: filter type '%s' is not supported", type);