November 27, 2006

1991 Mathematics Subject Classification. 42A38, 42B08, 42B15.The first author was supported by the Bulgarian Ministry of Science and Education under Project MM–1402/2004. The second author was supported in part by the National Science Foundation under Grant DMS-0201669 .
<ph f="cmbx">Reconstruction of a polynomial from its Radon projections</ph>

Borislav Bojanov and Yuan Xu

Department of Mathematics, University of Sofia, Blvd. James Boucher 5, 1164 Sofia, Bulgaria E-mail address : boris@fmi.uni-sofia.bg Department of Mathematics, University of Oregon, Eugene, Oregon 97403-1222.
E-mail address : yuan@math.uoregon.edu

1 Introduction

Let f   be a function defined on the unit disk B 2   on the plane. A Radon projection of f   is the integral of f   over a line segment inside B 2   . More precisely, for any given pair ( θ , t )   of a real number t ∈ [ − 1 , 1 ]   and any angle θ   , let I ( θ , t )   denote the line segment inside the unit disk B 2   , where the line passes through the point ( t cos θ , t sin θ )   and is perpendicular to the vector ( cos θ , sin θ )   . Then
ℛ s θ ( f ; t ) : = ∫ s I ( θ , t ) f ( x , y ) d x d y (1.1)
defines the Radon projection of f   on the line segment I ( θ , t )   .
The Radon transform f ↦ { ℛ s θ ( f ; t ) }   associates to f   a family of univariate functions of t   , parameterized by θ   . The problem of reconstructing f   from full or partial knowledge of { ℛ s θ ( f ; t ) }   has been studied by many authors. It plays an essential role in the computer tomography. It is well-known that the set of Radon projections { ℛ s θ ( f ; t ) : 0 ≤ t ≤ 1 , 0 ≤ θ ≤ 2 π }   determines f   completely. Furthermore, it is known that if f   has compact support in B 2   , then f   is uniquely determined by any infinite set of Radon projections [15] . In practice, however, the data set is usually finite. Thus, the main problem is to obtain a good approximation to the function from a large collection of its Radon projections. See [12] for the background and detail discussions. Using a finite data set, one may determine a polynomial whose Radon projections agree with the given data. Such an approach has been considered in [9, 10, 11] in connection with computer tomography.
The present paper concerns with the problem of whether a polynomial of degree n   in two variables can be uniquely determined from a set of ( n + 1 ) ( n + 2 ) / 2   distinct Radon projections. The solution depends on the arrangement of the lines on which the projections are taken. One solution was given early in [11] , in which the Radon projections are taken over all possible line segments [ x s i , x s j ]   joining x s i   and x s j   , where x s 0 , x s 1 , … , x s n + 1   are n + 2   equally spaced points on the boundary of the unit disk B 2   . In [6] (see also [4] ), a more general result was proved, in which x s 0 , x s 1 , … , x s n + 1   are n + 2   distinct points on the boundary of a convex set and the result was extended in [7] to higher dimensions. Recently this question was considered in [1] , in which the lines are taken in ( n + 1 )   directions, with one line in the first direction, two lines in the second direction, etc. This set has the drawback of lacking symmetry.
The main result of the present paper shows that a set of Radon projections taken over 2 ⌊ n / 2 ⌋ + 1   parallel lines on each of the 2 ⌊ ( n + 1 ) / 2 ⌋ + 1   equidistant directions can determine a polynomial of degree n   uniquely. The set of the corresponding line segments possesses a rotation symmetry. If n = 2 m   is a fixed positive integer, then our method requires the Radon data in 2 m + 1   directions, which are taken to be in equally spaced angles along the unit circle, and in each direction we take m + 1   parallel lines, associated to t s 0 , t s 1 , … , t s m   . We prove that for almost all choices of { t s k }   a polynomial P s n   of degree n = 2 m   in two variables is uniquely determined by the set of its Radon projections over these line segments; more precisely, for any given f   , P s n   is uniquely determined such that
ℛ s φ s j ( P s n ; t s k ) = ℛ s φ s j ( f ; t s k ) , 0 ≤ j ≤ 2 m , 0 ≤ k ≤ m . (1.2)
The similar construction works for n = 2 m − 1   . Several examples of the sets of points { t s k }   that define regular interpolation are given. Furthermore, the polynomial P s n   can be easily computed. Thus, the result appears to offer a simple way to recover a polynomial from its Radon data.
The number of conditions  1.2 is equal to ( n + 1 ) ( n + 2 ) / 2   , which is the dimension of the space of bivariate polynomials of total degree at most n   . Hence,  1.2 can be interpreted as an interpolation problem by polynomials based on line integrals. As in the case of interpolation based on points, this bivariate problem is not always regular. The only general configurations in the literature that hold for every integer n   are that of Hakopian [7] mentioned above and the non-symmetric configuration given recently in [1] . The result of the present paper provides another family of examples. The proof is based on an observation that recovering polynomials from their Radon projections can be reduced, using a formula in [11] , to a family of univariate interpolation problems that uses certain special classes of algebraic polynomials. In the case of [1] , the strategy of solving the interpolation problem resembles the approach that uses the Bezout theorem for pointwise interpolation.
The approach in the present paper follows the strategy of [2] , which uses equally spaced points on circles for pointwise interpolation, that allows one to step beyond the limitation of the Bezout Theorem. The paper is organized as follows. In the following section we prove the existence and the uniqueness of the reconstructing problem. In Section 3 we show how to construct the polynomial and give an outline of the algorithm.

2 The existence and the uniqueness of the solution

2.1 Preliminary

