demo_comp.sci

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

function []=demo_comp()
clear;
 
//! clc ignored
// clc
 ;
disp('This is a demo for reconstructing impedance changes of the form a+bi');
disp('********************************************************************');
disp(mtlb_sprintf('%s','\n'));
// 
mtlb_load('datacom','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*(2));
// 
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('datacom','gnd_ind','elec','no_pl','protocol','zc','sym');
// 
 
//! 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'));
// 
//Set the tolerance for the forward solver
tol = 0.000001;
// 
xpause(1000*(3));
// 
mat_ref = (1+0.5*sqrt(-1))*ones(828,1);
//%%%%% 
//Precalculated Jacobians are based on this.
// 
 
//! 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
 
//! mtlb(sym) can be replaced by sym() or sym whether sym is an m-file or not
[Eref,ppr] = fem_master_full(mtlb(vtx),mtlb(%simp),mat_ref,mtlb(gnd_ind),mtlb(elec),mtlb(zc),mtlb(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 complex changes ...');
disp(mtlb_sprintf('%s','\n'));
// 
xpause(1000*(2));
// 
// 
mat = mat_ref;
// 
mtlb_load('datacom','A','B');
//figure; [mat,grp] = set_inho(srf,simp,vtx,mat_ref,1.2+0.4i); 
 
//! mtlb(%i) can be replaced by %i() or %i whether %i is an m-file or not
sA = 1.2+0.4*mtlb(%i);
// A local complex change or
//sA = 1.2 + 0.5i; %just a local conductivity change
mtlb_sprintf('This one  %f + %f i',real(sA),imag(sA));
 
//! mtlb(A) can be replaced by A() or A whether A is an m-file or not
mat(mtlb(A),1) = sA(:);
// 
xset('window',max(winsid()+1))
max(winsid());
%v$2 = int((1-1)/2);%v$1 = 1-1-2*%v$2
xsetech([%v$1/2,%v$2/1,1/2,1/1]);
 
//! 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(real(mat),real(mat_ref),mtlb(vtx),mtlb(%simp));
xtitle('REAL',' ',' ');
 
//!! Unknown variable camlight ,the original calling sequence is used
camlight('lefth');
 
//!! Unknown variable lighting ,the original calling sequence is used
lighting('flat');
// 
%v$2 = int((2-1)/2);%v$1 = 2-1-2*%v$2
xsetech([%v$1/2,%v$2/1,1/2,1/1]);
 
//! 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(imag(mat),imag(mat_ref),mtlb(vtx),mtlb(%simp));
xtitle('IMAGINARY',' ',' ');
 
//!! Unknown variable camlight ,the original calling sequence is used
camlight('lefth');
 
//!! Unknown variable lighting ,the original calling sequence is used
lighting('flat');
// 
xpause(1000*(2));
xdel();
// 
//figure; [mat,grp] = set_inho(srf,simp,vtx,mat_ref,0.8+0.6i);       
//sB = 0.8 + 0.6i; % A local complex change
 
//! mtlb(%i) can be replaced by %i() or %i whether %i is an m-file or not
sB = 1+0.3*mtlb(%i);
// or a local permittivity change
mtlb_sprintf('And this one at  %f + %f i',real(sB),imag(sB));
 
//! mtlb(B) can be replaced by B() or B whether B is an m-file or not
mat(mtlb(B),1) = sB(:);
// 
xset('window',max(winsid()+1))
max(winsid());
%v$2 = int((1-1)/2);%v$1 = 1-1-2*%v$2
xsetech([%v$1/2,%v$2/1,1/2,1/1]);
 
//! 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(real(mat),real(mat_ref),mtlb(vtx),mtlb(%simp));
xtitle('REAL',' ',' ');
 
//!! Unknown variable camlight ,the original calling sequence is used
camlight('lefth');
 
//!! Unknown variable lighting ,the original calling sequence is used
lighting('flat');
// 
%v$2 = int((2-1)/2);%v$1 = 2-1-2*%v$2
xsetech([%v$1/2,%v$2/1,1/2,1/1]);
 
//! 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(imag(mat),imag(mat_ref),mtlb(vtx),mtlb(%simp));
xtitle('IMAGINARY',' ',' ');
 
//!! Unknown variable camlight ,the original calling sequence is used
camlight('lefth');
 
//!! Unknown variable lighting ,the original calling sequence is used
lighting('flat');
// 
xpause(1000*(2));
xdel();
// 
//[mat,grp] = set_inho(srf,simp,vtx,mat,0.9+0.41i); 
disp(mtlb_sprintf('%s','\n'));
// 
 
//! 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
 
//! mtlb(sym) can be replaced by sym() or sym whether sym is an m-file or not
[En,ppn] = fem_master_full(mtlb(vtx),mtlb(%simp),mat,mtlb(gnd_ind),mtlb(elec),mtlb(zc),mtlb(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
// 
// 
disp('Measurements calculated');
disp(mtlb_sprintf('%s','\n'));
// 
xpause(1000*(1));
// 
dva = refH-voltageH;
dvr = real(dva);
dvi = imag(dva);
// 
disp('Blended with Gaussian noise ...');
 
//! mtlb_mean(dvr) may be replaced by 
//!   mean(dvr) if dvris a vector
//!   mean(dvr,1) if dvris a matrix
dc = mtlb_mean(dvr);
//DC component of the noise
noi = dc ./ 7*ones(max(size(dvr)),1)+dc*rand(max(size(dvr)),1,'n');
//Add the AC component
dvrG = dvr+noi;
// 
 
//! mtlb_mean(dvi) may be replaced by 
//!   mean(dvi) if dviis a vector
//!   mean(dvi,1) if dviis a matrix
dc = mtlb_mean(dvi);
//DC component of the noise
noi = dc ./ 7*ones(max(size(dvi)),1)+dc*rand(max(size(dvi)),1,'n');
//Add the AC component
dviG = dvi+noi;
// 
// 
//disp('Retreiving the Jacobians ...')
mtlb_load('datacom','Jrr','Jri','Jir','Jii');
disp('Calculating Jacobians');
disp('Please wait, this may take some time ...');
//[IntGrad] = integrofgrad(vtx,simp,mat_ref);
//[v_f] = m_3d_fields(vtx,size(elec,1),indH,Eref,tol,gnd_ind);
// 
// 
//[Jrr,Jri,Jir,Jii] = ....
//             jacobian_3d_comp(I,elec,vtx,simp,gnd_ind,mat_ref,2,zc,IntGrad,v_f,dfr,tol,sym);
disp('Jacobians caclulated');
// 
disp('Computing priors');
// 
 
//! mtlb(%simp) can be replaced by %simp() or %simp whether %simp is an m-file or not
sm1 = iso_f_smooth(mtlb(%simp),2);
// 
disp('Calculating the inverse solution, this may take some time');
disp(mtlb_sprintf('%s','\n'));