Complete Sets of Radiating and Nonradiating Parts of a Source and Their Fields with Applications in Inverse Scattering Limited-Angle Problems

Many algorithms applied in inverse scattering problems use source-field systems instead of the direct computation of the unknown scatterer. It is well known that the resulting source problem does not have a unique solution, since certain parts of the source totally vanish outside of the reconstruction area. This paper provides for the two-dimensional case special sets of functions, which include all radiating and all nonradiating parts of the source. These sets are used to solve an acoustic inverse problem in two steps. The problem under discussion consists of determining an inhomogeneous obstacle supported in a part of a disc, from data, known for a subset of a two-dimensional circle. In a first step, the radiating parts are computed by solving a linear problem. The second step is nonlinear and consists of determining the nonradiating parts.


INTRODUCTION
Medical imaging involves the solution of mathematical inverse problems. This means that cause (the properties of living tissue) is inferred from effect (the observed signal). In the case of ultrasound tomography, the probe consists of ultrasonic pressure waves, and echoes inside the tissue show the internal structure.
The properties of the tissue inside a body are described by the function where c is the local speed of sound. Outside the body the speed c = c 0 is constant. The function a comprises the attenuation of a wave. The wavenumber k depends on the frequency ω of the ultrasound scanner. They are simply related by the equation Outside the body the function f becomes zero. Therefore, let us assume that the support of f is a subset of Ω = {x ∈ R 2 : |x| < R}.
Furthermore, the incident wave u i , which is generated by the ultrasound scanner, solves the homogeneous wave equation The scattered field u s , which is generated by the scatterer f and the incident wave u i , is described by the Lippmann-Schwinger integral equation where G is the Green function associated with the Helmholtz operator ( + k 2 ) and the Sommerfeld radiation condition in R 2 , lim r→∞ r 1/2 ∂u s ∂r − iku s = 0, r = x .
The total field u, is determined by the sum of the incident wave and the scattered field. Equation (4) is valid for all x ∈ R 2 . In the 2 International Journal of Biomedical Imaging 2-dimensional case, where H o is the Bessel function of the third kind of order 0 (see, e.g., [1]). In order to formulate the inverse problem, we assume that the data g := u s | Γ is known for a subset of the boundary Γ ⊂ ∂Ω of the reconstruction area. If the operator A is defined as then the inverse problem consists in solving the integral equation which is a nonlinear problem, since f and u s are unknown in the interior of the domain Ω. To get rid of the nonlinearity, this problem is reformulated as a source problem by defining Now the problem reads as Denote by B the operator B : L 2 (Ω) −→ L 2 (Ω), then u s = BΦ (13) is the field generated by the source Φ.
It is well known that the solution of the source problem (11) is not unique due to the existence of some nonradiating parts of the function Φ which completely vanish outside of the region Ω. In Sections 2 and 4, complete sets of functions of the radiating and the nonradiating parts of the source Φ are constructed. This decomposition already has successfully been used in [2] to solve the inverse scattering problem, where the authors discretise the operator A and computed the singular value decomposition numerically, that contains the radiating and the nonradiating parts, respectively.
To get uniqueness of the reconstruction problem, multiple experiments have to be performed using pairwise different incident waves u m i (see [3]). Then the corresponding sources are related by which follows from (10). It has been shown, for example, in [4,5], that there exists a unique solution f using incident plane waves u α i (x) = e ik α,x if data are collected for every direction α of the unit sphere S 1 . This proves the existence of a unique solution for every Φ α for all α ∈ S 1 in this case.
In Section 2, complete sets of radiating and nonradiating sources for the case Γ = ∂Ω are constructed. The corresponding fields to these sets are computed in Section 3. In Section 4, the limited angle problem Γ ⊂ ∂Ω is considered. Complete sets of radiating and nonradiating sources are developed by using the singular value decomposition (SVD) of the operator A. Algorithms and strategies, which use the decomposition of the sources to solve the inverse problem, are given in Section 5. In the last section, numerical results are presented.

