#### Abstract

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.

#### 1. 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.

#### 2. 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.

##### 2.1. Flat Scenery or Flat Landscape

Let us consider a source emitting isotropically and uniformly bursts of signals (wave-packets). 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 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 . At time , the return signal at is made of all the reflected energies at points situated on a circle centered at of radius . If is the emitted energy flux density at , then the received reflected signal from the two-dimensional landscape is the integral of on the circle . 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 , whose center is at abscissa on the orthogonal projection of the line trajectory of on the plane and radius Here is the integration measure of . The totality of 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.

##### 2.2. Hilly Landscape

A hilly landscape cast in three dimensions consists of a smooth surface in , on which there exists a reflectivity function , where 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 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 . Since has three variables, to generate the required data for reconstructing , one let move on a two-dimensional 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 , 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 . 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 . Later on to extend the integration range on , 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.

#### 3. The Spherical Radon Transform 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 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 Euclidean space with Cartesian coordinates .(ii)Image or Radon space is Euclidean space with Cartesian coordinates .(iii)Fourier space, dual to initial of object space, is 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 , , where is a -tuple in . Hankel transform, , of order and dimension with a -tuple in . The index indicates which of the variables of a function on is subjected to Fourier transform and leading to which of the dual variables in the Fourier transformed function. For example, 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 , where the index here has the same meaning as in Fourier transform. For example, the Hankel of order is applied to the third variable of , leaving the two first variables unaffected, (iii)The Radon transform on spheres , its inverse , and its adjoint will be given by their explicit expression in Sections 3.1, 3.2, and 3.3.

*Properties*. Are as follows:(i) and commute as they act in different variables of .(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

Let 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 which maps an integrable function in into its projections along such a sphere is given by An explicit form of in terms of spherical coordinates is where . 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 which enables us to derive our reconstruction algorithm in Section 4.

##### 3.2. Adjoint Transform

The dual operator is defined by the same kernel but the integration is carried on functions defined on . Let be an integrable function on , and then

Proposition 1. *The three-dimensional Fourier transform of is given by **where .*

*Proof. *From (9) and using Fubini’s theorem, the Fourier transform of is We proceed to the change of variables with in (11). Since , the right hand side becomes This equation consists of three integration parts. Let us examine each step separately. (i)The angular integration can be calculated in terms of the Bessel function using formula 4.624 of [12] (ii)The -integral becomes then and can be interpreted as a Hankel transform of order on the function , (iii)The last integration over is a two-dimensional Fourier transform on which completes the proof using .

##### 3.3. Inverse Radon Transform

We now provide the inversion formulae for . 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 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 , as This integral has to be understood as a Hankel transform of order with respect to a new variable , where where is the unit Heaviside step function. Taking the inverse Hankel transform, we finally get the result of the theorem.

Corollary 3. *Using Proposition 1, the three-dimensional Fourier transform of may be expressed in terms of as *

*Proof. *Since the multiplication by commutes with , (17) can be rewritten as By the same way, (10) becomes The proof is completed by taking .

#### 4. 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 and a Hankel transform whereas the second one uses the dual operator. Since we scan all the frequency range, the division by 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, , and of which corresponds to Riesz potential or order 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 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 three-dimensional approach for the Radon transform in many fields (see [14–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 when is continuous and satisfies plays so the role of a regularization parameter. We note the inner product of the -norm. Let the linear mapping be generated by linear functional as In order to approximate with , the linear functional will satisfy Then the linear mapping is now defined by elements in as where are solutions of the auxiliary problem and are called approximate reconstruction kernels by analogy with , the reconstruction kernel associated with the problem , defined as Putting for the sake of readibility and and omitting the indice in the inner product, , Corollary 3 provides an efficient way to derive the expression of in terms of , where . 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 and its dual enables us to compute the approximate reconstruction kernel for few points in data space. Indeed, the integral kernel of in (8) satisfies Therefore, does not require a computation for all but only for 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.

#### 5. 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 , 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: The parameters are as follows: (i),(ii),(iii),(iv),(v).

Figure 4 displays the slice 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 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].

**(a)**

**(b)**

**(a) With noise-free data**

**(b) With speckle noise**

#### 6. 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.

#### Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.