Let B 2 = { ( x , y ) : x 2 + y 2 ≤ 1 }   denote the unit disk on the plane. Let θ   be the angle in the polar coordinates x = r cos θ , y = r sin θ , r ≥ 0 , 0 ≤ θ ≤ 2 π .   A line ℓ   whose slope is − cot θ   is defined by the equation ℓ ( x , y ) : = x cos θ + y sin θ − t = 0 ,   where t   is a real number. Clearly the line ℓ   passes through the point ( t cos θ , t sin θ )   and is perpendicular to the vector ξ = ( cos θ , sin θ )   . Since, for a fixed t   , the line corresponding to ( θ , t )   coincides with the line corresponding to ( π + θ , − t )   , we could assume θ ∈ [ 0 , π )   . Alternatively, we could assume that θ ∈ [ 0 , 2 π )   and t ≥ 0   .
We will also use ℓ   to denote the set of points on the line ℓ   and introduce the notation I ( θ , t ) = ℓ ∩ B 2 , 0 ≤ θ < π , − 1 ≤ t ≤ 1 ,   to denote the line segment of ℓ   inside B 2   . The points on I ( θ , t )   can be represented as follows:
x = t cos θ − s sin θ , y = t sin θ + s cos θ ,   for s ∈ [ − 1 − t 2 , 1 − t 2 ]   .
The Radon projection of a function f   in the direction θ   with a parameter t ∈ [ − 1 , 1 ]   is denoted by ℛ s θ ( f ; t )   ,
ℛ s θ ( f ; t ) : = ∫ s I ( θ , t ) f ( x , y ) d x d y
= ∫ s − 1 − t 2 1 − t 2 f ( t cos θ − s sin θ , t sin θ + s cos θ ) d s .
In the literature it is also called an X   -ray. It is clear that ℛ s θ ( f ; t ) = ℛ s π + θ ( f ; − t )   so that we can assume, for the sake of definiteness, either 0 ≤ θ < π   or 0 ≤ t ≤ 1   .

2.2 Polynomial bases

