demo_real.m~

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

clear; 
clc;
disp('This is a demo for reconstructing real impedance changes')
disp(sprintf('\n'))

load datareal srf vtx simp;

trimesh(srf,vtx(:,1),vtx(:,2),vtx(:,3));
axis image;
set(gcf,'Colormap',[0 0 0]);
hold on;

disp('This is a cylindrical mesh, wait to attach the electrodes')
disp(sprintf('\n'))

pause(4);

load datacom sels;

  for u=1:size(sels)
      paint_electrodes(sels(u),srf,vtx);
  end
  
hidden off;
load datareal gnd_ind elec zc sym;

%Change sym to '{n}' for Cholesky instead of PCG in the forward solver
sym = '{n}';

load datareal protocol no_pl;

[I,Ib] = set_3d_currents(protocol,elec,vtx,gnd_ind,no_pl);

clc;
disp('Adjacent current patterns selected. Calculating reference measurements')
disp(sprintf('\n'))

pause(2);

mat_ref = ones(828,1); %%%%%%

%Set the tolerance for the forward solver
tol = 1e-4;

[Eref,ppr] = fem_master_full(vtx,simp,mat_ref,gnd_ind,elec,zc,sym);
[Vref] = forward_solver(vtx,Eref,I,tol,ppr);
[refH,refV,indH,indV,dfr]=get_3d_meas(elec,vtx,Vref,Ib,no_pl);
dfr = dfr(1:2:end); %Taking just the horrizontal measurements

close;
clc;
disp('Now let us create a couple of changes ...')
disp(sprintf('\n'))

pause(2);


mat=mat_ref;

load datareal A B
%figure; [mat,grp] = set_inho(srf,simp,vtx,mat_ref,1.1); 
sA = 1.2;
sprintf ('This one  %f',sA)
mat(A) = sA;

figure; 
trimesh(srf,vtx(:,1),vtx(:,2),vtx(:,3));
axis image;
set(gcf,'Colormap',[0 0 0]);
hidden off;
hold on;
repaint_inho(mat,mat_ref,vtx,simp); 
camlight lefth;
lighting flat;
drawnow;

pause(4);
close;
        
sB = 0.9;
mat(B) = sB;

figure; 
trimesh(srf,vtx(:,1),vtx(:,2),vtx(:,3));
axis image;
set(gcf,'Colormap',[0 0 0]);
hidden off;
hold on;
repaint_inho(mat,mat_ref,vtx,simp); 
camlight lefth;
lighting flat;

sprintf('And this one at  %f',sB)
disp(sprintf('\n'))

drawnow;
pause(4);
close;
        
%[mat,grp] = set_inho(srf,simp,vtx,mat,0.999); 
disp('Simulating measurements based on the CEM ...')


[En,ppn] = fem_master_full(vtx,simp,mat,gnd_ind,elec,zc,sym);
[Vn] = forward_solver(vtx,En,I,tol,ppn);
[voltageH,voltageV,indH,indV,dfv]=get_3d_meas(elec,vtx,Vn,Ib,no_pl);
dfv = dfv(1:2:end);

if size(dfr)~= size(dfv)
   error('Mismatched measurements')
end


dva = refH-voltageH;

disp('Measurements infused with Gaussian noise ...')


dc = mean(dva); %DC component of the noise
noi = dc./7 * ones(length(dva),1) + dc * randn(length(dva),1); %Add the AC component
dvaG = dva + noi;


disp('Retreiving the Jacobian and regularisation matrix ...')

load datareal J
load datareal sm1 %The smoothing prior 


% Or otherwise 
disp('Calculating the Jacobian')
disp('Please wait, this may take some time ..')
disp(sprintf('\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(sprintf('\n'))

tfac = 1e-8;

sol = (J'*J +  tfac*sm1'*sm1)\J' * dvaG;

load datareal fc;

h1 = figure;
set(h1,'NumberTitle','off');
set(h1,'Name','Simulated inhomogeneities');
trimesh(srf,vtx(:,1),vtx(:,2),vtx(:,3));
axis image;
set(gcf,'Colormap',[0 0 0]);
hidden off;
hold on;
repaint_inho(mat,mat_ref,vtx,simp); 
camlight lefth;
lighting flat;

h2 = figure;
set(h2,'NumberTitle','off');
set(h2,'Name','Reconstructed conductivity distribution');
disp('And its reconstructed image, planes from top to bottom')
subplot(2,3,1); [fc] = slicer_plot(2.63,sol,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=2.60'); 
subplot(2,3,2); [fc] = slicer_plot(2.10,sol,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=2.10'); 
subplot(2,3,3); [fc] = slicer_plot(1.72,sol,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=1.70'); 
subplot(2,3,4); [fc] = slicer_plot(1.10,sol,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=1.10'); 
subplot(2,3,5); [fc] = slicer_plot(0.83,sol,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=0.80');
subplot(2,3,6); [fc] = slicer_plot(0.10,sol,vtx,simp,fc); view(2); grid; colorbar; axis off; title('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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%