jacobian_3d_comp.sci

来自「用来实现三维阻抗及光学断层成像重建的matlab程序」· SCI 代码 · 共 135 行

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)
Jrr=[];Jri=[];Jir=[];Jii=[];
//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);
 
 
//! mtlb_sum(df) may be replaced by 
//!     sum(df) if dfis a vector
//!     sum(df,1) if dfis a matrix
if mtlb_sum(df)~=size(v_f,2) then
  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))'*matrix(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))'*matrix(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))'*matrix(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))'*matrix(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
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%