Last modified by smsiltan@helsinki_fi on 2024/02/13 07:34
Show last authors
1 (% class="p1" %)
2 == EIT with the D-bar method: discontinuous heart-and-lungs phantom ==
3
4 (% class="p1" %)
5 This page contains the computational Matlab files related to the book
6 [[//Linear and Nonlinear Inverse Problems with Practical Applications// >>url:http://www.ec-securehost.com/SIAM/CS10.html||shape="rect"]]
7 written by **Jennifer Mueller** and **Samuli Siltanen** and published by SIAM in 2012.
8
9 You can order the book at the [[SIAM webshop>>url:http://www.ec-securehost.com/SIAM/CS10.html||shape="rect"]].
10
11 [[Go to master page>>url:http://wiki.helsinki.fi/display/mathstatHenkilokunta/Inverse+Problems+Book+Page||shape="rect"]]
12
13 This page is related to Example 2, a discontinuous heart-and-lungs phantom.
14
15 ----
16
17 == Introduction ==
18
19 The D-bar method is a reconstruction method for the nonlinear inverse conductivity problem arising from
20 Electrical Impedance Tomography. This page contains Matlab routines implementing the D-bar method
21 for a discontinuous heart-and-lungs phantom.
22
23 The motivation for this example comes from the following experiment:
24
25 [[image:attach:ukko.png]]
26
27 We aim to image the heart and lungs of a patient using a row of electrodes arranged as in the picture above.
28 We think of maintaining a given voltage potential at each electrode and measuring the resulting electric current
29 through the electrodes. This measurement is repeated for several voltage patterns. Assuming that there are 32
30 electrodes in total, we can use at most 31 linearly independent voltage patterns since one of the electrodes is
31 considered as the ground potential to which the other potentials are compared to. We use trigonometric voltage
32 patterns (examples shown below on left), approximated by continuous sine curves (below right) at the boundary:
33
34 [[image:attach:trigpatterns.png]]
35
36 Furthermore, we approximate the above three-dimensional (3D) situation with a two-dimensional computational 
37 model. Our virtual patient will be modelled by a two-dimensional (2D) disc with varying electrical conductivity 
38 inside. This 2D approximation, while obviously incorrect, works surprisingly well in practice.
39
40 This page is for explaining how to simulate the above kind of voltage-to-current data and how to recover the inner 
41 conductivity from the boundary measurements using the D-bar method.
42
43 Please download the Matlab routines below to your working directory and run them in the order they appear. Note 
44 that you will need the PDE Toolbox for running some of the routines.
45
46 ----
47
48 == Definition of the conductivity ==
49
50 Our simulated conductivity is defined in the file [[attach:heartNlungs.m]], and you can plot it using the routine
51 [[attach:heartNlungs_plot.m]]. The result should look something like this:
52
53 [[image:attach:heartNlungs2D_mod.png]]
54
55
56
57 Here the background conductivity is 1, the heart is filled with blood and has higher conductivity 2, and the lungs 
58 are filled with air and have lower conductivity 0.5.
59
60 Wonder about the white color at the bottom of the colorbar above? It is for is for creating the white background 
61 in Matlab. Check out [[this page>>url:http://blogs.helsinki.fi/smsiltan/2012/05/10/displaying-image-data-for-comparison/||shape="rect"]] to learn more.
62
63 ----
64
65 == Simulation of EIT data ==
66
67 We simulate by first computing the current-to-voltage map using the Finite Element Method (FEM) and then inverting.
68 In what follows, the voltage-to-current map is called //DN map// for Dirichlet-to-Neumann, and the current-to-voltage
69 map is called //ND map// for Neumann-to-Dirichlet.
70
71 FEM is an efficient method for solving elliptic partial differential equations. The solution of the conductivity equation 
72 is given as a linear combination of basis functions that are piecewise linear in a triangular mesh. The mesh is constructed 
73 by the routine [[attach:mesh_comp.m]] containing the parameter Nrefine. Here are plots of the triangle meshes with Nrefine=0 (left)
74 and Nrefine=2 (right):
75
76 [[image:attach:mesh_both.png]]
77
78 In practice it is a good idea to use Nrefine=5 for accurate results. The file //mesh_comp.m// saves the mesh to a file called 
79 'data/mesh.mat'. (Note that //mesh_comp.m// creates a subdirectory called 'data'. If the subdirectory 'data' already exists, 
80 Matlab will show a warning. However, you don't need to care about the warning.)
81
82 The simulation of the continuum-model ND map involves solving Neumann problems of the form
83
84 [[image:attach:NeumannProblem.png]]
85
86 containing a Fourier basis function defined by
87
88 [[image:attach:Fourierbasis.png]]
89
90 Both the Neumann condition and the conductivity given to Matlab's PDE toolbox routines as user-defined functions.
91 Let us see how to do that. First of all, the user-defined routines need to define the functions at appropriate points related
92 to the FEM mesh. Here are some important points:
93
94 [[image:attach:meshes.png]]
95
96 The conductivity is constant inside each triangle, so it is naturally specified at the centers of triangles showin in (a).
97 This is implemented in the routine [[attach:FEMconductivity.m]]. The solution computed by FEM is represented as a function
98 which is linear inside each triangle and thus completely determined by its values at the vertices. Therefore, the gradient
99 of the solution is constant inside each triangle, and it is appropriate to specify the Neumann data at the midpoints (c) of
100 boundary segments. This is done in routine [[attach:BoundaryData.m]].(% style="font-size: 10.0pt;line-height: 13.0pt;" %) Finally, we need to find the trace of the solution for 
101 constructing the ND map; this can be done simply by picking the values of the FEM solution at the vertices (b) located 
102 at the boundary.
103
104 The routine [[attach:ND_comp.m]] (% style="font-size: 10.0pt;line-height: 13.0pt;" %)computes by FEM a matrix approximation to the ND map and saves it to data/ND.mat.
105
106 The DN map is the inverse of the ND map in the subspace of functions defined at the boundary and having mean zero. 
107 We can construct the DN matrix by first inverting the DN matrix and then adding a zero row and a zero column for 
108 dealing with the constant basis function. This is done in routine [[attach:ex2DN_comp.m]], which also computes the DN matrix
109 analytically for the constant conductivity 1. The results are saved to file data/ex2DN.mat.
110
111 ----
112
113 == Computation of the scattering transform via the boundary integral equation ==
114
115 We need to construct a grid in the complex k-plane for the evaluation of the scattering transform. Download and run
116 this file: [[attach:ex2Kvec_comp.m]]. You can view the grid using the routine [[attach:ex2Kvec_plot.m]]. You should see something like this:
117
118 [[image:attach:ex2Kvec.png]]
119
120 \\
121
122 The boundary integral equation involves a single-layer potential operator Sk with Faddeev Green's function Gk as kernel.
123 Now Gk(z)=G0(z)+H(z,k), where G0 is the standard logarithmic Green's function for the Laplacian in dimension two
124 and H(.,k) is a harmonic function with fixed k. We write Sk=S0+Hk, where S0 is the standard single-layer operator
125 (known analytically for the disc domain considered here) and Hk is an integral operator with H(.,k) as kernel. We need to
126 compute matrix approximations to the operators Hk. This is done by routine [[attach:ex2Hk_comp.m]], which needs the files [[attach:g1.m]]
127 and [[attach:H1tilde.m]](% style="font-size: 10.0pt;line-height: 13.0pt;" %). The matrices are saved to files data/ex2Hk_1.mat, data/ex2Hk_2.mat, ...
128
129
130
131 Note that the computation of the Hk matrices may take a long time to finish (with my MacBook pro it takes 10 hours). 
132 However, this computation needs to be done only once for a given k-grid. If the k-grid is not changed, the once-computed 
133 Hk matrices can be just loaded from the disc.
134
135 You can avoid the computation by downloading the file [[data.zip>>url:http://www.siltanen-research.net/data.zip||shape="rect"]]. (% style="font-size: 10.0pt;line-height: 13.0pt;" %)Of course, if you modify //ex2Kvec_comp.m// to get a different
136 k-grid, you need to actually run //ex2Hk_comp.m//.
137
138 (% style="font-size: 10.0pt;line-height: 13.0pt;" %)
139
140 Once //ex2Hk_comp.m// has been run (or the zip file downloaded and decompressed), we can solve the boundary integral 
141 equation for the traces of the complex geometric optics solutions. This is done by the routine(% style="font-size: 10.0pt;line-height: 13.0pt;" %) [[attach:ex2psi_BIE_comp.m]], 
142 and the result is saved to a file in the subdirectory 'data'. The next step is to evaluate the scattering transform by integration 
143 over the boundary using the routine [[attach:ex2tBIE_comp.m]].
144
145 We are now ready to plot the scattering transform. Run [[attach:ex2tBIE_plot.m]], and you should see something like this:
146
147 [[image:attach:ex2scatLSBIE_mod.png]]
148
149 Above is shown real part (left) and imaginary part (right) of the scattering transform in the disc of radius 7.
150
151 ----
152
153 == Reconstruction from scattering transform using the D-bar method ==
154
155 Now we are ready to reconstruct the conductivity. Download the routines [[attach:GV_grids.m]], [[attach:GV_project.m]], [[attach:GV_prolong.m]]
156 and [[attach:DB_oper.m]], and run [[attach:ex2tBIErecon_comp.m]].
157
158 You can look at the reconstruction using the routine [[attach:ex2tBIErecon_plot.m]]. You should see something like this:
159
160 [[image:attach:ex2recon_R4.png]]
161
162 The left image shows the original conductivity distribution, and the right image shows the reconstruction using R=4. 
163 The colors (and thus numerical values) are directly comparable. Note that the reconstruction underestimates the value 
164 at the "heart."
165
166 You can try a bigger frequency cutoff radius as well. Open the file //ex2tBIErecon_comp.m// and change R=4 to R=6, 
167 say. You will need more accuracy, so change M = 8 into M=9 as well. The result looks like this:
168
169 [[image:attach:ex2recon_R6.png]]
170
171 In the above picture, the colors (and thus numerical values) are directly comparable with each other. (However,
172 the colors in the previous R=4 reconstruction image are //not// comparable to the R=6 image.) Note that this
173 time the reconstruction //overestimates// the conductivity of the "heart;" this is a typical feature of the truncated
174 nonlinear Fourier transform.
175
176 \\