function [ imOut ] = Restoration_Discrete_TD( imIn, perc, kIni, kEps, Tol)
% Restoration_Discrete_TD restores an image using the Discrete
% Topological Derivative algorithm. The main idea behind this algorithm is
% to compute the topological derivative for an appropriate functional (the energy
% norm) and a perturbation given by the introduction of cracks between pixels.
% This derivative is used as an indicator function to find the best pixels
% to introduce the cracks that, in the presence of diffusion, will most remove noise
% preserving relevant image characteristics (see references in 'references.txt' file).
% The associated parameters are:
% imIn = Image to be restored.
% perc = (Optional) Percentage of pixels to put cracks in.
% kIni = (Optional) Conductivity in the pixels that dont have cracks.
% kEps = (Optional) Conductivity over the crack (for numerical reasosn must
% greater than zero).
% Tol = (Optional) A small value, used to stop the algorithm.
%
%
% LNCC - Laboratrio Nacional de Computao Cientfica.
% Petrpolis, RJ, Brazil, March 2007.
%
% Permission to use for evaluation purposes is granted provided that
% proper acknowledgments are given. For a commercial licence, contact
% nacho@lncc.br or feij@lncc.br.
%
% This software comes with no warranty, expressed or implied. By way of
% example, but not limitation, we make no representations of warranties
% of merchantability or fitness for any particular purpose or that the
% use of the software components or documentation will not infringe any
% patents, copyrights, trademarks, or other rights.
%
% The remaining files are copyright Ignacio Larrabide. Permission is
% granted to use the material for noncommercial and research purposes.
%
disp('Image restoration based on the discrete topological derivative.');
% im image to be processed.
im = double(imIn);
[m,n] = size(im);
% Default values for the optional parameters.
if ~exist('perc')
perc = 0.6;
end
if ~exist('kIni')
kIni = 1;
end
if ~exist('kEps')
kEps = 1e-2;
end
if ~exist('Tol')
Tol = 1e-3;
end
% To simplify the boundary treatment, the image is elarged by one pixel.
imTmp = zeros(m+2,n+2);
imTmp(2:m+1,1) = im(1:m,1);
imTmp(2:m+1,m+2) = im(1:m,n);
imTmp(m+2,n+2) = im(m,n);
imTmp(m+2,1) = im(m,1);
imTmp(1,n+2) = im(1,n);
imTmp(1,1) = im(1,1);
imTmp(1,2:n+1) = im(1,1:n);
imTmp(m+2,2:n+1) = im(m,1:n);
imTmp(2:m+1,2:n+1) = im;
% im is now resized and as boundary condition it repreats the last pixels
% value.
im = imTmp;
% ims im_star es the image that we are going to work with.
ims = im;
% dt is the time step for the evolution algorithm.
dt = 0.15;
% 16 different possibilities for the diffusion coeficient.
kTest = [kIni ,kEps ,kIni ,kEps ;
kEps ,kIni ,kIni ,kEps ;
kEps ,kIni ,kEps ,kIni ;
kIni ,kEps ,kEps ,kIni ;
kIni ,kIni ,kIni ,kEps ;
kIni ,kIni ,kEps ,kIni ;
kIni ,kEps ,kIni ,kIni ;
kEps ,kIni ,kIni ,kIni ;
kIni ,kEps ,kEps ,kEps ;
kEps ,kIni ,kEps ,kEps ;
kEps ,kEps ,kIni ,kEps ;
kEps ,kEps ,kEps ,kIni ;
kIni ,kIni ,kEps ,kEps ;
kEps ,kEps ,kIni ,kIni ;
kIni ,kIni ,kIni ,kIni ]';
% Number of possibilities for the diffusion coeficient.
Ncases = size(kTest,2);
% Initialize the topological derivative.
DT = zeros(Ncases,m+2,n+2);
%
% +-------+
% | im | j-1 -> neighbor 1
% | i-1 | j+1 -> neighbor 2
% +-------+-------+-------+ i-1 -> neighbor 3
% | im | im | im | i+1 -> neighbor 4
% | j-1 | | j+1 |
% +-------+-------+-------+
% | im |
% | i+1 |
% +-------+
%
% kRec is created bigger than the original image in orther to be able to
% use the same indexes that in the image itself.
kRec = kIni * ones(4,m+1,n+1);
% Compute the cost function in the first iteration.
cost = 0;
for j = 2:n+1,
for i = 2:m+1,
cost = cost + ( kRec(1,i,j) * (im(i,j-1)-im(i,j))^2 + kRec(2,i,j) * (im(i,j)-im(i,j+1))^2 ...
+ kRec(3,i,j) * (im(i-1,j)-im(i,j))^2 + kRec(4,i,j) * (im(i,j)-im(i+1,j))^2 );
end
end
cost = cost/(m*n);
Costs = [inf, cost];
it = 1;
disp('Iterative process started...');
% The algorithm stops when in two consecutive iterations the cost function
% decreased less that tol*[cost_function_value_at_the_begining]
while (abs(Costs(it)-Costs(it+1))>(Costs(2)*Tol))
% The easy way to set Newmann boundary conditions
imTmp = im;
imTmp(2:m+1,1) = im(2:m+1,2);
imTmp(2:m+1,n+2) = im(2:m+1,n+1);
imTmp(1,2:n+1) = im(2,2:n+1);
imTmp(m+2,2:n+1) = im(m+1,2:n+1);
im = imTmp;
% One step evolution with the same values of kRec.
for j = 2:n+1,
for i = 2:m+1,
g1 = im(i,j-1) - im(i,j);
g2 = im(i,j+1) - im(i,j);
g3 = im(i-1,j) - im(i,j);
g4 = im(i+1,j) - im(i,j);
imTmp(i,j) = im(i,j) + dt * (kRec(1,i,j) * g1 + kRec(2,i,j) * g2 + kRec(3,i,j) * g3 + kRec(4,i,j) * g4);
end
end
% Compute the topological derivativa as the total variation for each kTest
% with respect to kRec.
for j = 2:n+1,
for i = 2:m+1,
% Pre-calculate the gradient with each neigbour for the (new) evolved (t+1)
% image (to be used in the sensitivity computation)
gn1 = imTmp(i,j-1) - imTmp(i,j);
gn2 = imTmp(i,j+1) - imTmp(i,j);
gn3 = imTmp(i-1,j) - imTmp(i,j);
gn4 = imTmp(i+1,j) - imTmp(i,j);
% Pre-calculate the gradient with each neigbour for the (old) non evolved (t+1)
% image (to be used in the sensitivity computation)
go1 = im(i,j-1) - im(i,j);
go2 = im(i,j+1) - im(i,j);
go3 = im(i-1,j) - im(i,j);
go4 = im(i+1,j) - im(i,j);
for e = 1:Ncases,
% Evolve one iteration with the diffusion coeficient e, the gradient with
% each neighbour (old) are used here.
imEps(e,i,j) = im(i,j) + dt * (kTest(1,e)*go1 + kTest(2,e)*go2 + kTest(3,e)*go3 + kTest(4,e)*go4);
% Pre-calculate the gradient with each neigbour for the current diffusion
% coefficient (t+e, to be used in the sensitivity)
ge1 = imTmp(i,j-1) - imEps(e,i,j);
ge2 = imTmp(i,j+1) - imEps(e,i,j);
ge3 = imTmp(i-1,j) - imEps(e,i,j);
ge4 = imTmp(i+1,j) - imEps(e,i,j);
% Compute the sensitivity.
DT(e,i,j) = (kTest(1,e) * (ge1)^2 + kTest(2,e) * (ge2)^2 + kTest(3,e) * (ge3)^2 + kTest(4,e) * (ge4)^2) ...
- (kRec(1,i,j) * (gn1)^2 + kRec(2,i,j) * (gn2)^2 + kRec(3,i,j) * (gn3)^2 + kRec(4,i,j) * (gn4)^2);
end
end
end
MinDT = min(min(min(DT)));
% Reshape the topological derivative values.
reshDT = reshape(DT,Ncases,(m+2)*(n+2));
% Fixed poitn algorithm used to select the pixels to put cracks.
[minDTval,minDTcls] = min(reshDT,[],1);
Nchange = round(n*m*perc);
[ValOrd,IndOrd] = sort(minDTval);
% Take the diffusion coefficient kRec that produced the smaller value for each
% position.
for i = 1:Nchange,
IndOrd(i);
[indI,indJ] = ind2sub([m+2,n+2],IndOrd(i));
kRec(1,indI,indJ) = kTest(1,minDTcls(IndOrd(i)));
kRec(2,indI,indJ) = kTest(2,minDTcls(IndOrd(i)));
kRec(3,indI,indJ) = kTest(3,minDTcls(IndOrd(i)));
kRec(4,indI,indJ) = kTest(4,minDTcls(IndOrd(i)));
end
imOut = im;
% One step can be made, coltume the new processed image using kRec.
for j = 2:n+1,
for i = 2:m+1,
g1 = im(i,j-1) - im(i,j);
g2 = im(i,j+1) - im(i,j);
g3 = im(i-1,j) - im(i,j);
g4 = im(i+1,j) - im(i,j);
imOut(i,j) = im(i,j) + dt * (kRec(1,i,j) * g1 + kRec(2,i,j) * g2 + kRec(3,i,j) * g3 + kRec(4,i,j) * g4);
end
end
im = imOut;
% Compute the cost function value after the new image is obtained.
cost = 0;
for j = 2:n+1,
for i = 2:m+1,
cost = cost + (kRec(1,i,j) * (im(i,j-1)-im(i,j))^2 + kRec(2,i,j) * (im(i,j)-im(i,j+1))^2 ...
+ kRec(3,i,j) * (im(i-1,j)-im(i,j))^2 + kRec(4,i,j) * (im(i,j)-im(i+1,j))^2);
end
end
cost = cost/(m*n);
Costs = [Costs cost];
if ~mod(it,5)
disp(sprintf('Currently at iteration %d.',it));
end
it = it + 1;
end
% Return to the original size of the image.
imOut = imOut(2:m+1,2:n+1);
imOut = uint8(imOut);
disp(sprintf('Done in %d iterations.',it - 1));
