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 if 0 % Very old code
0080 trimesh(srf,vtx(:,1),vtx(:,2),vtx(:,3));
0081 axis image;
0082 set(gcf,'Colormap',[0 0 0]);
0083 hidden off;
0084 hold on;
0085 repaint_inho(real(mat),real(mat_ref),vtx,simp);
0086 else
0087    fmdl= eidors_obj('fwd_model','', ...
0088           'nodes',vtx,'elems',simp);
0089    img = eidors_obj('image','', ...
0090           'fwd_model',fmdl,'elem_data',real(mat));
0091    show_fem(img);
0092 end
0093 title('REAL');
0094 camlight left;
0095 lighting flat;
0096 
0097 subplot(1,2,2);
0098 if 0
0099     trimesh(srf,vtx(:,1),vtx(:,2),vtx(:,3));
0100     axis image;
0101     set(gcf,'Colormap',[0 0 0]);
0102     hidden off;
0103     hold on;
0104     repaint_inho(imag(mat),imag(mat_ref),vtx,simp);
0105 else
0106    img = eidors_obj('image','', ...
0107           'fwd_model',fmdl, ...
0108           'elem_data',imag(mat-mat_ref));
0109    show_fem(img);
0110 end
0111  title('IMAGINARY');
0112 camlight left;
0113 lighting flat;
0114 
0115 pause(2);
0116 close;
0117         
0118 %figure; [mat,grp] = set_inho(srf,simp,vtx,mat_ref,0.8-0.6i);
0119 %sB = 0.8 - 0.6i; % A local complex change
0120 sB = 1 + 0.75i; % or a local permittivity change
0121 sprintf('and this one at  %f + %f i',real(sB), imag(sB))
0122 mat(B) = sB;
0123 
0124 figure; subplot(1,2,1); 
0125 if 0
0126     trimesh(srf,vtx(:,1),vtx(:,2),vtx(:,3));
0127     axis image;
0128     set(gcf,'Colormap',[0 0 0]);
0129     hidden off;
0130     hold on;
0131     repaint_inho(real(mat),real(mat_ref),vtx,simp);
0132 else
0133    img = eidors_obj('image','', ...
0134           'fwd_model',fmdl, ...
0135           'elem_data',real(mat-mat_ref));
0136    show_fem(img);
0137 end
0138  title('REAL');
0139 camlight left;
0140 lighting flat;
0141 
0142 subplot(1,2,2); 
0143 if 0
0144     trimesh(srf,vtx(:,1),vtx(:,2),vtx(:,3));
0145     axis image;
0146     set(gcf,'Colormap',[0 0 0]);
0147     hidden off;
0148     hold on;
0149     repaint_inho(imag(mat),imag(mat_ref),vtx,simp);
0150 else
0151    img = eidors_obj('image','', ...
0152           'fwd_model',fmdl, ...
0153           'elem_data',imag(mat-mat_ref));
0154    show_fem(img);
0155 end
0156  title('IMAGINARY');
0157 camlight left;
0158 lighting flat;
0159 
0160 pause(2);
0161 close;
0162         
0163 %[mat,grp] = set_inho(srf,simp,vtx,mat,0.9+0.41i);
0164 disp(sprintf('\n'))
0165 
0166 [En,D,Ela,ppn] = fem_master_full(vtx,simp,mat,gnd_ind,elec,zc,sym);
0167 [Vn] = left_divide(En,I,tol,ppn);
0168 [voltageH,voltageV,indH,indV,dfv]=get_3d_meas(elec,vtx,Vn,Ib,no_pl);
0169 dfv = dfv(1:2:end);
0170 
0171 if size(dfr)~= size(dfv)
0172    error('Mismatched measurements')
0173 end
0174 
0175 disp('Measurements calculated')
0176 disp(sprintf('\n'))
0177 
0178 dva = voltageH-refH;
0179 
0180 disp('Measurements blended with Gaussian noise ...')
0181 dc = mean(real(dva)); %DC component of the noise
0182 dvrG = dc./7 * ones(length(dva),1) + dc * randn(length(dva),1); %Add the AC component
0183 
0184 dc = mean(imag(dva)); %DC component of the noise
0185 dviG = dc./7 * ones(length(dva),1) + dc * randn(length(dva),1); %Add the AC component
0186 
0187 dat = (real(dva) + dvrG) + (imag(dva) + dviG)*i;
0188 %dat = dva;
0189 
0190 disp('Calculating measurement fields')
0191 [v_f] = m_3d_fields(vtx,size(elec,1),indH,Eref,tol,gnd_ind);
0192 
0193 disp('Calculating a single complex jacobian')
0194 [J] = jacobian_3d(I,elec,vtx,simp,gnd_ind,mat_ref,zc,v_f,dfv,tol,sym);
0195 
0196 
0197 disp('Computing a smooting prior')
0198 [Reg] = iso_f_smooth(simp,vtx,1,2);
0199 
0200 disp('Calculating a linearised step inverse solution')
0201 disp(sprintf('\n'))
0202 tfac = 1e-7;
0203 sol = (J.'*J +  tfac*Reg'*Reg)\J.' * dat;
0204 sreal = real(sol);
0205 simag = imag(sol);
0206 
0207 truereal = real(mat-mat_ref);
0208 trueimag = imag(mat-mat_ref);
0209 
0210 v = version;
0211 
0212 if str2num(v(1)) > 5
0213 
0214 h01 = figure;
0215 set(h01,'NumberTitle','off');
0216 set(h01,'Name','True conductivity distribution');
0217 subplot(2,3,1); [fc] = slicer_plot_n(2.63,truereal,vtx,simp); view(2); grid; colorbar; axis off; title('z=2.63'); 
0218 %Calculates also fc. Just once!
0219 subplot(2,3,2); [fc] = slicer_plot_n(2.10,truereal,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=2.10'); 
0220 subplot(2,3,3); [fc] = slicer_plot_n(1.72,truereal,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=1.72'); 
0221 subplot(2,3,4); [fc] = slicer_plot_n(1.10,truereal,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=1.10'); 
0222 subplot(2,3,5); [fc] = slicer_plot_n(0.83,truereal,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=0.83');
0223 subplot(2,3,6); [fc] = slicer_plot_n(0.10,truereal,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=0.10');
0224 
0225 
0226 h02 = figure;
0227 set(h02,'NumberTitle','off');
0228 set(h02,'Name','True scaled permittivity distribution');
0229 subplot(2,3,1); [fc] = slicer_plot_n(2.63,trueimag,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=2.60'); 
0230 subplot(2,3,2); [fc] = slicer_plot_n(2.10,trueimag,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=2.10'); 
0231 subplot(2,3,3); [fc] = slicer_plot_n(1.72,trueimag,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=1.70'); 
0232 subplot(2,3,4); [fc] = slicer_plot_n(1.10,trueimag,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=1.10'); 
0233 subplot(2,3,5); [fc] = slicer_plot_n(0.83,trueimag,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=0.80');
0234 subplot(2,3,6); [fc] = slicer_plot_n(0.10,trueimag,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=0.10');
0235 
0236 
0237 
0238 h2 = figure;
0239 set(h2,'NumberTitle','off');
0240 set(h2,'Name','Reconstructed conductivity distribution');
0241 subplot(2,3,1); [fc] = slicer_plot_n(2.63,sreal,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=2.63'); 
0242 subplot(2,3,2); [fc] = slicer_plot_n(2.10,sreal,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=2.10'); 
0243 subplot(2,3,3); [fc] = slicer_plot_n(1.72,sreal,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=1.72'); 
0244 subplot(2,3,4); [fc] = slicer_plot_n(1.10,sreal,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=1.10'); 
0245 subplot(2,3,5); [fc] = slicer_plot_n(0.83,sreal,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=0.83');
0246 subplot(2,3,6); [fc] = slicer_plot_n(0.10,sreal,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=0.10');
0247 
0248 h3 = figure;
0249 set(h3,'NumberTitle','off');
0250 set(h3,'Name','Reconstructed scaled permittivity distribution');
0251 subplot(2,3,1); [fc] = slicer_plot_n(2.63,simag,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=2.60'); 
0252 subplot(2,3,2); [fc] = slicer_plot_n(2.10,simag,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=2.10'); 
0253 subplot(2,3,3); [fc] = slicer_plot_n(1.72,simag,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=1.70'); 
0254 subplot(2,3,4); [fc] = slicer_plot_n(1.10,simag,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=1.10'); 
0255 subplot(2,3,5); [fc] = slicer_plot_n(0.83,simag,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=0.80');
0256 subplot(2,3,6); [fc] = slicer_plot_n(0.10,simag,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=0.10');
0257 
0258 else
0259    
0260    disp('Change the plotting command above to slicer_plot or upgrate to MATLAB 6')
0261    disp('See demo_real for more details')
0262    
0263 end
0264 
0265 disp('Done')
0266 
0267 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0268 % This is part of the EIDORS suite.
0269 % Copyright (c) N. Polydorides 2003
0270 % Copying permitted under terms of GNU GPL
0271 % See enclosed file gpl.html for details.
0272 % EIDORS 3D version 2.0
0273 % MATLAB version 6.1 R12
0274 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%