demo_real.sci

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

function []=demo_real()
clear;
 
//! clc ignored
// clc
 ;
disp('This is a demo for reconstructing real impedance changes');
disp(mtlb_sprintf('%s','\n'));
// 
mtlb_load('datareal','srf','vtx','simp');
// 
 
//! mtlb(srf) can be replaced by srf() or srf whether srf is an m-file or not
 
//!! Unknown variable vtx ,the original calling sequence is used
 
//!! Unknown variable vtx ,the original calling sequence is used
 
//!! Unknown variable vtx ,the original calling sequence is used
 
//!! Unknown variable trimesh ,the original calling sequence is used
trimesh(mtlb(srf),vtx(:,1),vtx(:,2),vtx(:,3));
 
//!! Unknown variable 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]);
mtlb_hold('on');
// 
disp('This is a cylindrical mesh, wait to attach the electrodes');
disp(mtlb_sprintf('%s','\n'));
// 
xpause(1000*(4));
// 
mtlb_load('datacom','sels');
// 
 
//! mtlb(sels) can be replaced by sels() or sels whether sels is an m-file or not
for u = 1:size(mtlb(sels))
   
  //!! Unknown variable sels ,the original calling sequence is used
   
  //! mtlb(srf) can be replaced by srf() or srf whether srf is an m-file or not
   
  //! mtlb(vtx) can be replaced by vtx() or vtx whether vtx is an m-file or not
  paint_electrodes(sels(u),mtlb(srf),mtlb(vtx));
end
// 
 
//!! Unknown variable hidden ,the original calling sequence is used
hidden('off');
mtlb_load('datareal','gnd_ind','elec','zc','sym');
// 
//Change sym to '{n}' for Cholesky instead of PCG in the forward solver
sym = '{n}';
// 
mtlb_load('datareal','protocol','no_pl');
// 
 
//! mtlb(protocol) can be replaced by protocol() or protocol whether protocol is an m-file or not
 
//! mtlb(elec) can be replaced by elec() or elec whether elec is an m-file or not
 
//! mtlb(vtx) can be replaced by vtx() or vtx whether vtx is an m-file or not
 
//! mtlb(gnd_ind) can be replaced by gnd_ind() or gnd_ind whether gnd_ind is an m-file or not
 
//! mtlb(no_pl) can be replaced by no_pl() or no_pl whether no_pl is an m-file or not
[I,Ib] = set_3d_currents(mtlb(protocol),mtlb(elec),mtlb(vtx),mtlb(gnd_ind),mtlb(no_pl));
// 
 
//! clc ignored
// clc
 ;
disp('Adjacent current patterns selected. Calculating reference measurements');
disp(mtlb_sprintf('%s','\n'));
// 
xpause(1000*(2));
// 
mat_ref = ones(828,1);
//%%%%%
// 
//Set the tolerance for the forward solver
tol = 0.0001;
// 
 
//! mtlb(vtx) can be replaced by vtx() or vtx whether vtx is an m-file or not
 
//! mtlb(%simp) can be replaced by %simp() or %simp whether %simp is an m-file or not
 
//! mtlb(gnd_ind) can be replaced by gnd_ind() or gnd_ind whether gnd_ind is an m-file or not
 
//! mtlb(elec) can be replaced by elec() or elec whether elec is an m-file or not
 
//! mtlb(zc) can be replaced by zc() or zc whether zc is an m-file or not
[Eref,ppr] = fem_master_full(mtlb(vtx),mtlb(%simp),mat_ref,mtlb(gnd_ind),mtlb(elec),mtlb(zc),sym);
 
//! mtlb(vtx) can be replaced by vtx() or vtx whether vtx is an m-file or not
Vref = forward_solver(mtlb(vtx),Eref,I,tol,ppr);
 
//! mtlb(elec) can be replaced by elec() or elec whether elec is an m-file or not
 
//! mtlb(vtx) can be replaced by vtx() or vtx whether vtx is an m-file or not
 
//! mtlb(no_pl) can be replaced by no_pl() or no_pl whether no_pl is an m-file or not
[refH,refV,indH,indV,dfr] = get_3d_meas(mtlb(elec),mtlb(vtx),Vref,Ib,mtlb(no_pl));
 
//! mtlb(end) can be replaced by end() or end whether end is an m-file or not
dfr = dfr(1:2:mtlb(end));
//Taking just the horrizontal measurements
// 
xdel();
 
//! clc ignored
// clc
 ;