Let Π s n 2   denote the space of polynomials of total degree n   in two variables, which has dimension dim Π s n 2 = ( n + 1 ) ( n + 2 ) / 2   . If P ∈ Π s n 2   then P ( x ) = ∑ s k = 0 n ∑ s j = 0 k c s k , j x j y k − j , x = ( x , y ) .   Let U s k   denote the Chebyshev polynomial of the second kind, U s k ( x ) = sin ( k + 1 ) θ sin θ , x = cos θ .   For ξ = ( cos θ , sin θ )   and x = ( x , y )   , the ridge polynomial U s k ( θ ; ⋅ )   is defined by U s k ( θ ; x ) : = U s k ( 〈 x , ξ 〉 ) = U s k ( x cos θ + y sin θ ) .   Clearly U s k   is an element of Π s k 2   and it is constant on every line that is perpendicular to ξ   . The Radon projection of this function in any direction can be easily computed.
This is a result due to Marr [11] and it plays a central role in our discussion below.
Lemma 2.1. For each t ∈ ( − 1 , 1 )   , 0 ≤ θ , φ ≤ 2 π   , ℛ s φ ( U s k ( θ ; ⋅ ) ; t ) = 2 k + 1 1 − t 2 U s k ( t ) U s k ( cos ( φ − θ ) ) .  
The following useful relation follows easily from Marr's formula ([10] ),
1 π ∫ s B 2 U s k ( θ ; x ) U s k ( φ ; x ) d x = 1 k + 1 U s k ( cos ( φ − θ ) ) . (2.1)
Let V s n   denote the space of orthogonal polynomials of degree n   on B 2   with respect to the unit weight function; that is, P ∈ V s n   if P   is of degree n   and ∫ s B 2 P ( x ) Q ( x ) d x = 0 , f o r a l l Q ∈ Π s n − 1 2 .   For special choices of θ   and φ   , the equation  2.1 becomes an orthogonal relation.
Indeed, for k ∈ N   , let ξ s j , k = ( cos θ s j , k , sin θ s j , k ) , θ s j , k : = j π k + 1 , 0 ≤ j ≤ k .   For fixed k   , the points cos θ s j , k   , 1 ≤ j ≤ k   , are zeros of U s k   . The ridge polynomials U s k ( θ ; x )   have the remarkable orthogonal property that U s k ( θ ; ⋅ ) ∈ V s k   . Since the dimension of V s k   is k + 1   , it follows from  2.1 that the set P s k : = { U s k ( θ s j , k ; x ) : 0 ≤ j ≤ k }   is a basis for V s k   . In particular, this shows that the set { P s k : 0 ≤ k ≤ n }   is a basis for Π s n 2   . Together with Lemma  2.1 , this proves the following result:
Lemma 2.2. Every polynomial P s n ∈ Π s n 2   can be written uniquely as
P s n ( x ) = ∑ s k = 0 n ∑ s j = 0 k c s j , k U s k ( θ s j , k ; x ) . (2.2)
Furthermore, for each φ   and t   ,
ℛ s φ ( P s n ; t ) = 1 − t 2 ∑ s k = 0 n 2 k + 1 U s k ( t ) ∑ s j = 0 k c s j , k U s k ( cos ( φ − θ s j , k ) ) . (2.3)
There are several other explicit orthogonal bases for V s n   ; see, for example, [5] . If { Q s j k : 0 ≤ k ≤ n }   is an orthogonal basis of V s n   , then it was shown in [16] that
U s k ( θ s j , k ; x ) = 1 k + 1 ∑ s l = 0 k Q s l k ( x , y ) Q s l k ( cos θ s j , k , sin θ s j , k ) . (2.4)
One explicit basis of V s n   , denoted by Q s j , i ( x )   , is given in terms of polar coordinates as follows (cf. [5] ),
Q s j , 1 ( x , y ) = h s n − 2 j , 1 P s j ( 0 , n − 2 j ) ( 2 r 2 − 1 ) r n − 2 j cos ( n − 2 j ) θ , 0 ≤ 2 j ≤ n ,
Q s j , 2 ( x , y ) = h s n − 2 j , 2 P s j ( 0 , n − 2 j ) ( 2 r 2 − 1 ) r n − 2 j sin ( n − 2 j ) θ , 0 ≤ 2 j ≤ n − 1 ,
where P s j ( α , β )   denotes the usual Jacobi polynomial, h s 0 , 1 = 1 / ( n + 1 ) )   and h s j , 1 = h s j , 2 = 1 / ( 2 n + 2 )   for 1 ≤ 2 j ≤ n − 1   . Using this basis and restricting x   to the boundary of B 2   , the equation  2.4 becomes
U s 2 k ( cos ( φ − θ s j , 2 k ) ) = 1 + 2 ∑ s l = 1 k ( cos 2 l θ s j , 2 k cos 2 l φ + sin 2 l θ s j , 2 k sin 2 l φ ) , (2.5)
U s 2 k − 1 ( cos ( φ − θ s j , 2 k − 1 ) ) = 2 ∑ s l = 1 k ( cos ( 2 l − 1 ) θ s j , 2 k − 1 cos ( 2 l − 1 ) φ (2.6)
+ sin ( 2 l − 1 ) θ s j , 2 k − 1 sin ( 2 l − 1 ) φ ) . (2.7)
The above equations also follow from the elementary trigonometric identities 1 + 2 ∑ s j = 1 k cos 2 j θ = sin ( 2 k + 1 ) θ sin θ a n d 2 ∑ s j = 1 k cos ( 2 j − 1 ) θ = sin 2 k θ sin θ   with θ   replaced by φ − θ s j , 2 k   and φ − θ s j , 2 k − 1   , respectively. These relations allow us to rewrite  2.3 in a form that is easier to work with. In the following we use ⌊ x ⌋   to denote the integer part of x   .
Lemma 2.3. Let P s n   be given as in  2.2 . Then
ℛ s φ ( P s n ; t ) 1 − t 2 = ∑ s l = 0 ⌊ n / 2 ⌋ U s 2 l ( t ) [ a s 0 , 2 l + 2 ∑ s j = 1 l ( a s j , 2 l cos 2 j φ + b s j , 2 l sin 2 j φ ) ] (2.8)
+ ∑ s l = 1 ⌊ ( n + 1 ) / 2 ⌋ U s 2 l − 1 ( t ) [ 2 ∑ s j = 1 l ( a s j , 2 l − 1 cos ( 2 j − 1 ) φ + b s j , 2 l − 1 sin ( 2 j − 1 ) φ ) ] ,
in which the coefficients c s j , k   in  2.2 and a s j , k   , b s j , k   are related by
c s j , 2 l = 1 2 a s 0 , 2 l + ∑ s k = 1 l ( a s k , 2 l cos 2 k θ s j , 2 l + b s k , 2 l sin 2 k θ s j , 2 l ) (2.9)
for 0 ≤ j ≤ 2 l   , 0 ≤ l ≤ ⌊ n / 2 ⌋   , and
c s j , 2 l − 1 = ∑ s k = 1 l ( a s k , 2 l − 1 cos ( 2 k − 1 ) θ s j , 2 l − 1 + b s k , 2 l − 1 sin ( 2 k − 1 ) θ s j , 2 l − 1 ) (2.10)
for 0 ≤ j ≤ 2 l − 1   , 1 ≤ l ≤ ⌊ ( n + 1 ) / 2 ⌋   .
  • Proof. We only prove the case of n = 2 m   . Comparing  2.3 and  2.8 shows that
    1 2 l + 1 ∑ s j = 0 2 l c s j , 2 l U s 2 l ( cos ( φ − θ s j , 2 l ) ) = 1 2 a s 0 , 2 l + ∑ s j = 1 l ( a s j , 2 l cos 2 j φ + b s j , 2 l sin 2 j φ ) ,
    1 2 l ∑ s j = 0 2 l − 1 c s j , 2 l − 1 U s 2 l − 1 ( cos ( φ − θ s j , 2 l − 1 ) )
    = ∑ s j = 1 l ( a s j , 2 l − 1 cos ( 2 j − 1 ) φ + b s j , 2 l − 1 sin ( 2 j − 1 ) φ ) .
    Using the equations in  2.5 we can rewrite the left hand side in the form of the right hand side, which leads to the equations:
    a s j , 2 l = 2 2 l + 1 ∑ s k = 0 2 l c s k , 2 l cos 2 j θ s k , 2 l , b s j , 2 l = 2 2 l + 1 ∑ s k = 0 2 l c s k , 2 l sin 2 j θ s k , 2 l   for 1 ≤ j ≤ l   and 2 a s 0 , 2 l   satisfies the above formula with j = 0   , and
    a s j , 2 l − 1 = 1 l ∑ s k = 0 2 l − 1 c s k , 2 l − 1 cos ( 2 j − 1 ) θ s k , 2 l − 1 ,
    b s j , 2 l − 1 = 1 l ∑ s k = 0 2 l − 1 c s k , 2 l − 1 sin ( 2 j − 1 ) θ s k , 2 l − 1
    for 1 ≤ j ≤ l   . For fixed l   , the above linear relations can be written in the matrix form. After a proper normalization, the coefficient matrix turns out to be an orthogonal matrix. For example, for the even indices, we have 1 2 l + 1 [ 1 1 ⋅ ⋅ ⋅ 1 2 2 cos 2 π 2 l + 1 ⋅ ⋅ ⋅ 2 cos 2 l π 2 l + 1 0 2 sin 2 π 2 l + 1 ⋅ ⋅ ⋅ 2 sin 2 l π 2 l + 1 . . . . . . ⋅ ⋅ ⋅ . . . 2 2 cos 2 l 2 π 2 l + 1 ⋅ ⋅ ⋅ 2 cos 2 l 2 l π 2 l + 1 0 2 sin 2 l 2 π 2 l + 1 ⋅ ⋅ ⋅ 2 sin 2 l 2 l π 2 l + 1 ] ⋅ [ c s 0 , 2 l c s 1 , 2 l . . . c s 2 l , 2 l ] = 1 2 l + 1 [ a s 0 , 2 l a s 1 , 2 l 2 b s 1 , 2 l 2 . . . a s l , 2 l 2 b s l , 2 l 2 ] .   The ( l + 1 ) × ( l + 1 )   matrix in the left hand side, including the factor 1 / 2 l + 1   , is orthogonal. Hence, the linear system of equations can be easily reversed, which leads to the stated relations  2.9 and  2.10 . The formulas can also be verified directly by inserting  2.9 and  2.10 into the equations of a s j , k   and b s j , k   and using the well-known trigonometric identities such as 1 m + 1 ∑ s k = 0 m cos 2 j θ s k , m + i 1 m + 1 ∑ s k = 0 m sin 2 j θ s k , m = 1 m + 1 ∑ s k = 0 m e 2 i j θ s k , m = δ s 0 , j   for 0 ≤ j ≤ m   , where δ s 0 , j = 1   if j ≠ 0   and δ s 0 , j = 0   otherwise, and the elementary trigonometric identities such as 2 cos θ cos φ = cos ( θ + φ ) + cos ( θ − φ )   . □
