demo_real.sci~

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

//!  vtx() if vtx is  an m-file
Vn = forward_solver(mtlb(vtx),En,I,tol,ppn);
 
//! mtlb(elec) can be replaced by
//!  elec if elec is not an m-file
//!  elec() if elec is  an m-file
 
//! mtlb(vtx) can be replaced by
//!  vtx if vtx is not an m-file
//!  vtx() if vtx is  an m-file
 
//! mtlb(no_pl) can be replaced by
//!  no_pl if no_pl is not an m-file
//!  no_pl() if no_pl is  an m-file
[voltageH,voltageV,indH,indV,dfv] = get_3d_meas(mtlb(elec),mtlb(vtx),Vn,Ib,mtlb(no_pl));
 
//! mtlb(end) can be replaced by
//!  end if end is not an m-file
//!  end() if end is  an m-file
dfv = dfv(1:2:mtlb(end));
 
if size(dfr)~=size(dfv) then
  error('Mismatched measurements');
end
 
 
dva = refH-voltageH;
 
disp('Measurements infused with Gaussian noise ...');
 
 
 
//! mtlb_mean(dva) may be replaced by 
//!   mean(dva) if dvais a vector
//!   mean(dva,1) if dvais a matrix
dc = mtlb_mean(dva);
//DC component of the noise
noi = dc ./ 7*ones(max(size(dva)),1)+dc*rand(max(size(dva)),1,'n');
//Add the AC component
dvaG = dva+noi;
 
 
disp('Retreiving the Jacobian and regularisation matrix ...');
 
mtlb_load('datareal','J','');
mtlb_load('datareal','sm1','');
//The smoothing prior 
 
 
// Or else 
disp('Calculating the Jacobian');
disp('Please wait, this may take some time ..');
disp(mtlb_sprintf('%s','\n'));
// [v_f] = m_3d_fields(vtx,32,indH,Eref,tol,gnd_ind);
// [IntGrad] = integrofgrad(vtx,simp,mat_ref);
// [J] = jacobian_3d(I,elec,vtx,simp,gnd_ind,mat_ref,zc,IntGrad,v_f,dfr,tol,sym);
// and the smoothing regulariser
// [sm1] = iso_f_smooth(simp,4);
 
disp('Calculating a linear inverse solution, this may take some time');
disp(mtlb_sprintf('%s','\n'));
 
tfac = 0.00000001;
 
 
//! mtlb(J) can be replaced by
//!  J if J is not an m-file
//!  J() if J is  an m-file
 
//! mtlb(J) can be replaced by
//!  J if J is not an m-file
//!  J() if J is  an m-file
 
//! mtlb(sm1) can be replaced by
//!  sm1 if sm1 is not an m-file
//!  sm1() if sm1 is  an m-file
 
//! mtlb(sm1) can be replaced by
//!  sm1 if sm1 is not an m-file
//!  sm1() if sm1 is  an m-file
 
//! mtlb(J) can be replaced by
//!  J if J is not an m-file
//!  J() if J is  an m-file
sol = (mtlb(J)'*mtlb(J)+tfac*mtlb(sm1)'*mtlb(sm1))\mtlb(J)'*dvaG;
 
mtlb_load('datareal','fc');
 
xset('window',max(winsid()+1))
h1 = max(winsid());
 
//!! Scilab set function is a partial emulation of the Matlab one
set(h1,'NumberTitle','off');
 
//!! Scilab set function is a partial emulation of the Matlab one
set(h1,'Name','Simulated inhomogeneities');
 
//! mtlb(srf) can be replaced by
//!  srf if srf is not an m-file
//!  srf() if srf is  an m-file
 
//!! Unknown function vtx ,the original calling sequence is used
 
//!! Unknown function vtx ,the original calling sequence is used
 
//!! Unknown function vtx ,the original calling sequence is used
 
//!! Unknown function trimesh ,the original calling sequence is used
trimesh(mtlb(srf),vtx(:,1),vtx(:,2),vtx(:,3));
 
//!! Unknown function axis ,the original calling sequence is used
axis('image');
 
//!! Scilab set function is a partial emulation of the Matlab one
set(xget('window'),'Colormap',[0,0,0]);
 
//!! Unknown function hidden ,the original calling sequence is used
hidden('off');
mtlb_hold('on');
 
//! mtlb(vtx) can be replaced by
//!  vtx if vtx is not an m-file
//!  vtx() if vtx is  an m-file
 
//! mtlb(%simp) can be replaced by
//!  %simp if %simp is not an m-file
//!  %simp() if %simp is  an m-file
repaint_inho(mat,mat_ref,mtlb(vtx),mtlb(%simp));
 
//!! Unknown function camlight ,the original calling sequence is used
camlight('lefth');
 
//!! Unknown function lighting ,the original calling sequence is used
lighting('flat');
 
xset('window',max(winsid()+1))
h2 = max(winsid());
 
//!! Scilab set function is a partial emulation of the Matlab one
set(h2,'NumberTitle','off');
 
//!! Scilab set function is a partial emulation of the Matlab one
set(h2,'Name','Reconstructed conductivity distribution');
disp('And its reconstructed image, planes from top to bottom');
%v$2 = int((1-1)/3);%v$1 = 1-1-3*%v$2
xsetech([%v$1/3,%v$2/2,1/3,1/2]);
 
