function [ imOut , kRec ] = Restoration_Continuum_TD( imIn, perc, procType, kIni, kEps)
% Restoration_Continuum_TD restores an image using the Continuum
% Topological Derivative algorithm. The main idea behind this algorithm is
% to compute the topological derivative for an appropriate functional (the gradient
% norm) and a perturbation given by the introduction of a straight crack in
% the domain. This derivative is used as an indicator function to find the
% best places to introduce cracks in the domain that will 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.
% procType = (Optional) Indicates if the resulting image will be processed
% with isotropic or anisotropic diffusion. Values possible values are 'isotropic' or
% 'anisotropic'.
% kIni = (Optional) Initial conductivity for the non-perturbed problem.
% kEps = (Optional) Conductivity over the crack (for numerical reasosn must
% greater than zero).
%
%
% 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 continuum topological derivative.');
imIn = double(imIn);
[n,m] = size(imIn);
% Default values for the optional parameters.
if ~exist('perc')
perc = 0.7;
end
if ~exist('kIni')
kIni = 0.5;
end
if ~exist('kEps')
kEps = 1e-8;
end
if ~exist('procType')
procType = 'isotropic';
end
% Configuration of the numerical solver.
tol = 1e-4;
maxiter = 100;
% Number of elements that will present cracks.
Nchange = round((n-1)*(m-1)*perc);
v = reshape(imIn, n*m, 1);
% Difussion tensor
kRec = kIni * [ones(n*m,1), zeros(n*m,1), zeros(n*m,1), ones(n*m,1)];
if strcmp(procType, 'isotropic')
disp('Isotropic analysis started...');
[DT, IEN] = AnaliseIso(v, kIni, n, m);
% Fixed point algorithm
[ValOrd,IndOrd] = sort(DT);
I = IndOrd(1:Nchange);
kRec(I,1) = kEps;
kRec(I,4) = kEps;
elseif strcmp(procType, 'anisotropic')
disp('Anisotropic analysis started...');
[DT, IEN, Vec] = AnaliseAniso(v, kIni, n, m);
% Fixed point algorithm
[ValOrd,IndOrd] = sort(DT);
I = IndOrd(1:Nchange);
for i = I
vT = Vec(:,1,i);
vN = Vec(:,2,i);
kTmp = kEps * (vN * vN') + kIni * (vT * vT');
kRec(i,1) = kTmp(1);
kRec(i,2) = kTmp(2);
kRec(i,3) = kTmp(3);
kRec(i,4) = kTmp(4);
end
end
disp('done. Proceeding with restoration...');
% Obtain Final K and M matrices.
KG = MountKG(n, m, IEN, kRec);
MG = MountMG(n, m, IEN);
KFinal = KG + MG;
% Obtain final forcing term.
Mv = MulKX(MG, v);
% Solve the final problem (i.e., restore the image).
u = pcg(KFinal, Mv, tol, maxiter, [], [], v);
% Set the output in the final image form
imOut = uint8(reshape(u, n, m));
disp('Done!');
% End of main algorithm
function [DT, IEN] = AnaliseIso(v, kIni, n, m)
% AnaliseIso: Computes the topological derivative for the insertion of a
% small crack in each pixel.
% Difussion coeficient, initially constant for all pixels, must be
% positive.
kIni = abs(kIni);
% Total number of finite elements in the image. Mesh nodes are centered in
% the image pixels centers.
nel = (n-1)*(m-1);
% Symbolic elemental matrix.
IEN = MountIEN(n, m);
% Obtain the initial difusivity tensor for all the elements.
kElems = [ones(n*m,1), zeros(n*m,1), zeros(n*m,1), ones(n*m,1)];
kElems = kIni * kElems;
% Obtain the global K matrix.
KG = MountKG(n, m, IEN, kElems);
% Obtain the global M matrix.
MG = MountMG(n, m, IEN);
% Compute the forcing term.
Mv = MulKX(MG, v);
% Compute the Final Gloubal Matrix.
KFinal = KG + MG;
% Solve the direct problem.
tol = 1e-4; maxiter = 50;
u = pcg(KFinal, Mv, tol, maxiter);
% Obtain the forcing term for the adjoint equation.
KAdj = MountKG(n, m, IEN, kElems);
KAdju = MulKX(KAdj, u);
% Solve the adjoint equation.
p = -pcg(KFinal, KAdju, tol, maxiter);
% The topological derivative is computed as a function of the problem
% solution and the ajoint equation.
[uxy] = Grad(IEN,u);
[pxy] = Grad(IEN,p);
% Compute the topological derivative for each image element.
for i = 1:nel
% The topological derivative.
tensM = - (pxy(:,i)*uxy(:,i)'+uxy(:,i)*pxy(:,i)')/2 - uxy(:,i)*uxy(:,i)';
lambda = eig(tensM);
DT(i) = min(real(lambda));
end
function [DT, IEN, Vec] = AnaliseAniso(v, kIni, n, m)
% AnaliseAniso: Computes the topological derivative for the insertion of a
% small crack in each pixel. Vec are the direction normal and tangent to
% the optimal crack direction (anisotropic case).
% Difussion coeficient, initially constant for all pixels, must be
% positive.
kIni = abs(kIni);
% Total number of finite elements in the image. Mesh nodes are centered in
% the image pixels centers.
nel = (n-1)*(m-1);
% Symbolic elemental matrix.
IEN = MountIEN(n, m);
% Obtain the initial difusivity tensor for all the elements.
kElems = [ones(n*m,1), zeros(n*m,1), zeros(n*m,1), ones(n*m,1)];
kElems = kIni * kElems;
% Obtain the global K matrix.
KG = MountKG(n, m, IEN, kElems);
% Obtain the global M matrix.
MG = MountMG(n, m, IEN);
% Compute the forcing term.
Mv = MulKX(MG, v);
% Compute the Final Gloubal Matrix.
KFinal = KG + MG;
% Solve the direct problem.
tol = 1e-4; maxiter = 50;
u = pcg(KFinal, Mv, tol, maxiter);
% Obtain the forcing term for the adjoint equation.
KAdj = MountKG(n, m, IEN, kElems);
KAdju = MulKX(KAdj, u);
% Solve the adjoint equation.
p = -pcg(KFinal, KAdju, tol, maxiter);
% The topological derivative is computed as a function of the problem
% solution and the ajoint equation.
[uxy] = Grad(IEN,u);
[pxy] = Grad(IEN,p);
% Compute the topological derivative for each image element and the
% difusivity tensor considering the normal and tangent directions of the
% crack:
% M = -[\nabla u (x) \nabla u + k(\nabla p (x)_s \nabla u)]
Vec = zeros(2,2,nel);
for i = 1:nel
tensM = - (kIni*(pxy(:,i)*uxy(:,i)'+uxy(:,i)*pxy(:,i)')/2) - (uxy(:,i)*uxy(:,i)');
[av, lambda] = eig(tensM);
[DT(i), ind] = min([real(lambda(1,1)), real(lambda(2,2))]);
if ind == 1
Vec(:,1,i) = av(:,1)/norm(av(:,1));
Vec(:,2,i) = -av(:,2)/norm(av(:,2));
else
Vec(:,1,i) = -av(:,2)/norm(av(:,2));
Vec(:,2,i) = av(:,1)/norm(av(:,1));
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SOME SIMPLE FEM TOOLS %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ K ] = MountKG( n, m, IEN, kElem)
nel = (n-1)*(m-1);
K = sparse(n*m,n*m);
for i=1:nel
Kel = MountKel(kElem(i,:));
K(IEN(i,:),IEN(i,:)) = K(IEN(i,:),IEN(i,:)) + Kel;
end
function [ kel ] = MountKel( c )
c11 = c(1); c12 = c(2); c21 = c(3); c22 = c(4);
k1 = (1/12) * (4 * c11 + 3 * c12 + 3 * c21 + 4 * c22);
k2 = (1/12) * (2 * c11 - 3 * c12 + 3 * c21 - 4 * c22);
k3 = (1/12) * (-2* c11 - 3 * c12 - 3 * c21 - 2 * c22);
k4 = (1/12) * (-4* c11 + 3 * c12 - 3 * c21 + 2 * c22);
kel = [k1, k2, k3, k4;
k2, k1, k4, k3;
k3, k4, k1, k2;
k4, k3, k2, k1];
function [ M ] = MountMG( n, m, IEN)
Mel = [1/9,1/18,1/36,1/18;
1/18,1/9,1/18,1/36;
1/36,1/18,1/9,1/18;
1/18,1/36,1/18,1/9];
nel = (n-1)*(m-1);
M = sparse(n*m,n*m);
for i=1:nel
M(IEN(i,:),IEN(i,:)) = M(IEN(i,:),IEN(i,:)) + Mel;
end
function IEN = MountIEN(n,m)
IEN = zeros((m-1)*(n-1),4);
for i=1:m-1
for j=1:n-1
IENaux = [idx(j+1,i+1,n),idx(j+1,i,n),idx(j,i,n),idx(j,i+1,n)];
nel = (n-1)*(i-1)+j;
IEN(nel,:) = IENaux;
end
end
function uxy = Grad(IEN,u)
for i = 1:size(IEN, 1)
uloc = u(IEN(i,:));
uxy(1,i) = ((uloc(2)-uloc(1))+(uloc(3)-uloc(4)))/2;
uxy(2,i) = ((uloc(4)-uloc(1))+(uloc(3)-uloc(2)))/2;
end
function id = idx(j,i,n)
id = n*(i-1) + j;
function [ Ku ] = MulKX( K, u )
Ku = K * u;