# Approximate Image Reconstruction in Landscape Reflection Imaging

Authors: Rémi Régnier; Gaël Rigaud; Maï K. Nguyen
Journal: Mathematical Problems in Engineering
 (2015)

---

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


---

## Body

## 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. Section2 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 sourceS emitting isotropically and uniformly bursts of signals (wave-packets). This source is at first motionless and is maintained at a constant height h; see Figure 1. Below is a planar landscape represented by a reflectivity function f(x,y), which gives the percentage of signal energy sent back by reflection at point (x,y) of the plane. Let us assume that signals travel at a constant velocity c along all rectilinear trajectories in air and that S 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 S at time t=0 and will expand spherically around S at a distance cτ from S at time t=τ. It is clear that at time t=h/c, the signal burst hits the floor plane at site H and the reflected signal will be detected along the same propagation path at S at time 2h/c. At time t>h/c, the return signal at S is made of all the reflected energies at points M situated on a circle C(t) centered at H of radius (ct)2+h2. If I0 is the emitted energy flux density at S, then the received reflected signal from the two-dimensional landscape is the integral of I0f(x,y) on the circle C(t). 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 t. Collecting all such integrals will not be sufficient to find f(x,y), because it is a function of two variables (x,y). 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 S on its trajectory. One may take the simplest trajectory possible: a straight line parallel to the landscape plane at height h. In this situation and with all the stated assumptions, the reflected signal flux density is given by the integral of f(x,y) on the circle C(ξ,t), whose center is at abscissa ξ on the orthogonal projection of the line trajectory of S on the plane and radius ct(1)RCfξ,t=I0∫x,y∈Cξ,tdsfx,y.Here ds is the integration measure of C(ξ,t). The totality of RCf(ξ,t) for the unknown f(x,y) is what is called the Radon transform of f(x,y) 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 y-even part of f(x,y) can be reconstructed and a special scanning mode is required to obtain the full f(x,y).Figure 1
Flat landscape-scenery imaging.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 inR3, on which there exists a reflectivity function fs(u,v), where (u,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 f(x,y,z) 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 S after a time t following the emission of the signal burst by S is the integral of f(x,y,z) on the sphere S(t) of center S and radius ct. In fact the integral is practically nonzero only on one side of the intersection of the surface landscape and the sphere S(t). Since f(x,y,z) has three variables, to generate the required data for reconstructing f(x,y,z), one let S 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 f(x,y,z) on the sphere S(ξ,η,t), where (ξ,η) are the Cartesian coordinates of S on a plane above the landscape to be imaged. This will naturally introduce the concept of a Radon transform of f(x,y,z) on spheres centered on a plane. Such Radon transforms have been studied by many authors who have already derived an inversion formula: (2)RSfξ,η,t=I0∫x,y,z∈Sξ,η,tdafx,y,z,where da is the integration measure on the sphere S(ξ,η,t). A theoretical justification of these considerations may be found in [7] and an example of experimental realization in [8].Figure 2
Hilly landscape-scenery imaging.Hereh is the height of S, which is so chosen so as to have the landscape above the xOy plane. Therefore, we assume the function of reflectivity to be supported in the space z∈(0,h). Later on to extend the integration range on z∈R, we will use f(x,y,z)=f(x,y-z).


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


## 2.1. Flat Scenery or Flat Landscape

Let us consider a sourceS emitting isotropically and uniformly bursts of signals (wave-packets). This source is at first motionless and is maintained at a constant height h; see Figure 1. Below is a planar landscape represented by a reflectivity function f(x,y), which gives the percentage of signal energy sent back by reflection at point (x,y) of the plane. Let us assume that signals travel at a constant velocity c along all rectilinear trajectories in air and that S 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 S at time t=0 and will expand spherically around S at a distance cτ from S at time t=τ. It is clear that at time t=h/c, the signal burst hits the floor plane at site H and the reflected signal will be detected along the same propagation path at S at time 2h/c. At time t>h/c, the return signal at S is made of all the reflected energies at points M situated on a circle C(t) centered at H of radius (ct)2+h2. If I0 is the emitted energy flux density at S, then the received reflected signal from the two-dimensional landscape is the integral of I0f(x,y) on the circle C(t). 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 t. Collecting all such integrals will not be sufficient to find f(x,y), because it is a function of two variables (x,y). 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 S on its trajectory. One may take the simplest trajectory possible: a straight line parallel to the landscape plane at height h. In this situation and with all the stated assumptions, the reflected signal flux density is given by the integral of f(x,y) on the circle C(ξ,t), whose center is at abscissa ξ on the orthogonal projection of the line trajectory of S on the plane and radius ct(1)RCfξ,t=I0∫x,y∈Cξ,tdsfx,y.Here ds is the integration measure of C(ξ,t). The totality of RCf(ξ,t) for the unknown f(x,y) is what is called the Radon transform of f(x,y) 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 y-even part of f(x,y) can be reconstructed and a special scanning mode is required to obtain the full f(x,y).Figure 1
Flat landscape-scenery imaging.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 inR3, on which there exists a reflectivity function fs(u,v), where (u,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 f(x,y,z) 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 S after a time t following the emission of the signal burst by S is the integral of f(x,y,z) on the sphere S(t) of center S and radius ct. In fact the integral is practically nonzero only on one side of the intersection of the surface landscape and the sphere S(t). Since f(x,y,z) has three variables, to generate the required data for reconstructing f(x,y,z), one let S 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 f(x,y,z) on the sphere S(ξ,η,t), where (ξ,η) are the Cartesian coordinates of S on a plane above the landscape to be imaged. This will naturally introduce the concept of a Radon transform of f(x,y,z) on spheres centered on a plane. Such Radon transforms have been studied by many authors who have already derived an inversion formula: (2)RSfξ,η,t=I0∫x,y,z∈Sξ,η,tdafx,y,z,where da is the integration measure on the sphere S(ξ,η,t). A theoretical justification of these considerations may be found in [7] and an example of experimental realization in [8].Figure 2
Hilly landscape-scenery imaging.Hereh is the height of S, which is so chosen so as to have the landscape above the xOy plane. Therefore, we assume the function of reflectivity to be supported in the space z∈(0,h). Later on to extend the integration range on z∈R, we will use f(x,y,z)=f(x,y-z).


## 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 TransformRS 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 Section4 in order to develop a reconstruction algorithm. The constant factor I0 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 isR3 Euclidean space with Cartesian coordinates (x,y,z).(ii)
Image or Radon space isR3 Euclidean space with Cartesian coordinates (ξ,η,ρ).(iii)
Fourier space, dual to initial of object space, isR3 Euclidean space with Cartesian coordinates (p,q,k).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)
n-dimensional Fourier transform F(j), n={1,2,3}, where j is a n-tuple in {x,y,z}. Hankel transform, Hν(j), of order ν and dimension n={1,2,3} with j a n-tuple in {x,y,z}. The index j indicates which of the variables (x,y,z) of a function on R3 is subjected to Fourier transform and leading to which of the dual variables (p,q,k) in the Fourier transformed function. For example, F(x,y)f(p,q,z) is the Fourier transform of f(x,y,z) on the first two variables (x,y) of f(x,y,z), and it is a function of (p,q,z): (3)Fx,yfp,q,z=∫R2dxdye-2iπpx+qyfx,y,z,where the third variablez remains unaffected.(ii)
Hankel transform of orderνHν(j), where the index j here has the same meaning as in Fourier transform.For example, the Hankel of order1/2 is applied to the third variable z of f(x,y,z), leaving the two first variables (x,y) unaffected, (4)H1/2zfx,y,γ=∫0∞dzzγJ1/2zγfx,y,z.(iii)
The Radon transform on spheresRS, its inverse RS-1, and its adjoint RS† will be given by their explicit expression in Sections 3.1, 3.2, and 3.3.Properties. Are as follows:(i)
F(x,y) and Hν(z) commute as they act in different variables of R3.(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 ofRS

Lets=(ξ,η,h) be the Cartesian coordinates of the center S of a sphere located at a fixed height h and denote by ρ its radius. The associated spherical Radon transform RS which maps an integrable function f(x,y,z) in R3 into its projections along such a sphere is given by (5)RSfξ,η,ρ=∫n∈S2fs+ρndn.An explicit form of RSf(ξ,η,ρ) in terms of spherical coordinates is (6)RSfξ,η,ρ=ρ2∫02πdϕ∫0πsin⁡θdθfξ+ρsin⁡θcos⁡ϕ,η+ρsin⁡θsin⁡ϕ,h+ρcos⁡θ,where n=(sin⁡θcos⁡ϕ,sin⁡θsin⁡ϕ,cos⁡θ). This expression can be interpreted as an integral transform with a delta function kernel of the form (7)RSfξ,η,ρ=∫R3dxdydzKSξ,ηρ∣x,y,zfx,y,z,with (8)KSξ,η,ρ∣x,y,z=δρ-x-ξ2+y-η2+z-h2.Such an integral differs from original results [5, 10, 11] since here the sources of the spheres are located at a fixed height h 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 ofRS which enables us to derive our reconstruction algorithm in Section 4.


### 3.2. Adjoint TransformRS†

The dual operator is defined by the same kernel but the integration is carried on functions defined on(ξ,η,ρ)∈E=R2×R+. Let h(ξ,η,ρ) be an integrable function on E, and then (9)
R
S
†
h
x
,
y
,
z
=
∫
E
d
ξ
d
η
d
ρ
δ
ρ
-
x
-
ξ
2
+
y
-
η
2
+
z
-
h
2
h
ξ
,
η
,
ρ
.Proposition 1.
The three-dimensional Fourier transform ofRS†h(x,y,z) is given by (10)Fx,y,zRS†hp,q,k=2πe-2iπkhp2+q2+k2·H1/2ρFξ,ηh_p,q,2πp2+q2+k2,where h_(ξ,η,ρ)=ρh(ξ,η,ρ).Proof.
From (9) and using Fubini’s theorem, the Fourier transform of R† is (11)
F
x
,
y
,
z
R
S
†
h
p
,
q
,
k
=
∫
E
∫
R
3
d
ξ
d
η
d
ρ
d
x
d
y
d
z
e
-
2
i
π
p
x
+
q
y
+
k
z
K
S
ξ
,
η
,
ρ
∣
x
,
y
,
z
h
ξ
,
η
,
ρ
.We proceed to the change of variables (12)x=ξ+ρsin⁡θcos⁡ϕ,y=η+ρsin⁡θsin⁡ϕ,z=h+ρcos⁡θ,with dxdydz=ρ2sin⁡θdρdθdϕ in (11). Since (x-ξ)2+(y-η)2+(z-h)2=ρ, the right hand side becomes (13)∫R2e-2iπpξ+qηdξdη∫R3ρ2sin⁡θdρdθdϕe-2iπpx+qy+kze-2iπρpsin⁡θcos⁡ϕ+qsin⁡θsin⁡ϕ+kcos⁡θhξ,η,ρ.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 functionJ1/2(x) using formula 4.624 of [12] (14)∫S2sin⁡θdθdϕe-2iπρpsin⁡θcos⁡ϕ+qsin⁡θsin⁡ϕ+kcos⁡θ=2ρπρp2+q2+k21/4J1/22πρp2+q2+k2.(ii)
Theρ-integral becomes then (15)∫0∞dρρJ1/22πρp2+q2+k2ρhξ,η,ρand can be interpreted as a Hankel transform of order1/2 on the function h_(ξ,η,ρ)=ρh(ξ,η,ρ), (16)12π1p2+q2+k21/4·H1/2ρh_ξ,η,2πp2+q2+k2.(iii)
The last integration overdξdη is a two-dimensional Fourier transform on h_(ξ,η,2πp2+q2+k2) which completes the proof using F(ξ,η)H1/2(ρ)=H1/2(ρ)F(ξ,η).


### 3.3. Inverse Radon TransformRS-1

We now provide the inversion formulae forRS. 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 transformF(x,y,z)f(p,q,k) of the sought function f(x,y,z) can be obtained as (17)Fx,y,zfp,q,k=π2kcos⁡2πkhH1/2ρFξ,η·RSfρp,q,2πp2+q2+k2.Proof.
Using the same procedure compared to that for the previous proof, (6) can be rewritten in terms of Fourier transforms, assuming F(x,y,z)f(p,q,k)=F(x,y,z)f(p,q,-k), as (18)Fξ,ηRSfp,q,ρ=4πρ3/2∫0∞dkcos⁡2πkh·J1/22πρp2+q2+k2p2+q2+k21/4Fx,y,zfp,q,k.This integral has to be understood as a Hankel transform of order 1/2 with respect to a new variable t=2πp2+q2+k2, (19)Fξ,ηRSfξ,η,ρρ=∫0∞dttρJ1/2tρgp,q,t=H1/2tgp,q,ρ,where (20)gp,q,t=22πYt-2πp2+q2·Fx,y,zfp,q,t24π2-p2+q2·cos⁡ht2-4π2p2+q2t2-4π2p2+q2,where Y(k) is the unit Heaviside step function. Taking the inverse Hankel transform, we finally get the result of the theorem.Corollary 3.
Using Proposition1, the three-dimensional Fourier transform of f(x,y,z) may be expressed in terms of RS† as (21)Fx,y,zfp,q,k=k21+itan⁡2πhk·p2+q2+k2Fx,y,zRS†·RSfρ2p,q,k.Proof.
Since the multiplication by1/ρ commutes with F(ξ,η), (17) can be rewritten as (22)Fx,y,zfp,q,k=π2kcos⁡2πkhH1/2ρFξ,η1ρ·RSfp,q,2πp2+q2+k2.By the same way, (10) becomes (23)Fx,y,zRS†hp,q,k=2πe-2iπkhp2+q2+k2·H1/2ρFξ,ηρhp,q,2πp2+q2+k2.The proof is completed by taking ρh(ξ,η,ρ)=1/ρRSf(ξ,η,ρ).


## 3.1. Definition ofRS

Lets=(ξ,η,h) be the Cartesian coordinates of the center S of a sphere located at a fixed height h and denote by ρ its radius. The associated spherical Radon transform RS which maps an integrable function f(x,y,z) in R3 into its projections along such a sphere is given by (5)RSfξ,η,ρ=∫n∈S2fs+ρndn.An explicit form of RSf(ξ,η,ρ) in terms of spherical coordinates is (6)RSfξ,η,ρ=ρ2∫02πdϕ∫0πsin⁡θdθfξ+ρsin⁡θcos⁡ϕ,η+ρsin⁡θsin⁡ϕ,h+ρcos⁡θ,where n=(sin⁡θcos⁡ϕ,sin⁡θsin⁡ϕ,cos⁡θ). This expression can be interpreted as an integral transform with a delta function kernel of the form (7)RSfξ,η,ρ=∫R3dxdydzKSξ,ηρ∣x,y,zfx,y,z,with (8)KSξ,η,ρ∣x,y,z=δρ-x-ξ2+y-η2+z-h2.Such an integral differs from original results [5, 10, 11] since here the sources of the spheres are located at a fixed height h 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 ofRS which enables us to derive our reconstruction algorithm in Section 4.


## 3.2. Adjoint TransformRS†

The dual operator is defined by the same kernel but the integration is carried on functions defined on(ξ,η,ρ)∈E=R2×R+. Let h(ξ,η,ρ) be an integrable function on E, and then (9)
R
S
†
h
x
,
y
,
z
=
∫
E
d
ξ
d
η
d
ρ
δ
ρ
-
x
-
ξ
2
+
y
-
η
2
+
z
-
h
2
h
ξ
,
η
,
ρ
.Proposition 1.
The three-dimensional Fourier transform ofRS†h(x,y,z) is given by (10)Fx,y,zRS†hp,q,k=2πe-2iπkhp2+q2+k2·H1/2ρFξ,ηh_p,q,2πp2+q2+k2,where h_(ξ,η,ρ)=ρh(ξ,η,ρ).Proof.
From (9) and using Fubini’s theorem, the Fourier transform of R† is (11)
F
x
,
y
,
z
R
S
†
h
p
,
q
,
k
=
∫
E
∫
R
3
d
ξ
d
η
d
ρ
d
x
d
y
d
z
e
-
2
i
π
p
x
+
q
y
+
k
z
K
S
ξ
,
η
,
ρ
∣
x
,
y
,
z
h
ξ
,
η
,
ρ
.We proceed to the change of variables (12)x=ξ+ρsin⁡θcos⁡ϕ,y=η+ρsin⁡θsin⁡ϕ,z=h+ρcos⁡θ,with dxdydz=ρ2sin⁡θdρdθdϕ in (11). Since (x-ξ)2+(y-η)2+(z-h)2=ρ, the right hand side becomes (13)∫R2e-2iπpξ+qηdξdη∫R3ρ2sin⁡θdρdθdϕe-2iπpx+qy+kze-2iπρpsin⁡θcos⁡ϕ+qsin⁡θsin⁡ϕ+kcos⁡θhξ,η,ρ.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 functionJ1/2(x) using formula 4.624 of [12] (14)∫S2sin⁡θdθdϕe-2iπρpsin⁡θcos⁡ϕ+qsin⁡θsin⁡ϕ+kcos⁡θ=2ρπρp2+q2+k21/4J1/22πρp2+q2+k2.(ii)
Theρ-integral becomes then (15)∫0∞dρρJ1/22πρp2+q2+k2ρhξ,η,ρand can be interpreted as a Hankel transform of order1/2 on the function h_(ξ,η,ρ)=ρh(ξ,η,ρ), (16)12π1p2+q2+k21/4·H1/2ρh_ξ,η,2πp2+q2+k2.(iii)
The last integration overdξdη is a two-dimensional Fourier transform on h_(ξ,η,2πp2+q2+k2) which completes the proof using F(ξ,η)H1/2(ρ)=H1/2(ρ)F(ξ,η).


## 3.3. Inverse Radon TransformRS-1

We now provide the inversion formulae forRS. 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 transformF(x,y,z)f(p,q,k) of the sought function f(x,y,z) can be obtained as (17)Fx,y,zfp,q,k=π2kcos⁡2πkhH1/2ρFξ,η·RSfρp,q,2πp2+q2+k2.Proof.
Using the same procedure compared to that for the previous proof, (6) can be rewritten in terms of Fourier transforms, assuming F(x,y,z)f(p,q,k)=F(x,y,z)f(p,q,-k), as (18)Fξ,ηRSfp,q,ρ=4πρ3/2∫0∞dkcos⁡2πkh·J1/22πρp2+q2+k2p2+q2+k21/4Fx,y,zfp,q,k.This integral has to be understood as a Hankel transform of order 1/2 with respect to a new variable t=2πp2+q2+k2, (19)Fξ,ηRSfξ,η,ρρ=∫0∞dttρJ1/2tρgp,q,t=H1/2tgp,q,ρ,where (20)gp,q,t=22πYt-2πp2+q2·Fx,y,zfp,q,t24π2-p2+q2·cos⁡ht2-4π2p2+q2t2-4π2p2+q2,where Y(k) is the unit Heaviside step function. Taking the inverse Hankel transform, we finally get the result of the theorem.Corollary 3.
Using Proposition1, the three-dimensional Fourier transform of f(x,y,z) may be expressed in terms of RS† as (21)Fx,y,zfp,q,k=k21+itan⁡2πhk·p2+q2+k2Fx,y,zRS†·RSfρ2p,q,k.Proof.
Since the multiplication by1/ρ commutes with F(ξ,η), (17) can be rewritten as (22)Fx,y,zfp,q,k=π2kcos⁡2πkhH1/2ρFξ,η1ρ·RSfp,q,2πp2+q2+k2.By the same way, (10) becomes (23)Fx,y,zRS†hp,q,k=2πe-2iπkhp2+q2+k2·H1/2ρFξ,ηρhp,q,2πp2+q2+k2.The proof is completed by taking ρh(ξ,η,ρ)=1/ρRSf(ξ,η,ρ).


## 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 cos⁡(2πkh) 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 f.We note the presence of the norm of the frequency vector,p2+q2+k2, and of k 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 RS† 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 n-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 callSγ a regularization method of the linear, continuous, and injective mapping RS:L2(R3)→L2(E) when Sγ:L2(E)→L2(R3) is continuous and satisfies (24)limγ→0⁡SγRSf=f,∀f∈L2R3.γ plays so the role of a regularization parameter. We note ·,·L2(·) the inner product of the L2(·)-norm. Let the linear mapping SγRS be generated by linear functional δ(x,y,z)γ∈L2(R3) as (25)fγx,y,z=SγAfx,y,z=δx,y,zγ,fL2R3.In order to approximate f with fγ, the linear functional δ(x,y,z)γ will satisfy (26)limγ→0⁡δx,y,zγ,fL2R3=fx,y,z.Then the linear mapping Sγ is now defined by elements in L2(E) as (27)Sγgx,y,z=ψx,y,zγ,gL2E,where ψ(x,y,z)γ are solutions of the auxiliary problem (28)RS†ψx,y,zγ=δx,y,zγand are called approximate reconstruction kernels by analogy with ψ(x,y,z), the reconstruction kernel associated with the problem g=RSf, defined as (29)fx,y,z=∫Eψx,y,zξ,η,ρgξ,η,ρdξdηdρ=g,ψx,y,zL2E.Putting for the sake of readibility p=(p,q,k) and x=(x,y,z) and omitting the indice in the inner product, L2(R3), Corollary 3 provides an efficient way to derive the expression of ψ(x)γ in terms of δ(x)γ, (30)fγ=f,δx′γ=wpFxRS†RSfρ2p,Fxδx′γp=RS†RSfρ2x,FpwpFxδx′γx=RSf,ρ-2RSFpwpFxδx′γL2E=RSf,ψx′γL2E,where w(p)=1/2(1+itan⁡2πhk)|k|p2+q2+k2. Finally we get (31)ψxγ=ρ-2RSFpwpFxδxγ.This expression has the following advantages: (i)
It can be precomputed.(ii)
It is independent of the datag and hence is not affected by noise.(iii)
The shift invariance ofRS and its dual enables us to compute the approximate reconstruction kernel for few points in data space. Indeed, the integral kernel of RS in (8) satisfies (32)KSξ,η,ρ∣x-x′,y-y′,z=KSξ+x′,η+y′,ρ∣x,y,z.Therefore,ψ(x,y,z)γ(ξ,η,ρ) does not require a computation for all (x,y,z,ξ,η,ρ)∈R5×R+ but only for (z,ξ,η,ρ)∈R3×R+ which reduces substantially its computation time and size.The second step is to choose an appropriate mollifierδ(x,y,z)γ. 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.Table 1
Some approximate Dirac delta functions and their Fourier transforms.
δ
γ
r
F
r
δ
γ
p
(
2
π
)
(
-
3
/
2
)
γ
-
3
e
-
|
r
|
2
/
2
γ
2  (Gaussian function)
1
2
π
γ
2
e
-
(
γ
2
|
p
|
2
/
2
)  (Gaussian function)
χ
0,1
r
/
γ
1
-
|
r
|
2
/
γ
2 χ[0,1]: characteristic function
-
4
π
J
1
γ
p
|
p
|  (J1x: 1st kind of Bessel function of order 1)
γ
π
3
sin
⁡
γ
|
r
|
γ
|
r
| (cardinal sine function)
γ
2
π
3
Π
-
γ
-
1
,
γ
-
1
p Π-γ-1,γ-1x: gate function


## 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 of64×64×64, as shown in Figure 3.Figure 3
Three-dimensional hilly object landscape.Before proceeding the reconstruction method, one needs to produce the corresponding data. For this purpose, the implementation is performed using the alternative forward formula:(33)
R
S
ξ
,
η
,
ρ
=
∫
0
π
∫
-
ρ
ρ
f
ξ
+
ρ
2
-
z
2
cos
⁡
θ
,
η
+
ρ
2
-
z
2
sin
⁡
θ
,
h
-
z
ρ
d
z
d
θ
.The parameters are as follows: (i)
h=N=64,(ii)
ξ=η=-2N,-2N+1,…,2N,(iii)
ρ=0,1/2,…,3N,(iv)
θ=0,2π/(5N),…,2π,(v)
z=0,1,…,h-1.Figure4 displays the slice η=0 of the obtained measurement. Linear interpolations and trapezoidal rules were used to compute the integral along spheres.Figure 4
Aspect of Radon data from an object withξ on vertical axis, η=0 and ρ on horizontal axis.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 commandimnoise (.,‘speckle’, 0.1). The term w(p,q,k) 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].Figure 5
(a) Mean quadratic relative error with respect toγ. (b) Profile of the computed approximate reconstruction kernel ψ(0,0,h)γ(0,0,ρ) with γ=3.
(a)
(b)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
Reconstructed 3-dimensional hilly object landscape using the proposed method.
(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.



---

*Source: 268295-2015-10-19.xml*