jacobian_3d_comp.m

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

function [Jrr,Jri,Jir,Jii] = jacobian_3d_comp(I,elec,vtx,simp,gnd_ind,mat_ref,no_pl,zc,IntGrad,v_f,df,tol,sym);
%function [Jrr,Jri,Jir,Jii] = jacobian_3d_comp(I,elec,vtx,simp,gnd_ind,mat_ref,no_pl,zc,IntGrad,v_f,df,tol,sym);
%
%This function calculates the Jacobian matrices (with respect to conductivity) for the complex EIT system.
%
%
%
%I        = The currents used
%elec     = the electrodes matrix
%vtx      = The vertices matrix
%simp     = The simplices matrix
%gnd_ind  = The ground index (node)
%mat_ref  = The reference conductivity vector (complex)
%no_pl    = The number of electrode planes 
%zc       = The electrode contact impedance vector
%IntGrad  = The integrals of the gradients 
%v_f      = The measurement fields
%df       = Measurements per current pattern as used in v_f 
%Jrr      = The part of Jacobian which corresponds to a real dV after a real dsigma 
%Jri      = The part of Jacobian which corresponds to a real dV after an imaginary dsigma
%Jir      = The part of Jacobian which corresponds to an imaginary dV after a real dsigma
%Jii      = The part of Jacobian which corresponds to an imaginary dV after an imaginary dsigma


[vr,vc] = size(vtx);
[sr,sc] = size(simp);

el_no = size(elec,1);

if sum(df)~= size(v_f,2);
   error('Mismatched data input');
end

[E,pp] = fem_master_full(vtx,simp,mat_ref,gnd_ind,elec,zc,sym);

[V] = forward_solver(vtx,E,I,tol,pp);
 
%Select the part referring to the interior nodes
V = V(1:vr,:);
v_f = v_f(1:vr,:);


VI = imag(V);
VR = real(V);

for k=1:size(simp,1)
   
 Jcol =[];
   
   for p=1:size(V,2) 
      
      for m=1:df(p) 
         
     	JJ = real(v_f(:,sum(df(1:p-1))+m))'* reshape(IntGrad(:,k),vr,vr)*VR(:,p);     
		Jcol = [Jcol;JJ];
   	end
	end
   Jrr(:,k) = Jcol;
end

disp('Partition Jrr done')%Real delta V <-- Real delta S

for k=1:size(simp,1)
   
 Jcol =[];
   
   for p=1:size(V,2) 
      
      for m=1:df(p) 
         
      JJ = real(v_f(:,sum(df(1:p-1))+m))'* reshape(IntGrad(:,k),vr,vr)*VI(:,p);     
		Jcol = [Jcol;JJ];
   	end
	end
   Jri(:,k) = Jcol;
   
end

disp('Partition Jri done')%Real delta V <-- Imaginary delta S 


for k=1:size(simp,1)
   
 Jcol =[];
   
   for p=1:size(V,2) 
      
      for m=1:df(p) 
         
      JJ = imag(v_f(:,sum(df(1:p-1))+m))'* reshape(IntGrad(:,k),vr,vr)*VR(:,p);   
		Jcol = [Jcol;JJ];
   	end
	end
   Jir(:,k) = Jcol;
   
end

disp('Partition Jir done')%Imaginary delta V <-- Real delta S


for k=1:size(simp,1)
   
 Jcol =[];
   
   for p=1:size(V,2) 
      
      for m=1:df(p) 
         
      JJ = imag(v_f(:,sum(df(1:p-1))+m))'* reshape(IntGrad(:,k),vr,vr)*VI(:,p);   
		Jcol = [Jcol;JJ];
   	end
	end
   Jii(:,k) = Jcol;
end

disp('Partition Jii done')%Imaginary delta V <-- Imaginary delta S

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This is part of the EIDORS suite.
% Copyright (c) N. Polydorides and and W.R.B. Lionheart 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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%