octave_packages/image-1.0.15/imsmooth.m

   1 ## Copyright (C) 2007  Soren 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 3, or (at your option)
   6 ## any later version.
   7 ##
   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.
  12 ##
  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/>.
  15
  16 ## -*- texinfo -*-
  17 ## @deftypefn {Function File} @var{J} = imsmooth(@var{I}, @var{name}, @var{options})
  18 ## Smooth the given image using several different algorithms.
  19 ##
  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
  24 ##
  25 ## @table @asis
  26 ## @item  "Gaussian"
  27 ## Isotropic Gaussian smoothing. This is the default.
  28 ## @item  "Average"
  29 ## Smoothing using a rectangular averaging linear filter.
  30 ## @item  "Disk"
  31 ## Smoothing using a circular averaging linear filter.
  32 ## @item  "Median"
  33 ## Median filtering.
  34 ## @item  "Bilateral"
  35 ## Gaussian bilateral filtering.
  36 ## @item  "Perona & Malik"
  37 ## @itemx "Perona and Malik"
  38 ## @itemx "P&M"
  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.
  42 ## @end table
  43 ##
  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.
  47 ##
  48 ## @strong{Isotropic Gaussian smoothing}
  49 ##
  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.
  53 ##
  54 ## The image is extrapolated symmetrically before the convolution operation.
  55 ##
  56 ## @strong{Rectangular averaging linear filter}
  57 ##
  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}
  62 ## filter is used.
  63 ##
  64 ## The image is extrapolated symmetrically before the convolution operation.
  65 ##
  66 ## @strong{Circular averaging linear filter}
  67 ##
  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}.
  71 ##
  72 ## The image is extrapolated symmetrically before the convolution operation.
  73 ##
  74 ## @strong{Median filtering}
  75 ##
  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}.
  80 ##
  81 ## The image is extrapolated symmetrically before the filtering is performed.
  82 ##
  83 ## @strong{Gaussian bilateral filtering}
  84 ##
  85 ## The image is smoothed using Gaussian bilateral filtering as described by
  86 ## Tomasi and Manduchi [2]. The filtering result is computed as
  87 ## @example
  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))
  89 ##                  x   y
  90 ## @end example
  91 ## where @code{k} a normalisation variable, and
  92 ## @example
  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),
  96 ## @end example
  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).
 105 ##
 106 ## The image is extrapolated symmetrically before the filtering is performed.
 107 ##
 108 ## @strong{Perona and Malik}
 109 ##
 110 ## The image is smoothed using nonlinear isotropic diffusion as described by Perona and
 111 ## Malik [1]. The algorithm iteratively updates the image using
 112 ##
 113 ## @example
 114 ## I += lambda * (g(dN).*dN + g(dS).*dS + g(dE).*dE + g(dW).*dW)
 115 ## @end example
 116 ##
 117 ## @noindent
 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.
 121 ##
 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.
 125 ##
 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.
 130 ##
 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.
 141 ##
 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.
 145 ##
 146 ## @example
 147 ## @var{g} = @@(@var{d}) exp(-(@var{d}./50).^2);
 148 ## @var{J} = imsmooth(@var{I}, "p&m", 25, 0.25, @var{g});
 149 ## @end example
 150 ##
 151 ## @strong{Custom Gaussian Smoothing}
 152 ##
 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.
 159 ##
 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.
 163 ##
 164 ## @example
 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});
 167 ## @end example
 168 ##
 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.
 173 ##
 174 ## @strong{References}
 175 ##
 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.
 180 ##
 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.
 184 ##
 185 ## @seealso{imfilter, fspecial}
 186 ## @end deftypefn
 187
 188 ## TODO: Implement Joachim Weickert's anisotropic diffusion (it's soo cool)
 189
 190 function J = imsmooth(I, name = "Gaussian", varargin)
 191   ## Check inputs
 192   if (nargin == 0)
 193     print_usage();
 194   endif
 195   if (!ismatrix(I))
 196     error("imsmooth: first input argument must be an image");
 197   endif
 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");
 201   endif
 202   if (nargin == 2 && isscalar (name))
 203     varargin {1} = name;
 204     name = "Gaussian";
 205   endif
 206   if (!ischar(name))
 207     error("imsmooth: second input must be a string");
 208   endif
 209   len = length(varargin);
 210
 211   ## Save information for later
 212   C = class(I);
 213
 214   ## Take action depending on 'name'
 215   switch (lower(name))
 216     ##############################
 217     ###   Gaussian smoothing   ###
 218     ##############################
 219     case "gaussian"
 220       ## Check input
 221       s = 0.5;
 222       if (len > 0)
 223         if (isscalar(varargin{1}) && varargin{1} > 0)
 224           s = varargin{1};
 225         else
 226           error("imsmooth: third input argument must be a positive scalar when performing Gaussian smoothing");
 227         endif
 228       endif
 229       ## Compute filter
 230       h = ceil(3*s);
 231       f = exp( (-(-h:h).^2)./(2*s^2) ); f /= sum(f);
 232       ## Pad image
 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");
 237       endfor
 238
 239     ############################
 240     ###   Square averaging   ###
 241     ############################
 242     case "average"
 243       ## Check input
 244       s = [3, 3];
 245       if (len > 0)
 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))
 249           s = varargin{1};
 250         else
 251           error("imsmooth: third input argument must be a positive scalar or two-vector when performing averaging");
 252         endif
 253       endif
 254       ## Compute filter
 255       f2 = ones(1,s(1))/s(1);
 256       f1 = ones(1,s(2))/s(2);
 257       ## Pad image
 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");
 262       endfor
 263
 264     ##############################
 265     ###   Circular averaging   ###
 266     ##############################
 267     case "disk"
 268       ## Check input
 269       r = 5;
 270       if (len > 0)
 271         if (isscalar(varargin{1}) && varargin{1} > 0)
 272           r = varargin{1};
 273         else
 274           error("imsmooth: third input argument must be a positive scalar when performing averaging");
 275         endif
 276       endif
 277       ## Compute filter
 278       f = fspecial("disk", r);
 279       ## Pad image
 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");
 284       endfor
 285
 286     ############################
 287     ###   Median Filtering   ###
 288     ############################
 289     case "median"
 290       ## Check input
 291       s = [3, 3];
 292       if (len > 0)
 293         opt = varargin{1};
 294         if (isscalar(opt) && opt > 0)
 295           s = [opt, opt];
 296         elseif (isvector(opt) && numel(opt) == 2 && all(opt>0))
 297           s = opt;
 298         else
 299           error("imsmooth: third input argument must be a positive scalar or two-vector");
 300         endif
 301         s = round(s); # just in case the use supplies non-integers.
 302       endif
 303       ## Perform the filtering
 304       for i = imchannels:-1:1
 305         J(:,:,i) = medfilt2(I(:,:,i), s, "symmetric");
 306       endfor
 307
 308     ###############################
 309     ###   Bilateral Filtering   ###
 310     ###############################
 311     case "bilateral"
 312       ## Check input
 313       if (len > 0 && !isempty(varargin{1}))
 314         if (isscalar(varargin{1}) && varargin{1} > 0)
 315           sigma_d = varargin{1};
 316         else
 317           error("imsmooth: spread of closeness function must be a positive scalar");
 318         endif
 319       else
 320         sigma_d = 2;
 321       endif
 322       if (len > 1 && !isempty(varargin{2}))
 323         if (isscalar(varargin{2}) && varargin{2} > 0)
 324           sigma_r = varargin{2};
 325         else
 326           error("imsmooth: spread of similarity function must be a positive scalar");
 327         endif
 328       else
 329         sigma_r = 10/255;
 330         if (isinteger(I)), sigma_r *= intmax(C); endif
 331       endif
 332       ## Pad image
 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);
 337
 338     ############################
 339     ###   Perona and Malik   ###
 340     ############################
 341     case {"perona & malik", "perona and malik", "p&m"}
 342       ## Check input
 343       K = 25;
 344       method1 = @(d) exp(-(d./K).^2);
 345       method2 = @(d) 1./(1 + (d./K).^2);
 346       method = method1;
 347       lambda = 0.25;
 348       iter = 10;
 349       if (len > 0 && !isempty(varargin{1}))
 350         if (isscalar(varargin{1}) && varargin{1} > 0)
 351           iter = varargin{1};
 352         else
 353           error("imsmooth: number of iterations must be a positive scalar");
 354         endif
 355       endif
 356       if (len > 1 && !isempty(varargin{2}))
 357         if (isscalar(varargin{2}) && varargin{2} > 0)
 358           lambda = varargin{2};
 359         else
 360           error("imsmooth: fourth input argument must be a scalar when using 'Perona & Malik'");
 361         endif
 362       endif
 363       if (len > 2 && !isempty(varargin{3}))
 364         fail = false;
 365         if (ischar(varargin{3}))
 366           if (strcmpi(varargin{3}, "method1"))
 367             method = method1;
 368           elseif (strcmpi(varargin{3}, "method2"))
 369             method = method2;
 370           else
 371             fail = true;
 372           endif
 373         elseif (strcmp(typeinfo(varargin{3}), "function handle"))
 374           method = varargin{3};
 375         else
 376           fail = true;
 377         endif
 378         if (fail)
 379           error("imsmooth: fifth input argument must be a function handle or the string 'method1' or 'method2' when using 'Perona & Malik'");
 380         endif
 381       endif
 382       ## Perform the filtering
 383       I = double(I);
 384       for i = imchannels:-1:1
 385         J(:,:,i) = pm(I(:,:,i), iter, lambda, method);
 386       endfor
 387
 388     #####################################
 389     ###   Custom Gaussian Smoothing   ###
 390     #####################################
 391     case "custom gaussian"
 392       ## Check input
 393       if (length (varargin) < 2)
 394         error ("imsmooth: not enough input arguments");
 395       elseif (length (varargin) == 2)
 396         varargin {3} = 0; # default theta value
 397       endif
 398       arg_names = {"third", "fourth", "fifth"};
 399       for k = 1:3
 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});
 406           endif
 407         else
 408           error ("imsmooth: %s input argument must be a scalar or a matrix", arg_names {k});
 409         endif
 410         if (!strcmp (class (varargin {k}), "double"))
 411           error ("imsmooth: %s input argument must be of class 'double'", arg_names {k});
 412         endif
 413       endfor
 414
 415       ## Perform the smoothing
 416       for i = imchannels:-1:1
 417         J(:,:,i) = __custom_gaussian_smoothing__ (I(:,:,i), varargin {:});
 418       endfor
 419
 420     ######################################
 421     ###   Mean Shift Based Smoothing   ###
 422     ######################################
 423     # NOT YET IMPLEMENTED
 424     #case "mean shift"
 425     #  J = mean_shift(I, varargin{:});
 426
 427     #############################
 428     ###   Unknown filtering   ###
 429     #############################
 430     otherwise
 431       error("imsmooth: unsupported smoothing type '%s'", name);
 432   endswitch
 433
 434   ## Cast the result to the same class as the input
 435   J = cast(J, C);
 436 endfunction
 437
 438 ## Perona and Malik for gray-scale images
 439 function J = pm(I, iter, lambda, g)
 440   ## Initialisation
 441   [imrows, imcols] = size(I);
 442   J = I;
 443
 444   for i = 1:iter
 445     ## Pad image
 446     padded = impad(J, 1, 1, "replicate");
 447
 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;
 453
 454     gN = g(dN);
 455     gS = g(dS);
 456     gE = g(dE);
 457     gW = g(dW);
 458
 459     ## Update
 460     J += lambda*(gN.*dN + gS.*dS + gE.*dE + gW.*dW);
 461   endfor
 462 endfunction
 463
 464 ## Mean Shift smoothing for gray-scale images
 465 ## XXX: This function doesn't work!!
 466 #{
 467 function J = mean_shift(I, s1, s2)
 468   sz = [size(I,2), size(I,1)];
 469   ## Mean Shift
 470   [x, y] = meshgrid(1:sz(1), 1:sz(2));
 471   f = ones(s1);
 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 );
 478
 479   ms = sub2ind(sz, ms_y, ms_x);
 480   %for i = 1:10
 481   i = 0;
 482   while (true)
 483     disp(++i)
 484     ms(ms) = ms;
 485     #new_ms = ms(ms);
 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)
 491     %    ms(j) = ms(ms(j));
 492     %  endif
 493     %endfor
 494   endwhile
 495   %endfor
 496
 497   ## Compute result
 498   J = I(ms);
 499 endfunction
 500 #}