COMPLETE SETS OF RADIATING AND NONRADIATING SOURCES
In this section, it is assumed that the data are given on the complete boundary, that is, Γ = ∂Ω. The general case Γ ⊂ ∂Ω will be discussed in Section 4. Some efforts were made in [6] to construct complete sets of nonradiating sources. The authors used the differential equation which is equivalent to the integral equation (4). They determined a basis of functions for radial symmetric nonradiating sources. As an alternative approach in the current paper the integral equation is considered to get a complete set of nonradiating sources. From (15) and the definition of the source (10), it follows that a field u s generated by a source Φ ∈ L 2 (Ω), which is zero at the boundary, completely vanishes outside of Ω. Therefore, AΦ = 0 is a necessary and sufficient condition for Φ to be a nonradiating source. Now the problem is to determine a complete set of the null space N(A), which contains all nonradiating sources. On the other hand, a basis of N(A) ⊥ represents a complete set of the radiating parts of the source.
A useful tool for investigating an operator is the singular value decomposition (SVD) (v m , u m , σ m ), which was determined in [7] in the case Γ = ∂Ω. The singular functions v m and u m are given up to constants by The functions J m denote the Bessel functions of the first kind of order m, and Y m are the spherical harmonics of degree |m|, which are in the 2-dimensional case given by Y m (ω(ϕ)) = e imϕ , where ω(ϕ) are vectors of the unit sphere S 1 depending on the angle ϕ. The singular functions v m generate an orthogonal basis of the space N(A) ⊥ and hence {v m } m∈Z is a complete set of radiating sources. To complete this set to a base of L 2 (Ω) the null space N(A), which includes all nonradiating functions, has to be characterized. Assuming that a nonradiating source can be separated in a part q, which only depends on the radius, and a part Y p , depending only on E. Wallacher and A. K. Louis 3 the direction of a point x = rω ∈ Ω, Φ(rω) = q(r)Y p (ω), then the function Φ has to be orthogonal to all functions v m ∈ N(A) ⊥ . That means Introducing polar coordinates one can deduce Finally, the orthogonality property of the spherical harmonics leads to the conclusion that has to be fulfilled for Φ to be a nonradiating source. The problem changes into finding sets of functions {q p,n } n∈N , that provide in combination with the function J p (k·) a basis of the Hilbert space L 2 ([0, R], w −1 ) with the weight w(r) = r. A helpful tool to construct such systems is given by a special case of Lommel's theorem (see, e.g., [8]), which is summarized in Appendix 7. Now it is possible to get an orthogonal basis of L 2 (Ω), which includes all radiating parts {v m } m∈Z . The existence of the zeros, claimed in the theorem below, follows from Appendix 7 1-4. From formula 5 in Appendix 7 one can deduce that the number of zeros is countable.
Furthermore, let λ m,n (n ∈ N) be the positive zeros of a m J m (r)+ b m rJ m (r) (in the case where a + |m| < 0, any of the two zeros in {ix | x ∈ R + } are to be added). After permutation let λ m,o = kR. Then the functions build an orthogonal basis of L 2 (Ω) and especially Proof. It follows from Lommel's theorem (see Appendix 7) and the orthogonality property (19) of the spherical har-monics that the functions w m,n are orthogonal on L 2 (Ω). The fact that λ m,o = kR is a zero of a m J m (r) + b m rJ m (r) follows from the definition of a m and b m . Therefore, the functions w m,o (rω) are equal to J m (kr)Y m (ω) and thus the radiating parts are all contained in the basis.
To complete the proof, it has to be shown that for any m ∈ Z the set {q m,n (r)} n∈N := {J m ((λ m,n /R)r)} n∈N is a complete basis for L 2 ([0, R], w −1 ). For the case b m = 0 this was shown in [9]. In a similar way it can also be proved for the case b m = 1. A detailed proof is given in [10]. Now let g ∈ L 2 (Ω) be any function, that is orthogonal to every function w m,n . A consequence of the completeness of the spherical harmonics {Y m } m∈Z on the space L 2 (S 1 ) and the completeness of the functions {q m,n } n∈N in L 2 ([0, R], w −1 ) is that g has to be zero, which finishes the proof.
Sometimes it is better to deal with normalized functions. An approximation can always be computed by numerical integration. But in this case one can also do this analytically, which follows from Theorem 2 below.
Proof. Applying Lommel's theorem (see Appendix 7) and the orthogonality properties (19) of the spherical harmonics leads directly to the norms of w m,n .
In the next theorem the main results of this section are resumed.
Proof. Let w m,n denote the functions of Theorem 1. Because of the singular value decomposition (16) of the operator ,it is obvious that (i) holds. Again it is a consequence of Theorem 1 that (ii) is also true.

