1 ## Copyright (C) 2007 Soren 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 3, or (at your option)
8 ## This program is distributed in the hope that it will be useful, but
9 ## WITHOUT ANY WARRANTY; without even the implied warranty of
10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 ## General Public License for more details.
13 ## You should have received a copy of the GNU General Public License
14 ## along with this file. If not, see <http://www.gnu.org/licenses/>.
17 ## @deftypefn {Function File} @var{J} = imsmooth(@var{I}, @var{name}, @var{options})
18 ## Smooth the given image using several different algorithms.
20 ## The first input argument @var{I} is the image to be smoothed. If it is an RGB
21 ## image, each color plane is treated separately.
22 ## The variable @var{name} must be a string that determines which algorithm will
23 ## be used in the smoothing. It can be any of the following strings
27 ## Isotropic Gaussian smoothing. This is the default.
29 ## Smoothing using a rectangular averaging linear filter.
31 ## Smoothing using a circular averaging linear filter.
35 ## Gaussian bilateral filtering.
36 ## @item "Perona & Malik"
37 ## @itemx "Perona and Malik"
39 ## Smoothing using nonlinear isotropic diffusion as described by Perona and Malik.
40 ## @item "Custom Gaussian"
41 ## Gaussian smoothing with a spatially varying covariance matrix.
44 ## In all algorithms the computation is done in double precision floating point
45 ## numbers, but the result has the same type as the input. Also, the size of the
46 ## smoothed image is the same as the input image.
48 ## @strong{Isotropic Gaussian smoothing}
50 ## The image is convolved with a Gaussian filter with spread @var{sigma}.
51 ## By default @var{sigma} is @math{0.5}, but this can be changed. If the third
52 ## input argument is a scalar it is used as the filter spread.
54 ## The image is extrapolated symmetrically before the convolution operation.
56 ## @strong{Rectangular averaging linear filter}
58 ## The image is convolved with @var{N} by @var{M} rectangular averaging filter.
59 ## By default a 3 by 3 filter is used, but this can e changed. If the third
60 ## input argument is a scalar @var{N} a @var{N} by @var{N} filter is used. If the third
61 ## input argument is a two-vector @code{[@var{N}, @var{M}]} a @var{N} by @var{M}
64 ## The image is extrapolated symmetrically before the convolution operation.
66 ## @strong{Circular averaging linear filter}
68 ## The image is convolved with circular averaging filter. By default the filter
69 ## has a radius of 5, but this can e changed. If the third input argument is a
70 ## scalar @var{r} the radius will be @var{r}.
72 ## The image is extrapolated symmetrically before the convolution operation.
74 ## @strong{Median filtering}
76 ## Each pixel is replaced with the median of the pixels in the local area. By
77 ## default, this area is 3 by 3, but this can be changed. If the third input
78 ## argument is a scalar @var{N} the area will be @var{N} by @var{N}, and if it's
79 ## a two-vector [@var{N}, @var{M}] the area will be @var{N} by @var{M}.
81 ## The image is extrapolated symmetrically before the filtering is performed.
83 ## @strong{Gaussian bilateral filtering}
85 ## The image is smoothed using Gaussian bilateral filtering as described by
86 ## Tomasi and Manduchi [2]. The filtering result is computed as
88 ## @var{J}(x0, y0) = k * SUM SUM @var{I}(x,y) * w(x, y, x0, y0, @var{I}(x0,y0), @var{I}(x,y))
91 ## where @code{k} a normalisation variable, and
93 ## w(x, y, x0, y0, @var{I}(x0,y0), @var{I}(x,y))
94 ## = exp(-0.5*d([x0,y0],[x,y])^2/@var{sigma_d}^2)
95 ## * exp(-0.5*d(@var{I}(x0,y0),@var{I}(x,y))^2/@var{sigma_r}^2),
97 ## with @code{d} being the Euclidian distance function. The two paramteres
98 ## @var{sigma_d} and @var{sigma_r} control the amount of smoothing. @var{sigma_d}
99 ## is the size of the spatial smoothing filter, while @var{sigma_r} is the size
100 ## of the range filter. When @var{sigma_r} is large the filter behaves almost
101 ## like the isotropic Gaussian filter with spread @var{sigma_d}, and when it is
102 ## small edges are preserved better. By default @var{sigma_d} is 2, and @var{sigma_r}
103 ## is @math{10/255} for floating points images (with integer images this is
104 ## multiplied with the maximal possible value representable by the integer class).
106 ## The image is extrapolated symmetrically before the filtering is performed.
108 ## @strong{Perona and Malik}
110 ## The image is smoothed using nonlinear isotropic diffusion as described by Perona and
111 ## Malik [1]. The algorithm iteratively updates the image using
114 ## I += lambda * (g(dN).*dN + g(dS).*dS + g(dE).*dE + g(dW).*dW)
118 ## where @code{dN} is the spatial derivative of the image in the North direction,
119 ## and so forth. The function @var{g} determines the behaviour of the diffusion.
120 ## If @math{g(x) = 1} this is standard isotropic diffusion.
122 ## The above update equation is repeated @var{iter} times, which by default is 10
123 ## times. If the third input argument is a positive scalar, that number of updates
124 ## will be performed.
126 ## The update parameter @var{lambda} affects how much smoothing happens in each
127 ## iteration. The algorithm can only be proved stable is @var{lambda} is between
128 ## 0 and 0.25, and by default it is 0.25. If the fourth input argument is given
129 ## this parameter can be changed.
131 ## The function @var{g} in the update equation determines the type of the result.
132 ## By default @code{@var{g}(@var{d}) = exp(-(@var{d}./@var{K}).^2)} where @var{K} = 25.
133 ## This choice gives privileges to high-contrast edges over low-contrast ones.
134 ## An alternative is to set @code{@var{g}(@var{d}) = 1./(1 + (@var{d}./@var{K}).^2)},
135 ## which gives privileges to wide regions over smaller ones. The choice of @var{g}
136 ## can be controlled through the fifth input argument. If it is the string
137 ## @code{"method1"}, the first mentioned function is used, and if it is @var{"method2"}
138 ## the second one is used. The argument can also be a function handle, in which case
139 ## the given function is used. It should be noted that for stability reasons,
140 ## @var{g} should return values between 0 and 1.
142 ## The following example shows how to set
143 ## @code{@var{g}(@var{d}) = exp(-(@var{d}./@var{K}).^2)} where @var{K} = 50.
144 ## The update will be repeated 25 times, with @var{lambda} = 0.25.
147 ## @var{g} = @@(@var{d}) exp(-(@var{d}./50).^2);
148 ## @var{J} = imsmooth(@var{I}, "p&m", 25, 0.25, @var{g});
151 ## @strong{Custom Gaussian Smoothing}
153 ## The image is smoothed using a Gaussian filter with a spatially varying covariance
154 ## matrix. The third and fourth input arguments contain the Eigenvalues of the
155 ## covariance matrix, while the fifth contains the rotation of the Gaussian.
156 ## These arguments can be matrices of the same size as the input image, or scalars.
157 ## In the last case the scalar is used in all pixels. If the rotation is not given
158 ## it defaults to zero.
160 ## The following example shows how to increase the size of an Gaussian
161 ## filter, such that it is small near the upper right corner of the image, and
162 ## large near the lower left corner.
165 ## [@var{lambda1}, @var{lambda2}] = meshgrid (linspace (0, 25, columns (@var{I})), linspace (0, 25, rows (@var{I})));
166 ## @var{J} = imsmooth (@var{I}, "Custom Gaussian", @var{lambda1}, @var{lambda2});
169 ## The implementation uses an elliptic filter, where only neighbouring pixels
170 ## with a Mahalanobis' distance to the current pixel that is less than 3 are
171 ## used to compute the response. The response is computed using double precision
172 ## floating points, but the result is of the same class as the input image.
174 ## @strong{References}
176 ## [1] P. Perona and J. Malik,
177 ## "Scale-space and edge detection using anisotropic diffusion",
178 ## IEEE Transactions on Pattern Analysis and Machine Intelligence,
179 ## 12(7):629-639, 1990.
181 ## [2] C. Tomasi and R. Manduchi,
182 ## "Bilateral Filtering for Gray and Color Images",
183 ## Proceedings of the 1998 IEEE International Conference on Computer Vision, Bombay, India.
185 ## @seealso{imfilter, fspecial}
188 ## TODO: Implement Joachim Weickert's anisotropic diffusion (it's soo cool)
190 function J = imsmooth(I, name = "Gaussian", varargin)
196 error("imsmooth: first input argument must be an image");
198 [imrows, imcols, imchannels, tmp] = size(I);
199 if ((imchannels != 1 && imchannels != 3) || tmp != 1)
200 error("imsmooth: first input argument must be an image");
202 if (nargin == 2 && isscalar (name))
207 error("imsmooth: second input must be a string");
209 len = length(varargin);
211 ## Save information for later
214 ## Take action depending on 'name'
216 ##############################
217 ### Gaussian smoothing ###
218 ##############################
223 if (isscalar(varargin{1}) && varargin{1} > 0)
226 error("imsmooth: third input argument must be a positive scalar when performing Gaussian smoothing");
231 f = exp( (-(-h:h).^2)./(2*s^2) ); f /= sum(f);
233 I = double(impad(I, h, h, "symmetric"));
234 ## Perform the filtering
235 for i = imchannels:-1:1
236 J(:,:,i) = conv2(f, f, I(:,:,i), "valid");
239 ############################
240 ### Square averaging ###
241 ############################
246 if (isscalar(varargin{1}) && varargin{1} > 0)
247 s = [varargin{1}, varargin{1}];
248 elseif (isvector(varargin{1}) && length(varargin{1}) == 2 && all(varargin{1} > 0))
251 error("imsmooth: third input argument must be a positive scalar or two-vector when performing averaging");
255 f2 = ones(1,s(1))/s(1);
256 f1 = ones(1,s(2))/s(2);
258 I = impad(double(I), floor([s(2), s(2)-1]/2), floor([s(1), s(1)-1]/2), "symmetric");
259 ## Perform the filtering
260 for i = imchannels:-1:1
261 J(:,:,i) = conv2(f1, f2, I(:,:,i), "valid");
264 ##############################
265 ### Circular averaging ###
266 ##############################
271 if (isscalar(varargin{1}) && varargin{1} > 0)
274 error("imsmooth: third input argument must be a positive scalar when performing averaging");
278 f = fspecial("disk", r);
280 I = impad(double(I), r, r, "symmetric");
281 ## Perform the filtering
282 for i = imchannels:-1:1
283 J(:,:,i) = conv2(I(:,:,i), f, "valid");
286 ############################
287 ### Median Filtering ###
288 ############################
294 if (isscalar(opt) && opt > 0)
296 elseif (isvector(opt) && numel(opt) == 2 && all(opt>0))
299 error("imsmooth: third input argument must be a positive scalar or two-vector");
301 s = round(s); # just in case the use supplies non-integers.
303 ## Perform the filtering
304 for i = imchannels:-1:1
305 J(:,:,i) = medfilt2(I(:,:,i), s, "symmetric");
308 ###############################
309 ### Bilateral Filtering ###
310 ###############################
313 if (len > 0 && !isempty(varargin{1}))
314 if (isscalar(varargin{1}) && varargin{1} > 0)
315 sigma_d = varargin{1};
317 error("imsmooth: spread of closeness function must be a positive scalar");
322 if (len > 1 && !isempty(varargin{2}))
323 if (isscalar(varargin{2}) && varargin{2} > 0)
324 sigma_r = varargin{2};
326 error("imsmooth: spread of similarity function must be a positive scalar");
330 if (isinteger(I)), sigma_r *= intmax(C); endif
333 s = max([round(3*sigma_d),1]);
334 I = impad(I, s, s, "symmetric");
335 ## Perform the filtering
336 J = __bilateral__(I, sigma_d, sigma_r);
338 ############################
339 ### Perona and Malik ###
340 ############################
341 case {"perona & malik", "perona and malik", "p&m"}
344 method1 = @(d) exp(-(d./K).^2);
345 method2 = @(d) 1./(1 + (d./K).^2);
349 if (len > 0 && !isempty(varargin{1}))
350 if (isscalar(varargin{1}) && varargin{1} > 0)
353 error("imsmooth: number of iterations must be a positive scalar");
356 if (len > 1 && !isempty(varargin{2}))
357 if (isscalar(varargin{2}) && varargin{2} > 0)
358 lambda = varargin{2};
360 error("imsmooth: fourth input argument must be a scalar when using 'Perona & Malik'");
363 if (len > 2 && !isempty(varargin{3}))
365 if (ischar(varargin{3}))
366 if (strcmpi(varargin{3}, "method1"))
368 elseif (strcmpi(varargin{3}, "method2"))
373 elseif (strcmp(typeinfo(varargin{3}), "function handle"))
374 method = varargin{3};
379 error("imsmooth: fifth input argument must be a function handle or the string 'method1' or 'method2' when using 'Perona & Malik'");
382 ## Perform the filtering
384 for i = imchannels:-1:1
385 J(:,:,i) = pm(I(:,:,i), iter, lambda, method);
388 #####################################
389 ### Custom Gaussian Smoothing ###
390 #####################################
391 case "custom gaussian"
393 if (length (varargin) < 2)
394 error ("imsmooth: not enough input arguments");
395 elseif (length (varargin) == 2)
396 varargin {3} = 0; # default theta value
398 arg_names = {"third", "fourth", "fifth"};
400 if (isscalar (varargin {k}))
401 varargin {k} = repmat (varargin {k}, imrows, imcols);
402 elseif (ismatrix (varargin {k}) && ndims (varargin {k}) == 2)
403 if (rows (varargin {k}) != imrows || columns (varargin {k}) != imcols)
404 error (["imsmooth: %s input argument must have same number of rows "
405 "and columns as the input image"], arg_names {k});
408 error ("imsmooth: %s input argument must be a scalar or a matrix", arg_names {k});
410 if (!strcmp (class (varargin {k}), "double"))
411 error ("imsmooth: %s input argument must be of class 'double'", arg_names {k});
415 ## Perform the smoothing
416 for i = imchannels:-1:1
417 J(:,:,i) = __custom_gaussian_smoothing__ (I(:,:,i), varargin {:});
420 ######################################
421 ### Mean Shift Based Smoothing ###
422 ######################################
423 # NOT YET IMPLEMENTED
425 # J = mean_shift(I, varargin{:});
427 #############################
428 ### Unknown filtering ###
429 #############################
431 error("imsmooth: unsupported smoothing type '%s'", name);
434 ## Cast the result to the same class as the input
438 ## Perona and Malik for gray-scale images
439 function J = pm(I, iter, lambda, g)
441 [imrows, imcols] = size(I);
446 padded = impad(J, 1, 1, "replicate");
448 ## Spatial derivatives
449 dN = padded(1:imrows, 2:imcols+1) - J;
450 dS = padded(3:imrows+2, 2:imcols+1) - J;
451 dE = padded(2:imrows+1, 3:imcols+2) - J;
452 dW = padded(2:imrows+1, 1:imcols) - J;
460 J += lambda*(gN.*dN + gS.*dS + gE.*dE + gW.*dW);
464 ## Mean Shift smoothing for gray-scale images
465 ## XXX: This function doesn't work!!
467 function J = mean_shift(I, s1, s2)
468 sz = [size(I,2), size(I,1)];
470 [x, y] = meshgrid(1:sz(1), 1:sz(2));
472 tmp = conv2(ones(sz(2), sz(1)), f, "same"); # We use normalised convolution to handle the border
473 m00 = conv2(I, f, "same")./tmp;
474 m10 = conv2(I.*x, f, "same")./tmp;
475 m01 = conv2(I.*y, f, "same")./tmp;
476 ms_x = round( m10./m00 ); # konverter ms_x og ms_y til linære indices og arbejd med dem!
477 ms_y = round( m01./m00 );
479 ms = sub2ind(sz, ms_y, ms_x);
486 if (i >200), break; endif
487 #idx = ( abs(I(ms)-I(new_ms)) < s2 );
488 #ms(idx) = new_ms(idx);
489 %for j = 1:length(ms)
490 % if (abs(I(ms(j))-I(ms(ms(j)))) < s2)