IDENTIFICATION OF SOUND-SOFT 3D OBSTACLES FROM PHASELESS DATA

The inverse problem for time-harmonic acoustic wave scattering to recover a sound-soft obstacle from a given incident field and the far field pattern of the scattered field is considered. We split this problem into two subproblems; first to reconstruct the shape from the modulus of the data and this is followed by employing the full far field pattern in a few measurement points to find the location of the obstacle. We extend a nonlinear integral equation approach for shape reconstruction from the modulus of the far field data [6] to the three-dimensional case. It is known, see [13], that the location of the obstacle cannot be reconstructed from only the modulus of the far field pattern since it is invariant under translations. However, employing the underlying invariance relation and using only few far field measurements in the backscattering direction we propose a novel approach for the localization of the obstacle. The efficient implementation of the method is described and the feasibility of the approach is illustrated by numerical examples.


Introduction.
In practical applications such as nondestructive testing, radar, sonar or medical imaging, the inverse obstacle scattering problem for acoustic waves occurs for frequencies in the resonance region, that is, for scatterers and wave numbers k such that the wavelength 2π/k is of a comparable size to the diameter of the scatterer.This paper is devoted to a method for shape reconstruction of threedimensional obstacles from phaseless data.Although the problem of phase retrieval is extensively studied for high frequencies only a few results are available for intermediate frequencies in the resonance region.The idea of the method originates from the combination of two papers by Kress and Rundell [13,14] and it was applied already for the inverse scattering problem in two dimensions, see [6].
We proceed by formulating the inverse problem in its mathematical context.Given an obstacle D, i.e., a bounded domain D ⊂ IR 3 with C 2 boundary Γ such that Γ can be bijectively mapped onto a sphere, consider the scattering of a plane wave u i (x) = e ikx• d with wave number k > 0 and a unit vector d describing the direction of propagation.The direct obstacle scattering problem consists of finding the total field u = u i + u s as a solution to the Helmholtz equation (1) ∆u + k 2 u = 0 in IR 3 \ D satisfying the Dirichlet boundary condition (sound-soft) (2) u = 0 on Γ such that the scattered wave u s fulfills the Sommerfeld radiation condition (3) lim r→∞ r ∂u s ∂r − iku s = 0, r = |x|, uniformly with respect to all directions.This condition ensures uniqueness of the solution to the direct scattering problem and guarantees that the scattered wave is outgoing.Furthermore, the radiation condition is equivalent to an asymptotic behavior of the form uniformly in all directions, with the far field pattern u ∞ defined on the unit sphere Ω in IR 3 .
In this paper we investigate the following three inverse problems: Inverse Problem (IP1): Given the far field pattern u ∞ on Ω for one incident wave u i , determine the shape and the location of the boundary surface Γ of the scatterer D. Inverse Problem (IP2): Given the modulus of far field pattern |u ∞ | on Ω for one incident plane wave u i , determine the shape of the boundary surface Γ of the scatterer D. Inverse Problem (IP3): Given the far field pattern u ∞ in a few measurement points on Ω for one incident plane wave u i and knowing the shape of the surface Γ, determine the location of the scatterer D.
The inverse problem (IP1) is well studied, in particular due to Colton and Sleeman [2] and Gintides [4] it is known that a sound-soft obstacle contained in a ball of radius R can be uniquely identified from the knowledge of the far field pattern for one incident plane wave provided kR ≤ 4.49.
The solution to (IP2) is not unique since the modulus of u ∞ is invariant under translations [13], that is, for the shifted domain D ℓ := {x + ℓ : x ∈ D} with ℓ ∈ IR 3 a fixed vector, the far field pattern u ℓ,∞ satisfies the relation This ambiguity cannot be remedied by using finitely many incident waves with different wave numbers or different incident directions [13].The relation (4) also holds in the case of the Neumann and impedance boundary conditions, see [15].Unfortunately, even for an obstacle with a fixed location currently there are no uniqueness results available for reconstructing its shape from the modulus of far field data, except for the uniqueness result by Liu and Zhang, see [16], for a soundsoft ball centered at the origin.The plan of the paper is as follows.In the next section a system of nonlinear integral equations equivalent to (IP2) is derived and an iterative scheme is suggested.In Section 3 we discuss the injectivity of the linearized operator in the proposed method.Section 4 is devoted to presenting the numerical solution of the integral equations and in Section 5 implementation details are described.A novel approach for recovering the location of the obstacle, that is solving (IP3), is introduced in Section 6.Finally, in Section 7 we present numerical examples and conclude by discussing some closely related methods for obstacle reconstruction from full far field data.

