demo_comp.m~

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

clear;
clc;
disp('This is a demo for reconstructing impedance changes of the form a+bi')
disp('********************************************************************')
disp(sprintf('\n'))

load datacom 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(2);

load datacom sels;

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

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

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

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

pause(3);

mat_ref = (1+0.5*sqrt(-1))*ones(828,1); %%%%%% 
%Precalculated Jacobians are based on this.

[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 complex changes ...')
disp(sprintf('\n'))

pause(2);


mat=mat_ref;

load datacom A B
%figure; [mat,grp] = set_inho(srf,simp,vtx,mat_ref,1.2+0.4i); 
sA = 1.2 + 0.4i; % A local complex change or
%sA = 1.2 + 0.5i; %just a local conductivity change
sprintf ('This one  %f + %f i',real(sA), imag(sA))
mat(A) = sA;

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

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

pause(2);
close;
        
%figure; [mat,grp] = set_inho(srf,simp,vtx,mat_ref,0.8+0.6i);       
%sB = 0.8 + 0.6i; % A local complex change
sB = 1 + 0.3i; % or a local permittivity change
sprintf('And this one at  %f + %f i',real(sB), imag(sB))
mat(B) = sB;

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

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

pause(2);
close;
        
%[mat,grp] = set_inho(srf,simp,vtx,mat,0.9+0.41i); 
disp(sprintf('\n'))

[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


disp('Measurements calculated')
disp(sprintf('\n'))

pause(1);

dva = refH-voltageH;
dvr=real(dva);
dvi=imag(dva);

disp('Blended with Gaussian noise ...')
dc = mean(dvr); %DC component of the noise
noi = dc./7 * ones(length(dvr),1) + dc * randn(length(dvr),1); %Add the AC component
dvrG = dvr + noi;

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


%disp('Retreiving the Jacobians ...')
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')

[sm1] = iso_f_smooth(simp,2);

disp('Calculating the inverse solution, this may take some time')
disp(sprintf('\n'))

tfac = 1e-7;

MR = [sm1, zeros(size(sm1)); zeros(size(sm1)), sm1];

MJ = [Jrr Jri; Jir -Jii];

dat = [dvrG; dviG];

%load datacom MJ MR

sol = (MJ'*MJ +  tfac*MR)\MJ' * dat;

sreal = sol(1:size(simp,1));
simag = sol(size(simp,1)+1:end);

load datacom fc;

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

subplot(1,2,2); trimesh(srf,vtx(:,1),vtx(:,2),vtx(:,3));
axis image;
set(gcf,'Colormap',[0 0 0]);
hidden off;
hold on;
repaint_inho(imag(mat),imag(mat_ref),vtx,simp); title('Original permittivity distribution');
camlight left;
lighting flat;


disp('And their reconstructed images')
h2 = figure;
set(h2,'NumberTitle','off');
set(h2,'Name','Reconstructed conductivity distribution');
subplot(2,3,1); [fc] = slicer_plot(2.63,sreal,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=2.60'); 
subplot(2,3,2); [fc] = slicer_plot(2.10,sreal,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=2.10'); 
subplot(2,3,3); [fc] = slicer_plot(1.72,sreal,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=1.70'); 
subplot(2,3,4); [fc] = slicer_plot(1.10,sreal,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=1.10'); 
subplot(2,3,5); [fc] = slicer_plot(0.83,sreal,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=0.80');
subplot(2,3,6); [fc] = slicer_plot(0.10,sreal,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=0.10');

h3 = figure;
set(h3,'NumberTitle','off');
set(h3,'Name','Reconstructed scaled permittivity distribution');
subplot(2,3,1); [fc] = slicer_plot(2.63,simag,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=2.60'); 
subplot(2,3,2); [fc] = slicer_plot(2.10,simag,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=2.10'); 
subplot(2,3,3); [fc] = slicer_plot(1.72,simag,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=1.70'); 
subplot(2,3,4); [fc] = slicer_plot(1.10,simag,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=1.10'); 
subplot(2,3,5); [fc] = slicer_plot(0.83,simag,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=0.80');
subplot(2,3,6); [fc] = slicer_plot(0.10,simag,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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%