Programming of numerical algorithms and useful tools

Let us start with a simple OCTAVE code to generate the butterfly shown in figure 5.1 or figure 6.1, which is one of the standard examples for a dynamical system introduced by Lorenz in 1963. It is very popular as a study case for data assimilation algorithms and consists of three coupled ordinary differential equations

$$\begin{eqnarray}&&\dot{x}=\sigma (y-x)\end{eqnarray} \tag{ 6.1.1 }$$

$$\begin{eqnarray}&&\dot{y}=x(\rho -z)-y\end{eqnarray} \tag{ 6.1.2 }$$

$$\begin{eqnarray}&&\dot{z}={xy}-\beta z,\end{eqnarray} \tag{ 6.1.3 }$$

with constants $\sigma ,\rho ,\beta ,$ known as the Prandtl number, the Rayleigh number and a non-dimensional wave number, respectively. The classical values $\sigma =10,$ $\beta =8/3$ and $\rho =28$ have been suggested by Lorenz.

The state of the system is given by $\varphi ={(x,y,z)}^{T}\in {{\mathbb{R}}}^{3}$ and the above system can be summarized into the form

$$\begin{eqnarray}&&\dot{\varphi }(t)=F(t,\varphi ),\qquad \quad \varphi (0)={\varphi }_{0}\end{eqnarray} \tag{ 6.1.4 }$$

with

$$\begin{eqnarray}&&F(t,\varphi ):=\left(\begin{array}{c}\sigma (y-x)\\ x(\rho -z)-y\\ {xy}-\beta z\end{array}\right).\end{eqnarray} \tag{ 6.1.5 }$$

The Runge–Kutta scheme of 4th-order (see for example [2]) to calculate a numerical solution to the system (6.1.1)–(6.1.3) with initial values x₀, y₀ and z₀ iteratively calculates

$$\begin{eqnarray}&&{{\bf{k}}}_{1}=F({t}_{k},{\varphi }_{k})\end{eqnarray} \tag{ 6.1.6 }$$

$$\begin{eqnarray}&&{{\bf{k}}}_{2}=F\left({t}_{k}+\frac{1}{2}h,{\varphi }_{k}+\frac{1}{2}h{{\bf{k}}}_{1}\right)\end{eqnarray} \tag{ 6.1.7 }$$

$$\begin{eqnarray}&&{{\bf{k}}}_{3}=F\left({t}_{k}+\frac{1}{2}h,{\varphi }_{k}+\frac{1}{2}h{{\bf{k}}}_{2}\right)\end{eqnarray} \tag{ 6.1.8 }$$

$$\begin{eqnarray}&&{{\bf{k}}}_{4}=F\left({t}_{k}+h,{\varphi }_{k}+h{{\bf{k}}}_{3}\right)\end{eqnarray} \tag{ 6.1.9 }$$

$$\begin{eqnarray}&&{\varphi }_{k+1}={\varphi }_{k}+\frac{1}{6}h({{\bf{k}}}_{1}+2{{\bf{k}}}_{2}+2{{\bf{k}}}_{3}+{{\bf{k}}}_{4}).\end{eqnarray} \tag{ 6.1.10 }$$

Code 6.1.1.   The following OCTAVE or MATLAB function simulates the solution to the ordinary differential equation (6.1.1)–(6.1.3) by a 4th-order Runge–Kutta scheme.

Zoom In Zoom Out Reset image size

Here, we use the function sim_06_1_2_F as given by the following code.

Code 6.1.2.   Forcing term for the Lorenz 1963 dynamical system.

MATLAB or OCTAVE allows both functions and scripts. Functions hide the values of variables and allow clean calls of functionality. For code development, often simple scripting is helpful at the beginning. In any case, parameters need to be set and functions need to be called. We recommend carrying this out using well-defined scripts as in the following code 6.1.3, since this allows a clear debugging of how values are set.

Code 6.1.3.   A simple script sim_06_1_3_butterfly.m to call the Runge–Kutta solution of the Lorenz 1963 system to generate the 'butterfly' shown in figure 5.1.

In the next section we will show how data assimilation methods can easily be applied to the Lorenz 1963 model.

For data assimilation algorithms, we need to be able to apply the model to some initial state and run it some time T into the future. Here, we realize this model by a routine called Lorenz63. It is based on the Runge–Kutta solution of the Lorenz 1963 model described in section 6.1. We have set up this model routine such that we can generate an ensemble of $L\in {\mathbb{N}}$ states by running the model for an ensemble of L initial conditions.

