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

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

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

Content of output folder

Output: Figures

A table of possible plotted and saved figures and their number.

Fig. no.Content
01-10 Preliminary definitions...
expSetup.m
01source and measurement points plot...
setContrast.m
02predefined contrast (real part)
03predefined contrast (imag part)
04-10currently unused
11-20 Reconstruction
11Subplots (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
12discrepancy and error (as in subplot)
13Tikhonov functional (as in subplots): M1, M2, MT
14reconstructed contrast (real part)
15reconstructed contrast (imag part)
16difference of reconstructed and true contrast (abs) (switched off by if 0)
17-20currently unused
21-40 Additional for 3D reconstruction (plotAndSaveFigures.m)
(contrast in a sectional plane through the scatterer)
21-23sectional plane through the predefined contrast (real part) (1st, 2nd, 3rd direction)
24-26sectional plane through the predefined contrast (imaginary part) (1st, 2nd, 3rd direction)
27-29sectional plane through the reconstructed contrast (real part) (1st, 2nd, 3rd direction)
30-32sectional plane through the reconstructed contrast (imaginary part) (1st, 2nd, 3rd direction)
33-40currently unused
41-50 PDA specific additional reconstruction plots
over inner iterations (pdaPlot.m)
41Parts of minimization functional: fd, fs, fg, fp
42Error in pda and relLinDisInPda
43pda lin. disc. and nonlin. disc.
44parts of min. func. M1 = F(Kh), M2 = G(h) (switched off by if 0)
45ThetaiOut (inner tolerance principle) (switched off by if 0)
46tau, sigma values (switched off by if 0, because currently sigma and tau are not saved)
47-50currently 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).

Usage of advanced starts

We demontrate it in case of various $\alpha$ 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 $\alpha$ = 100, 500, and 1000 are stored inside.

Code: Part 1: Initialization, set input and start diary

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:

% 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

% 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

% 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 $u^\mathrm{i}: \bf{R}^d \to \bf{C}$ in dimension $d = 2, 3$ with time-dependence $\exp(-\mathrm{i} \omega t)$, where $\omega > 0$ is the angular frequency.

When we omit the time-dependence, the incident field satisfies the Helmholtz equation

$$ \Delta u^\mathrm{i}(x) + k^2 u^\mathrm{i}(x) = 0, \quad x \in \bf{R}^d $$

with constant wave number $k>0$.

The scattering object is described by a refractive index function $n: \bf{R}^d \to \bf{C}$. This function equals $1$ outside a bounded domain.

When the incident wave interacts with a scattering object a scattered wave $u^\mathrm{s}$ is generated. The total wave

$$ (1) \quad u^\mathrm{t} := u^\mathrm{i} + u^\mathrm{s} $$

admits the Helmholtz equation

$$ (2) \quad \Delta u^\mathrm{t} + k^2 n^2 u^\mathrm{t} = 0 \quad \textrm{in } \bf{R}^d $$.

The scattered wave $u^\mathrm{s}$ satisfies Sommerfeld's radiation condition

$$ (3) \quad \lim\limits_{|x|\to\infty} |x|^{(d-1)/2}\left(
\frac{\partial}{\partial |x|} - \mathrm{i} k\right) u^\mathrm{s}(x) = 0, $$

uniformly in all directions $\hat x = x/|x| \in \bf{S} := \{ y \in \bf{R}^d: \, |y| =1\}$, where $|y|:= \sqrt{y_1^2 + \ldots + y_d^2}$.

The direct scattering problem is to find a solution $u$ 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

$$ q := n^2-1 \quad \textrm{in } \bf{R}^d. $$

The inverse scattering problem

Keep in mind that $u^\mathrm{s} = u^\mathrm{t}-u^\mathrm{i}$. If we know the incident wave $u^\mathrm{i}$ and measure the total wave $u^\mathrm{t}$, we can compute the scattered wave $u^\mathrm{s}$.

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 $q$, that minimizes the Tikhonov functional of the non-linearized problem, i.e.:

$$ \min_{q \in X}
     \underbrace{\frac{1}{2}\|\mathcal{F}(q)-F_\mathrm{meas}^\delta\|_\mathrm{F}^2}_{\mathrm{discrepancy}}
   + \underbrace{\alpha \|q\|_\mathrm{1}}_{\mathrm{sparsity\ penalty}}
   + \underbrace{\beta \| \nabla q \|_\mathrm{1}}_{\mathrm{total\ variation\ penalty}}
   + \underbrace{
         \delta_{[a,b]}( \mathrm{Re}(q) ) +
         \delta_{[c,d]}(\mathrm{Im}(q) )}_{\mathrm{penalty\ for\ physical\ bounds}
         }.
   $$

... What we do is to minimize the Tikhonov functional of the linearized problem, i.e. (this is done in pda.html)

$$ \min_{h \in X}
     \underbrace{\frac{1}{2}\|\mathcal{F}'(q)[h]+\mathcal{F}(q) - F_\mathrm{meas}^\delta\|_\mathrm{F}^2}_{\mathrm{discrepancy\ (linearized\ problem)}}
   + \underbrace{\alpha \|q+h\|_\mathrm{1}}_{\mathrm{sparsity\ penalty}}
   + \underbrace{\beta \| \nabla (q+h) \|_\mathrm{1}}_{\mathrm{total\ variation\ penalty}}
   + \underbrace{
         \delta_{[a,b]}( \mathrm{Re}(q+h) ) +
         \delta_{[c,d]}(\mathrm{Im}(q+h) )}_{\mathrm{penalty\ for\ physical\ bounds}
         }
   $$

and then update $q := q + h$ (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 $q$.

Authors

The authors of IPscatt are:

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

See Also