inverse_solver

function [solf,solp] = inverse_solver(I,voltage,tol,mat_ref,vtx,simp,elec,no_pl,zc,perm_sym,gnd_ind,tfac,Reg,it);

Calculates a Newton non-linear inverse solution by iteration.



solf    = The non-linear 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

0001 function [solf,solp] = inverse_solver(I,voltage,tol,mat_ref,vtx,simp,elec,no_pl,zc,perm_sym,gnd_ind,tfac,Reg,it);
0002 %function [solf,solp] = inverse_solver(I,voltage,tol,mat_ref,vtx,simp,elec,no_pl,zc,perm_sym,gnd_ind,tfac,Reg,it);
0003 %
0004 %Calculates a Newton non-linear inverse solution by iteration.
0005 %
0006 %
0007 %
0008 %solf    = The non-linear inverse solution
0009 %mat_ref = Initial guess on the solution
0010 %I       = The current patterns
0011 %elec    = The electrode faces
0012 %zc      = The contact impedance vector
0013 %voltage = The measurements
0014 %tfac    = The regularisation parameter
0015 %Reg     = The regularisation matrix
0016 %it      = Number of iterations
0017 %vtx     = The vertices matrix
0018 %simp    = The simplices matrix
0019 %gnd_ind = The ground index (node/electrode)
0020 %no_pl   = The number of planes
0021 
0022 
0023  tol = 1e-4; %Inverse calculations error tolerance. Change accodringly
0024 
0025  sol_upd = mat_ref; %Initial estimate - homogeneous background
0026  
0027  solp = zeros(size(simp,1),1); %Total change (plotting purposes)
0028  
0029  Ib = I(size(vtx,1)+1:end,:);
0030  
0031  el_no = size(elec,1);
0032  
0033  if it==1
0034      disp('A linear solution will be calculated')
0035  end
0036  
0037 
0038 for i=1:it
0039  
0040  [E,D,Ela,pp] = fem_master_full(vtx,simp,sol_upd,gnd_ind,elec,zc,perm_sym);
0041  
0042  if i==1
0043    %sprintf('Current fields for iteration %d',i)
0044    [V] = forward_solver(E,I,tol,pp);
0045    [viH,viV,indH,indV,df] = get_3d_meas(elec,vtx,V,Ib,no_pl);
0046    dfv = df(1:2:end);
0047    vi = viH; 
0048    %sprintf('Measurement fields for iteration %d',i)
0049    [v_f] = m_3d_fields(vtx,el_no,indH,E,tol,gnd_ind);
0050 else
0051    %sprintf('Current fields for iteration %d',i)
0052    [V] = forward_solver(E,I,tol,pp,V);
0053    [viH,viV,indH,indV,df] = get_3d_meas(elec,vtx,V,Ib,no_pl);
0054    dfv = df(1:2:end);
0055    vi = viH; 
0056    %sprintf('Measurement fields for iteration %d',i)
0057    [v_f] = m_3d_fields(vtx,el_no,indH,E,tol,gnd_ind,v_f);
0058 end
0059 
0060 [J] = jacobian_3d(I,elec,vtx,simp,gnd_ind,sol_upd,zc,v_f,dfv,tol,perm_sym);    
0061 
0062 sol = (J.'*J + tfac*Reg.'*Reg)\ (J.' * (voltage - vi));
0063 %sol = pinv([J;sgrt(tfac)*Reg],tol) * [sqrt(tfac)*Reg*sol_upd; (voltage - vi)];
0064 sol_upd = sol_upd + sol;
0065 solp = solp + sol;
0066 
0067 h1 = figure;
0068 set(h1,'NumberTitle','off');
0069 set(h1,'Name','Reconstructed conductivity distribution');
0070 subplot(2,3,1); [fc] = slicer_plot_n(2.63,sol_upd,vtx,simp); view(2); grid; colorbar; axis off; title('z=2.63'); 
0071 subplot(2,3,2); [fc] = slicer_plot_n(2.10,sol_upd,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=2.10'); 
0072 subplot(2,3,3); [fc] = slicer_plot_n(1.72,sol_upd,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=1.72'); 
0073 subplot(2,3,4); [fc] = slicer_plot_n(1.10,sol_upd,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=1.10'); 
0074 subplot(2,3,5); [fc] = slicer_plot_n(0.83,sol_upd,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=0.83');
0075 subplot(2,3,6); [fc] = slicer_plot_n(0.10,sol_upd,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=0.10');
0076 drawnow;
0077 
0078 
0079 sprintf('Error norm at iteration %d is %f',i,norm(voltage - vi))
0080 
0081 end %for it
0082 
0083 solf = sol_upd;
0084 
0085 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0086 % This is part of the EIDORS suite.
0087 % Copyright (c) N. Polydorides 2003
0088 % Copying permitted under terms of GNU GPL
0089 % See enclosed file gpl.html for details.
0090 % EIDORS 3D version 2.0
0091 % MATLAB version 6.1 R12
0092 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%