//! mtlb(vtx) can be replaced by
//!  vtx if vtx is not an m-file
//!  vtx() if vtx is  an m-file
 
//! mtlb(%simp) can be replaced by
//!  %simp if %simp is not an m-file
//!  %simp() if %simp is  an m-file
 
//! mtlb(fc) can be replaced by
//!  fc if fc is not an m-file
//!  fc() if fc is  an m-file
fc = slicer_plot(2.63,sol,mtlb(vtx),mtlb(%simp),mtlb(fc));
 
//!! Unknown function view ,the original calling sequence is used
view(2);
 
//!! Unknown function grid ,the original calling sequence is used
grid();
 
//!! Unknown function colorbar ,the original calling sequence is used
colorbar();
 
//!! Unknown function axis ,the original calling sequence is used
axis('off');
xtitle('z=2.60',' ',' ');
%v$2 = int((2-1)/3);%v$1 = 2-1-3*%v$2
xsetech([%v$1/3,%v$2/2,1/3,1/2]);
 
//! mtlb(vtx) can be replaced by
//!  vtx if vtx is not an m-file
//!  vtx() if vtx is  an m-file
 
//! mtlb(%simp) can be replaced by
//!  %simp if %simp is not an m-file
//!  %simp() if %simp is  an m-file
fc = slicer_plot(2.1,sol,mtlb(vtx),mtlb(%simp),fc);
 
//!! Unknown function view ,the original calling sequence is used
view(2);
 
//!! Unknown function grid ,the original calling sequence is used
grid();
 
//!! Unknown function colorbar ,the original calling sequence is used
colorbar();
 
//!! Unknown function axis ,the original calling sequence is used
axis('off');
xtitle('z=2.10',' ',' ');
%v$2 = int((3-1)/3);%v$1 = 3-1-3*%v$2
xsetech([%v$1/3,%v$2/2,1/3,1/2]);
 
//! mtlb(vtx) can be replaced by
//!  vtx if vtx is not an m-file
//!  vtx() if vtx is  an m-file
 
//! mtlb(%simp) can be replaced by
//!  %simp if %simp is not an m-file
//!  %simp() if %simp is  an m-file
fc = slicer_plot(1.72,sol,mtlb(vtx),mtlb(%simp),fc);
 
//!! Unknown function view ,the original calling sequence is used
view(2);
 
//!! Unknown function grid ,the original calling sequence is used
grid();
 
//!! Unknown function colorbar ,the original calling sequence is used
colorbar();
 
//!! Unknown function axis ,the original calling sequence is used
axis('off');
xtitle('z=1.70',' ',' ');
%v$2 = int((4-1)/3);%v$1 = 4-1-3*%v$2
xsetech([%v$1/3,%v$2/2,1/3,1/2]);
 
//! mtlb(vtx) can be replaced by
//!  vtx if vtx is not an m-file
//!  vtx() if vtx is  an m-file
 
//! mtlb(%simp) can be replaced by
//!  %simp if %simp is not an m-file
//!  %simp() if %simp is  an m-file
fc = slicer_plot(1.1,sol,mtlb(vtx),mtlb(%simp),fc);
 
//!! Unknown function view ,the original calling sequence is used
view(2);
 
//!! Unknown function grid ,the original calling sequence is used
grid();
 
//!! Unknown function colorbar ,the original calling sequence is used
colorbar();
 
//!! Unknown function axis ,the original calling sequence is used
axis('off');
xtitle('z=1.10',' ',' ');
%v$2 = int((5-1)/3);%v$1 = 5-1-3*%v$2
xsetech([%v$1/3,%v$2/2,1/3,1/2]);
 
//! mtlb(vtx) can be replaced by
//!  vtx if vtx is not an m-file
//!  vtx() if vtx is  an m-file
 
//! mtlb(%simp) can be replaced by
//!  %simp if %simp is not an m-file
//!  %simp() if %simp is  an m-file
fc = slicer_plot(0.83,sol,mtlb(vtx),mtlb(%simp),fc);
 
//!! Unknown function view ,the original calling sequence is used
view(2);
 
//!! Unknown function grid ,the original calling sequence is used
grid();
 
//!! Unknown function colorbar ,the original calling sequence is used
colorbar();
 
//!! Unknown function axis ,the original calling sequence is used
axis('off');
xtitle('z=0.80',' ',' ');
%v$2 = int((6-1)/3);%v$1 = 6-1-3*%v$2
xsetech([%v$1/3,%v$2/2,1/3,1/2]);
 
//! mtlb(vtx) can be replaced by
//!  vtx if vtx is not an m-file
//!  vtx() if vtx is  an m-file
 
//! mtlb(%simp) can be replaced by
//!  %simp if %simp is not an m-file
//!  %simp() if %simp is  an m-file
fc = slicer_plot(0.1,sol,mtlb(vtx),mtlb(%simp),fc);
 
//!! Unknown function view ,the original calling sequence is used
view(2);
 
//!! Unknown function grid ,the original calling sequence is used
grid();
 
//!! Unknown function colorbar ,the original calling sequence is used
colorbar();
 
//!! Unknown function axis ,the original calling sequence is used
axis('off');
xtitle('z=0.10',' ',' ');
disp('Done');
 
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// 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
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%