disp('Now let us create a couple of changes ...');
disp(mtlb_sprintf('%s','\n'));
// 
xpause(1000*(2));
// 
// 
mat = mat_ref;
// 
mtlb_load('datareal','A','B');
//figure; [mat,grp] = set_inho(srf,simp,vtx,mat_ref,1.1); 
sA = 1.2;
mtlb_sprintf('This one  %f',sA);
 
//! mtlb(A) can be replaced by A() or A whether A is an m-file or not
mat(mtlb(A)) = sA
// 
xset('window',max(winsid()+1))
max(winsid());
 
//! mtlb(srf) can be replaced by srf() or srf whether srf is an m-file or not
 
//!! Unknown variable vtx ,the original calling sequence is used
 
//!! Unknown variable vtx ,the original calling sequence is used
 
//!! Unknown variable vtx ,the original calling sequence is used
 
//!! Unknown variable trimesh ,the original calling sequence is used
trimesh(mtlb(srf),vtx(:,1),vtx(:,2),vtx(:,3));
 
//!! Unknown variable 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 variable hidden ,the original calling sequence is used
hidden('off');
mtlb_hold('on');
 
//! mtlb(vtx) can be replaced by vtx() or vtx whether vtx is an m-file or not
 
//! mtlb(%simp) can be replaced by %simp() or %simp whether %simp is an m-file or not
repaint_inho(mat,mat_ref,mtlb(vtx),mtlb(%simp));
 
//!! Unknown variable camlight ,the original calling sequence is used
camlight('lefth');
 
//!! Unknown variable lighting ,the original calling sequence is used
lighting('flat');
 
//!! drawnow ignored
// drawnow
 ;
// 
xpause(1000*(4));
xdel();
// 
sB = 0.9;
 
//! mtlb(B) can be replaced by B() or B whether B is an m-file or not
mat(mtlb(B)) = sB
// 
xset('window',max(winsid()+1))
max(winsid());
 
//! mtlb(srf) can be replaced by srf() or srf whether srf is an m-file or not
 
//!! Unknown variable vtx ,the original calling sequence is used
 
//!! Unknown variable vtx ,the original calling sequence is used
 
//!! Unknown variable vtx ,the original calling sequence is used
 
//!! Unknown variable trimesh ,the original calling sequence is used
trimesh(mtlb(srf),vtx(:,1),vtx(:,2),vtx(:,3));
 
//!! Unknown variable 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 variable hidden ,the original calling sequence is used
hidden('off');
mtlb_hold('on');
 
//! mtlb(vtx) can be replaced by vtx() or vtx whether vtx is an m-file or not
 
//! mtlb(%simp) can be replaced by %simp() or %simp whether %simp is an m-file or not
repaint_inho(mat,mat_ref,mtlb(vtx),mtlb(%simp));
 
//!! Unknown variable camlight ,the original calling sequence is used
camlight('lefth');
 
//!! Unknown variable lighting ,the original calling sequence is used
lighting('flat');
// 
mtlb_sprintf('And this one at  %f',sB);
disp(mtlb_sprintf('%s','\n'));
// 
 
//!! drawnow ignored
// drawnow
 ;
xpause(1000*(4));
xdel();
// 
//[mat,grp] = set_inho(srf,simp,vtx,mat,0.999); 
disp('Simulating measurements based on the CEM ...');
// 
// 
 
//! mtlb(vtx) can be replaced by vtx() or vtx whether vtx is an m-file or not
 
//! mtlb(%simp) can be replaced by %simp() or %simp whether %simp is an m-file or not
 
//! mtlb(gnd_ind) can be replaced by gnd_ind() or gnd_ind whether gnd_ind is an m-file or not
 
//! mtlb(elec) can be replaced by elec() or elec whether elec is an m-file or not
 
//! mtlb(zc) can be replaced by zc() or zc whether zc is an m-file or not
[En,ppn] = fem_master_full(mtlb(vtx),mtlb(%simp),mat,mtlb(gnd_ind),mtlb(elec),mtlb(zc),sym);
 
//! mtlb(vtx) can be replaced by vtx() or vtx whether vtx is an m-file or not
Vn = forward_solver(mtlb(vtx),En,I,tol,ppn);
 
//! mtlb(elec) can be replaced by elec() or elec whether elec is an m-file or not
 
//! mtlb(vtx) can be replaced by vtx() or vtx whether vtx is an m-file or not
 
//! mtlb(no_pl) can be replaced by no_pl() or no_pl whether no_pl is an m-file or not
[voltageH,voltageV,indH,indV,dfv] = get_3d_meas(mtlb(elec),mtlb(vtx),Vn,Ib,mtlb(no_pl));
 
