demo_complex

This demo function shows how the EIT problem can be formulated in a complex
form. The complex measurements along with a single complex Jacobian are then 
fed into the reconstruction process. This is a different forulation from the 
demo_comp.m function.

0001 %This demo function shows how the EIT problem can be formulated in a complex
0002 %form. The complex measurements along with a single complex Jacobian are then
0003 %fed into the reconstruction process. This is a different forulation from the
0004 %demo_comp.m function.
0005 
0006 warning off;
0007 disp('This is a demo for reconstructing admittivity changes of the form a + bi')
0008 
0009 load datacom srf vtx simp;
0010 %srf : the boundary surfaces (triangles)
0011 %vtx : the vertices of the model (co-ordinates of the nodes)
0012 %simp: the simplices of the model (connectivity in tetrahedral)
0013 
0014 trimesh(srf,vtx(:,1),vtx(:,2),vtx(:,3));
0015 axis image;
0016 set(gcf,'Colormap',[0 0 0]);
0017 hold on;
0018 
0019 disp('This is a cylindrical mesh with homogeneous distribution 1 + 0.5i')
0020 disp('Wait to attach the electrodes')
0021 disp(sprintf('\n'))
0022 
0023 pause(2);
0024 
0025 load datacom sels;
0026 % FIXME: this needs to use the new show_fem functions
0027 %sels :Index in srf matrix denoting the faces to be assigned as electrodes
0028 % for u=1:size(sels)
0029 %     paint_electrodes(sels(u),srf,vtx);
0030 % end
0031 hidden off;
0032   
0033 load datacom gnd_ind elec no_pl protocol zc sym;
0034 %gnd_ind : The ground index here a node, can also be an electrode
0035 %elec : The electrodes matrix.
0036 %np_pl : Number of electrode planes (in planar arrangements)
0037 %protocol : Adjacent or Opposite or Customized.
0038 %zc : Contact impedances of the electrodes
0039 %sym : Boolean entry for efficient forward computations
0040 %Direct solvers :'{y}' / Iterative : '{n}'
0041 %sym='{y}';
0042 
0043 [I,Ib] = set_3d_currents(protocol,elec,vtx,gnd_ind,no_pl);
0044 
0045 disp('Adjacent current patterns selected')
0046 disp('Calculating reference measurements')
0047 disp(sprintf('\n'))
0048 
0049 %Set the tolerance for the forward solver
0050 tol = 1e-5;
0051 
0052 pause(2);
0053 
0054 mat_ref = (1+0.5i)*ones(828,1); %%%%%%
0055 %Jacobians will be calculated based on this
0056 
0057 [Eref,D,Ela,ppr] = fem_master_full(vtx,simp,mat_ref,gnd_ind,elec,zc,sym);
0058 [Vref] = left_divide(Eref,I,tol,ppr);
0059 [refH,refV,indH,indV,dfr]=get_3d_meas(elec,vtx,Vref,Ib,no_pl);
0060 dfr = dfr(1:2:end); %Taking just the horrizontal measurements
0061 
0062 close;
0063 disp('Allow a couple of complex changes ...')
0064 disp(sprintf('\n'))
0065 pause(2);
0066 
0067 
0068 mat=mat_ref;
0069 
0070 load datacom A B %Indices of the elements to represent the inhomogeneities
0071 %figure; [mat,grp] = set_inho(srf,simp,vtx,mat_ref,1.2-0.4i);
0072 sA = 1.2 + 0.4i; % A local complex change or
0073 %sA = 1.2 + 0.5i; %just a local conductivity change
0074 sprintf ('This one  at %f + %f i',real(sA), imag(sA))
0075 mat(A) = sA;
0076 
0077 figure; 
0078 subplot(1,2,1);
0079 trimesh(srf,vtx(:,1),vtx(:,2),vtx(:,3));
0080 axis image;
0081 set(gcf,'Colormap',[0 0 0]);
0082 hidden off;
0083 hold on;
0084 repaint_inho(real(mat),real(mat_ref),vtx,simp); title('REAL');
0085 camlight left;
0086 lighting flat;
0087 
0088 subplot(1,2,2);
0089 trimesh(srf,vtx(:,1),vtx(:,2),vtx(:,3));
0090 axis image;
0091 set(gcf,'Colormap',[0 0 0]);
0092 hidden off;
0093 hold on;
0094 repaint_inho(imag(mat),imag(mat_ref),vtx,simp); title('IMAGINARY');
0095 camlight left;
0096 lighting flat;
0097 
0098 pause(2);
0099 close;
0100         
0101 %figure; [mat,grp] = set_inho(srf,simp,vtx,mat_ref,0.8-0.6i);
0102 %sB = 0.8 - 0.6i; % A local complex change
0103 sB = 1 + 0.75i; % or a local permittivity change
0104 sprintf('and this one at  %f + %f i',real(sB), imag(sB))
0105 mat(B) = sB;
0106 
0107 figure; subplot(1,2,1); 
0108 trimesh(srf,vtx(:,1),vtx(:,2),vtx(:,3));
0109 axis image;
0110 set(gcf,'Colormap',[0 0 0]);
0111 hidden off;
0112 hold on;
0113 repaint_inho(real(mat),real(mat_ref),vtx,simp); title('REAL');
0114 camlight left;
0115 lighting flat;
0116 
0117 subplot(1,2,2); 
0118 trimesh(srf,vtx(:,1),vtx(:,2),vtx(:,3));
0119 axis image;
0120 set(gcf,'Colormap',[0 0 0]);
0121 hidden off;
0122 hold on;
0123 repaint_inho(imag(mat),imag(mat_ref),vtx,simp); title('IMAGINARY');
0124 camlight left;
0125 lighting flat;
0126 
0127 pause(2);
0128 close;
0129         
0130 %[mat,grp] = set_inho(srf,simp,vtx,mat,0.9+0.41i);
0131 disp(sprintf('\n'))
0132 
0133 [En,D,Ela,ppn] = fem_master_full(vtx,simp,mat,gnd_ind,elec,zc,sym);
0134 [Vn] = left_divide(En,I,tol,ppn);
0135 [voltageH,voltageV,indH,indV,dfv]=get_3d_meas(elec,vtx,Vn,Ib,no_pl);
0136 dfv = dfv(1:2:end);
0137 
0138 if size(dfr)~= size(dfv)
0139    error('Mismatched measurements')
0140 end
0141 
0142 disp('Measurements calculated')
0143 disp(sprintf('\n'))
0144 
0145 dva = voltageH-refH;
0146 
0147 disp('Measurements blended with Gaussian noise ...')
0148 dc = mean(real(dva)); %DC component of the noise
0149 dvrG = dc./7 * ones(length(dva),1) + dc * randn(length(dva),1); %Add the AC component
0150 
0151 dc = mean(imag(dva)); %DC component of the noise
0152 dviG = dc./7 * ones(length(dva),1) + dc * randn(length(dva),1); %Add the AC component
0153 
0154 dat = (real(dva) + dvrG) + (imag(dva) + dviG)*i;
0155 %dat = dva;
0156 
0157 disp('Calculating measurement fields')
0158 [v_f] = m_3d_fields(vtx,size(elec,1),indH,Eref,tol,gnd_ind);
0159 
0160 disp('Calculating a single complex jacobian')
0161 [J] = jacobian_3d(I,elec,vtx,simp,gnd_ind,mat_ref,zc,v_f,dfv,tol,sym);
0162 
0163 
0164 disp('Computing a smooting prior')
0165 [Reg] = iso_f_smooth(simp,vtx,1,2);
0166 
0167 disp('Calculating a linearised step inverse solution')
0168 disp(sprintf('\n'))
0169 tfac = 1e-7;
0170 sol = (J.'*J +  tfac*Reg'*Reg)\J.' * dat;
0171 sreal = real(sol);
0172 simag = imag(sol);
0173 
0174 truereal = real(mat-mat_ref);
0175 trueimag = imag(mat-mat_ref);
0176 
0177 v = version;
0178 
0179 if str2num(v(1)) > 5
0180 
0181 h01 = figure;
0182 set(h01,'NumberTitle','off');
0183 set(h01,'Name','True conductivity distribution');
0184 subplot(2,3,1); [fc] = slicer_plot_n(2.63,truereal,vtx,simp); view(2); grid; colorbar; axis off; title('z=2.63'); 
0185 %Calculates also fc. Just once!
0186 subplot(2,3,2); [fc] = slicer_plot_n(2.10,truereal,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=2.10'); 
0187 subplot(2,3,3); [fc] = slicer_plot_n(1.72,truereal,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=1.72'); 
0188 subplot(2,3,4); [fc] = slicer_plot_n(1.10,truereal,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=1.10'); 
0189 subplot(2,3,5); [fc] = slicer_plot_n(0.83,truereal,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=0.83');
0190 subplot(2,3,6); [fc] = slicer_plot_n(0.10,truereal,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=0.10');
0191 
0192 
0193 h02 = figure;
0194 set(h02,'NumberTitle','off');
0195 set(h02,'Name','True scaled permittivity distribution');
0196 subplot(2,3,1); [fc] = slicer_plot_n(2.63,trueimag,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=2.60'); 
0197 subplot(2,3,2); [fc] = slicer_plot_n(2.10,trueimag,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=2.10'); 
0198 subplot(2,3,3); [fc] = slicer_plot_n(1.72,trueimag,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=1.70'); 
0199 subplot(2,3,4); [fc] = slicer_plot_n(1.10,trueimag,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=1.10'); 
0200 subplot(2,3,5); [fc] = slicer_plot_n(0.83,trueimag,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=0.80');
0201 subplot(2,3,6); [fc] = slicer_plot_n(0.10,trueimag,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=0.10');
0202 
0203 
0204 
0205 h2 = figure;
0206 set(h2,'NumberTitle','off');
0207 set(h2,'Name','Reconstructed conductivity distribution');
0208 subplot(2,3,1); [fc] = slicer_plot_n(2.63,sreal,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=2.63'); 
0209 subplot(2,3,2); [fc] = slicer_plot_n(2.10,sreal,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=2.10'); 
0210 subplot(2,3,3); [fc] = slicer_plot_n(1.72,sreal,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=1.72'); 
0211 subplot(2,3,4); [fc] = slicer_plot_n(1.10,sreal,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=1.10'); 
0212 subplot(2,3,5); [fc] = slicer_plot_n(0.83,sreal,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=0.83');
0213 subplot(2,3,6); [fc] = slicer_plot_n(0.10,sreal,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=0.10');
0214 
0215 h3 = figure;
0216 set(h3,'NumberTitle','off');
0217 set(h3,'Name','Reconstructed scaled permittivity distribution');
0218 subplot(2,3,1); [fc] = slicer_plot_n(2.63,simag,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=2.60'); 
0219 subplot(2,3,2); [fc] = slicer_plot_n(2.10,simag,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=2.10'); 
0220 subplot(2,3,3); [fc] = slicer_plot_n(1.72,simag,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=1.70'); 
0221 subplot(2,3,4); [fc] = slicer_plot_n(1.10,simag,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=1.10'); 
0222 subplot(2,3,5); [fc] = slicer_plot_n(0.83,simag,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=0.80');
0223 subplot(2,3,6); [fc] = slicer_plot_n(0.10,simag,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=0.10');
0224 
0225 else
0226    
0227    disp('Change the plotting command above to slicer_plot or upgrate to MATLAB 6')
0228    disp('See demo_real for more details')
0229    
0230 end
0231 
0232 disp('Done')
0233 
0234 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0235 % This is part of the EIDORS suite.
0236 % Copyright (c) N. Polydorides 2003
0237 % Copying permitted under terms of GNU GPL
0238 % See enclosed file gpl.html for details.
0239 % EIDORS 3D version 2.0
0240 % MATLAB version 6.1 R12
0241 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%