2.
Iterative scheme based on nonlinear integral equations.For the scattering of an entire field u i from a sound-soft obstacle D by Huygens' principle, see Theorem 3.12 in [1], we have ( 5) and the far field pattern of the scattered field u s is given by where ν is the outward unit normal vector to Γ and is the fundamental solution to the Helmholtz equation in IR 3 .
We introduce the single-layer operator S Γ : L 2 (Γ) → L 2 (Γ) and the far field operator S Γ,∞ : The sound-soft boundary condition (2), ( 5) and ( 6) imply the equivalence of the inverse obstacle scattering problem (IP2) and the two-by-two system of integral equations ( 7) where h = ∂u ∂ν Γ , in the sense that a solution to (IP2) satisfies the system and vice versa, see [6].
For the numerical solution of the integral equations ( 7)-( 8) we assume that the surface Γ is C 2 -smooth, homeomorphic to the unit sphere Ω and has a star-shaped representation Γ := Γ r = {r(x)x : x ∈ Ω} with r(x) > 0 for x ∈ Ω.With the aid of the substitution where J r = r r 2 + |∇ s r| 2 is the Jacobian of the transformation r and ∇ s denotes the surface gradient, we replace the integral operators and the right-hand sides of the system ( 7)-( 8) by their parameterized form Hence, the parametrized form of the nonlinear system ( 7)-( 8) reads (12) S r υ = w r , and Three closely related approaches for solving the latter system for the unknown r, that is, for the reconstruction of the unknown scatterer have been presented in the literature, see the survey [10] and the references therein.A first method was suggested by Sleeman [19,11] for sound-soft obstacle reconstruction from the full far field in two dimensions.It consists in solving the mildly ill-posed linear equation (12) for υ and then linearizing (13) to find r.
Reversing the roles of ( 12) and ( 13) a second approach is obtained.In a slight modification, termed a hybrid method, it was investigated in a series of papers by Kress and Serranho, among those also the three-dimensional problem was considered in [18].
In this paper we concentrate on a third approach that was suggested by Kress and Rundell [14] for the Laplace equation and then extended to different scattering problems including, in particular, the shape reconstruction from only the modulus of the far field pattern [6].This method is based on a simultaneous linearization of both equations with respect to both unknowns.In the first step of the iterative procedure we make an initial guess for the surface Γ parametrized by r and find the density υ from (12).Then given a current approximation for r and υ we look for density and surface updates ϑ and q that satisfy the linearized system (14) S r υ + S r ϑ + S ′ r [υ]q = w r + w ′ r q and (15) With these, the approximations for the density and the surface are updated via υ = υ + ϑ and r = r + q, respectively.As a stopping criterion for the iterative procedure we choose (16) ǫ for some sufficiently small parameter τ > 0 depending on the noise level.
In the iterative scheme ( 14)- (15) the Fréchet derivatives of the operators S r and S r,∞ and the right-hand side w r are required.They can be obtained via differentiating their kernels with respect to the radial function r, see [17], and have the form To solve the linear integral equations ( 14)-( 15) numerically we introduce spherical coordinates on Ω, i.e. (19) p(θ, φ) := (sin θ cos φ, sin θ sin φ, cos θ), θ ∈ [0, π], φ ∈ [0, 2π], and choose appropriate approximation spaces for the density υ and the radial function r.Since the spherical harmonics Y l,j (θ, φ) = c j l P |j| l (cos θ) e ij φ , j = −l, . . ., l, l = 0, 1, 2, . . ., with coefficients form a complete orthonormal system in L 2 (Ω) they provide a natural choice for the approximation spaces for the unknown density.Here P |j| l are the associated Legendre functions of degree l and order |j|.Since the radial function r is real valued we modify the basis functions into (20) Thus, the density approximation reads and for the surface approximation we employ 3. On noninjectivity of the linearized operator.Due to [6, Theorem 3.1] the first iteration of our procedure is equivalent to one step of the method suggested in [13], therefore we can reduce the analysis for the injectivity of the operator in the linearized equations ( 14) and ( 15) at the exact solution to the latter case.The inverse obstacle scattering problem (IP2) can be formulated in terms of the operator F : r → u ∞ that maps the radial function r onto the far field pattern and is defined via the solution to the direct problem ( 1)-( 3) with a fixed incident wave.(IP2) can be solved via Newton's method, i.e. by solving the linearized equation It has been proved [12,17] that the operator F is Fréchet differentiable with the derivative given by where v ∞ is the far field pattern of the solution to the Helmholtz equation in the exterior of Γ satisfying the Sommerfeld radiation condition and the Dirichlet boundary condition with x q (x) = q(x) x and x is defined through the relation r(x) x = x for x ∈ Γ.We also recall [13, Lemma 1] which states that the solution w to the Helmholtz equation given by and has far field pattern where ℓ ∈ IR 3 is a fixed vector.Let us consider the special case when D is the unit ball.Choosing ℓ = (1, 0, 0) in ( 26) we obtain ℓ • ν = sin θ cos φ and hence ν • x q = ℓ • ν for q(θ, φ) = sin θ cos φ.Then from ( 24), (26) via uniqueness for the direct Dirichlet problem we have that Repeating the same considerations for ℓ = (0, 1, 0) and ℓ = (0, 0, 1) we observe that the relation (28) also holds for q(θ, φ) = sin θ sin φ and q(θ, φ) = cos θ, respectively.
Summarizing the facts listed above, since according to (28) the expression F r F ′ [r]q is purely imaginary we have established the following lemma.
Lemma 3.1.The columns in the matrix of the discretized version of ( 14)- (15) corresponding to the coefficients of Y IR 1,j with j = −1, 0, 1 for the radial function q vanish.
The vanishing columns in the matrix correspond to the fact that the location of the obstacle cannot be recovered from the modulus and the algorithm is modified to find only (N + 1) 2 − 3 coefficients in (22), that is, to reconstruct the shape.The three coefficients of Y IR 1,j for j = −1, 0, 1 are kept fixed.We also want to emphasize on the fact that the method suggested in this paper for (IP2) as opposed to (23) does not require the solution of the forward problem in each iteration step.Furthermore, the Fréchet derivatives with respect to Γ of the boundary integral operators involved in our method can be explicitly characterized as integral operators, see ( 17) and (18), which reduces the computational costs.
4. Numerical solution of the integral equations.For the numerical solution of the surface integral equations ( 14)-( 15) we use Wienert's method [20] which is based on spherical harmonics and on the transformation of the boundary surface to a sphere.For the numerical integration of a continuous function over the unit sphere we apply the so called Gauss trapezoidal product rule where p is the parameterization of the unit sphere Ω given by (19).The quadrature points are given by (30) with the zeros z s ′ of the Legendre polynomials P n ′ +1 of degree n ′ + 1 and the corresponding weights are This quadrature rule is derived by introducing a suitable projection operator L n ′ onto the space of spherical harmonics of order less than or equal to n ′ + 1 and integrating the projected function L n ′ f exactly observing the characteristics of the trapezoidal and the Gauss-Legendre rules, see [1].This way the quadrature rule (29) inherits its name and by construction it is exact for spherical harmonics of order less than or equal to n ′ + 1.Some of the integral operators which we need to treat have weakly singular kernels, e.g. S r and S ′ r .To cope with this difficulty we will introduce a coordinate system change in Ω and transfer the singularity to the north pole n = (0, 0, 1).The quadrature rule for such an integrand has the form This rule is based on the orthonormality of spherical harmonics and the relation In the rest of this section we will show that the quadrature rules stated above can be applied for the integral operators involved in the solution of the inverse problem.It is easy to see that the kernels of S r,∞ and S ′ r,∞ are smooth.The kernel of the integral operator S r we split into a weakly singular and a continuous part, i.e. with and Here, by the parameter r in the definition of the functions we indicate the dependence of the kernels on the unknown surface.
In order to move the singularity in the first integral to the north pole we consider an orthogonal linear transformation T x such that T x x = n, see [20,1,3].We also introduce an induced transformation T x as x ŷ), f ∈ C(Ω), and its bivariate analogue The orthogonality of T −1 x yields the identity In addition to this, we introduce a function R by which is smooth as a function with the first component kept fixed at the north pole, i.e.R(n, •) : η → R(n, η) for η ∈ Ω, [1,5].Finally, the integral operator S r given by ( 10) can be represented in the form with the kernel H 1 (x, ŷ; r) = R(x, ŷ)H 1 (x, ŷ; r) that is smooth in the sense of the mapping η → T xH 1 (n, η; r), see [5,Lemma 4.6].Moreover, the singularity is canceled out by ds(η) = sin θ dθ dφ.Hence the quadrature rules (31) and ( 29) can be applied to the first and the second integrals in (32), respectively.
Here, we note that in the numerical evaluation of the integral operator (32) it is sufficient to only apply the rotation T x to the singular part since the object function υ in both integrals is approximated by (21) and, hence, the discretized version of the operator acts on the vector of its coefficients.