//! mtlb(end) can be replaced by end() or end whether end is an m-file or not
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() or J whether J is an m-file or not
 
//! mtlb(J) can be replaced by J() or J whether J is an m-file or not
 
//! mtlb(sm1) can be replaced by sm1() or sm1 whether sm1 is an m-file or not
 
//! mtlb(sm1) can be replaced by sm1() or sm1 whether sm1 is an m-file or not
 
//! mtlb(J) can be replaced by J() or J whether J is an m-file or not
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() or srf whether srf is an m-file or not
 
//!! Unknown variable vtx ,the original calling sequence is used
 
//!! Unknown variable vtx ,the original calling sequence is used
 
//!! Unknown variable vtx ,the original calling sequence is used
 
//!! Unknown variable trimesh ,the original calling sequence is used
trimesh(mtlb(srf),vtx(:,1),vtx(:,2),vtx(:,3));
 
//!! Unknown variable 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 variable hidden ,the original calling sequence is used
hidden('off');
mtlb_hold('on');
 
//! mtlb(vtx) can be replaced by vtx() or vtx whether vtx is an m-file or not
 
//! mtlb(%simp) can be replaced by %simp() or %simp whether %simp is an m-file or not
repaint_inho(mat,mat_ref,mtlb(vtx),mtlb(%simp));
 
//!! Unknown variable camlight ,the original calling sequence is used
camlight('lefth');
 
//!! Unknown variable 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() or vtx whether vtx is an m-file or not
 
//! mtlb(%simp) can be replaced by %simp() or %simp whether %simp is an m-file or not
 
//! mtlb(fc) can be replaced by fc() or fc whether fc is an m-file or not
fc = slicer_plot(2.63,sol,mtlb(vtx),mtlb(%simp),mtlb(fc));
 
//!! Unknown variable view ,the original calling sequence is used
view(2);
 
//!! Unknown variable grid ,the original calling sequence is used
grid();
 
//!! Unknown variable colorbar ,the original calling sequence is used
colorbar();
 
//!! Unknown variable 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() or vtx whether vtx is an m-file or not
 
//! mtlb(%simp) can be replaced by %simp() or %simp whether %simp is an m-file or not
fc = slicer_plot(2.1,sol,mtlb(vtx),mtlb(%simp),fc);
 
//!! Unknown variable view ,the original calling sequence is used
view(2);
 
//!! Unknown variable grid ,the original calling sequence is used
grid();
 
//!! Unknown variable colorbar ,the original calling sequence is used
colorbar();
 
//!! Unknown variable 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() or vtx whether vtx is an m-file or not
 
//! mtlb(%simp) can be replaced by %simp() or %simp whether %simp is an m-file or not
fc = slicer_plot(1.72,sol,mtlb(vtx),mtlb(%simp),fc);
 
//!! Unknown variable view ,the original calling sequence is used
view(2);
 
//!! Unknown variable grid ,the original calling sequence is used
grid();
 
//!! Unknown variable colorbar ,the original calling sequence is used
colorbar();
 
//!! Unknown variable 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() or vtx whether vtx is an m-file or not
 
//! mtlb(%simp) can be replaced by %simp() or %simp whether %simp is an m-file or not
fc = slicer_plot(1.1,sol,mtlb(vtx),mtlb(%simp),fc);
 
//!! Unknown variable view ,the original calling sequence is used
view(2);
 
//!! Unknown variable grid ,the original calling sequence is used
grid();
 
//!! Unknown variable colorbar ,the original calling sequence is used
colorbar();
 
//!! Unknown variable 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() or vtx whether vtx is an m-file or not
 
//! mtlb(%simp) can be replaced by %simp() or %simp whether %simp is an m-file or not
fc = slicer_plot(0.83,sol,mtlb(vtx),mtlb(%simp),fc);
 
//!! Unknown variable view ,the original calling sequence is used
view(2);
 
//!! Unknown variable grid ,the original calling sequence is used
grid();
 
//!! Unknown variable colorbar ,the original calling sequence is used
colorbar();
 
//!! Unknown variable 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() or vtx whether vtx is an m-file or not
 
//! mtlb(%simp) can be replaced by %simp() or %simp whether %simp is an m-file or not
fc = slicer_plot(0.1,sol,mtlb(vtx),mtlb(%simp),fc);
 
//!! Unknown variable view ,the original calling sequence is used
view(2);
 
//!! Unknown variable grid ,the original calling sequence is used
grid();
 
//!! Unknown variable colorbar ,the original calling sequence is used
colorbar();
 
//!! Unknown variable 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
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%