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: smooth and radial case ==
3
4 (% class="p1" %)
5 (% style="font-size: 10.0pt;line-height: 13.0pt;" %)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||style="font-size: 10.0pt;line-height: 13.0pt;" shape="rect"]](% style="font-size: 10.0pt;line-height: 13.0pt;" %)
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 a smooth and rotationally symmetric conductivity.
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 smooth and rotationally symmetric conductivity.
22
23 Note carefully that although we use the rotational symmetry of the conductivity to speed up some computations, 
24 the reconstruction process including the solution of the D-bar equation is two-dimensional, not one-dimensional.
25
26 Please download the Matlab routines below to your working directory and run them in the order they appear.
27
28 ----
29
30 == Definition of the example ==
31
32 The first example concerns a rotationally symmetric and smooth conductivity that equals one near the unit circle.
33 Outside the unit disc the conductivity has value 1.
34
35 The following file defines a rotationally symmetric and smooth conductivity in the unit disc: [[attach:sigma.m]].
36
37 Furthermore, this file implements the Schrödinger potential related to the conductivity: [[attach:poten.m]].
38
39 The Laplace operator appearing in the definition of the potential is implemented by finite differences in poten.m.
40 The following routines plot the conductivity and the potential, respectively: [[attach:sigma_plot.m]], [[attach:poten_plot.m]].
41 (% style="font-size: 10.0pt;line-height: 13.0pt;" %)Please run the plot commands before continuing to make sure that everything is working properly. 
42 You should see something like this:
43
44 (% style="font-size: 10.0pt;line-height: 13.0pt;" %)[[image:attach:sigma.png]]
45
46 (% style="font-size: 10.0pt;line-height: 13.0pt;" %)[[image:attach:poten.png]]
47
48 (% style="font-size: 10.0pt;line-height: 13.0pt;" %)
49
50 ----
51
52 (% style="font-size: 10.0pt;line-height: 13.0pt;" %)
53
54 == Computation of the scattering transform via the Lippmann-Schwinger equation ==
55
56 (% style="font-size: 10.0pt;line-height: 13.0pt;" %)Next we define a set of points in the k-plane for evaluating the scattering transform t(k).
57
58 Because of the rotational symmetry of this example, it is enough to choose k-values along the positive real axis: 
59 the scattering transform is known in this case to be rotationally symmetric and real-valued. Why? See [[proof>>attach:radial_t.pdf]].
60
61 So please download and run this file: [[attach:kvec_comp.m]].
62
63 The above file //kvec_comp.m// defines a set of k-points and saves them to a file called 'data/kvec.mat'. 
64 (Note that //kvec_comp.m// creates a subdirectory called 'data'. If you already created it before, Matlab will show 
65 a warning. However, you don't need to care about the warning.)
66
67 When running the example for the first time you might just use the file kvec_comp.m as it is. 
68 Later you might want to modify it to choose a different set of k-values.
69
70 Now that we have decided on the k-points, it's time to evaluate the scattering transform. Here we do it first by 
71 'cheating', or by knowing the actual conductivity, because then there are no ill-posed steps involved. 
72 Later we will compute the scattering transform also honestly from (simulated) EIT measurements. 
73 This is the file that evaluates the scattering transform t(k) at the k-points: [[attach:tLS_comp.m]].
74
75 The result will be saved to a file called 'data/tLS.mat'. Here 'LS' refers to the use of the Lippmann-Schwinger
76 equation in the computation of the complex geometric optics solutions. //tLS_comp.m// needs these files:
77 \\[[attach:green_faddeev.m]], [[attach:g1.m]], [[attach:GV_grids.m]], [[attach:GV_LS.m]], [[attach:GVLS_solve.m]], [[attach:GV_project.m]], [[attach:GV_prolong.m]]
78
79 In the names of the above files, 'GV' refers to Gennadi Vainikko, a numerical analyst who invented 
80 the periodization-based algorithm used here for solving Lippmann-Schwinger type equations.
81
82 It is interesting to compare the linear and nonlinear Fourier transform. To that end, we compute the Fourier 
83 transform of the potential q using the routines [[attach:Fq_comp.m]] and [[attach:gaussint.m]].
84
85 After running the files tLS_comp.m and Fq_comp.m, you can plot the results using [[attach:tLS_plot.m]].
86 You should see something like this:
87
88 [[image:attach:tLS_Fq.png]]
89
90 In the above image, the scattering transform is not quite approaching the value t(0)=0 as k tends to zero, 
91 although we know from theoretical results that this should be the case. Why such an error at zero? 
92 This is because we only used the value M=7 for constructing the grid, which was then of size 128x128. 
93 The Faddeev fundamental solution has a log(|k|) singularity at the origin, and the Lippmann-Schwinger 
94 type approach has always difficulties near the origin. You can either use as big M value as your patience 
95 and computer memory allows, or you can use the boundary integral equation approach below to compute 
96 t(k) for k near zero.
97
98 The rule of thumb is: The Lippmann-Schwinger approach for computing t(k) is good for k somewhat 
99 away from the origin, and the boundary integral approach for computing t(k) is good (only) for k near zero.
100
101 ----
102
103 == Simulation of EIT data ==
104
105 Since the conductivity is rotationally symmetric, the Dirichlet-to-Neumann map can be approximated by
106 a diagonal matrix in the Fourier basis. Why? See [[proof>>attach:eigenvals.pdf]].
107
108 The diagonal elements have been precomputed and are simply given as a data file [[attach:DN1eigs.mat]].
109
110 The DN matrix is constructed by the routine [[attach:DN_comp.m]].
111
112 For details of the computation of the matrix elements, see
113 **Mueller J L and Siltanen S 2003**,
114 //Direct reconstructions of conductivities from boundary measurements//,
115 SIAM Journal of Scientific Computation **24**(4), pp. 1232-1266. [[PDF (617 KB)>>url:http://www.siltanen-research.net/SIAM_EIT.pdf||shape="rect"]]
116
117 ----
118
119 == Computation of the scattering transform via the boundary integral equation ==
120
121 We need to build matrices for the single layer operators S_k parametrized by the complex number k.
122 This is done by the routines [[attach:Hk_comp.m]] and [[attach:H1.m]].
123
124 Here //Hk_comp.m// computes the matrices and //H1.m// is an auxiliary function. Note that the order Ntrig 
125 of trigonometric approximation (in other words, the number of basis functions used) has been chosen 
126 in the routine //DN_comp.m// above and saved to disc for later reference. The routine Hk_comp.m loads 
127 Ntrig from file.
128
129 Running the routine //Hk_comp.m// is computationally the most demanding task on this page. 
130 Do not be surprised if it takes 10 minutes or more.
131
132 Note that we save time by not running //Hk_comp.m// for too large values of |k|. Namely, the ill-posedness 
133 of the EIT problem has the effect that the boundary integral equation cannot be solved for |k| exceeding 
134 a certain threshold value R; for k values satisfying |k|>R the computation will produce numerical garbage.
135
136 Once Hk_comp.m has been run, we can solve the boundary integral equation for the traces of the 
137 complex geometric optics solutions. This is done by the routine [[attach:psi_BIE_comp.m]], and the result is saved
138 to a file in the subdirectory 'data'. The next step is to evaluate the scattering transform by integration over
139 the boundary using the routine [[attach:tBIE_comp.m]].
140
141 We are now ready to plot the result and compare it to the scattering transform computed using the
142 Lippmann-Schwinger equation approach. Run the file [[attach:t_plot.m]], and you should see something like this:
143
144 [[image:attach:scatLSBIE.png]]
145
146 Note how the magenta line (computed using the boundary integral equation) diverges for |k|>8.
147 However, the values of the magenta line for k near zero are very accurate.
148
149 ----
150
151 == Reconstruction from scattering transform using the D-bar method ==
152
153 Now we are ready to reconstruct the conductivity. Download the routines [[attach:tBIErecon_comp.m]] and [[attach:DB_oper.m]],
154 and run //tBIErecon_comp.m//.\\
155
156 You can look at the reconstruction using the routine [[attach:recon_plot.m]]. You should see something like this:
157
158 [[image:attach:recon.png]]
159
160 Here we used R=4, so the reconstruction is not very close to the original. Try setting M=8 and R=7 in 
161 //tBIErecon_comp.m// and see what //recon_plot.m// produces then.
162
163 Note carefully that although we used the rotational symmetry of the conductivity to speed up some of the above 
164 computations, the solution of the D-bar equation is a two-dimensional process, not one-dimensional. 
165 We did choose the reconstruction points along the positive x1-axis for convenience, but any planar point x 
166 could be chosen. Also, note that the reconstruction at one x point is completely independent from the reconstruction 
167 at another point, so the D-bar method allows region-of-interest imaging and trivial parallellization.
168
169 You can experiment with the truncation radius R. When you use tBIE computed using the boundary integral equation, 
170 you can take R up to 7 with no problems. However, when R is so large that the bad-quality parts of the above magenta 
171 plot are being used, the reconstruction will be bad.
172
173 You can take the experiment further by using tLS instead of tBIE; then you can push the reconstruction to higher 
174 values of R. Also, you can replace the low-quality tLS values near k=0 with the higher-quality tBIE.
175
176 \\
177
178 (% style="font-size: 10.0pt;line-height: 13.0pt;" %)