To reconstruct the polynomial, we will determine the coefficients a s j , k   and b s j , k   and use them in  2.9 and  2.10 to determine c s j , k   in  2.2 . The following representation is useful.
Lemma 2.4. Let n = 2 m   or n = 2 m − 1   . Let P s n   be given as in  2.2 .
Then
ℛ s φ ( P s n ; t ) 1 − t 2 = A s 0 ( t ) + ∑ s j = 1 n [ A s j ( t ) cos j φ + B s j ( t ) sin j φ ] , (2.11)
where A s 0 ( t ) = ∑ s l = 0 ⌊ n / 2 ⌋ a s 0 , 2 l U s 2 l ( t ) ,   and for 1 ≤ j ≤ m   ,
A s 2 j ( t ) = 2 ∑ s l = j ⌊ n / 2 ⌋ a s j , 2 l U s 2 l ( t ) , B s 2 j ( t ) = 2 ∑ s l = j ⌊ n / 2 ⌋ b s j , 2 l U s 2 l ( t ) , (2.12)
A s 2 j − 1 ( t ) = 2 ∑ s l = j ⌊ n / 2 ⌋ a s j , 2 l − 1 U s 2 l − 1 ( t ) , B s 2 j − 1 ( t ) = 2 ∑ s l = j ⌊ n / 2 ⌋ b s j , 2 l − 1 U s 2 l − 1 ( t ) . (2.13)
  • Proof. This comes from changing the order of summations in  2.8 . □

2.3 Existence and uniqueness of the solution

