Approximate Image Reconstruction in Landscape Reflection Imaging

Simple reflection imaging of landscape (scenery or extended objects) poses the inverse problem of reconstructing the landscape reflectivity function from its integrals on some particular family of spheres. Such data acquisition is encoded in the framework of a Radon transform on this family of spheres. In spite of the existence of an exact inversion formula, the numerical landscape reflectivity function reconstitution is best obtained with an approximate but judiciously chosen reconstruction kernel. We describe the working of this reflection imaging modality and its theoretical handling, introduce an efficient and stable image reconstruction algorithm, and present simulation results to prove the validity of this choice as well as to demonstrate the feasibility of this imaging process.


Introduction
Imaging science is a rapidly developing field in all areas of human activity ranging from medical diagnostics to industrial nondestructive evaluation.In the last several decades, it has expanded vigorously in environmental/navigational surveillance, national security monitoring, weather forecast, hazard assessment, and so forth.It plays an essential role in remote sensing by providing information about objects or areas from a distance, typically from aircraft or satellites.By collecting data across a wide range of the electromagnetic spectrum at small spectral resolution (5-15 nm) and high spatial resolution (1-5 m), it allows detailed spectral signatures to be identified for different imaged materials.
So the aim is to obtain rapidly accurate images of large areas (or landscapes/sceneries) of the earth surface.Two main technologies have been conceived to this end: aerial or satellite photography (including television imaging) [1] and radar imaging, in particular the so-called Synthetic Aperture Radar Imaging (SAR) [2], which has the advantage of being weather independent.
The aim of this paper is not to focus on neither aerial photography nor SAR imaging and discuss their specific functioning problems.Its object is to single out an imaging concept based on the phenomena of wave reflection on more or less opaque objects and the registering of reflected wave energy by a single detector.It turns out that the appropriate mathematical description for this imaging modality is an integral transform, which is a generalization of the Radon transform, popular in medicine and industrial control.The crucial point is to show that imaging with this principle is viable and exploitable in practice.
The paper is organized as follows.Section 2 is devoted to defining active reflection imaging.It discusses the way information is recorded and used to produce images.Section 3 reviews the main mathematical tool which supports this active reflection imaging: the Radon transform on spheres centered on a plane.The next point is the derivation of an approximate reconstruction formula for the reflectivity function in Section 4. Numerical simulations and comparison comments on the results obtained by the exact and approximate reconstruction formulas are given in Section 5. A conclusion closes the papers with perspectives on possible future research directions.

Reflection Imaging
Reflection imaging is the simplest way to acquire an image of an extended object or a landscape of macroscopic dimensions.Viewing an object under the illumination of a light is the simplest example of reflection imaging.More generally if a signal (or wave pulse) is sent through space and gets reflected by an opaque surface before being recorded by an apparatus, information on the presence of this surface may be obtained provided that the signal propagation properties are known.A scanning of the object by a large number of such signals and their detection may allow us to obtain some image of this object.This is the principle of reflection imaging.
One important mode of reflection imaging known to everyone is human vision.Natural light reflected on the surface of objects passes through the eye aperture and gets projected on the retina.This way of producing an image is also used in photographic and television cameras.However there are limitations to this modality of reflection imaging when large objects or sceneries are to be imaged due to weather conditions (precipitation, fog, and clouds), variations of radiation intensity, and distance related blurring [3].Therefore it may be useful to seek an alternative reflection imaging principle.In this paper, we discuss a simple way of acquiring reflection data to obtain the image of an object.
Concretely, we will be concerned here in particular by large objects such as a landscape.We will consider first the case of flat scenery or landscape before going to the case of a hilly or structural landscape when no aperture (with retina or photographic film) is used.The first case is meant to facilitate the understanding of the mechanism of reflection imaging but is not a topic of main interest in itself.

