Inversion of the elliptical Radon transform arising in migration imaging using the regular Radon transform
UM Energy Institute, University of Michigan, Ann Arbor, MI 48109, U.S.
Abstract
In recent years, many types of elliptical Radon transforms that integrate functions over various sets of ellipses/ellipsoids have been considered, relating to studies in bistatic synthetic aperture radar, ultrasound reflection tomography, radio tomography, and migration imaging. In this article, we consider the transform that integrates a given function in over a set of ellipses (when ) or ellipsoids of rotation (when ) with foci restricted to a hyperplane. We show a relation between this elliptical Radon transform and the regular Radon transform, and provide the inversion formula for the elliptical Radon transform using this relation. Numerical simulations are performed to demonstrate the suggested algorithms in two dimensions, and these simulations are also provided in this article.
Keywords: Radon transform; tomography; elliptical; migration; seismic reflection
AMS Subject Classification: 44A12; 65R10; 92C55
1 Introduction
Radontype transforms that assign to a given function its integrals over various sets of ellipses/ellipsoids arise in migration imaging under an assumption that the medium is acoustic and homogeneous. The aim of migration is to construct an image of the inside of the earth from seismic reflections recorded at its surface [9, 25]. A graphical approach called classical migration was developed systematically by Hagedoorn [11]. Classical migration had been abandoned after the waveequation method was introduced by Claerbout [5] in 1971. Gazdag and Sguazzero pointed out that the classical migration procedures that existed at that time were not based on a completely sound theory [9]. However, a correct construction for the waveequation method was often difficult to find because the experiment data did not fit into a single wave equation. To adapt the diffraction stack to borehole seismic experiments, a new approach to seismic migration was found. This approach gave classical migration a sound theory. After that discovery, classical migration has attracted many researchers in the field. The underlying idea is that seismic data in the far field can be regarded as if the data are coming from integrals of the earth’s acoustic scattering potential over surfaces determined by the velocity model [16]. These Radontype transforms relate to migration imaging as well as Bistatic Synthetic Aperture Radar (BiSAR) [2, 6, 13, 14, 32], Ultrasound Reflection Tomography (URT) [1, 10, 26], and radio tomography [27, 28, 29].
Because of these applications, there have been several papers devoted to the topic of elliptical Radon transforms. The family of ellipses with one focus fixed at the origin and the other one moving along a given line was considered in [14]. In the same paper, the family of ellipses with a fixed focal distance was also studied. The authors of [1, 10] dealt with the case of circular acquisition, when the two foci of ellipses with a given focal distance are located on a given circle. A family of ellipses with two moving foci was also handled in [6]. Radio tomography, which uses a wireless network of radio transmitters and receivers to image the distribution of attenuation within the network, was discussed in [27, 28, 29]. They approximated the obtained signal by the volume integral of the attenuation over this ellipsoid. One work [17, 18] derived two inversion formulas of this volume integral of the attenuation over this ellipsoid and studied its properties. Many works found an approximate inversion for elliptical Radon transforms.
Here we consider migration imaging and introduce a new type of an elliptical Radon transform obtained by restricting the position of the source and receiver in migration imaging. We find an explicit inversion formula for this elliptical Radon transform arising in migration imaging which is the line/area integral of the function over the ellipse/ellipsoid with foci restricted to a hyperplane.
The rest of this paper is organized as follows. The problem of interest is stated precisely and the elliptical Radon transform is formulated in section 2. In section 3, we show how to reduce the elliptical Radon transform to the regular Radon transform. The numerical simulation to demonstrate the suggested 2dimensional algorithm is presented in section 4.
2 Formulation of the problem
Let and represent 3dimensional source and receiver positions, respectively. For fixed points , an isochron surface ^{2}^{2}2Hagedoorn called this set a surface of maximum concavity (see [11]). is a surface consisting of image points associated by the travel time function , which gives the travel time between two points and with a known velocity (see [16]).
Mathematically, can be described as the set of the image points satisfying the constraint that the total travel time from the source though the image point to the receiver is constant and equal to . The isochron surface can be represented as
Seismic experiments yield data which are functions of the source position , receiver position , and time . Assuming an object function on , the data is modelled by the integral of over , i.e.,
The two detectors have nonzero sizes and time is also passing, so it is reasonable to assume that what we measure is the “average”concentrated near the location of the detectors and nearly zero sufficiently far away from the detectorsover a small region of space and a small time interval preceding the time . In mathematical terms, our data can be written as
where is the Dirac delta function. (Most works [1, 2, 6, 10, 13, 14] dealing with an elliptical Radon transform consider the arc length measure instead.) Miller, Oristaglio, and Beylkin suggested this model and approximately recovered the object function by appropriately weighted back projection of the data in [16]. However, we present the explicit formula for reconstructing using the specific minimum data of by changing the variables of positions of the two detectors.
The total dimension of the data is 7. To reduce the overdeterminacy, we assume that and are located on a hyperplane, say and also that the difference vector is parallel to a given line, say the axis. (In 2 dimensions, the difference vector is automatically parallel to the axis because the hyperplane is the line.) Hence, for any , we write and . Also, assuming a velocity is set to a constant 1, becomes the distance between and . Since satisfies , a point can be described by the ellipsoid equation
Let us choose a constant between 0 and 1. We consider the fixed value of satisfying the equation , which is a new additional requirement. Then a point can be also described by the equation
If we set , then our data become
where . Here, we restrict the positions of the source and receiver in migration imaging, say, and for . In general, reducing the space of positions for two devices is more useful and economical. If the function is odd with respect to , then is equal to zero. We thus assume the function to be even with respect to : where . We call an elliptical Radon transform, since this is the surface area integral of over the set of these ellipsoids. We generalize this dimension setup to dimension and define a more general form of an elliptical Radon transform.
Definition 2.1.
Let be given numbers, let denote the diagonal matrix with diagonal entries , and let be a locally integrable function on , even with respect to . The elliptical Radon transforms of is defined by
Note that we do not need the condition any more. If all , , are equal to 1, the elliptical Radon transform becomes the spherical Radon transform with the centers of the sphere of integration located on the hyperplane, a wellstudied transform (see [3, 4, 7, 8, 15, 12, 19, 20, 21, 22, 23, 24, 30, 31]).
3 Inversion formula
Here we assume that the object function does not touch the detectors; that is, the support of does not intersect the hyperplane where the two detectors are located.
To obtain an inversion formula for the elliptical Radon transform, we manipulate . By changing the variables , we can write
(1) 
where and . Here is the area measure of the unit sphere.
Let and . We define a map by
where is the diagonal matrix.
Proposition 3.1.