We now state the Radon projections from which the polynomial P s n   can be uniquely recovered. These are given by
ℛ s φ s j , m ( P s n ; t s k ) , 0 ≤ k ≤ ⌊ n / 2 ⌋ , 0 ≤ j ≤ 2 m , m = ⌊ ( n + 1 ) / 2 ⌋ , (2.14)
which consists of total ( ⌊ n / 2 ⌋ + 1 ) ( ( 2 ⌊ ( n + 1 ) / 2 ⌋ + 1 ) = ( n + 1 ) ( n + 2 ) / 2   many projections, the same as the dimension of Π s n 2   . The angles are chosen to be equidistant,
Θ s m : = { φ s j , m = 2 j π 2 m + 1 : 0 ≤ j ≤ 2 m } . (2.15)
The reason that we choose the equidistant angles lies in the lemma below, which plays a key role in our study. Such a lemma appears first in [2] and has been used in [3] and [17] .
Lemma 2.5. For φ ∈ Θ s m   and P s n ∈ Π s n 2   with n = 2 m   or n = 2 m − 1   ,
ℛ s φ ( P s n ; t ) 1 − t 2 = A s 0 ( t )
+ ∑ s j = 1 m [ ( A s j ( t ) + A s 2 m − j + 1 ( t ) ) cos j φ + ( B s j ( t ) − B s 2 m − j + 1 ( t ) ) sin j φ ] ,
where we assume that A s 2 m = B s 2 m = 0   if n = 2 m − 1   .
  • Proof. Using the expression in the previous lemma, the proof follows from the fact that cos ( 2 m − j + 1 ) φ = cos j φ a n d sin ( 2 m − j + 1 ) φ = − sin j φ   for φ ∈ Θ s m   . □
For the expression in this lemma, the variable t   is in ( − 1 , 1 )   . We can also state the result for φ ∈ Θ ~ s m = { ( 2 j + 1 ) π / ( 2 m + 1 ) : 0 ≤ j ≤ 2 m }   , for which the sign before A s 2 m − j + 1   and B s 2 m − j + 1   will be reversed. It is easy to see that Θ ~ s m = Θ s m + π   modulus 2 π   . The two expressions are consistent with the fact that ℛ s θ ( f ; t ) = ℛ s π + θ ( f ; − t )   .
We are now ready to prove our main result, which shows that the polynomial P s n   can be uniquely determined by the data  2.14 . We need the following function classes: For n = 2 m   define X s j ( t ) = { U s 2 m ( t ) , U s 2 m − 2 ( t ) , … , U s 2 j ( t ) , U s 2 m − 1 ( t ) , U s 2 m − 3 ( t ) , … , U s 2 m − 2 j + 1 ( t ) }   for 1 ≤ j ≤ m − 1   . We also consider X s j   as a column vector in R m + 1   and regard Ξ s j ( t ) = [ X s j ( t s 0 ) , X s j ( t s 1 ) , … , X s j ( t s m ) ]   as an ( m + 1 ) × ( m + 1 )   matrix. A set of points t = { t s 0 , t s 1 , … , t s m }   in ( − 1 , 1 )   is said to be asymmetric if t   and − t   do not both belong to t   .
Theorem 2.6. Let t : = { t s 0 , t s 1 , … , t s m }   be a set of asymmetric distinct points in ( − 1 , 1 )   such that the matrices Ξ s 1 ( t ) , Ξ s 2 ( t ) , … , Ξ s m − 1 ( t )   are all nonsingular, then the polynomial P ∈ Π s 2 m 2   is uniquely determined by the set  2.14 of its Radon projections.
  • Proof. In order to prove that there is a unique solution, it is sufficient to show that if ℛ s φ s j , m ( P ; t s k ) = 0   for 0 ≤ j ≤ 2 m   and 0 ≤ k ≤ m   , then P ( x ) ≡ 0   .
    Setting its left hand side zero, the expression in Lemma  2.5 shows that A s 0 ( t s k ) + ∑ s j = 1 m [ ( A s j ( t s k ) + A s 2 m − j + 1 ( t s k ) ) cos j φ + ( B s j ( t s k ) − B s 2 m − j + 1 ( t s k ) ) sin j φ ] = 0   for 0 ≤ k ≤ m   and φ ∈ Θ s m   . The left hand side is a trigonometric polynomial of degree m   and it vanishes on 2 m + 1   points, hence, it follows from the uniqueness of the trigonometric interpolation that the coefficients have to be zero; that is,
    A s 0 ( t s k ) = 0 , A s j ( t s k ) + A s 2 m − j + 1 ( t s k ) = 0 , (2.16)
    B s j ( t s k ) − B s 2 m − j + 1 ( t s k ) = 0 , 1 ≤ j ≤ m (2.17)
    for 0 ≤ k ≤ m   . Since A s j   and A s 2 m − j + 1   have different parity, the coefficients in A s j ( t ) + A s 2 m − j + 1 ( t )   will not combine and there are exactly m + 1   terms in this polynomial, written as a linear combination of the Chebyshev polynomials { U s k }   . Let Y s 0 ( t ) = { U s 2 m ( t ) , U s 2 m − 2 ( t ) , … , U s 2 ( t ) , U s 0 ( t ) }   and let
    Y s 2 j ( t ) = { U s 2 m ( t ) , U s 2 m − 2 ( t ) , … , U s 2 j ( t ) , U s 2 m − 1 ( t ) , U s 2 m − 3 ( t ) , … , U s 2 m − 2 j + 1 ( t ) } ,
    Y s 2 j − 1 ( t ) = { U s 2 m − 1 ( t ) , U s 2 m − 3 ( t ) , … , U s 2 j + 1 ( t ) , U s 2 m ( t ) , U s 2 m − 2 ( t ) , … , U s 2 m − 2 j ( t ) } .
    Furthermore, let Ξ s Y s j   denote the matrix Ξ s Y s j = [ Y s j ( t s 0 ) , Y s j ( t s 1 ) , … , Y s j ( t s m ) ]   .
    The coefficients of the linear systems of equations in  2.16 are Ξ s Y s 0   , Ξ s Y s 2 j   for 1 ≤ j ≤ m / 2   , and Ξ s Y s 2 j − 1   for 1 ≤ j ≤ ( m + 1 ) / 2   . Hence, in order to prove that  2.16 implies A s 0 ( t ) ≡ 0   , A s j ( t ) = B s j ( t ) ≡ 0   for 1 ≤ j ≤ m   , we need to show that these matrices are all invertible. However, it is easy to see that we have the relation Ξ s Y s 2 ( m − j ) − 1 = Ξ s Y s 2 j   , which shows, in particular, that Ξ s Y s m = Ξ s Y s m − 1   . Thus, we can deal with Ξ s Y s 2 j   for 0 ≤ j ≤ m − 1   .
    Furthermore, Y s 0   contains only even polynomials U s 2 l   , 0 ≤ l ≤ m   . Using the notation s = t 2   , it follows that U s 2 l ( s )   is a polynomial of degree l   , so that the matrix Ξ s Y s 0   is always invertible as t   is a set of asymmetric distinct points. Thus, since Ξ s Y s 2 j = Ξ s j ( t )   , our assumption on t   implies that all matrices are invertible. □
A similar theorem holds in the case of n = 2 m − 1   . First we need to define for n = 2 m − 1   , X s j ( t ) = { U s 2 m − 2 ( t ) , U s 2 m − 4 ( t ) , … , U s 2 j ( t ) , U s 2 m − 1 ( t ) , U s 2 m − 3 ( t ) , … , U s 2 m − 2 j + 1 ( t ) }   for 1 ≤ j ≤ m − 1   and X s m ( t ) = { U s 2 m − 1 ( t ) , U s 2 m − 3 ( t ) , … , U s 3 ( t ) , U s 1 ( t ) } .   We will keep the notion Ξ s j ( t )   for the matrix built upon from the columns of X s j   , which is an m × m   matrix.
Theorem 2.7. Let t = { t s 0 , t s 1 , … , t s m − 1 }   be a set of asymmetric distinct points in ( − 1 , 0 ) ∪ ( 0 , 1 )   such that the matrices Ξ s 1 ( t ) , Ξ s 2 ( t ) , … , Ξ s m ( t )   are all nonsingular, then the polynomial P ∈ Π s 2 m − 1 2   is uniquely determined by the set  2.14 of its Radon projections.
That the points t s 0 , t s 1 , … , t s m   are all in ( − 1 , 0 ) ∪ ( 0 , 1 )   means that none of the points can be zero. In fact, since U s 2 k − 1   is an odd polynomial, Ξ s m ( t )   contains only odd polynomials, which vanish at the origin. Otherwise, the proof of this theorem is similar to that of the previous theorem. In the case of n = 2 m   , the set X s 0 = { U s 2 m , U s 2 m − 2 , … , U s 0 }   does not appear in the conditions of the theorem since the interpolation at distinct points in [ 0 , 1 )   by X s 0   is always regular.
Such a reduction of conditions does not appears to happen for n = 2 m − 1   .
For each fixed j   , the determinant of Ξ s j ( t )   can be consider as a polynomial function in t = ( t s 0 , t s 1 , … , t s m )   , so that det Ξ s j ( t ) = 0   defines a hypersurface of dimension m   . Hence, the determinant is not zero for almost all choices of t ∈ R m + 1   .
Evidently, the same holds true for m + 1   matrices. Hence, we have the following corollary:
Corollary 2.8. For n = 2 m   or 2 m − 1   and for almost all choices of distinct points t s 0 , t s 1 , … , t s m   in ( − 1 , 1 )   for n = 2 m   or in ( − 1 , 0 ) ∪ ( 0 , 1 )   for n = 2 m − 1   , the polynomial P s n ∈ Π s n 2   is uniquely determined by the set  2.14 of its Radon projections.

2.4 Choices of t s k  

We have shown that almost all choices of t = { t s 0 , t s 1 , … , t s m }   will lead to the unique solution of recovering the polynomial. One may ask the question if any choice of t ⊂ ( − 1 , 1 ) m + 1   will work. The answer, however, is negative as the following example shows.
Example: m = 2   . By Theorem  2.6 we only have to choose t s 0 , t s 1 , t s 2   in ( 0 , 1 )   such that the set X s 2 = { U s 2 , U s 3 , U s 4 }   has nonsingular determinant Ξ s 2 ( t )   .
Recall that U s 2 ( t ) = 4 t 2 − 1 , U s 3 ( t ) = 8 t 3 − 4 t , U s 4 ( t ) = 16 t 4 − 12 t 2 + 1 .   We can compute the determinant Ξ s 2 ( t )   explicitly and the result is
det Ξ s 2 ( t ) = 32 ∏ s 1 ≤ i < j ≤ 3 ( t s i − t s j )
× [ 8 t s 1 2 t s 2 2 t s 3 2 + 4 t s 1 t s 2 t s 3 ( t s 1 + t s 2 + t s 3 ) + ( 2 t s 1 2 − 1 ) ( 2 t s 2 2 − 1 ) ( 2 t s 3 2 − 1 ) ] .
Since t s i ∈ ( 0 , 1 )   , the first two terms in the square bracket are positive, while the third term could be negative. This shows, in particular, that det Ξ s 2 ( t ) ≠ 0   if one of the t s i   is 2 / 2   or if two of t s 1 , t s 2 , t s 3   are less than 2 / 2   and the other one is greater than 2 / 2   . However, det Ξ s 2 ( t )   can be zero for some choices of t s 0 , t s 1 , t s 2   .
On the positive side, we give two results that provides sets of points that will ensure the uniqueness of the reconstruction. The first result uses the following theorem due to Obrechkoff [13] .
Lemma 2.9. Let d μ   be a nonnegative weight function on an interval and let P s 0 , P s 1 , P s 2 , … ,   be the orthogonal polynomials with respect to d μ   . Denote by α s n   the largest zero of P s n   . Then the number of zeros of the polynomial a s 0 P s 0 + a s 1 P s 2 + … + a s n P s n   , where a s 0 , a s 1 , … , a s n   are real numbers, in the interval ( α s n , + ∞ )   is at most the number of sign changes in the sequence of the coefficients a s 0 , a s 1 , … , a s n   .
The Chebyshev polynomials U s 0 , U s 1 , U s 2 , …   are orthogonal with respect to the weight function 1 − x 2   on [ − 1 , 1 ]   , so that the above lemma implies the following:
Proposition 2.10. Let t = { t s 0 , t s 1 , … , t s m }   be a set of numbers that satisfies
cos π 2 m + 1 < t s 0 < t s 1 < … < t s m < 1 . (2.18)
Then the matrices Ξ s 1 ( t ) , … , Ξ s m − 1 ( t )   are all nonsingular.
  • Proof. If one of the matrices, say Ξ s j ( t )   , were singular, there would be a nonzero polynomial P ∈ span X s j   that vanishes on t   . This means that the number of zeros of P   in ( α s n , ∞ )   would be m + 1   . However, each set Ξ s j   has cardinality m + 1   so that the number of sign changes in the sequence of the coefficients of P   is at most m   . This is a contradiction to the conclusion of the lemma. □
The set of points given in this proposition ensures the uniqueness of determining a polynomial by the set of its Radon projections. The condition  2.18 implies that all points are clustered toward the one end of the interval [ − 1 , 1 ]   . This appears to be neither practical nor a good choice for computation. In fact, the numerical test indicates that some of the matrices Ξ s j   tend to have very large condition numbers.
Our second positive result is more interesting. Here the points { t s k }   are based on the zeros of the Chebyshev polynomials U s 2 m   . It is well-known that these zeros are given by η s j , 2 m : = cos j π 2 m + 1 , j = 1 , 2 , … , 2 m .   We state the result for n = 2 m   first.
Theorem 2.11. Let t s 0   be any point in ( − 1 , 1 )   such that U s 2 m ( t s 0 ) ≠ 0   .
Let t s j = η s 2 j , 2 m   for j = 1 , 2 , … , m   . Then the matrices Ξ s 1 ( t ) , Ξ s 2 ( t ) , … , Ξ s m − 1 ( t )   in Theorem 2.6 are all nonsingular. Consequently, the polynomial P ∈ Π s 2 m 2   is uniquely determined by the set  2.14 of its Radon projection.
  • Proof. It is easy to see that the set t   is asymmetric. Hence, according to Theorem 2.6, we only need to show that Ξ s j ( t )   is invertible for each j = 1 , 2 , … , m − 1   .
    Using the explicit expression of U s 2 m ( t )   , it is easy to see that
    U s 2 m − 2 j ( t s k ) = U s 2 m − 2 j ( η s 2 k , 2 m ) = − U s 2 j − 1 ( η s 2 k , 2 m ) = − U s 2 j − 1 ( t s k ) . (2.19)
    Indeed, since sin ( 2 m + 1 ) η s 2 k , 2 m = 0   and cos ( 2 m + 1 ) η s 2 k , 2 m = 1   , it follows from the addition formula that sin ( 2 m − 2 j + 1 ) η s 2 k , 2 m = − sin 2 j η s 2 k , 2 m   , from which the equation follows.
    Since U s 2 m ( t s 1 ) = … = U s 2 m ( t s m ) = 0   , we evidently have det Ξ s j ( t ) = U s 2 m ( t s 0 ) det Ξ ~ s j ( t ) ,   where Ξ ~ s j ( t )   is a submatrix of Ξ s j ( t )   with the column X s j ( t s 0 )   and the row containing U s 2 m ( t s j )   of Ξ s j ( t )   removed. Using the equation  2.19 , we can replace all rows of Ξ ~ s j ( t )   that contain even Chebyshev polynomials by rows that contain odd Chebyshev polynomials. More precisely, we replace the rows ( U s 2 i ( t s 1 ) , … , U s 2 i ( t s m ) )   for i = j , j + 1 , … , m − 1   by rows ( U s 2 m − 2 i − 1 ( t s 1 ) , … , U s 2 m − 2 i − 1 ( t s m ) )   , respectively. The new matrix has all elements given in terms of the Chebyshev polynomials of the odd degree; hence, det Ξ ^ s j ( t ) = ɛ det ( U s 2 i − 1 ( t s k ) ) s i = 1 , k = 1 m , m   where ɛ = ± 1   . Assuming now that this determinant is zero. Then there exists a non-zero polynomial, Q   , of the form Q ( x ) = ∑ s i = 1 m b s i U s 2 i − 1 ( x ) ,   which vanishes at t s 1 , … , t s m   . Since Q   is odd, it vanishes also at − t s 1 , … , − t s m   .
    Observe that the sets ( − t s 1 , … , − t s m )   and ( t s 1 , … , t s m )   do not overlap. We conclude that the polynomial Q   of degree 2 m − 1   vanishes at 2 m   distinct points and, consequently, it vanishes identically, a contradiction. □
A similar theorem holds for n = 2 m − 1   , which we state below.
Theorem 2.12. Let t s j = η s 2 j , 2 m   for j = 1 , 2 , … , m   . Then the matrices Ξ s 1 ( t ) , Ξ s 2 ( t ) , … , Ξ s m ( t )   in Theorem 2.7 are all nonsingular. Consequently, the polynomial P ∈ Π s 2 m − 1 2   is uniquely determined by the set  2.14 of its Radon projection.
The proof is similar to that of Theorem  2.11 and it uses the equation U s 2 m − 2 j ( η s 2 k − 1 , 2 m ) = U s 2 j − 1 ( η s 2 k − 1 , 2 m ) ,   which can be verified using elementary trigonometric identities.
We have conducted some numerical tests for identifying other sets of t s k   for which the matrices Ξ s 1 ( t ) , … , Ξ s m − 1 ( t )   are nonsingular, so that the reconstruction of a polynomial from the set  2.14 of its Radon projections is unique. For n ≤ 20   , it turns out that both the equidistant points in ( 0 , 1 )   and the Chebyshev points in ( 0 , 1 )   work out. See, however, the discussion at the end of the next section.

3 Reconstruction of polynomials

In this section we set up the algorithm that can be used to compute P s n   from the Radon projections. We only consider the case n = 2 m   . Let γ s j , k   be the data
γ s j , k = ℛ s φ s j ( f ; t s k ) / 1 − t s k 2 , 0 ≤ j ≤ 2 m , 0 ≤ k ≤ m , (3.1)
where φ s j , m   are given as in  2.15 and t s k   are distinct numbers in [ 0 , 1 )   , chosen in advance, such that the linear systems of equations  3.2 and  3.3 below have unique solutions for all j   .
Recall that the Lagrange interpolation by trigonometric polynomials is used in the proof of Theorem  2.6 . The Lagrange interpolation based on the points in  2.15 is given explicitly by (see, for example, [18] ) L s n f ( φ ) = ∑ s j = 0 2 m f ( φ s j , m ) ℓ s j ( φ ) , ℓ s j ( φ ) = sin ( m + 1 2 ) ( φ − φ s j , m ) ( 2 m + 1 ) sin 1 2 ( φ − φ s j , m ) ,   where we assume that f   is the function being interpolated; that is, L s n f ( φ s j , m ) = f ( φ s j , m )   for 0 ≤ j ≤ 2 m   . Using the well-known formula 1 + 2 cos φ + … + 2 cos m φ = sin ( m + 1 2 ) φ sin φ 2 ,   we can write L s n f   in the standard form of a trigonometric polynomial, L s n f ( φ ) = m s 0 ( f ) + 2 ∑ s j = 1 ( m s j C ( f ) cos j φ + m s j S ( f ) sin j φ )   where
m s j C ( f ) = 1 2 m + 1 ∑ s l = 0 2 m f ( φ s l , m ) cos j φ s l , m
m s j S ( f ) = 1 2 m + 1 ∑ s l = 0 2 m f ( φ s l , m ) sin j φ s l , m .
For each k   , 0 ≤ k ≤ m   , we will use the above formulas with f ( φ s j , m ) = γ s j , k   and we shall write m s j , k C = m s j C ( f )   and m s j , k S = m s j S ( f )   for this particular f   .
This will allow us to determine the values A s 0 ( t s k )   , A s j ( t s k ) + A s 2 m − j + 1 ( t s k )   and B s j ( t s k ) − B s 2 m − j + 1 ( t s k )   for 1 ≤ j ≤ m   .
The next step is to fix j   and use the values of A s 0 ( t s k )   , or A s j ( t s k ) + A s 2 m − j + 1 ( t s k )   , or B s j ( t s k ) − B s 2 m − j + 1 ( t s k )   for k = 0 , 1 , … , m   to determine the coefficients of A s j   and B s j   . For this purpose we consider the following linear systems of equations: For 1 ≤ j ≤ m / 2   ,
∑ s l = j m d s j , 2 l U s 2 l ( t s k ) + ∑ s l = m − j + 1 m d s m − j + 1 , 2 l U s 2 l − 1 ( t s k ) = m s 2 j , k , 0 ≤ k ≤ m , (3.2)
and for 1 ≤ j ≤ ( m + 1 ) / 2   ,
∑ s l = j m d s j , 2 l − 1 U s 2 l − 1 ( t s k ) + ∑ s l = m − j + 1 m d s m − j + 1 , 2 l U s 2 l ( t s k ) = m s 2 j − 1 , k , 0 ≤ k ≤ m . (3.3)
The quantities m s j , k   will be specified later. It is easy to see that the coefficient matrix of these systems of equations are Ξ s j   in Theorem  2.6 . There are total m   system of linear equations, each of size ( m + 1 ) × ( m + 1 )   . Solving these equations with proper m s j , k   will determine the coefficients of A s j   and B s j   .
The last step is to use  2.9 and  2.10 to get c s j , k   , which are the coefficients of P s n   in  2.2 .
Algorithm: Step 1. For 1 ≤ l ≤ m   and 0 ≤ k ≤ m   , compute
m s l , k C = 1 2 m + 1 ∑ s j = 0 2 m γ s j , k cos l θ s j a n d m s l , k S = 1 2 m + 1 ∑ s j = 0 2 m γ s j , k sin l θ s j . (3.4)
Step 2. Solve the systems of equations  3.2 and  3.3 to get the coefficients a s j , k   and b s k , j   in  2.8 :
Case 1. Solve the m   systems of  3.2 and  3.3 for m s j , k = m s j , k C   to get a s j , k : = d s j , k , 1 ≤ j ≤ m 0 ≤ k ≤ 2 m .   Case 2. Solve the m   systems of  3.2 and  3.3 for m s j , k = m s j , k S   to get
b s j , 2 l : = d s j , 2 l , 1 ≤ j ≤ m / 2 , a n d b s j , 2 l : = − d s j , 2 l , m + 2 2 ≤ j ≤ m ,
b s j , 2 l − 1 : = d s j , 2 l − 1 , 1 ≤ j ≤ m + 1 2 , a n d b s j , 2 l − 1 : = − d s j , 2 l − 1 , m + 1 2 ≤ j ≤ m ,
where 1 ≤ l ≤ m   .
Step 3. Substitute the output a s j , k   and b s j , k   in Step 2 into  2.9 and  2.10 to get c s j , k   . Then the polynomial P   is given by  2.2 .
The output of this algorithm is the polynomial P   that satisfies  1.2 . The main computation appears to be the Step 2. Initial numerical experiments indicate that the linear systems of equations  3.2 and  3.3 are ill-conditioned. However, among the set of points that we tested, the largest condition number of the matrices corresponding to the equidistant points t s k = ( k + 1 ) / ( m + 2 )   , 0 ≤ k ≤ m   , of ( 0 , 1 )   is smaller than that corresponding to the Chebyshev points t s k = cos ( k + 1 ) π / ( 2 m + 4 )   , 0 ≤ k ≤ m   or the points in Theorem  2.11 .
As pointed out by one referee, it is perhaps not surprising that the matrices of the linear systems  3.2 and  3.3 are ill-conditioned. Solving these linear systems means inverting the Radon transform, while it is known that the Radon reconstruction problem is ill posed. Thus, for m   large, the algorithm may not be useful for practical computation. On the other hand, our result shows that it is possible to determine a polynomial of degree n   uniquely from a set of ( n + 1 ) ( n + 2 ) / 2   Radon data that concurs with the parallel geometry. This appears to be of independent interest.
References

  1. B. Bojanov and I. K. Georgieva, Interpolation by bivariate polynomials based on Radon projections, Studia Math, 162 (2004), 141 160.
  2. B. Bojanov and Yuan Xu, On a Hermite interpolation by polynomials of two variables, SIAM J. Numer. Anal. 39 (2002), 1780-1793.
  3. B. Bojanov and Yuan Xu, On polynomial interpolation of two variables, J. Approx. Theory, 120 (2003), 267-282.
  4. A.S. Cavaretta, C.A. Micchelli and A. Sharma, Multivariate interpolation and the Radon transform, Part I, Math. Z. 174 (1980), 263–279; Part II, in Quantitive Approximation, DeVore and K. Schemer eds., Acad. Press, New York, 1980, 49–62.
  5. C. F. Dunkl and Yuan Xu, Orthogonal polynomials of several variables, Cambridge Univ. Press, 2001.
  6. H. Hakopian, Multivariate divided differences and multivariate interpolation of Lagrange and Hermite type, J. Approx. Theory, 34 (1982), 286-305.
  7. H. Hakopian, Multivariate spline functions, B-spline basis and polynomial interpolation II, Studia Math., 79 (1984), 91-102.
  8. H. Hakopian and S. Ismaeil, On a bivariate interpolation problem, J. Approx. Theory, 116 (2002), 76-99.
  9. C. Hamaker and D. C. Solmon, The angles between the null spaces of X   -rays, J. Math. Anal. Appl. 62 (1978), 1-23.
  10. B. Logan and L. Shepp, Optimal reconstruction of a function from its projections, Duke Math. J. 42 (1975), 645-659.
  11. R. Marr, On the reconstruction of a function on a circular domain from a sampling of its line integrals, J. Math. Anal. Appl., 45 (1974), 357-374.
  12. F. Natterer, The mathetmaics of computerized tomography, Reprint of the 1986 original. Classics in Applied Mathematics, 32. SIAM, Philadelphia, PA, 2001.
  13. N. Obrechkoff, Zeros of Polynomials, BAN, Sofia, 1963 (in Bulgarian); English transl.: Bulgarian Academic Monographs 7, Marin Drinov Academic Publishing House, Sofia, 2003.
  14. T. J. Rivlin, Chebyshev Polynomials: from Approximation Theory to Algebra and Number Theory, 2nd ed., Wiley, New York, 1990.
  15. D. C. Solmon, The X-ray transform, J. Math. Anal. Appl. 56 (1976), 61-83.
  16. Yuan Xu, Funk-Hecke formula for orthogonal polynomials on spheres and on balls, Bull. London Math. Soc. 32 (2000), 447-457.
  17. Yuan Xu, Polynomial interpolation on the unit sphere, SIAM J. Numer. Anal. 41 (2003), p. 751-766.
  18. A. Zygmund, Trigonometric series, Cambridge University Press, Cambridge, 1959.

Department of Mathematics, University of Sofia, Blvd. James Boucher 5, 1164 Sofia, Bulgaria E-mail address : boris@fmi.uni-sofia.bg Department of Mathematics, University of Oregon, Eugene, Oregon 97403-1222.
E-mail address : yuan@math.uoregon.edu