THE CORRESPONDING FIELDS AND PROPERTIES
Let w m,n be the functions of the orthogonal base of Theorem 1, then they are associated to the radiating and nonradiating parts of a source as described in Theorem 3. To develop an algorithm that solves the inverse problem, the 4 International Journal of Biomedical Imaging corresponding fields u m,n of the sources w m,n have to be computed. They are related by (13) u m,n = Bw m,n .

Theorem 4.
The fields u m,n corresponding to the sources w m,n are given by where Proof. Using the expansion of Hankel function (see Appendix 7 formula 10) and Lommel's theorem (see Appendix 7) leads to the statement given in Theorem 4.
A more detailed examination shows that the fields u m,n corresponding to the nonradiating sources w m,n (n ≥ 1) can be expressed by the functions w m,n themselves in the following way. where are constants.
Proof. Let I 1 m,n (r) and I 2 m,n (r) be the functions due to Theorem 4. From Theorem 1, we obtain where Applying formulas 9 and 11 of Appendix 7 leads to Inserting T m and S m yields E. Wallacher and A. K. Louis

5
From this we deduce that Finally, the properties of the source-field systems (u m,n , w m,n ) are summarized in the next theorem. Theorem 6. Let u m,n and w m,n be the fields and sources defined as above. Then, (v) the space of the radiating sources is orthogonal to the space of nonradiating sources.
These properties were proved in [11] for any complexvalued or vector source-field system described by a linear scalar or vector partial differential equation. In this case, they follow immediately from the outlines made before and the fact that the functions w m,n are eigenfunctions of the Laplace operator with eigenvalues −(λ m,n /R) 2 .

THE LIMITED-ANGLE PROBLEM
The singular value decomposition (SVD) is a powerful tool to analyze a problem. The range and the orthogonal complement of the nullspace of an operator is given by the singular functions. The behavior of the singular values characterizes the ill-posedness of the problem. A classification is, for example, outlined in [12]. The goal of this section is to derive the SVD of the operator A for the general case Γ ⊂ ∂Ω. To handle this problem, it is assumed that Γ is parameterized in the following way: where η is an arbitrary function defined on the sphere The measurement setup is shown in Figure 1. In this case, the center of the detector is η(α) and the length of the detector is 2β. The operator under consideration depends on the direction of the incident wave as well as of the length of the detector. The function η only ensures that the center of the detector is independent of the direction of the incident wave u α i . For instance, one can imagine that for applications the detector can always be installed in the opposite of the incident wave. That means η(α) = α. On the other hand, it is possible that the detector is fixed somewhere and so η(α) ≡ c becomes a constant. For the sake of simplicity, we replace in the following η(α) by α because the arguments in the prove are the same: To get the SVD of the operator A := A α,β , the eigenfunctions and eigenvalues of the operator A * A have to be computed. Therefore, a suitable representation of this operator is needed. In the following, α is considered as the same time as a vector on the sphere S 1 as well as the corresponding angle.
Proof. We have the identity International Journal of Biomedical Imaging Thus, we obtain Using polar coordinates x = r x ω x and the series representation of the Hankel function (see Appendix 7 formula 10), one can deduce Uniform convergence leads finally to Some useful properties of the Sinc-matrix (S(m, n, β)) m,n∈Z are summarized in the following lemma. Lemma 1. By using the notation above, the following statements hold: is an isomorphism. One can conclude that Hence C is well defined. From (v) it follows immediately that C is selfadjoint. It is well known that A is a compact operator and thus A * A is compact. From (49) one can deduce that C is also compact. Now it is possible to compute the SVD of the operator A. which is contradiction to u = v p for all p ∈ Z. Finally, the equation Av p = σ p u p has to hold for every p ∈ Z. By using polar coordinates and the series representation of the Hankel function (see formula 10 in Appendix 7), one can conclude There are similarities to incomplete data problems in Xray computerized tomography (CT), where the Sinc-matrix also played an important role for the SVD (see [13]). In CT it is well known that there exist a unique solution of the limited angle problem. Immediately the question arises, if there is also a unique solution in the inverse limited angle scattering Proof. Let It is obvious that Then, x can be expanded in a series of functions from V 2 , that is, Let P be the orthogonal projection onto V 1 . We have the identity The vectors (v β p,m ) m∈Z ∈ l 2 build an orthonormal basis and the conclusion is Hence V 1 = V 2 and Theorem 3 completes the proof.
A direct consequence of the theorem above is that the data g α,β = u α s | Γα,β can be extended to the entire boundary ∂Ω. Note that the minimum energy solution Φ min of the problem A α,β Φ = g α,β can be expressed with the help of the SVD by (63) Theorem 10. Let α, β be fixed and g α,β = u α s | Γα,β . Then where Φ min is the minimum energy solution of A α,β Φ = g α,β .
Proof. The minimum energy solution is the projection of the source Φ α onto N ⊥ (A α,β ). But from Theorems 3 and 9 it follows that N(A α,β ) ⊥ is independent on β. Furthermore, we have N(A α,β ) ⊥ = N(B| ∂Ω ) ⊥ which proves the statement.