Flat Scenery or Flat Landscape.
Let us consider a source  emitting isotropically and uniformly bursts of signals (wavepackets).This source is at first motionless and is maintained at a constant height ℎ; see Figure 1.Below is a planar landscape represented by a reflectivity function (, ), which gives the percentage of signal energy sent back by reflection at point (, ) of the plane.Let us assume that signals travel at a constant velocity  along all rectilinear trajectories in air and that  can simultaneously register returning signal energies, and this is advocated, for example, in [4], and also called monostatic mode in radar technology.An emitted burst of signals emerges from  at time  = 0 and will expand spherically around  at a distance  from  at time  = .It is clear that at time  = ℎ/, the signal burst hits the floor plane at site  and the reflected signal will be detected along the same propagation path at  at time 2ℎ/.At time  > ℎ/, the return signal at  is made of all the reflected energies at points  situated on a circle C () centered at  of radius √ () 2 + ℎ 2 .If  0 is the emitted energy flux density at , then the received reflected signal from the two-dimensional landscape is the integral of  0 (, ) on the circle C () .To keep the discussion simple, we have neglected signal spreading and attenuation along the propagation direction.So such circle integral of the reflectivity function is for the moment only a function of one variable .Collecting all such integrals will not be sufficient to find (, ), because it is a function of two variables (, ).
To overcome this problem of insufficient data, we can move the point source along a given trajectory (or curve) which will introduce a second variable: the curvilinear abscissa  of  on its trajectory.One may take the simplest trajectory possible: a straight line parallel to the landscape plane at height ℎ.In this situation and with all the stated assumptions, the reflected signal flux density is given by the integral of (, ) on the circle C (,) , whose center is at abscissa  on the orthogonal projection of the line trajectory of  on the plane and radius  R   (, ) =  0 ∫ (,)∈C (,)   (, ) . ( Here  is the integration measure of C (,) .The totality of R  (, ) for the unknown (, ) is what is called the Radon transform of (, ) on the family of circles centered on a straight line parallel to the landscape plane.This integral functional transform has been studied by many authors who have worked out the inverse transform; see, for example, [5].
In the inversion of this Radon transform on circles centered on a line, only -even part of (, ) can be reconstructed and a special scanning mode is required to obtain the full (, ).
We now extend these considerations to a hilly landscape in three dimensions.

Hilly Landscape.
A hilly landscape cast in three dimensions consists of a smooth surface in R 3 , on which there exists a reflectivity function   (, V), where (, V) are curvilinear coordinates of a point on this surface (Figure 2).We have assumed a very simple reflection mechanism in the sense that the reflected ray is always in the same direction as the incident ray so that detailed structure at the reflecting point is averaged out and represented by an effective reflectivity function (, , ) defined throughout space and possibly experiencing jump discontinuities at the separating surface with air.Following the arguments of the previous subsection, we see that the reflected signal at  after a time  following the emission of the signal burst by  is the integral of (, , ) on the sphere S () of center  and radius .In fact the integral is practically nonzero only on one side of the intersection of the surface landscape and the sphere S () .Since (, , ) has three variables, to generate the required data for reconstructing (, , ), one let  move on a twodimensional surface, which will be taken generally as a plane for convenience.The case of spheres centered on a sphere is discussed in [6].In this context, the reflection data is the integral of (, , ) on the sphere S (,,) , where (, ) are the Cartesian coordinates of  on a plane above the landscape to be imaged.This will naturally introduce the concept of a Radon transform of (, , ) on spheres centered on a plane.Such Radon transforms have been studied by many authors who have already derived an inversion formula: where  is the integration measure on the sphere S (,,) .A theoretical justification of these considerations may be found in [7] and an example of experimental realization in [8].
Here ℎ is the height of , which is so chosen so as to have the landscape above the  plane.Therefore, we assume the function of reflectivity to be supported in the space  ∈ (0, ℎ).Later on to extend the integration range on  ∈ R, we will use (, , ) = (,  − ).

2.3.
Applications.This simplified model of reflection imaging may be used in principle in many areas.One of the obvious areas is the domain of earth surveillance for environmental and meteorological purposes.In this respect it has the ingredients of Synthetic Aperture Radar (SAR) if it uses bursts of microwaves [9].However SAR can be achieved without collecting surface integral data and some of the problems of SAR are not addressed in this approach of Radon transform on spheres.Yet it may be used in any weakly reflecting medium in which strongly reflecting objects are placed assuming that signals can travel without distortion at a constant velocity [4].Of course this type of imaging can also serve as first approximation to a more involved one.
The problem addressed here consists of constructing an efficient computational algorithm based on exact inverse formulas of the Radon transform on some classes of spheres to obtain an image of reasonable quality.The idea put forward here is the use of Approximate Inverses (or mollifiers) to achieve this goal.

The Spherical Radon Transform R 𝑆 over Spheres Centered on a Plane
We now present the mathematical foundations for inverting the spherical Radon transform studied here.These results are then used in Section 4 in order to develop a reconstruction algorithm.The constant factor  0 in (2) will be discarded to simply the writing.The following notations and rules will be used.
Coordinates.Are as follows: (i) Initial or object space is R 3 Euclidean space with Cartesian coordinates (, , ).(ii) Image or Radon space is R 3 Euclidean space with Cartesian coordinates (, , ).(iii) Fourier space, dual to initial of object space, is R 3  Euclidean space with Cartesian coordinates (, , ).
Transforms.We are concerned here by integrable functions of three variables defined on various spaces of three dimensions.They are subjected to integral transforms, which may act at each step either on one, two, or three variables.An appropriate notation is introduced to clarify these actions when they are compounded.
(i) -dimensional Fourier transform F () ,  = {1, 2, 3}, where  is a -tuple in {, , }.Hankel transform, H () ] , of order ] and dimension  = {1, 2, 3} with  a -tuple in {, , }.The index  indicates which of the variables (, , ) of a function on R 3 is subjected to Fourier transform and leading to which of the dual variables (, , ) in the Fourier transformed function.For example, F (,) (, , ) is the Fourier transform of (, , ) on the first two variables (, ) of (, , ), and it is a function of (, , ): where the third variable  remains unaffected.(ii) Hankel transform of order ] H ()  ] , where the index  here has the same meaning as in Fourier transform.For example, the Hankel of order 1/2 is applied to the third variable  of (, , ), leaving the two first variables (, ) unaffected, Mathematical Problems in Engineering (iii) The Radon transform on spheres R  , its inverse R −1  , and its adjoint R †  will be given by their explicit expression in Sections 3.1, 3.2, and 3.3.

