1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71
|
function diff_im = anisodiffmoi(im, num_iter, delta_t, kappa)
% Convert input image to double.
im = double(im);
% PDE (partial differential equation) initial condition.
diff_im = im;
% Center pixel distances.
dx = 1;
dy = 1;
dd = sqrt(2);
% 2D convolution masks - finite differences.
hN = [0 1 0; 0 -1 0; 0 0 0];
hS = [0 0 0; 0 -1 0; 0 1 0];
hE = [0 0 0; 0 -1 1; 0 0 0];
hW = [0 0 0; 1 -1 0; 0 0 0];
hNE = [0 0 1; 0 -1 0; 0 0 0];
hSE = [0 0 0; 0 -1 0; 0 0 1];
hSW = [0 0 0; 0 -1 0; 1 0 0];
hNW = [1 0 0; 0 -1 0; 0 0 0];
% Anisotropic diffusion.
for t = 1:num_iter
% Finite differences. [imfilter(.,.,'conv') can be replaced by conv2(.,.,'same')]
nablaN = imfilter(diff_im,hN,'conv');
nablaS = imfilter(diff_im,hS,'conv');
nablaW = imfilter(diff_im,hW,'conv');
nablaE = imfilter(diff_im,hE,'conv');
nablaNE = imfilter(diff_im,hNE,'conv');
nablaSE = imfilter(diff_im,hSE,'conv');
nablaSW = imfilter(diff_im,hSW,'conv');
nablaNW = imfilter(diff_im,hNW,'conv');
A=[nablaN nablaS nablaW nablaE nablaNE nablaSE nablaSW nablaNW];
M= max(A);
nablaN1 = imfilter(nablaN,hN,'conv');
nablaS1 = imfilter(nablaS,hS,'conv');
nablaW1 = imfilter(nablaW,hW,'conv');
nablaE1 = imfilter(nablaE,hE,'conv');
nablaNE1 = imfilter(nablaNE,hNE,'conv');
nablaSE1 = imfilter(nablaSE,hSE,'conv');
nablaSW1 = imfilter(nablaSW,hSW,'conv');
nablaNW1 = imfilter(nablaNW,hNW,'conv');
cN = (kappa*exp(-nablaN/M))/(max(exp(nablaN1),1+nablaN));
cS = kappa*exp(-nablaS/M)/max(exp(nablaS1),1+nablaS);
cW = kappa*exp(-nablaW/M)/max(exp(nablaW1),1+nablaW);
cE = kappa*exp(-nablaE/M)/max(exp(nablaE1),1+nablaE);
cNE = kappa*exp(-nablaNE/M)/max(exp(nablaNE1),1+nablaNE);
cSE = kappa*exp(-nablaSE/M)/max(exp(nablaSE1),1+nablaSE);
cSW = kappa*exp(-nablaSW/M)/max(exp(nablaSW1),1+nablaSW);
cNW = kappa*exp(-nablaNW/M)/max(exp(nablaNW1),1+nablaNW);
% Discrete PDE solution.
diff_im = diff_im + ...
delta_t*(...
(1/(dy^2))*cN.*nablaN + (1/(dy^2))*cS.*nablaS + ...
(1/(dx^2))*cW.*nablaW + (1/(dx^2))*cE.*nablaE + ...
(1/(dd^2))*cNE.*nablaNE + (1/(dd^2))*cSE.*nablaSE + ...
(1/(dd^2))*cSW.*nablaSW + (1/(dd^2))*cNW.*nablaNW );
% Iteration warning.
fprintf('\rIteration %d\n',t);
end |
Partager