EIDORS: Electrical Impedance Tomography and Diffuse Optical Tomography Reconstruction Software 
EIDORS
(mirror) Main Documentation Tutorials Download Contrib Data GREIT Browse Docs Browse SVN News Mailing list (archive) FAQ Developer

EIDORS3D tutorialThis tutorial was created by Nick Polydorides, David Stephenson and Bill Lionheart for use with EIDORS3D (EIDORS version 2). This code will still function with EIDORS version 3.x, but this style of code is now deprecated.%Proposed solutions for the VCIPT Eidors "Hands On" Session %Banff, Canada: September 1st 2003 %Original Solutions: Nick Polydorides Dept. of Mathematics  UMIST %Edited by: David Stephenson Dept. of Chemical Engineering  UMIST %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Introduction: Setting Up Matlab And Run Demo's. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %These demo loads mesh geometry, attaches electrodes, sets the %background conductivity, sets the current injection protocol, adds an %inhomogeneity, calculates the forward and inverse solutions and %visualises them. %View the demonstration routine for real conductivity reconstruction by %copying the black text below into the Matlab command window and pressing %return. demo_real; %And for complex admittivity reconstruction, demo_complex; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Exercise 1: Loading mesh geometry, generate surface and plot. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Load the example mesh geometry (vertices matrix and simplices matrix) provided by the datareal MAT file, load datareal vtx simp; %Generate the outer surfaces of the mesh, [srf] = find_boundary(simp); %Now plot the outer surfaces, trimesh(srf,vtx(:,1),vtx(:,2),vtx(:,3)); colormap([0 0 0]); daspect([1 1 1]); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Exercise 2: Set electrode configuration. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Load the example electrode configuiration provided (then move to "Save") load datareal elec sels %  %  %  %Alternatively you may select your own electrode faces using, %  [elec_face,sels,cnts,VV] = set_electrodes_n(vtx,srf,[],[],[],[30 60]); %  %to select the first electrode face and then, %  [elec_face,sels,cnts,VV] = set_electrodes(vtx,srf,elec_face,sels,cnts,VV); %  %until you select all the electrode faces. %  %  %Once you have gathered the faces to be assigned as electrodes, %  %provided that you have the same N number of faces per electrode and these are indexed in %  %sequence inside "elec_face", then use something like %  elec = reshape(elec_face',(N*3),number of electrodes)'; %  %  %i.e for 2 planes of 8 electrodes, of 4 faces each use %  elec = reshape(elec_face',(4*3),16)'; % v %Save workspace as %name_work elec sels %i.e saving workspace with the new variables "elec" and "sels" defined %sels is the array holding the indices of the faces selected in %elec within the boundary surfaces matrix srf %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Exercise 3: Select the ground node (or electrode). %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %To select the centre node on the top surface of the cylinder as the ground point use, load datareal gnd_ind %(This loads the grounding node used by the datareal example) %Plot the position of the ground node trimesh(srf,vtx(:,1),vtx(:,2),vtx(:,3)); colormap([0 0 0]); daspect([1 1 1]); hold on; plot3(vtx(gnd_ind,1),vtx(gnd_ind,2),vtx(gnd_ind,3),'.', ... 'MarkerSize',20,'Color','g'); %Alternatively setting the second electrode elec(2,:) to ground gnd_ind = elec(2,:) ; %Save workspace as %name_work gnd_ind %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Exercise 4: Set background conductivity. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %For a homogeneous isotropic background (real) conductivity of 1 S/m mat_ref = ones(size(simp,1),1); %Save workspace as %name_work mat_ref %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Exercise 5: Set current injection pattern. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %For some "adjacent" current patterns use the function below [I,Ib]=set_3d_currents('{ad}',elec,vtx,gnd_ind,2); %Alternatively, for some custom made patterns use Ibt = ones(1,31); Ibr =  eye(31,31); Ib = [Ibt;Ibr]; I = [zeros(size(vtx,1),size(Ib,2)); Ib]; %And check that Sum I = 0 (current in = current out) by sum(I) %Save workspace as %name_work I Ib %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Exercise 6: Enter contact impedance, compute the system matrix, apply %reference conditions. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Enter a contact impedance vector, say of 20 Ohms zc = 20* ones(size(elec,1),1); %Compute the E_f matrix, the systrem matrix for %the complete electrode model [Ef] = bld_master_full(vtx,simp,mat_ref,elec,zc); %Check the size of E_f, and its sparsity size(Ef) figure(1); spy(Ef); %Then apply the reference conditions [Er] = ref_master(Ef,vtx,gnd_ind,0); %for the grounding node gnd_ind or % [Er] = ref_master(Ef,vtx,gnd_ind,1); %for the grounding electrode gnd_ind %and then compare numerical ranks before rank(full(Ef)) %and after grounding rank(full(Er)) %Er is also possitive definite, check its eigenvalues/singular values eigs(Er) svds(Ef) %For the gradient operator D get a percentage measure of its %sparsity as ((size(D,1)*size(D,2))  nnz(D)) / (size(D,1)*size(D,2))*100 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Exercise 7: The reference forward solution. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Recompute the system matrix based on mat_ref, call this matrix Eref [Eref,D,Ela,ppr] = fem_master_full(vtx,simp,mat_ref,gnd_ind,elec,zc,'{n}'); %and calculte the reference forward solution Vref [Vref] = forward_solver(Eref,I,1e5,ppr); %or use Vref = Eref\I; %Then animate the forward solution % potplot(vtx,srf,Vref); %Save workspace as %name_work zc Eref D Ela Vref %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Exercise 8: Calculate reference measurements, check nonlinearity and %place inhomogeneity. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Calculate some reference measurements refH frm Vref [refH,refV,indH,indV,dfr]=get_3d_meas(elec,vtx,Vref,Ib,2); %and their index dfj = dfr(1:2:end); %Check the nonlinearity of the problem (homogeneous isotropic background %conductivity now set to 3 S/m) mat_ref3 = 3*ones(size(simp,1),1); [Eref3,D,Ela,ppr3] = fem_master_full(vtx,simp,mat_ref3,gnd_ind,elec,zc,'{n}'); [Vref3] = forward_solver(Eref3,I,1e5,ppr); [refH3,refV3,indH,indV,dfr3]=get_3d_meas(elec,vtx,Vref3,Ib,2); %Then plot the difference in the data sets figure(1); plot(refH  refH3); %or alternatively use figure(1); plot(refH) hold on; plot(refH3,'r') %Create target inhomogeneity of magnitude 1.2 (remember that mat_ref = 1) %[mat,grp] = set_inho(srf,simp,vtx,mat_ref,1.2); mat = mat_ref; mat(A) = 1.2; %and maybe a second one as %[mat,grp] = set_inho(srf,simp,vtx,mat,0.85); mat(B) = 0.85; %The target conductivity distribution is now mat %Calculating its corresponding measurements voltageH [E,D,Ela,pp] = fem_master_full(vtx,simp,mat,gnd_ind,elec,zc,'{n}'); [V] = forward_solver(E,I,1e5,pp); [voltageH,voltageV,indH,indV,df]=get_3d_meas(elec,vtx,V,Ib,2); %Plot the perturbation on the data dv due to the introduced inhomogeneities dv = voltageH  refH; figure(1); plot(dv); %Save workspace as %name_work E D Ela refH voltageH dfj %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Exercise 9: Calculate measurement fields and compute the Jacobian. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %The measurement fields v_f for the measurement patterns indH [v_f] = m_3d_fields(vtx,size(elec,1),indH,Eref,1e5,gnd_ind); %Determine the integral of the gradients matrix for first order tetrahedral elements [IntGrad] = integrofgrad(vtx,simp,mat_ref); %Define the tolerance tol = 1e5 %Compute the Jacobian at the first homogeneous guess mat_ref [J] = jacobian_3d(I,elec,vtx,simp,gnd_ind,mat_ref,zc,v_f,dfj,tol,'{n}'); %Its condition number cond(J) %and its singular values s = svd(J); figure(1); semilogy(s,'.'); %Check that cond(J) > condest(Eref) %Are the measurements linearly independent? rank(J) %Save workspace as %name_work v_f J %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Exercise 10: Linear inverse solution 1. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %The Tikhonov regularized direct linear solver, %Add couple of smoothing regularization matrices %based on neighboring vertices [R1] = iso_f_smooth(simp,1); %and faces [R2] = iso_f_smooth(simp,3); %Let the relaxation parameter a = 1e7; %and gengerate some white Gaussian noise ns = 1e6*randn(size(refH)); %or maybe even more noise %The corresponding linear regularized solutions for R1 xt1 = mat_ref + (J.'*J + a*R1.'*R1)\J.' * ((voltageH  refH)+ns); %and R2 xt2 = mat_ref + (J.'*J + a*R2.'*R2)\J.' * ((voltageH  refH)+ns); %Write the solutions to a VTK (Visualization ToolKit) file for looking at later, writevtkfile('xt1',vtx,simp,xt1) %and writevtkfile('xt2',vtx,simp,xt2) %Save workspace as %name_work xt1 xt2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Exercise 11: Linear inverse solution 2. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %The linear Landweber iteration b = (voltageH  refH)+ns; x = zeros(size(mat_ref)); t = max(svd(J.'*J)); for i=1:1000 %usually a few thousands are just about enough xL = x + t*J.'*(b  J*x); end %and the preconditioned conjugate gradient (pcg) solution [xpcg,flag,relres,iter,resvec] = pcg(J.'*J,J.'*b,1e4,50); %Save workspace as %name_work xL xpcg %Write the solutions to a VTK (Visualization ToolKit) file using, writevtkfile('xL',vtx,simp,xL) %and writevtkfile('xpcg',vtx,simp,xpcg) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Exercise 12: Nonlinear inverse solution. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Newton's nonlinear iterative solver %For 5 iterations use [solf] = inv_sol(I,voltageH,1e5,mat_ref,vtx,simp,elec,2, ... zc,'{n}',gnd_ind,a,R1,5); %Save workspace as %name_work solf %Write the solutions to a VTK (Visualization ToolKit) file using, writevtkfile('solf',vtx,simp,solf) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Exercise 13: Complex admittivity computations. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Define a complex reference admittivity mat_refC = (1+0.5i)*ones(size(simp,1),1); %and compute the reference system matrix, forward solution, and complex data [ErefC,D,ElaC,ppr] = fem_master_full(vtx,simp,mat_refC,gnd_ind,elec,zc,'{n}'); [VrefC] = forward_solver(ErefC,I,1e5,pp); [refHC,refVC,indH,indV,df]=get_3d_meas(elec,vtx,VrefC,Ib,2); %complex measurement fields [v_fC] = m_3d_fields(vtx,size(elec,1),indH,ErefC,1e5,gnd_ind); %Determine the integral of the gradients matrix for first order tetrahedral elements [IntGradC] = integrofgrad(vtx,simp,mat_refC); %Define the tolerance tolC = 1e5 %complex Jacobian [JC] = jacobian_3d(I,elec,vtx,simp,gnd_ind,mat_refC,zc,IntGradC,v_fC,dfj,tolC,'{n}'); %Introduce some complex changes [matC,grp] = set_inho(srf,simp,vtx,mat_refC,1.1 + 0.6i); %and [matC,grp] = set_inho(srf,simp,vtx,matC,0.9 + 0.6i); %Compute the data with the inhomogeneity in place [EC,D,Ela,pp] = fem_master_full(vtx,simp,matC,gnd_ind,elec,zc,'{n}'); [VC] = forward_solver(EC,I,1e5,pp); [voltageHC,voltageVC,indH,indV,df]=get_3d_meas(elec,vtx,VC,Ib,2); %and a Tikhonov regularized solution xtC = mat_refC + (JC.'*JC + a*R2.'*R2)\JC.' * ((voltageHC  refHC)+ns); srC = real(matC); rrC = real(xtC); siC = imag(matC); riC = imag(xtC); %Plot the real part of the solution hsrC = figure(1); set(hsrC,'NumberTitle','off'); set(hsrC,'Name','Original conductivity distribution'); subplot(2,3,1); [fc] = slicer_plot_n(2.63,srC,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=2.63'); subplot(2,3,2); [fc] = slicer_plot_n(2.10,srC,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=2.10'); subplot(2,3,3); [fc] = slicer_plot_n(1.72,srC,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=1.72'); subplot(2,3,4); [fc] = slicer_plot_n(1.10,srC,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=1.10'); subplot(2,3,5); [fc] = slicer_plot_n(0.83,srC,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=0.83'); subplot(2,3,6); [fc] = slicer_plot_n(0.10,srC,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=0.10'); hrrC = figure(1); set(hrrC,'NumberTitle','off'); set(hrrC,'Name','Reconstructed conductivity distribution'); subplot(2,3,1); [fc] = slicer_plot_n(2.63,rrC,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=2.63'); subplot(2,3,2); [fc] = slicer_plot_n(2.10,rrC,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=2.10'); subplot(2,3,3); [fc] = slicer_plot_n(1.72,rrC,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=1.72'); subplot(2,3,4); [fc] = slicer_plot_n(1.10,rrC,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=1.10'); subplot(2,3,5); [fc] = slicer_plot_n(0.83,rrC,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=0.83'); subplot(2,3,6); [fc] = slicer_plot_n(0.10,rrC,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=0.10'); %and the imaginary part hsiC = figure(1); set(hsiC,'NumberTitle','off'); set(hsiC,'Name','Original permitivity distribution'); subplot(2,3,1); [fc] = slicer_plot_n(2.63,siC,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=2.63'); subplot(2,3,2); [fc] = slicer_plot_n(2.10,siC,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=2.10'); subplot(2,3,3); [fc] = slicer_plot_n(1.72,siC,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=1.72'); subplot(2,3,4); [fc] = slicer_plot_n(1.10,siC,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=1.10'); subplot(2,3,5); [fc] = slicer_plot_n(0.83,siC,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=0.83'); subplot(2,3,6); [fc] = slicer_plot_n(0.10,siC,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=0.10'); hriC = figure(1); set(hriC,'NumberTitle','off'); set(hriC,'Name','Reconstructed permitivity distribution'); subplot(2,3,1); [fc] = slicer_plot_n(2.63,riC,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=2.63'); subplot(2,3,2); [fc] = slicer_plot_n(2.10,riC,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=2.10'); subplot(2,3,3); [fc] = slicer_plot_n(1.72,riC,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=1.72'); subplot(2,3,4); [fc] = slicer_plot_n(1.10,riC,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=1.10'); subplot(2,3,5); [fc] = slicer_plot_n(0.83,riC,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=0.83'); subplot(2,3,6); [fc] = slicer_plot_n(0.10,riC,vtx,simp,fc); view(2); grid; colorbar; axis off; title('z=0.10'); %Save workspace as %name_work complex %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Exercise 14: Anisotropic computations. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %First create an anisotropic admittivity distribution gamma_xx = 1.0 + 0.5i; gamma_yy = 1.2 + 0.7i; gamma_zz = 0.9 + 0.3i; mat_ani = zeros(size(simp,1)*3,1); mat_ani(1:3:end) = gamma_xx; mat_ani(2:3:end) = gamma_yy; mat_ani(3:3:end) = gamma_zz; mat_ani = diag(mat_ani); %Then in bld_master replace the line which calculates Ela %.... %This is for the global conductance matrix (includes conductivity) %Ela = sparse((1:dimen*ns),(1:dimen*ns),kron((a.* Vols).',ones(1,dimen))); %with Ela1 = kron( Vols.',ones(1,dimen)); Ela2 = []; for i=1:ns Ela2 = sparse(blkdiag(Ela2,a(3*i2:3*i,3*i2:3*i))); end Ela = diag(Ela1) .* Ela2; %Save the modified bld_master function under a different name. %To generate the anisotropic Jacobian the procedure is even more trivial %All is needed is to comment the line where the x, the y and z components %of the sensitivity are added together in the isotropic case, i.e. %Jrow_u = Jrow_x3(1:3:end) + Jrow_x3(2:3:end) + Jrow_x3(3:3:end); %and use the anisotropic normalized (by the volume) admittivities Jrow = Jrow_x3 .* diag(Ela); 
Last Modified: $Date: 20170228 13:12:08 0500 (Tue, 28 Feb 2017) $ by $Author: aadler $