Properties. Are as follows:
(i) F (,) and H ()  ] commute as they act in different variables of R 3 .
(ii) The multiplication by an arbitrary function of a subset of variables disjoint of the set of variables of an integral transform (Fourier or Hankel) commutes with this integral transform.
3.1.Definition of R  .Let s = (,,ℎ) be the Cartesian coordinates of the center  of a sphere located at a fixed height ℎ and denote by  its radius.The associated spherical Radon transform R  which maps an integrable function (, , ) in R 3 into its projections along such a sphere is given by An explicit form of R  (, , ) in terms of spherical coordinates is where n = (sin  cos , sin  sin , cos ).This expression can be interpreted as an integral transform with a delta function kernel of the form with Such an integral differs from original results [5,10,11] since here the sources of the spheres are located at a fixed height ℎ whereas in former works the source was located on the plane crossing the origin.This choice appears more intuitive for SAR applications and makes new factors appear in the different formulae derived below.Although these factors do not change the inversion process, they are crucial in the computation of different formulae.Now we provide the important results of R  which enables us to derive our reconstruction algorithm in Section 4.
This equation consists of three integration parts.Let us examine each step separately.

Inverse Radon Transform R −1
.We now provide the inversion formulae for R  .The proofs correspond with [9] except for the parametrization of the source which changes the factors.We give then a very short version of the proofs.Theorem 2. The three-dimensional Fourier transform F (,,) (, , ) of the sought function (, , ) can be obtained as Proof.Using the same procedure compared to that for the previous proof, (6) can be rewritten in terms of Fourier transforms, assuming F (,,) (, , ) = F (,,) (, , −), as F (,,)  (, , ) . ( This integral has to be understood as a Hankel transform of order 1/2 with respect to a new variable  = 2√ 2 +  2 +  2 , where where () is the unit Heaviside step function.Taking the inverse Hankel transform, we finally get the result of the theorem.

Reconstruction Algorithm Based Approximate Inverse
We now address the problem of reconstruction for the proposed SAR modality using the Approximate Inverse [13].The last section provides essentially two reconstruction formulae given in Theorem 2 and Corollary 3. The first one involves a division by cos(2ℎ) and a Hankel transform whereas the second one uses the dual operator.Since we scan all the frequency range, the division by cos produces numerical stability; this is why we focus on the second solution for recovering .
We note the presence of the norm of the frequency vector, √ 2 +  2 +  2 , and of  which corresponds to Riesz potential or order −1 on the vector (, , ) and on the time/radius coordinate, , respectively.This step is well-known since it is involved in the Filtered Backprojection (FBP) algorithm used in conventional Computerized Tomography and requires an apodization to cut off high frequencies.
However an important difference with the FBP algorithm, here, is the position of the dual/backprojection operator.The backprojection operator which spreads the data on the corresponding manifolds (lines, spheres,. ..) is a smoothing and so regularizing operator.Since R †  is first applied on the data, its smoothing effect does not affect the filter part with Riesz potentials.
Because of the difficulty and unstability in computing directly this inversion, we propose to apply the Approximate Inverse as described below.This approach was proposed in the first place by Louis and now it uses for the threedimensional approach for the Radon transform in many fields (see [14][15][16][17]).It was successfully used in [11] for the -dimensional spherical Radon transform.But the designed reconstruction method was not suited for a 3D case essentially due to the computation time.We propose here to simplify this process to make it fit with our applications.
The Approximate Inverse method intends to reconstruct the solution in a more regular space in order to get rid of the singularities in the inversion process.Considering the studied inverse problem here, we call   a regularization method of the linear, continuous, and injective mapping plays so the role of a regularization parameter.We note ⟨⋅, ⋅⟩  2 (⋅) the inner product of the  2 (⋅)-norm.Let the linear mapping   R  be generated by linear functional In order to approximate  with   , the linear functional   (,,) will satisfy lim Then the linear mapping   is now defined by elements in  2 (E) as where   (,,) are solutions of the auxiliary problem and are called approximate reconstruction kernels by analogy with  (,,) , the reconstruction kernel associated with the problem  = R  , defined as Putting for the sake of readibility p = (, , ) and x = (, , ) and omitting the indice in the inner product,  2 (R 3 ) , Corollary 3 provides an efficient way to derive the expression of   (x) in terms of   (x) , where (p) = (1/2)(1 +  tan 2ℎ)||√ 2 +  2 +  2 .Finally we get This expression has the following advantages: (i) It can be precomputed.
(ii) It is independent of the data  and hence is not affected by noise.
(iii) The shift invariance of R  and its dual enables us to compute the approximate reconstruction kernel for few points in data space.Indeed, the integral kernel of R  in (8) satisfies Therefore,   (,,) (, , ) does not require a computation for all (, , , , , ) ∈ R 5 × R + but only for (, , , ) ∈ R 3 × R + which reduces substantially its computation time and size.
The second step is to choose an appropriate mollifier   (,,) .The mollifier and the associated approximate reconstruction kernel in [11] were designed according to the frequential properties of the inverse problem.But the designed process turns out to be incomputable in 3D.According to the theory of the Approximate Inverse, the choice of mollifier is motivated by its degree of smoothness, its invariance properties, and its properties of computation.Here we need a shift invariant function with a well-known Fourier transform and approximating the delta distribution.Table 1 presents three well-known functions satisfying these sought properties.With the problem of finite aperture of the source/detector system, measured data are incomplete and so the inversion process is severely ill-posed.For this reason, we choose the Gaussian mollifier which is the most smoothing function among these three mollifiers.In the next section, we present numerical simulations using the proposed approximate reconstruction kernel built with a Gaussian mollifier and compared to the straightforward implementation of the exact inversion formula, Corollary 3.

Simulation Results
In order to attest the efficiency and relevancy to use such a method for simple reflection imaging, we present now numerical simulations on a simple object representing a hilly landscape and built with exponential bumps on a grid of 64 × 64 × 64, as shown in Figure 3. Before proceeding the reconstruction method, one needs to produce the corresponding data.For this purpose, the implementation is performed using the alternative forward formula:  We reconstruct this object first by using the exact reconstruction formula and next with an approximate reconstruction kernel, obtained with a Gaussian function.Moreover the value of the regularization parameter  is fixed at  ≈ 3 when the normalized mean square error reaches its minimum value; see Figure 5(a).A profile of the computed approximate reconstruction kernel is depicted in Figure 5(b).Figure 6 consists of simulation results using the exact inversion formula based on Corollary 3 and using the approximate reconstruction kernel specified above with noise-free data and with speckle noise.This noise was produced with Matlab using the command imnoise (., 'speckle', 0.1).The term (, , ) involved in the exact inversion formula was apodized by a Gaussian function similarly with the mollifier chose for the Approximate Inverse method.The initial object and the associated reconstructions are represented in terms of height   slice by slice.In addition, Figure 7 displays the reconstructed 3D surfaces using the proposed method.All reconstructions were truncated by a "mask" function generated from the null set in the studied data.As expected the straightforward method succeeds in the recovery of the contours but to the price of a strong regularization.By contrast, the Approximate Inverse method appears more suited since the accuracy of the contours is higher even in presence of noise and the artifacts are far more reduced.These results reinforce the feasibility of such a technique when the Approximate Inverse method is used as pointed out in [11].

Conclusion
In this paper, we have shown in the context of reflection imaging using integral data that image reconstruction can be successfully performed using the method of approximate reconstruction integral operator.This method avoids the complicated numerical setting required by the exact reconstruction formula, which involves discontinuous operations and their composition.The computational algorithm may also include technical features such as finite aperture of the emitting antenna or the shape of the lobes of the emission pattern when signals are of electromagnetic nature.But in essence this scheme works also for pulses of waves in elastic matter as proposed by [4] on a smaller scale landscape.Thus we believe that this concept has a potential for development in the future.

Figure 4 :
Figure 4: Aspect of Radon data from an object with  on vertical axis,  = 0 and  on horizontal axis.

Figure 4
Figure4displays the slice  = 0 of the obtained measurement.Linear interpolations and trapezoidal rules were used to compute the integral along spheres.We reconstruct this object first by using the exact reconstruction formula and next with an approximate reconstruction kernel, obtained with a Gaussian function.Moreover the value of the regularization parameter  is fixed at  ≈ 3 when the normalized mean square error reaches its minimum value; see Figure5(a).A profile of the computed approximate reconstruction kernel is depicted in Figure5(b).Figure6consists of simulation results using the exact inversion formula based on Corollary 3 and using the approximate reconstruction kernel specified above with noise-free data and with speckle noise.This noise was produced with Matlab using the command imnoise (., 'speckle', 0.1).The term (, , ) involved in the exact inversion formula was apodized by a Gaussian function similarly with the mollifier chose for the Approximate Inverse method.The initial object and the associated reconstructions are represented in terms of height

Figure 6 :
Figure 6: From bottom up: column 1, level representation of landscape; column 2, its reconstruction by exact inverse formula; column 3, its reconstruction by approximate inversion kernel with  = 3; and column 4, its reconstruction by approximate inversion kernel with  = 4 and in presence of speckle noise.

Figure 7 :
Figure 7: Reconstructed 3-dimensional hilly object landscape using the proposed method.