Code 6.2.1.   The dynamics is generated by the model routine sim_06_2_1_Lorenz63.m.

The code here calls the routine sim_06_1_1_Lorenz63.m, which generates the full trajectory between time t = 0 and time t = T with a spacing h of the time steps and ${N}_{\text{Time}}$ such steps to reach T. This has been chosen just for simplicity, of course for other applications one would modify the original routine and not store all intermediate states.

Now, we first run the model to generate some kind of truth. This is usually denoted as a nature run. In practical applications the nature run is different from the model—as nature is different from our approximations to it. So it is important to thoroughly choose the parameters and set-up of the nature run to reflect this difference between the truth and the model which is to be used for the data assimilation. Here, we choose to modify the Lorenz 1963 parameter σ for the nature run.

Observations are modeled by the observation operator H. Here, we choose a linear operator, which observes either one component of the system, or the sum of two components, or the full state.

Code 6.2.2.   Code sim_06_2_2_Generate_Original_and_Data.m to generate the nature run and the measurement data.

The next step is running the data assimilation algorithm. Here, we start with three-dimensional variational assimilation (3D-VAR) or Tikhonov regularization. We need to define the covariance matrices, which here we choose by hand as diagonal matrices.

A key point here is to make sure that the model is different from the nature run model. Here, we modify the parameter σ. The initial state for our assimilation is chosen to be identical to the initial state of the nature run.

The assimilation cycle takes model runs and assimilation steps in turn. We carry out Nnat such cycling steps for all available data y_j at time t_j , $j=1,\ldots ,\;{\mathtt{\text{Nnat}}}.$ From the analysis xa at some time step ${t}_{j-1}$ we run the model forward and calculate the background xb at time t_j . Then, the analysis step following equation (5.2.14) with $\alpha =1$ is carried out and xa is calculated at time t_j . We store it in some vector and also calculate the error e_j at time t_j , i.e. the norm difference between the analysis and the truth. Now, the whole analysis cycle can be summarized into a short script.

Code 6.2.3.   Code sim_06_2_3_Carry_out_Assimilation.m to carry out the data assimilation by a 3D-VAR algorithm.

Finally, we display the nature run, the background states and the analysis states in typical three-dimensional graphics as shown in figure 6.2. The following code also generates a figure with the error curve e_j , $j=1,\ldots ,\;{\mathtt{\text{Nnat}}}.$

Zoom In Zoom Out Reset image size

Code 6.2.4.   Code sim_06_2_4_Plot_Original_and_Analysis.m to display the nature run, the background states and the analysis states.

Finally, we note that these scripts can be called by a control script to run them all one by one, just call sim_06_2_5_Control.m.

Code 6.2.5.   Script to control the nature run, the data assimilation cycle and visualization.

The goal of this section is to show how the ensemble Kalman square root filter (SRF) can be set up in a simple way for the Lorenz 1963 system to carry out the assimilation. We will also show that for some observation operators this ensemble data assimilation system is better than 3D-VAR.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The generation of the nature run and the simulation of measurements is described in code 6.2.2. Then, we need to generate an initial ensemble and carry out the data assimilation step in a cycled way as described by equation (5.2.14) with the covariance matrix B approximated according to equation (5.5.15) and the calculation of the analysis ensemble by the SRF as described in (5.5.22). The control code is given as follows, generating figure 6.3.

Code 6.3.1.   Script sim_06_3_1_control_EnKF.m to carry out the ensemble Kalman SRF assimilation cycle for the Lorenz 1963 system.

And the ensemble Kalman filter is carried out by code 6.3.2.

Code 6.3.2.   Script sim_06_3_2_EnKF_square_root_filter.m to carry out the ensemble Kalman SRF assimilation cycle for the Lorenz 1963 system.

We note that for a three-dimensional system such as Lorenz 1963 we do not need techniques such as localization. Also, we need to remark that we employed some inflation by a factor $\gamma =1.2$ in line 19 of code 6.3.2, i.e. we have pushed the spread of the analysis ensemble up by a multiplicative factor.

If we do not use some kind of inflation to take care of the model error, the spread of the background ensemble is only increased by the divergence of the model dynamics when different ensemble initial conditions are used. In this case, it will become too small and not properly take into account the difference between the truth (the nature run) and the model ensembles.

To learn more about regularization methods let us solve an integral equation of the first kind with some Gaussian kernel. Such equations can be employed as models for atmospheric temperature or humidity reconstruction in an atmospheric column from measured infrared or microwave radiation.