THE INVERSE PROBLEM
Many algorithms (see, e.g., [2,14]) use the definition of the sources to reconstruct the unknown scatterer f . The problem can be divided into two separated steps. In a first step, the linear integral equation (11) has to be solved. This leads to the radiating parts of the sources Φ. In a second nonlinear step, the nonradiating parts are to be computed.

The radiating parts of a source
Because the minimum energy solution of the problem (11) is the unique solution in N(A) ⊥ , this solution contains the radiating parts of the source Φ. If the source is split into radiating and nonradiating parts Φ = Φ R + Φ NR , the radiating part can be computed using the singular value decomposition (v m , u m , σ m ) of the operator A. More precisely, The operator A is a compact operator. Thus the problem is ill-posed in the following sense (see also [15]). Small errors in the data g lead to large errors in the reconstruction Φ R . Hence a regularization has to be performed. A summary of regularization methods is given in [12]. All the methods mentioned in this book are closely related to a general framework, the approximate inverse (see [16][17][18]), where the basic idea is to approximate the solution Φ R by smoothing here the mollifier e γ (x, ·) is an approximation of the deltadistribution (e.g., e γ (x, y) = (1/2πγ)e − x−y 2 /2γ 2 the Gaussian) and the reconstruction kernel ψ γ (x, ·) is the solution of A * ψ γ (x, ·) = e γ (x, ·). The constant γ denotes the regularization parameter, which depends on the noise level of the given data g. The advantages of this method are as follows.
(i) The reconstruction kernel ψ γ can be computed before the measurement process.
E. Wallacher and A. K. Louis

9
(ii) The kernel does not depend on the data, that are afflicted by errors. (iii) If the same problem is to be solved many times (like the case for every given data g α ), then the kernel ψ γ has to be computed only once. This fact yields a good performance.
Again we apply the SVD to get the kernel ψ For the numerical implementation one has to cut off the series. For γ = 0 which corresponds to the Deltadistribution and restricting the summation from m = 0 and up to m = M that approach provides the truncated SVD.
A problem of the approximate inverse is the large memory consumption, because for every reconstruction point x ∈ Ω and for every incident wave direction α the kernel ψ α,β γ (x, ·) has to be computed. To overcome this problem, we use some invariance properties of the kernel (see also [18]).

