% ANISODIFF - Anisotropic diffusion.%% Usage:%  diff = anisodiff(im, niter, kappa, lambda, option)%% Arguments:%         im     - input image%         niter  - number of iterations.%         kappa  - conduction coefficient 20-100 ?%         lambda - max value of .25 for stability%         option - 1 Perona Malik diffusion equation No 1%                  2 Perona Malik diffusion equation No 2%% Returns:%         diff   - diffused image.%% kappa controls conduction as a function of gradient.  If kappa is low% small intensity gradients are able to block conduction and hence diffusion% across step edges.  A large value reduces the influence of intensity% gradients on conduction.%% lambda controls speed of diffusion (you usually want it at a maximum of% 0.25)%% Diffusion equation 1 favours high contrast edges over low contrast ones.% Diffusion equation 2 favours wide regions over smaller ones.% Reference: % P. Perona and J. Malik. % Scale-space and edge detection using ansotropic diffusion.% IEEE Transactions on Pattern Analysis and Machine Intelligence, % 12(7):629-639, July 1990.%% Peter Kovesi  % School of Computer Science & Software Engineering% The University of Western Australia% pk @ csse uwa edu au% http://www.csse.uwa.edu.au%% June 2000  original version.       % March 2002 corrected diffusion eqn No 2.function diff = anisodiff(im, niter, kappa, lambda, option)if ndims(im)==3  error('Anisodiff only operates on 2D grey-scale images');endim = double(im);[rows,cols] = size(im);diff = im;  for i = 1:niter%  fprintf('\rIteration %d',i);  % Construct diffl which is the same as diff but  % has an extra padding of zeros around it.  diffl = zeros(rows+2, cols+2);  diffl(2:rows+1, 2:cols+1) = diff;  % North, South, East and West differences  deltaN = diffl(1:rows,2:cols+1)   - diff;  deltaS = diffl(3:rows+2,2:cols+1) - diff;  deltaE = diffl(2:rows+1,3:cols+2) - diff;  deltaW = diffl(2:rows+1,1:cols)   - diff;  % Conduction  if option == 1    cN = exp(-(deltaN/kappa).^2);    cS = exp(-(deltaS/kappa).^2);    cE = exp(-(deltaE/kappa).^2);    cW = exp(-(deltaW/kappa).^2);  elseif option == 2    cN = 1./(1 + (deltaN/kappa).^2);    cS = 1./(1 + (deltaS/kappa).^2);    cE = 1./(1 + (deltaE/kappa).^2);    cW = 1./(1 + (deltaW/kappa).^2);  end  diff = diff + lambda*(cN.*deltaN + cS.*deltaS + cE.*deltaE + cW.*deltaW);%  Uncomment the following to see a progression of images%  subplot(ceil(sqrt(niterations)),ceil(sqrt(niterations)), i)%  imagesc(diff), colormap(gray), axis imageend%fprintf('\n');