Let $k :[0,H]\times [0,H]\to {\mathbb{R}}$ for some maximal height $H\gt 0$ be a continuous kernel function, $f\in C([0,H])$ be some function representing observations and $\varphi \in C([0,H])$ be an unknown profile function. In atmospheric radiance retrievals, φ could be the temperature or humidity in an atmospheric column and $f(\nu )$ could be measurements of brightness temperature at particular frequencies ν. Here we consider f as a function of the maximum value of its kernel $k(\nu ,\cdot )$ in $[0,H]$—such that ν becomes a height variable and is no longer a frequency or some channel number. We search for a solution to

$$\begin{eqnarray}&&{\int }_{0}^{H}k(\nu ,s)\varphi (s) {\rm{d}}s=f(\nu ),\quad \nu \in [0,H]\end{eqnarray} \tag{ 6.4.1 }$$

by a collocation method as follows. Let ${\bf{x}}\in {{\mathbb{R}}}^{m}$ be a discretization of the function φ at quadrature points

$$\begin{eqnarray}&&{s}_{\xi }=\frac{\xi }{n}\cdot H\in [0,H],\quad \xi =1,\ldots ,\;n.\end{eqnarray} \tag{ 6.4.2 }$$

We replace the integral by a quadrature formula, which in the most simple case is the rectangular rule, evaluate the formula for ${\nu }_{\xi } := \xi /n\cdot H$ for $\xi =1,\ldots ,\;n,$ employ the definitions

$$\begin{eqnarray}&&{\bf{y}}:=\;{(f({\nu }_{1}),\ldots ,f({\nu }_{n}))}^{T}\in {{\mathbb{R}}}^{m},\quad {\gamma }_{\xi } := \frac{H}{n},\end{eqnarray} \tag{ 6.4.3 }$$

for the right-hand side at the points ${y}_{\xi },$ $\xi =1,\ldots ,\;n$ and with quadrature weights ${\gamma }_{\xi },$ to obtain the finite-dimensional discrete linear system

$$\begin{eqnarray}&&\sum _{\xi =1}^{n}k({s}_{j},{s}_{\xi }){\gamma }_{\xi }\ {{\bf{x}}}_{\xi }={{\bf{y}}}_{j},\quad j=1,\ldots ,\;n.\end{eqnarray} \tag{ 6.4.4 }$$

The integral operator $A\;:C([0,H])\to C([0,H])$ defined by

$$\begin{eqnarray}&&(A\varphi )(\nu ):={\int }_{0}^{H}k(\nu ,s)\varphi (s) {\rm{d}}s,\quad \nu \in [0,H],\end{eqnarray} \tag{ 6.4.5 }$$

is approximated by the matrix

$$\begin{eqnarray}&&{\bf{A}}:={\left(k({s}_{j},{s}_{\xi }){\gamma }_{\xi }\right)}_{j,\xi =1,\ldots ,n},\end{eqnarray} \tag{ 6.4.6 }$$

such that the operator equation (6.4.1) is approximated by

$$\begin{eqnarray}&&{\bf{Ax}}={\bf{y}}.\end{eqnarray} \tag{ 6.4.7 }$$

Let k be a Gaussian kernel, i.e.

$$\begin{eqnarray}&&k(\nu ,s)={{\rm{e}}}^{-\sigma | \nu -s{| }^{2}},\quad \nu ,s\in [0,H].\end{eqnarray} \tag{ 6.4.8 }$$

We define a true solution ${{\bf{x}}}_{\text{true}} :=\;{(\mathrm{cos}({s}_{\xi }))}_{\xi =1,\ldots ,n},$ data by ${\bf{y}}:={{\bf{Ax}}}_{\text{true}}+\delta$ where δ is some random error generated by a uniform distribution between ±0.01. The unregularized solution is calculated by $x={{\bf{A}}}^{-1}{\bf{y}}$ and the regularized solution via Tikhonov regularization (3.1.24) ${{\bf{x}}}_{\alpha }:=\;{(\alpha I+{{\bf{A}}}^{*}{\bf{A}})}^{-1}{{\bf{A}}}^{*}y.$ The OCTAVE code is given by code 6.4.1 and the results are displayed in figure 6.4.

Code 6.4.1.   Script sim_06_4_1_radiance_igl.m to demonstrate the implicit ill-posedness of radiance data retrieval and radiance assimilation.

