On the τ-p Transform and Seismic Data

In this work, an improved algebraic reconstruction for Radon transform of seismic type is proposed. Namely, the τ-p transform. The discrete formulations of this reconstruction problem are accomplished via the standard square function on the 2-D region [0,1]×[0,1], since its transform is known exactly. The interior-point algorithm approach with some MATLAB optimization functions is used for solving the system and recovering the image. The new algorithm produces accurate reconstructed images and is suitable for the case of limitation on the number of available projections.


Introduction
In this work, an improved algebraic reconstruction for Radon transform of seismic type is proposed.Namely, for the τ -p transform.Our approach relates this form of Radon transform to the traditional form that was originally proposed by Radon [1][2][3][4].Radon transform provides the mathematical framework for a large class of reconstruction problems, and there are several types of this transform that fits different applications, nature of data, or geometry of measurements.The different forms of this transform are involved with a variety of applications in seismic data, digital rock physics, radar imaging, medical modalities such as computed tomography, magnetic resonance, and others [5][6][7][8][9][10][11][12][13][14][15][16].Fundamentally, the Radon transform of a function f on 2-D region for example is the line integral of f along all possible lines, or some curves.Ultimately, of course, the aim is to obtain some structural information or some properties of the object in study.For example, Radon's original work (the traditional form) assumes that the integral is taken over a family of straight lines, and these lines are described using the normal form of line equation.This form found application in the field of computed tomography [2][3][4].In another context, the τ -p transform or slant stack, maps seismic data set from the time-offset distance (t, x) plane to the transform domain defined by the intercept time-slowness plane (τ, p), for instance [6,9,10,11].Another important form the Radon transform in seismic data processing is the parabolic Radon transform, where the integrals are taken over parabolas.[9][10][11].
There are several reconstruction algorithms for inverting Radon transform, [17][18][19][20][21][22][23][24][25][26][27][28].As said, the goal of this work is to introduce an improved algebraic reconstruction for the τ -p transform.In the literature there are some recent works that relate different types of Radon transformation to the traditional one.Our work is motivated by [9][10] but is different in several ways: integrating over lines only and a linear system with a precise coefficient matrix.In fact, we calculate the τ -p transform of the unit square function via a known formula for the traditional transform of the square function.This allows us to construct a linear system of the unknown (pixels).The 'interior-point' algorithm approach with some MATLAB optimization functions is used for solving the system and recovering the image.
The work is presented in the following order: in section 2, background material is reviewed.In section 3, our theory is presented, this includes the analytical formula for the traditional Radon transform and the τ -p of the square function; representation of an image using the square functions as well as an algebraic procedure for the reconstruction of an image from the τ-p transform.Experiments and demonstrations are presented in section 4, the conclusion is given in section 5.
In a different context, the slant stacking or the  −  transform of the function  is known as (3) This also can be written as,  ~(, ) = ∫ (, where  is the slope of the line  ,  is its intercept, shown in figure 1, and  is the delta function.Several fundamental properties of the transform can be derived from (1-4).In particular, the shifting rules are relevant to our work.Let  and  are related through the parameters x ' ,  ' : (5) and,

The Square Function
It is possible to derive the Radon transform analytically of some simple mathematical functions.In particular, the square function is relevant to this discussion.Define The image representation of this function is shown in figure 2a.The Radon transform of (, ) is calculated in [2], We can derive the  −  transform of (, ) by linking it to (8).It is possible to show that,

Representation of an Image Via the Square Function
The following image representation can be useful in two major computations: the first is the use of the analytical tools for direct manipulation of certain aspects of an image such as the Radon transform, that can detect any linear transformation applied on the image, the computation of image moments, projection moments, and many other important measures and features of a given image.The second, is the use of this representation in some tomography problems such as algebraic and iterative reconstruction algorithms.We now describe this method.Let  *4 (, ) = [ − ( − 1),  − ( − )], where  = 1, … , , and  = 1, … ,

A Reconstruction Method
Given  projections  ~(, ) of an image at the slopes  0,  !, … ,  9 .We build the system A =

Using Constrained Optimization
One way to handle the large system  =  from ( 16) is to minimize the function  16) are constructed.We may use more projections so that matrix  is not a square matrix where the system would be overdetermined.Let  * be the actual image vector;  * = [ 00 ,  0! , …  07 ,  !0 ,  !! … ,  !7 ,  70 ,  7! … ,  77 ] : as validating data, we find that  = ‖ * − ‖ !!≈ 10 $<= .This suggests that the matrix  is almost exact.In this way, all the lines intersect the image that is seated in the first quadrant.Using (9) that relates ,  to ,  we define all values of  and .Therefore, the number of intercepts () per projection = 2 = 512.In this way, the experiment involves  !equations against  !unknowns.In addition to the visual comparisons; we use the quantitative measure: Figure 4a shows the original image and figure 4b shows the reconstructed image.In this experiment,  = 3.14 × 10 $0! .

Example3
Figure 4a, shows a synthetic image that is generated randomly with  = 64, a total of 4096 unknowns.System ( 16) is constructed three different times; using different numbers of projections  ~(, ) with 2 lines ( values) per slope (rounded).Each experiment involves a number of equations as shown in table 1. Observe that, unlike the filter-back projection method, the algebraic approach does not require a full range of projections and works with some limitations of the number of available projections.
Figure 5a shows the original image, (b-d) are the output images using different numbers of equations, error D is shown in table 1.From these experiments we observe that, with the absence of noise in  ~(, ), the reconstruction is exact as in figures 4b,5d.

Conclusion
In this work we have seen an improved algebraic reconstruction for the τ -p transform.The τ -p transform of the unit square function is calculated via the known formula for the traditional transform of the square function.This allowed us to construct the linear system ( 14) of the unknowns (pixels) with an exact matrix of coefficients.The interior-point algorithm approach with some MATLAB optimization functions is used for solving the system and recovering the image.The approach is simple, practical, accurate, and is suitable for limited data tomography.

Figure 2 .
Figure 2. (a-b) Square and translated Square, (c) spatial coordinates and the pixel indexing.

Figure 5 .
Figure 5. (a) original image, (b-d) reconstructed images using different numbers of equations of system (16)