start
Main file of framework IPscatt in MATLAB, i.e. Computational Framework in MATLAB for the Inverse Medium Problem in Scattering.
Warning: start clears variables and stores figures and files in the folder output. The routine eval or feval is used in setInput and setContrast.
Contents
- Quick start
- Input "Arguments"
- Output "Arguments"
- Output: Figures
- Advanced starts
- Code: Part 1: Initialization, set input and start diary
- Code: Part 2: Set data
- Code: Part 3: Variational reconstruction
- Code: Part 4: End
- More About... Theory: The direct and inverse scattering from inhomogeneous media
- Authors
- References
- See Also
Quick start
Adviced usage of MATLAB
The toolbox IPscatt generates several figures and data, which are saved in the folder output. We recommend to start MATLAB with the option nodisplay to avoid annoying pop-ups of figures.
matlab -nodisplay
Simple Example
To run the computational framework, i.e. solve the direct scattering problem, add some noise to the data and start the reconstruction process you can just type in
start
The framework IPscatt will set some default parameters, do computation, and save the results in the folder output.
Change parameters
To compute with other parameters you can set them in a file, e.g. example, in the folder inseti and refer to this file by Parameter inseti, e.g.
inseti = 'example';
start
Attention in usage of the framework
- start saves figures and files in the folder output without demand. (Adapt setInput and set out to 0 or 1 instead of 2 to change this behaviour).
- start calls setInput, in which pre-defined variables (except closed, test, and inseti) are deleted.
- Therefore you must define the struct seti in such a file (and not in terminal).
- The name of a new file in inseti must differ from existing functions (because the file in inseti is called by eval).
- The code was written and tested with MATLAB in version R2016a.
Input "Arguments"
There are no classical input arguments because start is not a function.
From default deviating parameters in struct seti are set as described in the section Quick Start.
Some input Parameters are described in example.html.
Output "Arguments"
There are no classical output arguments because start is not a function.
Output results in folder
Output results are saved in a folder created after start like this
output/20160907T111711_example_test
Structure of foldername: seti.dirOutput + / + seti.dirDatetime + _ + seti.inseti + seti.dirSuffix
- seti.dirOutput is set to "output"
- seti.dirDatetime is date and time separated by the character "T"
- seti.inseti is parameter inseti set before start (otherwise seti.inseti is set to "noinseti")
- seti.dirSuffix can set in the corresponding file in inseti
Content of output folder
- figures with prefix fig_, see table in the next section (possible output formats are png, eps and fig, see example.html)
- diary.log: contains terminal output.
- inseti_...: contains a copy of the parameters set in the file committed by inseti.
- save_Fmeas.mat: contains exact data
and noisy simulated data
as FmeasExact and FmeasDelta.
- save_qROIcomp_iOutStop.mat: contains the reconstructed contrast qROIcomp and stop index iOutStop.
- save_dis.mat, save_err.mat, and save_dif.mat are saved if seti.savedata equals 1. They contain relative discrepancy, error, and difference of computed contrast to its predecessor.
Output: Figures
A table of possible plotted and saved figures and their number.
Fig. no. | Content |
01-10 Preliminary definitions... | |
expSetup.m | |
01 | source and measurement points plot... |
setContrast.m | |
02 | predefined contrast (real part) |
03 | predefined contrast (imag part) |
04-10 | currently unused |
11-20 Reconstruction | |
11 | Subplots (subplots.m )(1,1) ground truth contrast (real part) (1,2) reconstructed contrast (real part) (2,1) relative error and relative discrepancy (2,2) value of Tikhonov function (function to minimize) |
plotAndSaveFigures.m | |
12 | discrepancy and error (as in subplot) |
13 | Tikhonov functional (as in subplots): M1, M2, MT |
14 | reconstructed contrast (real part) |
15 | reconstructed contrast (imag part) |
16 | difference of reconstructed and true contrast (abs) (switched off by if 0 ) |
17-20 | currently unused |
21-40 Additional for 3D reconstruction (plotAndSaveFigures.m )(contrast in a sectional plane through the scatterer) | |
21-23 | sectional plane through the predefined contrast (real part) (1st, 2nd, 3rd direction) |
24-26 | sectional plane through the predefined contrast (imaginary part) (1st, 2nd, 3rd direction) |
27-29 | sectional plane through the reconstructed contrast (real part) (1st, 2nd, 3rd direction) |
30-32 | sectional plane through the reconstructed contrast (imaginary part) (1st, 2nd, 3rd direction) |
33-40 | currently unused |
41-50 PDA specific additional reconstruction plots over inner iterations ( pdaPlot.m ) | |
41 | Parts of minimization functional: fd, fs, fg, fp |
42 | Error in pda and relLinDisInPda |
43 | pda lin. disc. and nonlin. disc. |
44 | parts of min. func. M1 = F(Kh), M2 = G(h) (switched off by if 0 ) |
45 | ThetaiOut (inner tolerance principle) (switched off by if 0 ) |
46 | tau, sigma values (switched off by if 0 , because currently sigma and tau are not saved) |
47-50 | currently unused |
For Output figures see also:
Advanced starts
It is possible to run the computational framework with various input parameters. Three cases are prepared (others can add similarly).
- varalpha.m: various input parameters for the regularization parameter
(sparsity).
- varbeta.m: various input parameters for the regularization parameter
(total variation).
- vardelta.m: various input parameters for the relative noise level
.
- vartol.m: various input parameters for the tolerance in GMRES.
Usage of advanced starts
We demontrate it in case of various inputs.
In file varalpha:
alpha = [100; 500; 1000];
In terminal:
inseti = 'example'; % setting of seti.alpha inside is overwrittten by varalpha varalpha
A folder in output with suffix varalpha is created. Reconstruction results for = 100, 500, and 1000 are stored inside.
Code: Part 1: Initialization, set input and start diary
- init: init.html: Addpath and set parameter closed in init.m, which differs between closed and public version of code (new developments are not published yet).
- setInput: setInput.html: Creates folder for output; clears variables except of closed, test, and inseti; deals with various input parameters (varalpha, varbeta, vardelta).
- start diary log (output of terminal is stored in diary.log)
disp(' ') disp('################ -------- start of IPscatt ------- ################') disp(' '); disp(' Computational Framework in MATLAB for ') disp(' the Inverse Medium Problem in Scattering ') disp(' ') % Initialization disp(' ') disp('## init -- initialization -----------------------------------------') disp(' ') init; % Set input and make directories for output disp(' ') disp('## setInput -- set input and make directories for output ----------') disp(' ') setInput; % Diary log: start disp(' ') disp('## diary -- start diary log ---------------------------------------') disp(' ') diary(sprintf('%s/diary%s.log',seti.dirname,seti.fileSuffix)); % store all output in diary diary on; % later diary off % Set output depth out dispDepth = 4; % 0: nearly no display messages; 4: default; 5: too much information out = 2; % 0: no figures, 1: figures, 2: save figures or files
Code: Part 2: Set data
setData.m (setData.html) does the following:
- Load real-world data (if set).
- Set geometry and simulation (grid, kernel, experimental set-up, predefined contrast, general settings for figures)
- Compute exact and noisy data.
% Set data disp(' ') disp('## setData -- geometry and experimental set-up, noisy data -------') disp(' ') seti = setData(seti,dispDepth,out); % geometry and simulation and fresnel etc. is inside setData...
Code: Part 3: Variational reconstruction
- setRecon.m (setRecon.html): define functions for reconstruction.
- recon.m (recon.html): reconstruction process: minimize the Tikhonov functional of the linearized problem by primal-dual algorithm and update the contrast
(in public version). (The functionals are given in the section "More About...").
% Settings for variational reconstruction disp(' ') disp('## setRecon -- variational reconstruction -------------------------') disp(' ') seti = setRecon(seti,dispDepth); % Variational reconstruction (process) disp(' ') disp('## recon -- variational reconstruction -------------------------') disp(' ') seti = recon(seti,dispDepth,out);
Code: Part 4: End
- Define end of diray log.
- Save workspace (if seti.saveworkspace was set to 1)
% Diary log: end disp(' ') disp('## diary -- end diary log -----------------------------------------') disp(' ') diary off; % Save workspace seti = checkfield(seti,'saveworkspace',0); if seti.saveworkspace == 1 && out >= 2 disp(' ') disp('## save workspace ---------------------------------------------') disp(' ') save(sprintf('%s/workspace.mat',seti.dirname)); end
More About... Theory: The direct and inverse scattering from inhomogeneous media
We give a brief description of the direct and inverse scattering problem. A reference is e.g. Chapter 8 in [1].
The direct scattering problem
We consider a time-harmonic incident wave in dimension
with time-dependence
, where
is the angular frequency.
When we omit the time-dependence, the incident field satisfies the Helmholtz equation
with constant wave number .
The scattering object is described by a refractive index function . This function equals
outside a bounded domain.
When the incident wave interacts with a scattering object a scattered wave is generated. The total wave
admits the Helmholtz equation
.
The scattered wave satisfies Sommerfeld's radiation condition
uniformly in all directions , where
.
The direct scattering problem is to find a solution to equations (1)-(3).
We will not solve the above PDE. Instead we will solve the "Lippmann-Schwinger equation", which is a equivalent reformulation.
Further, instead of the refractive index we will consider the contrast
The inverse scattering problem
Keep in mind that . If we know the incident wave
and measure the total wave
, we can compute the scattered wave
.
Then the scattered wave is used as data for the inverse scattering problem: Reconstruct the contrast function q from several measurements of the scattered wave.
For this we will need multi-static measurements, i.e. several transmitters propagate incident waves one after another and the generated scattered fields are measured by several receivers. (In real-world total fields are measured and scattered fields are computed subtracting the corresponding incident fields.)
More theoretical background
The exact techniques how to tackle the direct and inverse scattering from inhomogeneous media are described in [2]. In the next section we have a brief look at the function we minimize for reconstruction.
For an overview of the software framework as well as its application areas, we refer to [3]. In addition, IPscatt is compared with existing software packages (solving the same problem) in [3].
Function to minimize...
... We would like to find a contrast , that minimizes the Tikhonov functional of the non-linearized problem, i.e.:
is the Frobenius norm.
- Actually all norms are weighted. (We ommit it to keep notation simple.)
is a indicator function, i.e. is 0, if all entries of the vector x are between
and
, and is
, otherwise.
... What we do is to minimize the Tikhonov functional of the linearized problem, i.e. (this is done in pda.html)
and then update (this is done in minPda.html).
Next, the Tikhonov functional of the linearized problem is minimized again. (This minimization is done by application of primal-dual algorithm, see [4], on this problem, see [2].)
The whole process minimizes the Tikhonov functional of the non-linearized inversion problem and leads to the seeked contrast .
Authors
The authors of IPscatt are:
- Florian Bürgel (Center for Industrial Mathematics, University of Bremen, Germany, fbuergel@uni-bremen.de.)
- Kamil S. Kazimierski (Institute of Mathematics and Scientific Computing, University of Graz, Austria, kazimier@uni-graz.at. The Institute of Mathematics and Scientific Computing is a member of NAWI Graz (www.nawigraz.at) and BioTechMed Graz (www.biotechmed.at).)
- Armin Lechleiter (Center for Industrial Mathematics, University of Bremen, Germany.)
Note that the source code of gmresKelley.html in the folder 3rdparty was not written by the authors and contains the corresponding license information.
References
- [1] David Colton and Rainer Kress. Inverse Acoustic and Electromagnetic Scattering Theory. Springer, New York, 2013. URL: https://doi.org/10.1007/978-1-4614-4942-3.
- [2] Florian Bürgel, Kamil S. Kazimierski, and Armin Lechleiter. A sparsity regularization and total variation based computational framework for the inverse medium problem in scattering. Journal of Computational Physics, 339:1-30, 2017. URL: https://doi.org/10.1016/j.jcp.2017.03.011.
- [3] Florian Bürgel, Kamil S. Kazimierski, and Armin Lechleiter. Algorithm 1001: IPscatt—A MATLAB Toolbox for the Inverse Medium Problem in Scattering. ACM Transactions on Mathematical Software, 45(4), Article 45, 20 pages, 2019. URL: https://doi.org/10.1145/3328525.
- [4] Antonin Chambolle and Thomas Pock. A first-order primal-dual algorithm for convex problems with applications to imaging. Journal of Mathematical Imaging and Vision, 40(1):120-145, 2011. URL: https://doi.org/10.1007/s10851-010-0251-1.