# 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:

- pdaPlot.html
- plotAndSaveFigures.html
- plotExpSetup.html
- plotPredefinedContrast.html
- subplots.html

## 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.

## See Also

- example.html

- varalpha.html
- varbeta.html
- varalphabeta.html
- vardelta.html
- vartol.html

- init.html
- setInput.html
- setData.html
- setRecon.html
- recon.html