inv_sol.m

用来实现三维阻抗及光学断层成像重建的matlab程序

function [solf,solp] = inv_sol(I,voltage,tol,mat_ref,vtx,simp,elec,no_pl,zc,sym,gnd_ind,tfac,Reg,it);
%function [solf,solp] = inv_sol(J,I,voltage,tol,mat_ref,vtx,simp,elec,no_pl,zc,sym,gnd_ind,tfac,Reg,it);
%
%Calculates a Newton non-linear inverse solution.
%
%
%
%solf    = The inverse solution
%mat_ref = Initial guess on the solution
%I       = The current patterns 
%elec    = The electrode faces
%zc      = The contact impedance vector
%voltage = The measurements
%tfac    = The regularisation parameter
%Reg     = The regularisation matrix 
%it      = Number of iterations
%vtx     = The vertices matrix
%simp    = The simplices matrix
%gnd_ind = The ground index (node/electrode)
%no_pl   = The number of planes


 tol = 1e-4; %Inverse calculations error tolerance. Change accodringly 

 sol_upd = mat_ref; %Initial estimate - homogeneous background
 
 solp = zeros(size(simp,1),1); %Total change (plotting purposes)
 
 Ib = I(size(vtx,1)+1:end,:);
 
 el_no = size(elec,1);
 
for i=1:it
 
 [E,pp] = fem_master_full(vtx,simp,sol_upd,gnd_ind,elec,zc,sym);
 
 [IntGrad] = integrofgrad(vtx,simp,sol_upd);
 
 if i==1
   sprintf('Current fields for iteration %d',i)
   [V] = forward_solver(vtx,E,I,tol,pp);
   [viH,viV,indH,indV,df] = get_3d_meas(elec,vtx,V,Ib,no_pl);
   dfv = df(1:2:end);
   vi = viH; 
   sprintf('Measurement fields for iteration %d',i)
   [v_f] = m_3d_fields(vtx,el_no,indH,E,tol,gnd_ind);
else
   sprintf('Current fields for iteration %d',i)
   [V] = forward_solver(vtx,E,I,tol,pp,V);
   [viH,viV,indH,indV,df] = get_3d_meas(elec,vtx,V,Ib,no_pl);
   dfv = df(1:2:end);
   vi = viH; 
   sprintf('Measurement fields for iteration %d',i)
   [v_f] = m_3d_fields(vtx,el_no,indH,E,tol,gnd_ind,v_f);
end

[J] = jacobian_3d(I,elec,vtx,simp,gnd_ind,sol_upd,zc,IntGrad,v_f,dfv,tol,sym);
    
    
sol = (J'*J + tfac*Reg'*Reg)\ (J' * (vi - voltage));% - tfac*Reg'*Reg*sol);
%sol = pinv([J;sgrt(tfac)*Reg],tol) * [sqrt(tfac)*Reg*sol_upd; (ref-voltage)];
sol_upd = sol_upd + sol;
solp = solp + sol;
   
sprintf('Error norm at iteration %d is %f',i,norm(vi-voltage))

end %for it

solf = mat_ref + sol;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This is part of the EIDORS suite.
% Copyright (c) N. Polydorides 2001
% Copying permitted under terms of GNU GPL
% See enclosed file gpl.html for details.
% EIDORS 3D version 1.0
% MATLAB version 5.3 R11
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%