The map is a bijection with the inverse map defined by

We have
where .
The map transforms an ellipsoid into a hyperplaneâ. Changing variables using this map plays a critical role in reducing the elliptical Radon transform to the regular Radon transform.
Proof.
We can easily check that for and for , so is a bijection.
Consider for with :
(2) 
∎
We define the function on by
where . This is equivalent for , to
where . By the evenness of , we have
(3) 
Let the regular Radon transform be defined by
where and for
This can be represented by
(4) 
for and . The following theorem shows the relation between and .
Theorem 3.2.
Let be even in and have compact support in . Then we have
(5) 
Again, .
Proof.
Combining two equations (1) and (3), we have
Let us consider the variable change
The Jacobian of this transformation is , so
Since has compact support in , we have
(6) 
We recognize the right hand side as the integral along the hyperplane perpendicular to
with (signed) distance
from the origin. In this case, the measure for the hyperplane becomes Setting
in equation (4) we have the desired formula from equation (6). ∎
Using the projection slice theorem for the regular Radon transform, we obtain an analog of the projection slice theorem for the elliptical Radon transform :
Theorem 3.3.
Let be even in and have compact support in . Then we have for ,
where is the dimensional Fourier transform of .
Proof.
The projection slice theorem implies
To use this theorem, we multiply equation (5) by and integrate from to infinity:
Since has compact support in , the plane perpendicular to with distance from the origin less than does not intersect the compact support of . Hence, is equal to
where in the second line, we used Theorem 3.2. Setting and completes our proof. ∎
Taking the inverse Fourier transform of a function and using equation (3), we have the following inversion formula:
Corollary 3.4.
Let be even in and have compact support in . Then can be reconstructed as follows:
(7) 
where
4 2dimensional numerical implementation
Here we discuss the results of 2dimensional numerical implementations.
In 2dimension, equation (5) becomes
(8) 
Setting and
in equation (8), we have
where
(While is used in Theorem 3.2, is used to utilize the builtin function in Matlab. Hence a small change is required.) Again, since has compact support in , the plane perpendicular to with distance from the origin less than does not intersect the compact support of . Therefore, we have
(9) 
First of all, and were set to be 0.8 and 1, respectively. In the experiments presented here we use the phantom shown in Figure 1 (a). The phantom, supported within the rectangle , is the sum of eight characteristic functions of disks. Notice that it has to be even with respect to the axis and there are four characteristic functions of disks centered at , , , and with radii 0.2 0.15 0.05 and 0.05, whose values are 1, 0.5, 1.5, and 2 above the axis. Hence we also include their reflection below the axis. (Actually, our phantom has support in
This implies that has support in the unit ball and this makes it sufficient to consider the range in .) The images are used in Figure 1. To reconstruct the image in Figure 1 (b), we have projections for and in equation (9). All projections are computed by numerical integration. After finding the function using the usual inversion code for the regular Radon transform, we obtain the function using equation (3). (When using the inversion code for the regular Radon transform, the builtin function “iradon” in Matlab was used. The function “iradon” is the inversion of the builtin function “radon” in Matlab which considers the number of the pixels where the line passes through. Thus when computing , we also considers the number of the pixels where the ellipse passes through. We used the default version of the function in which the filter, whose aim is to deemphasize high frequencies, is set to the RamLak filter and the interpolation is set to be linear.)
While Figure 1 (b) demonstrates the image reconstructed from the exact data, Figure 1 (c) shows the absolute value of the reconstruction from noisy data. The noise is modeled by normally distributed random numbers and this is scaled so that its norm was equal to 5 of the norm of the exact data. In Figure 1 (c) the noisy data is modeled by adding the noise values scaled to 5 of the norm of the exact data to the exact data. In Figure 2, the surface plots of Figure 1 (a) and (b) are provided. Another phantom and reconstruction are shown in Figure 3. The reconstructed two dimensional data sets consisting of 256256 projections using the implemented inversion formula have been computed in less than one second (around 0.3 second) on a Intel(R) Cor(TM) i53470 CPU @ 3.20 GHz.
5 Conclusion
This paper is devoted to the study of the elliptical Radon transform arising in seismic imaging. We provide an inversion formula for the elliptical Radon transform by reducing this transform to the regular Radon transform. Also, we demonstrate our algorithm by providing numerical simulations.
Acknowledgements
The first author thanks C. Jung for fruitful discussions. The first author was supported in part by Basic Science Research Program through the National Research Foundation of Korea(NRF) funded by the Ministry of science, ICT and future planning (2015R1C1A1A02037662)).
References
 [1] G. Ambartsoumian, J. Boman, V.P. Krishnan, and E.T. Quinto. Microlocal analysis of an ultrasound transform with circular source and receiver trajectories. Contemporary Mathematics, 598, 2013.
 [2] G. Ambartsoumian, R. Felea, V.P. Krishnan, C. Nolan, and E.T. Quinto. A class of singular Fourier integral operators in Synthetic Aperture Radar imaging. Journal of Functional Analysis, 264(1):246 – 269, 2013.
 [3] L. Andersson. On the determination of a function from spherical averages. SIAM Journal on Mathematical Analysis, 19(1):214–232, 1988.
 [4] A. Beltukov. Inversion of the spherical mean transform with sources on a hyperplane. arXiv:0919.1380v1, 2009.
 [5] J. Claerbout. Toward a unified theory of reflector mapping. GEOPHYSICS, 36(3):467–481, 1971.
 [6] J.D. Coker and A.H. Tewfik. Multistatic SAR image reconstruction based on an ellipticalgeometry Radon transform. In Waveform Diversity and Design Conference, pages 204 –208, June 2007.
 [7] R. Courant and D. Hilbert. Methods of Mathematical Physics. Number V. 2 in Methods of Mathematical Physics. Interscience Publishers, 1962.
 [8] J. Fawcett. Inversion of ndimensional spherical averages. SIAM Journal on Applied Mathematics, 45(2):336–341, 1985.
 [9] J. Gazdag and P. Sguazzero. Migration of seismic data. Proceedings of the IEEE, 72(10):1302–1315, 1984.
 [10] R. GouiaZarrad and G. Ambartsoumian. Approximate inversion algorithm of the elliptical Radon transform. In Mechatronics and its Applications (ISMA), 2012 8th International Symposium on, pages 1 –4, April 2012.
 [11] J. G. Hagedoorn. a process of seismic reflection interpretation*. Geophysical Prospecting, 2(2):85–127, 1954.
 [12] J. Klein. Inverting the spherical Radon transform for physically meaningful functions. arXiv:math/0307348, 2003.
 [13] V. Krishnan and E.T. Quinto. Microlocal aspects of common offset synthetic aperture radar imaging. Inverse Problems and Imaging, 5:659–674, 2011.
 [14] V.P. Krishnan, H. Levinson, and E.T. Quinto. Microlocal analysis of elliptical Radon transforms with foci on a line. In Irene Sabadini and Daniele C Struppa, editors, The mathematical legacy of Leon Ehrenpreis, volume 16 of Springer Proceedings in Mathematics, pages 163–182. Springer Milan, 2012.
 [15] M.M. Lavrentiev, V.G. Romanov, and V.G. Vasiliev. Multidimensional inverse problems for differential equations. Lecture notes in mathematics. SpringerVerlag, 1970.
 [16] D. Miller, M. Oristaglio, and G. Beylkin. A new slant on seismic imaging: Migration and integral geometry. Geophysics, 52(7):943–964, 1987.
 [17] S. Moon. Properties of some integral transforms arising in tomography. Dissertation, Texas A&M University, 2013.
 [18] S. Moon. On the determination of a function from an elliptical Radon transform. Journal of Mathematical Analysis and Applications, 416(2):724 – 734, 2014.
 [19] E.K. Narayanan and Rakesh. Spherical means with centers on a hyperplane in even dimensions. Inverse Problems, 26(3):035014, 2010.
 [20] F. Natterer and F. Wübbeling. Mathematical methods in image reconstruction. SIAM Monographs on mathematical modeling and computation. SIAM, Society of industrial and applied mathematics, Philadelphia (Pa.), 2001.
 [21] S. Nilsson. Application of fast backprojection techniques for some inverse problems of integral geometry. Linköping studies in science and technology: Dissertations. Department of Mathematics, Linköping University, 1997.
 [22] S. J. Norton. Reconstruction of a reflectivity field from line integrals over circular paths. The Journal of the Acoustical Society of America, 67(3):853–863, 1980.
 [23] V. P. Palamodov. Reconstructive integral geometry, volume 98. Springer, 2004.
 [24] N.T. Redding and G.N. Newsam. Inverting the circular Radon transform. DTSO Research Report DTSORu0211, August 2001.
 [25] E.A. Robinson. Migration of geophysical data. International Human Resources Development Corp., 1983.
 [26] S. Roy, V. P. Krishnan, P. Chandrashekar, and AS V. Murthy. An efficient numerical algorithm for the inversion of an integral transform arising in ultrasound imaging. Journal of Mathematical Imaging and Vision, pages 1–14, 2014.
 [27] J. Wilson and N. Patwari. Throughwall tracking using variancebased radio tomography networks. ArXiv eprints, September 2009.
 [28] J. Wilson and N. Patwari. Radio tomographic imaging with wireless networks. IEEE Transactions on Mobile Computing, 9(5):621 –632, May 2010.
 [29] J. Wilson, N. Patwari, and F.G. Vasquesz. Regularization methods for radio tomographic imaging. Virginia Tech Symposium on Wireless Personal Communications, June 2009.
 [30] C.E. Yarman and B. Yazici. Inversion of circular averages using the Funk transform. In Acoustics, Speech and Signal Processing, 2007. ICASSP 2007. IEEE International Conference on, volume 1, pages I–541–I–544, 2007.
 [31] C.E. Yarman and B. Yazici. Inversion of the circular averages transform using the Funk transform. Inverse Problems, 27(6):065001, 2011.
 [32] C.E. Yarman, B. Yazici, and M. Cheney. Bistatic synthetic aperture radar imaging for arbitrary flight trajectories. Image Processing, IEEE Transactions on, 17(1):84–93, Jan 2008.