adjoint_spin.sci

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

function [JTb]=adjoint_spin(vtx,%simp,elec,x,gnd_ind,zc,I,no_pl,Vmes)
JTb=[];
//function [JTb] = adjoint_spin(vtx,simp,elec,x,gnd_ind,zc,I,no_pl,Vmes);
// 
//The function calculates the product J'*b, i.e. Jacobian transpose times 
//a (measurements) vector b, using the adjoint sources formulation.
// 
// 
// 
//JTb     = J'*b
//vtx     = The vertices matrix
//simp    = The simplices matrix
//elec    = The electrodes matrix
//x       = The solution estimate based on which J is supposed to be caclulated
//gnd_ind = The grounf index
//zc      = The electrode contact impedance
//I       = The current patterns
//no_pl   = Number of electrode planes
//Vmes    = The measurements vector
 
 
spin = size(I,2);
vr = size(vtx,1);
 
//for i=1:it %Outter loop
 
ptr = 0;
 
//Update model based on last x 
[E,pp] = fem_master_full(vtx,%simp,x,gnd_ind,elec,zc,'{n}');
EB = eye(spin,spin).*.inv(E);
 
 
for j = 1:spin:size(I,2)
  //For each group of current patterns
   
  Ib = I(:,j:j+spin-1);
  IB = matrix(Ib,size(Ib,1)*size(Ib,2),1);
   
  VB = EB*IB;
   
   
  //The current fields for I(: ,spin) only
  VB1 = matrix(VB,size(vtx,1)+size(elec,1),spin);
  VB2 = VB1(1:size(vtx,1),:);
  //These are the standard current fields used in the calculation of the Jacobian
  Vjcf = matrix(VB2,size(vtx,1)*spin,1);
  //Electrode potentials removed
   
  volts = [];
  //Some simulated data based on x 
  ind = [];
  dfh = [];
   
   
  //! mtlb(end) can be replaced by end() or end whether end is an m-file or not
  [voltageH,voltageV,indH,indV,df] = get_3d_meas(elec,vtx,VB1,I(size(vtx,1)+1:mtlb(end),j:j+spin-1),no_pl);
  volts = [volts;voltageH];
  ind = [ind;indH];
   
  //! mtlb(end) can be replaced by end() or end whether end is an m-file or not
  dfh = [dfh;df(1:2:mtlb(end))];
   
   
  //Construct the adjoint operator
  IntGrad = integrofgrad(vtx,%simp,x);
   
   
  //The measurement fields for the measurements during I(:,spin) are active
  VmI = Vmes(ptr+1:ptr+sum(dfh));
  //Part of the measurements
  ptr = ptr+sum(dfh);
   
   
  //Set up the measurement patterns for all I(:,spin)
  Imb = [];
  MC = [];
  kap = 1;
   
  for ij = 1:size(ind,1)
    m_n = zeros(size(elec,1),1);
    m_n(ind(ij,1)) = 1
    m_n(ind(ij,2)) = -1
    MC = [MC,m_n];
    if ij==sum(dfh(1:kap)) then
      if ij~=size(ind,1) then
        kap = kap+1;
      end
       
      //! mtlb(end) can be replaced by end() or end whether end is an m-file or not
       
      //!! Unknown function blkdiag ,the original calling sequence is used
      Imb = sparse(blkdiag(Imb,MC(:,1+ij-dfh(kap):mtlb(end))));
    end
  end
   
   
  //These are measurement residuals of the j'th current pattern
  b = VmI;
   
  Im = Imb*b;
   
  //Reshape it in columns. 
  Im_col = matrix(Im,size(elec,1),spin);
  mul = 1;
   
  I_s = [];
  for t = 1:spin
    I_s = [I_s;zeros(size(vtx,1),1);Im_col(:,t)*mul];
  end
   
  Vmf = EB*I_s;
  Vmf1 = matrix(Vmf,size(vtx,1)+size(elec,1),spin);
  Vmf2 = Vmf1(1:size(vtx,1),:);
  Vjmf = matrix(Vmf2,size(vtx,1)*spin,1);
  // 
   
   
  JTb = [];
   
  for p = 1:size(%simp,1)
    JTb(1,p) = (Vjmf.')*(eye(spin,spin).*.matrix(IntGrad(:,p),size(vtx,1),size(vtx,1)))*Vjcf;
  end
end