Remark 2.
If β < π, that is, Γ is a real subset of ∂Ω, no invariance has been found. In fact it is possible to extent the data on the entire boundary by now with the help of Theorem 10, but this also needs a regularization because of the ill-posedness of the problem. For this problem all reconstruction kernels have to be computed for every α and every reconstruction point x. At least the integrals e γ (x, ·), w m,o which appear in the inversion formula for the minimum energy solution of the reconstruction kernel should be precomputed. This procedure still provides an acceptable performance of reconstructing the radiating parts of the source. In the following only the case β = π is considered.
Hence the kernel at any point x = re iϕx ∈ Ω is given by a rotation of the kernel at the point r on the positive real axis.
Proof. Let ψ α,π γ (r x e iϕx , ·) be the minimum energy solution of A * ψ γ (r x e iϕx , ·) = e γ (r x e iϕx , ·). Hence, (69) The inner products in (69) are given up to constants (see [7]) by With the approximate inverse it is possible to compute a smoothed solution of the radiating part of the source where D ϕ denotes the rotation operator with respect to the angle ϕ. This method was tested in a numerical test for a homogeneous scatterer f , which is given by a disc with radius ρ < R. In this case the source and the field can be computed analytically, (see [19]). Table 1 shows the errors between the exact radiating source and the source reconstructed by the approximate inverse. Uniform noise g where added to the exact data g. Different noise levels where added to the exact data. The regularization parameter γ was adjusted to the data error.
For this example α is fixed and β = π. The radius is R = 1.5 and the wavenumber k = 2π/0.9. The object f is a homogeneous disc with radius 0.5 and value 0.1.

10
International Journal of Biomedical Imaging

The nonradiating parts of the sources Φ α
Suppose that the data g α,β are given for pairwise different incident waves u α i . Then the radiating parts Φ α R of the sources Φ α can be computed, as it was shown in Section 5.1. The nonradiating parts and the corresponding fields u α NR can be represented by the sets given in Theorems 3 and 5, leading to a α m,n w m,n (rω), where the coefficients a α m,n are to be determined. To solve the problem numerically, the solution has to be calculated on a finite subspace. Two basic ideas which rely on [2] can be applied.
(1) Divide the problem into two parts. In the first step, the cost function has to be minimized by varying the scatterer f . In a second step, the cost functions for every given direction α ∈ S 1 have to be minimized by variation of the coefficients a α m,n of the nonradiating parts Φ α NR of the sources Φ α . Return to the first step until changes in the profile are less than some specified tolerance. The two steps include a linear problem.
(2) Another idea is to solve the problem in a simultaneous way. As a consequence of (14) we have that Φ α NR and Φ β NR are related by If the data (g αk ) P k=1 are given for pairwise different incident waves u αk i , then the functional has to be minimized, for example, with the Gauss-Newton algorithm.
Remark 3. The iteration can also be performed for pairwise different wavenumbers k, because the scatterer f does not depend on k (see (14)). Figure 2 shows in the first row a complex valued scatterer. The real part is given on the left. The middle row shows the reconstruction of the scatterer. The data are given accurately. An error of 5% is added on the data for the reconstruction in the last row. The parameters of the problem and the measuring instrument are selected as follows: radius R = 1, wavenumber k = 25, 128 measuring directions equidistantly distributed between 0 and 2π, and exactly as many measuring points per measuring direction. The scatterer f consists of 5 single figures. The radiating parts of the source can be computed by solving linear systems of equations. The approximate inverse showed its ability to regularize this ill-posed problems in a fast way. The main advantages are listed in Section 5.1. The high memory consumption for the general case Γ ⊂ ∂Ω is still a problem. But because the time for solving the linear problems is negligible compared to the time for solving the nonlinear problem, Remark 2 shows an acceptable way out of this problem.

NUMERICAL RESULTS
The nonradiating parts of the source were determined with the methods proposed in Section 5.2 and compared with one another. The quality of the reconstructions do not differ from one to the other. The linear algorithm needs more iterations than the simultaneous algorithm but less time per iteration. The speed strongly depends on the given parameters and the measuring setup. In fact, various approaches have been formulated, which try to solve the nonlinear inverse scattering problem with more or less approximations.
It is not the intention of the current paper to compare the performance of all reconstruction methods. It gives a better comprehension of the source-field system. Especially the solver for the linear problem has been improved by the approximate inverse.

CONCLUSION
In the first section of this paper, a special orthogonal basis of the space L 2 (Ω) for Ω a unit ball has been derived. Using this basis it is possible to separate the source in its radiating and nonradiating parts, respectively. The radiating parts of a source can be obtained by solving an ill-posed linear problem. In our paper, they are computed by means of the approximate inverse. A numerical example for a known scatterer and different noise levels shows the applicability of the method. The nonradiating parts can be calculated by nonlinear optimization. The theoretical results are verified in Section 6 by numerical experiments, indicating that the here presented method is applicable for the problem under consideration.

APPENDIX B
Some useful properties of Bessel functions are summarized in Table 2.