Figure 6.4 demonstrates the unregularized and regularized solutions to the integral equation (6.4.1), where we chose regularization parameter $\alpha ={10}^{-3}.$ The behavior of the solutions depends on the number of quadrature points chosen. The more quadrature points we choose, the higher the ill-posedness of the inverse problem, reflected in figure 6.4(a) with n = 24 quadrature points and (b) with n = 26 quadrature points.

The following chapter 7 will introduce inverse theory for neural fields. Here, we want to carry out the basic first steps to simulate such fields. The neural field equation in its simplest form is an integro-differential equation of the form

$$\begin{eqnarray}&&\tau \dot{u}(x,t)=-u(x,t)+{\int }_{{\rm{\Omega }}}w(x,y)f(u(y,t))\quad {\rm{d}}y\end{eqnarray} \tag{ 6.5.1 }$$

for $x\in {\rm{\Omega }},t\in [0,T]$ on some domain Ω in ${{\mathbb{R}}}^{m}$ for $m=1,2,3$ and time interval $[0,T],$ where τ is some parameter, w is a neural kernel and f is a sigmoidal function

$$\begin{eqnarray}&&f(s,\eta )=\frac{1}{1+{{\rm{e}}}^{-(s-\eta )}}.\end{eqnarray} \tag{ 6.5.2 }$$

As is well known, it is possible to numerically solve such an equation by explicit time iteration known as an Euler scheme defined by

$$\begin{eqnarray}&&\frac{\tau }{{h}_{t}}({{\bf{u}}}_{k+1}-{{\bf{u}}}_{k})=-{{\bf{u}}}_{k}+{\bf{W}}f({{\bf{u}}}_{k},\eta ),\quad k=0,1,2,\ldots ,\end{eqnarray} \tag{ 6.5.3 }$$

for the state vectors ${{\bf{u}}}_{\xi }\in {{\mathbb{R}}}^{m}$ at time ${t}_{\xi },$ where we have chosen a discretization

$$\begin{eqnarray}&&{({{\bf{u}}}_{k})}_{\xi }:=u({x}_{\xi },\quad {t}_{k}),\quad k=0,1,2,\ldots ,\; \xi =1,\ldots ,\; n\end{eqnarray} \tag{ 6.5.4 }$$

where ${x}_{\xi },$ $\xi =1,\ldots ,\; n$ with $n\in {\mathbb{N}}$ points is a discretization of the domain Ω and ${t}_{k}:={h}_{t}k,$ $k=0,1,2,\ldots$ with the time grid spacing ${h}_{t}\gt 0.$ The kernel ${\bf{W}}$ based on $w({x}_{\xi },{x}_{j})$ for $\xi ,j=1,\ldots ,\;n$ is defined as the matrix

$$\begin{eqnarray}&&{\bf{W}} :=\;{((w({x}_{\xi },{x}_{j}){s}_{j}))}_{\xi ,j=1,\ldots ,n}\end{eqnarray} \tag{ 6.5.5 }$$

with integration weight s_j at point x_j , $j=1,\ldots ,\;n.$ From (6.5.3) we obtain

$$\begin{eqnarray}&&{{\bf{u}}}_{k+1}={{\bf{u}}}_{k}+\frac{{h}_{t}}{\tau }\left(-{{\bf{u}}}_{k}+{\bf{W}}f({{\bf{u}}}_{k},\eta )\right),\quad k=0,1,2,\ldots .\end{eqnarray} \tag{ 6.5.6 }$$

For our example, we choose a kernel defined by

$$\begin{eqnarray}&&w(x,y)\;:= c{{\rm{e}}}^{-{\sigma }_{1}\cdot {({x}_{1}-{y}_{1})}^{2}-{\sigma }_{2}{({x}_{2}-{y}_{2})}^{2}}\cdot ({x}_{1}-{y}_{1}),\quad x,y\in {\rm{\Omega }},\end{eqnarray} \tag{ 6.5.7 }$$

which we denote as a directed Gaussian in direction x₁. We start some simulation in a rectangular neural patch Ω with the initial state given by a Gaussian

$$\begin{eqnarray}&&u(x,0):={{\rm{e}}}^{-\sigma {| x-{x}_{0}| }^{2}},\quad x\in {\rm{\Omega }},\end{eqnarray} \tag{ 6.5.8 }$$

with some ${x}_{0}\in {\rm{\Omega }}.$

For many examples we need to quickly generate points on a regular grid. This is carried out in an elementary way in lines 1–13 of code 6.5.1, where p1v and p2v contain the x₁ and x₂ coordinates of the calculation nodes on a regular grid on ${\rm{\Omega }}=[0,{a}_{1}]\times [0,{a}_{2}].$ This is explained in detail in figure 8.2 and section 8.1.