5.
Implementation.With the intention to obtain a fully discrete linear system we apply the quadrature rules ( 29) and (31) to the integral operators and evaluate them in the same points as that used for the quadrature rules.For the sake of brevity we introduce the following notations with the quadrature points θ s , θ s ′ , φ ρ , φ ρ ′ defined in (30) and the spherical parametrization function p given by (19).The linear transformation operator is defined via T p(θ,φ) = P (φ)Q(−θ)P (−φ) with 3 × 3 matrices P (ψ) and Q(ψ) corresponding to positive rotations by ψ around the x 3 -and x 2 -axis, i.e.
For the efficient implementation of the algorithm we use the representation of the rotated spherical harmonics derived in [3] (34) with is the normalized Jacobi polynomial, with In the case when m − j or −m − j is negative one can compute d jm (π/2) by using the symmetry relation From (34) for the rotated real valued spherical harmonics (20) via representing the imaginary and real parts by subtracting or adding the complex conjugate we obtain (35) We start with discussing the implementation of the discrete version of the operator S ′ r [υ] defined in (33) since it has the largest computational costs.The approximation is given by (36) As can be seen from the latter representation, in a naive implementation, one needs at least O((n ′ + 1) 6 × (N + 1) 2 ) operations with N ≤ n ≤ n ′ .In the implementation scheme which we will suggest further below the computational cost can be reduced to O((n ′ + 1) 6 × (N + 1)).
For each of the four terms in the summation of the expression (36) we will set up the matrix separately, denoting them D sρlj , A sρlj , F 1 sρlj and F 2 sρlj , respectively.Then the elements of the matrix S ′ r [ υ] are given by Recalling the definition of M 1 and keeping in mind that υ and r are given by ( 21) and (22) in each iteration step one can easily see that precomputing the values Y lj ρsρ ′ s ′ = Y l,j (ŷ ρ ′ s ′ ρs ) and Y IR lj ρsρ ′ s ′ = Y IR l,j (ŷ ρ ′ s ′ ρs ) before the first iteration significantly speeds up the algorithm.We also precompute Y lj ρs = Y l,j (x ρs ) and Y IR lj ρs = Y IR l,j (x ρs ) for the evaluation of the kernel M 2 .The coefficients of the 2(n ′ + 1) 2 × ((N + 1) 2 − 3) matrix S ′ r [υ] can be evaluated recursively through the scheme The matrix corresponding to the first line of (36) is built up by two steps where IR indicates that the function is real valued and it is defined analogously to (20) and The computationally most involved of all matrices is the one arising from the singular kernel, in particular, from the product of the function M 1 and the rotated spherical harmonics.To fill in the matrix A sρlj we perform the following three recursive steps To evaluate the matrices F 1 sρlj and F 2 sρlj which correspond to the continuous part of the kernel in the integral operator we proceed via The matrix for the discretized operator S ′ r,∞ [υ] is set up in the same way as F 2 sρlj .Let us consider the discrete single-layer operator The recursive scheme for its coefficients is given by The elements of ( S r,∞ ) sρlj are found by the same procedure as the elements of the latter matrix D sρlj .
The fully discrete system for ( 14)-( 15) with the unknown complex valued coefficients ϑ lj and the real valued coefficients q lj can be written in a brief form as Here A is a 2(n ′ +1) 2 × (n+1) 2 +(N +1) 2 −3 matrix.Since the system inherits the ill-posedness of the problem we incorporate for its numerical solution a Tikhonov regularization, i.e. instead of (37) we solve (38) where α it , β it are regularizations parameters, and the matrix I 1 corresponds to the Sobolev H 1 (Ω) penalty term.Since the associated Legendre functions are the canonical solutions of the Legendre equation we obtain a simple representation for its elements by 6. Recovering the location of the scatterer.In this section we introduce an algorithm for the numerical solution of the problem (IP3).Assume that we have already reconstructed the shape of the unknown obstacle from the modulus of the far field pattern by the method described in Section 2 and denote by Γ 0 = {r(x)x : x ∈ Ω} the surface of the scatterer centered at the origin and by Γ ℓ = {r(x)x + ℓ : x ∈ Ω} the surface of the scatterer shifted to the point ℓ.Recalling the relation (4) between the far fields produced by these scatterers we obtain a system of equations for recovering the unknown location ℓ.We note that the value u 0,∞ is available as a byproduct of the algorithm for the shape reconstruction and this spares us from solving an additional forward problem.Taking the logarithm on both sides of equation (39) we have However, the logarithm ln is a multi-valued function and we have to choose one branch of it, for example, the principal branch.This induces a limitation on the length of the vector ℓ which can be reconstructed.The right-hand side of ( 40 Then d − x can be represented as For an arbitrary vector ℓ we can estimate the inner product ( d − x) • ℓ by Finally, we obtain the condition | on the length of the vector ℓ that is controlling the unknown location.From this estimate we observe that for small values of α and β one can find the location of obstacles with larger distance from the origin.Hence, for the numerical solution of (40) we take the far field patterns u 0,∞ and u ℓ,∞ in several points in the vicinity of the backscattering direction.For the two dimensional case the estimate (41) is simplified to .
Concluding this section we would like to remark that although for the location reconstruction some a priori information is necessary it is not very restrictive as compared with the a priori information required by the Newton method for the inverse scattering problem with the full far field data.Furthermore, since the translation invariance (4) also holds in the case of sound-hard and impedance obstacles this algorithm for recovering the location carries over without any changes.
7. Numerical examples.In this section we illustrate the robustness of the proposed methods for the shape and the location reconstruction.To avoid committing an inverse crime we generate a far field pattern by solving the direct problem via a Galerkin method for the coupled single-and double-layer potential approach, see [3].
The synthetic data consists of 128 values of the modulus of the far field data for the unknown scatterer and 10 values of the full (amplitude and phase) far field data measured in the backscattering direction.Noise was then added pointwise, i.e. 1% of noise level means that for each data point with value u a value δu = 0.01ηu was added where η is a random number from the interval [−1, 1].In the case of the complex valued data, noise was added separately to the real and imaginary part.
For all examples the same parameters were set up.The number of discretization and quadratures points is equal to 128, i.e. n ′ = 7; the dimension of the approximation spaces for the density and the surface is n = 7 and N = 6, respectively.The wave number was chosen k = 1 and the incident direction d = (0, 1, 0) T , which is indicated by the arrow on the figures.As an initial guess we take a ball of radius 3Y 0,0 ≈ 0.85 centered at the origin, i.e. r 1j = 0 for j = −1, 0, 1.We note here that the coefficients r 1j for j = −1, 0, 1 will be kept fixed during the iterations due to Lemma 3.1.As long as βγ it ≤ 0.0001 we update the regularization parameters in each iteration step of (38) by α it = αγ it and β it = βγ it with α = 0.01, β = 0.5, γ = 2/3.The parameter τ for the stopping criterion ( 16) is chosen as τ = 0.01.
In the figures below the reconstructions of the shape and location of the surfaces are presented from data with a 2% noise level.In the captions of the figures we point out the number of iterations and the obtained relative error (16).
In the upper part of Figure 1 we present the reconstruction of the shape of a pinched ball parameterized by z(θ, φ) = 1.44 + 0.5 cos 2φ(cos 2θ − 1) p(θ, φ), θ ∈ [0, π], φ ∈ [0, 2π], from the modulus of the far field pattern, which is centered at the origin in the figure .The far field data are given for the obstacle which was shifted to the point ℓ † = (0.3, −0.5, 4) T .
Using the full far field data in 10 points in the backscattering direction via (40) we find the location ℓ of the obstacle and place the exact object to the point ℓ † − ℓ.
To illustrate the difference between the reconstructed and the exact surfaces we added to the plot their projections onto the x 1 x 2 , x 2 x 3 and x 1 x 3 planes, which are indicated as meshed and solid regions, respectively.
As one can see from Figure 1, recovering of the location is very accurate as well as the reconstruction of the surface in the illuminated part.In the shadow region of the surface the quality of the reconstructions deteriorates.
In the second example, Figure 2, we consider the reconstruction of a cushionshaped surface with the parametrization which was shifted to ℓ † = (−2, 0.3, 3) T .The notations are kept the same as for the first example.
For the third example we consider a surface which does not belong to the class of obstacles that are star-shaped with respect to the origin.The parametrization of this bean-like surface reads and it is located at the point ℓ † = (1, 1, −1) T .
In Figure 3 we observe that the approximate location is still very close to the exact location although the reconstruction of the surface is not as good as in the previous two examples.This was expected to happen since the surface is not starshaped with respect to the origin and has a large concave part and the data is not exact.
Summarizing, we conclude that the method for reconstructing the shape from only the modulus of the far field data provides very satisfactory results.Moreover the simple method for recovering the location turned out to be feasible, stable against noise in the data and to require only a few input data.
8. Iterative approaches to the numerical solution of (IP1).The final section we devote to some closely related methods for solving the inverse problem (IP1), that is, to reconstruct the obstacle from the full far field data.Via marginal modifications, analogous to two dimensions [9], the method suggested in Section 2 can be adapted to the full data case.The injectivity of the operator proved in [7] carries over without changes to the three dimensional case.
However, a second method as introduced by Sleeman [19,11] cannot be extended directly to three dimensions and requires more modifications.Let us consider the operator (11) but without the substitution (9), i.e. we introduce the modified operator as where ϕ(x) = h(r(x)x) and J r is the Jacobian of the transformation r.Analogously we modify the operator S r given by (10).
The Fréchet derivative of the operator S r,∞ is given by 4π Ω e −ik x•ŷ r(ŷ) ϕ(ŷ) J r (ŷ)q(ŷ) ds(ŷ) The iterative scheme consists in two steps; for the current approximation r find the density ϕ from the mildly ill-posed equation S r ϕ = w r , then solve the linearized equation (43) S r,∞ ϕ + S ′ r,∞ [ϕ]q = w ∞ for the update q and calculate new radial surface approximation as r = r + q.This procedure is repeated until a suitable stopping criterion is satisfied.Due to ill-posedness of (43) one has to apply a regularization method, e.g.Tikhonov regularization.
The shape reconstructions from the full far field pattern via the method with simultaneous linearization and the method (43) are of the comparable quality with those presented in Section 7 and the reconstruction from [18].Naturally, the method based on (43) can be extended to the inverse scattering problem with the modulus of the far field pattern as data, the detailed comparison of the approaches is deferred to a future research.
In the paper [8] the methods were investigated for general parameterizations.It was found out that the largest distance between the location of the exact obstacle and the initial guess is restricted by the condition that both objects should be contained in the disc of radius 2.4048/k.Furthermore, the reconstruction strongly depends on whether the initial guess is placed in the illuminated or the shadow region of the scatterer and uses the full far field pattern.The approach suggested in Section 6 allows to recover the obstacles with much weaker a priori information on the location (42) and requires only a few measurements of the full far field data.For the exact data it is enough to know the far field pattern in 2 or 3 appropriate measurement points under the condition (42) and (41) for two or three dimensions, correspondingly.
The methods discussed in this paper can be adopted to the inverse scattering problem for a three dimensional sound-hard obstacle by using the indirect potential approach, i.e. by representing the scattered field via a single-layer potential.In this case the kernels of the integral operators will contain at most a weak singularity as opposed to the direct approach which involves a hypersingular integral operator.
For phaseless near field data the approach suggested in Section 2 can be applied with minor changes and it does not suffer from the translation invariance.

) takes its values from the interval 1 k
[−π, π] and, hence, the maximum of the absolute value of inner product ( d − x) • ℓ should be less than or equal to π k .Let us consider a parametric representation of the vectors d and x in the form d =   sin θ cos φ sin θ sin φ cos θ   and x =   sin(θ + α) cos(φ + β) sin(θ + α) sin(φ + β) cos(θ + α)   .

Figure 1 .
Figure 1.Reconstruction of a pinched ball after 6 iterations with ǫ r = 0.0098 and the exact surface

Figure 2 .
Figure 2. Reconstruction of a cushion after 16 iterations with ǫ r = 0.0099 and the exact surface

Figure 3 .
Figure 3. Reconstruction of a bean after 26 iterations with ǫ r = 0.0156 and the exact surface x• dq(x).