Zoom In Zoom Out Reset image size

Based on these points, the n × n-matrix representing the neural kernel on the grid with n points is defined in lines 14–18, compare figure 6.5. Initialization of the neural field equation is carried out in lines 19–25. The Euler scheme (6.5.6) is performed by the loop of the lines 26–30.

Code 6.5.1.   Script sim_06_5_1_neural_simulation.m to calculate the temporal evolution of some neural activity function according to the neural field equation (6.5.1). Further, we show the code for the sigmoidal function f.m.

Zoom In Zoom Out Reset image size

We display the state evolution of the neural activity function $u(x,t)$ in figure 6.6 by showing the activity at three different time steps. The OCTAVE code to calculate these images is shown in code 6.5.2.

Code 6.5.2.   Script sim_06_5_1_show_graphics.m to display the temporal evolution (figure 6.6) of some neural activity function according to the neural field equation (6.5.1).

For many inversion tasks, working with masks and operators defined on two- or three-dimensional images is very advantageous. Here, we briefly provide a survey of such operators.

Zoom In Zoom Out Reset image size

Definition of masks. A straightforward realization of the masks is carried out as follows. We define

$$\begin{eqnarray}{\chi }_{G}(x):=\left\{\begin{array}{ll}1 & x\in G\\ 0 & {\rm{\text{otherwise}}}.\end{array}\right.\end{eqnarray} \tag{ 6.6.1 }$$

The discretized version of the mask in the form of a column vector is defined by

$$\begin{eqnarray}&&{{\boldsymbol{\chi }}}_{G}:={\left({\chi }_{G}({p}_{l})\right)}_{l=0,\ldots ,M}.\end{eqnarray} \tag{ 6.6.2 }$$

Examples are given in figure 6.7. We reorder the column vector according to (8.1.17) to obtain a matrix version of the mask ${{\bf{m}}}_{G}$ of dimension ${M}_{1}\times {M}_{2}.$

In the terminology of image processing the mask ${{\bf{m}}}_{G}$ is a binary digital image of dimension ${M}_{1}\times {M}_{2}.$ There is a large range of possible binary operations which can be employed for image processing and is very useful for the construction of scattering domains. Before we introduce several of these operations, we need to consider the construction of binary masks from boundary representations of domains.

Construction of masks. The construction of the masks can be a non-trivial task, depending on the particular form of the test domain G. If G is a ball with center z and radius R, then the mask can be easily constructed by testing the distance of a point $x\in { \mathcal G }$ to z. The same procedure can be applied to all star-like domains which are given in parametrized form.

For more general domains it is more difficult to decide whether a point $x\in { \mathcal G }$ is inside or outside G. We briefly mention two possibilities.

• (a)  
  It is possible to restrict oneʼs attention to a simple class of domains for which every point x outside G has a straight line connecting x with infinity. In this case we can test all lines passing from a point x as to whether they intersect $\partial G.$ If we find one which is not intersecting $\partial G$ we know that x is on the outside, otherwise it is in $\bar{G}.$
• (b)  
  In two dimensions, using the winding number of a curve, we can calculate a mask for some domain G with given boundary parametrization. In complex coordinates $z={z}_{1}+{\rm{i}}{z}_{2}$ with ${z}_{1},{z}_{2}\in {\mathbb{R}},$ it is carried out by the Cauchy integral

  $$\begin{eqnarray}&&f(z):=\frac{1}{2\pi i}{\int }_{\partial G}\frac{1}{y-z}{\rm{d}}y,\quad z\in {{\mathbb{R}}}^{2}.\end{eqnarray} \tag{ 6.6.3 }$$

  For our purposes transform (6.6.3) into the form

  $$\begin{eqnarray}&&f(z)=\frac{1}{2\pi }{\int }_{0}^{2\pi }\left\{\frac{({z}_{1}-{y}_{1})\cdot {{\rm{\Gamma }}}_{2}^{^{\prime} }(t)}{{| z-y| }^{2}}-\frac{({z}_{2}-{y}_{2})\cdot {{\rm{\Gamma }}}_{1}^{^{\prime} }(t)}{{| z-y| }^{2}}\right\}{\rm{d}}t\end{eqnarray} \tag{ 6.6.4 }$$

  with the tangential vector ${\rm{\Gamma }}(t)^{\prime}$ of the curve Γ parametrized over $[0,2\pi ].$ Figure 6.8 demonstrates the generation of masks via the Cauchy integral (6.6.4).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Masking operations. The elementary set theoretic operations such as intersections and unions of sets can easily be implemented via masking operations. Let G and H denote two domains in ${{\mathbb{R}}}^{2}$ with masks ${{\bf{m}}}_{G}$ and ${{\bf{m}}}_{H}.$ Then the intersection is carried out by

$$\begin{eqnarray}&&{{\bf{m}}}_{G\cap H}={{\bf{m}}}_{G}\odot {{\bf{m}}}_{H},\end{eqnarray} \tag{ 6.6.5 }$$

where the symbol $\odot$ denotes pointwise (element-wise) multiplication of matrices or vectors. Unions of G and H are calculated by

$$\begin{eqnarray}&&{{\bf{m}}}_{G\cup H}={\rm{sign}}\left({{\bf{m}}}_{G}+{{\bf{m}}}_{H}\right),\end{eqnarray} \tag{ 6.6.6 }$$

which means that we add the masks of G and H by regular matrix addition and then pointwise take the signum function. Both operations (6.6.5) and (6.6.6) can be generalized to products or sums with n terms, i.e. we have

$$\begin{eqnarray}&&{{\bf{m}}}_{{G}_{1}\cap \ldots \cap {G}_{m}}={{\bf{m}}}_{{G}_{1}}\odot \ldots \odot {{\bf{m}}}_{{G}_{m}},\end{eqnarray} \tag{ 6.6.7 }$$

$$\begin{eqnarray}&&{{\bf{m}}}_{{G}_{1}\cup \ldots \cup {G}_{m}}={\rm{sign}}\left({{\bf{m}}}_{{G}_{1}}+\cdots +{{\bf{m}}}_{{G}_{m}}\right).\end{eqnarray} \tag{ 6.6.8 }$$

Another useful operation is the following shrinking process. On a binary matrix ${\bf{m}}$ we define

$$\begin{eqnarray}&&S{{\bf{m}}}_{i,j} :=\;\mathrm{min}\left\{{{\bf{m}}}_{i,j},{{\bf{m}}}_{i+1,j},{{\bf{m}}}_{i-1,j},{{\bf{m}}}_{i,j+1},{{\bf{m}}}_{i,j-1}\right\},\end{eqnarray} \tag{ 6.6.9 }$$

i.e. the matrix element $S{{\bf{m}}}_{i,j}$ is one if and only if all the corresponding adjacent matrix elements of ${\bf{m}}$ are one. The application of the shrinking operator S removes all boundary elements of a domain G given by the binary matrix ${\bf{m}}.$ Analogously, we can define the thickening operator

$$\begin{eqnarray}&&T{{\bf{m}}}_{i,j} :=\;\mathrm{max}\left\{{{\bf{m}}}_{i,j},{{\bf{m}}}_{i+1,j},{{\bf{m}}}_{i-1,j},{{\bf{m}}}_{i,j+1},{{\bf{m}}}_{i,j-1}\right\}.\end{eqnarray} \tag{ 6.6.10 }$$

A smoothing operator of level ${\ell }$ can be defined by

$$\begin{eqnarray}{L}_{{\ell }}{{\bf{m}}}_{i,j} :=\;\left\{\begin{array}{ll}1, & \sum \left\{{{\bf{m}}}_{i,j},{{\bf{m}}}_{i+1,j},{{\bf{m}}}_{i-1,j},{{\bf{m}}}_{i,j+1},{{\bf{m}}}_{i,j-1}\right\}\gt {\ell }\\ 0, & {\rm{\text{otherwise}}}.\end{array}\right.\end{eqnarray} \tag{ 6.6.11 }$$

This operator removes pixels which have less or equal to ${\ell }$ neighbors of value 1.

The operators of shrinking, thickening and smoothing are usually employed for purposes of noise removal and smoothing. We employ them for our probe and sampling implementations. The smoothing effect is demonstrated in figure 6.9, where we apply the shrinking operation two times to a noisy mask and then thicken it by an application of the thickening operator.

Multiplication operators. Multiplication operators will be helpful for several inversion schemes. Given a function $f\in {L}^{2}(M)$ with some set M and a function $g\in C(M),$ let M_g be the operator which maps ${L}^{2}(M)$ into itself by multiplication of every function f by g, i.e.

$$\begin{eqnarray}&&({M}_{g}\quad f)(z):=g(z)\cdot f(z),\quad z\in M.\end{eqnarray} \tag{ 6.6.12 }$$