Next Article in Journal
Mean-Value-at-Risk Portfolio Optimization Based on Risk Tolerance Preferences and Asymmetric Volatility
Previous Article in Journal
Tomographic Reconstruction: General Approach to Fast Back-Projection Algorithms
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Chebyshev Interpolation Using Almost Equally Spaced Points and Applications in Emission Tomography

by
Vangelis Marinakis
1,
Athanassios S. Fokas
2,3,
George A. Kastis
3,4 and
Nicholas E. Protonotarios
3,*
1
Department of Civil Engineering, University of the Peloponnese, 26334 Patras, Greece
2
Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Cambridge CB3 0WA, UK
3
Mathematics Research Center, Academy of Athens, 11527 Athens, Greece
4
Institute of Nuclear and Radiological Science and Technology, Energy and Safety (INRASTES), National Center for Scientific Research “Demokritos”, 15310 Aghia Paraskevi, Greece
*
Author to whom correspondence should be addressed.
Mathematics 2023, 11(23), 4757; https://doi.org/10.3390/math11234757
Submission received: 27 October 2023 / Revised: 18 November 2023 / Accepted: 23 November 2023 / Published: 24 November 2023
(This article belongs to the Special Issue Advances in Inverse Problems and Imaging)

Abstract

:
Since their introduction, Chebyshev polynomials of the first kind have been extensively investigated, especially in the context of approximation and interpolation. Although standard interpolation methods usually employ equally spaced points, this is not the case in Chebyshev interpolation. Instead of equally spaced points along a line, Chebyshev interpolation involves the roots of Chebyshev polynomials, known as Chebyshev nodes, corresponding to equally spaced points along the unit semicircle. By reviewing prior research on the applications of Chebyshev interpolation, it becomes apparent that this interpolation is rather impractical for medical imaging. Especially in clinical positron emission tomography (PET) and in single-photon emission computerized tomography (SPECT), the so-called sinogram is always calculated at equally spaced points, since the detectors are almost always uniformly distributed. We have been able to overcome this difficulty as follows. Suppose that the function to be interpolated has compact support and is known at q equally spaced points in 1 , 1 . We extend the domain to a , a , a > 1 , and select a sufficiently large value of a, such that exactly q Chebyshev nodes are included in 1 , 1 , which are almost equally spaced. This construction provides a generalization of the concept of standard Chebyshev interpolation to almost equally spaced points. Our preliminary results indicate that our modification of the Chebyshev method provides comparable, or, in several cases including Runge’s phenomenon, superior interpolation over the standard Chebyshev interpolation. In terms of the L norm of the interpolation error, a decrease of up to 75% was observed. Furthermore, our approach opens the way for using Chebyshev polynomials in the solution of the inverse problems arising in PET and SPECT image reconstruction.

1. Introduction

Chebyshev polynomials of the first kind have been extensively employed in the context of polynomial interpolation. Since their introduction by the renowned Russian mathematician Pafnuty Chebyshev [1], these orthogonal polynomials have played a pivotal role in approximation theory [2], providing an invaluable tool in numerical analysis [3]. Furthermore, Chebyshev polynomials are widely used in several fields including scientific computing [4], matrix theory [5], integral transforms [6], and image processing [7], as well as engineering [8,9], machine learning [10], quantum computing [11], and medical imaging [12,13].
Chebyshev interpolation is a numerical method that involves the approximation of an arbitrary function via a polynomial that passes through a set of given data points [14]. In particular, it is known that Chebyshev polynomials of the first kind have proved particularly useful. In the interval [ 1 , 1 ] they are orthogonal, i.e., pairwise perpendicular, in the L 2 -sense. This specific property renders them particularly useful for the efficient calculation of the coefficients of the interpolating polynomial. In addition, Chebyshev polynomials of the first kind converge rapidly to the function being interpolated, especially near the endpoints of the interval; this fact is important when approximating functions with singularities [15]. Another benefit of Chebyshev interpolation is its ability to achieve high accuracy with relatively few data points, making it particularly useful in cases of limited data [16]. This is the result not only of the intrinsic properties of Chebyshev polynomials but also of the choice of interpolation points, referred to as nodes.
Although standard polynomial interpolation methods usually employ equally spaced points, Chebyshev interpolation, rather than equally spaced points along a line, instead uses the roots of Chebyshev polynomials, known as Chebyshev nodes, corresponding to equally spaced points along the unit semicircle (see Figure 1). It is well known (e.g., [16]) that Chebyshev nodes tend to minimize certain interpolation errors and are thus preferred over equidistant nodes, which may suffer from the so-called Runge phenomenon [17,18], manifested as large oscillations near the endpoints of the interpolation interval. Furthermore, Chebyshev polynomials are known to provide the best, in the L -norm sense, polynomial approximation to a function [19].
It should be noted that in medical imaging applications, common practice dictates that data should be stored at uniformly spaced intervals. Especially in clinical positron emission tomography (PET) and in single-photon emission computerized tomography (SPECT), the so-called sinogram data are always calculated at equally spaced points. Because of the many advantages of Chebyshev polynomials, despite this difficulty, several attempts have been made to incorporate Chebyshev interpolation in a tomographic setting employing Chebyshev nodes, both in PET [20,21,22,23] and in magnetic particle imaging (MPI) [24,25].
In this work, we attempt to overcome the difficulty imposed by the non-uniform nature of Chebyshev nodes as follows. Suppose that the function to be interpolated has compact support, as commonly attributed to emission tomography [26,27,28], and is known at q equally spaced points in 1 , 1 . We extend the domain to a , a , a > 1 , and select a sufficiently large value of a, such that exactly q Chebyshev nodes are included in 1 , 1 , which are almost equally spaced. In this way, we provide a generalization of the concept of standard Chebyshev interpolation to almost equally spaced points. Our preliminary results indicate that our modification of the Chebyshev interpolation provides comparable, or, in several cases including Runge’s phenomenon, superior interpolation over the standard Chebyshev and Lagrangian techniques. Furthermore, our approach opens the way for using Chebyshev polynomials in the solution of the inverse problems arising in PET and SPECT image reconstruction.

2. Chebyshev Polynomials of the First Kind

2.1. The Case of [−1, 1]

We consider Chebyshev polynomials of the first kind and of degree n N , denoted by T n ( x ) , namely
T n ( x ) = cos ( n arccos x ) , x [ 1 , 1 ] , n N .
It is straightforward to see that
T 0 ( x ) = 1 , a n d T 1 ( x ) = x .
Furthermore, it is well known [29] that the following recursive relation is valid for Chebyshev polynomials of the first kind:
T n + 1 ( x ) = 2 x T n ( x ) T n 1 ( x ) , n 1 .
It is worth noting that T n ( x ) is a polynomial of degree n,
T n ( x ) = i = 0 n a i x i , with a n = 2 n 1 ,
and that T ( x ) = 0 n form a basis for P n , i.e., the set of polynomials with degree up to n.
The zero set of the polynomials T n ( x ) is defined as follows:
x ker T n : = x 1 , 1 | T n ( x ) = 0 .
The definition of T n ( x ) given in Equation (1) implies that its roots, denoted by x k ( n ) , are located in
x k ( n ) = cos 2 k 1 n π 2 , k = 1 , 2 , , n .
Equation (5), along with the fact that Chebyshev polynomials of the first kind are polynomials of degree n, imply that T n ( x ) has n distinct roots, namely x k ( n ) k = 1 n , hence | ker T n | = n . Furthermore, Equation (3) implies that
T n ( x ) = 2 n 1 k = 1 n ( x x k ( n ) ) .
Let f ( x ) denote an arbitrary function on the interval [ 1 , 1 ] , and c j denote a set of n coefficients, for j = 0 , 1 , , n 1 , defined as follows:
c j = 2 n k = 1 n f ( x k ( n ) ) T j ( x k ( n ) ) = 2 n k = 1 n f ( x k ( n ) ) cos ( 2 k 1 ) j n π 2 ,
where x k ( n ) ker T n ( x ) , given in Equation (5). It follows that the approximation formula
f ( x ) 1 2 c 0 + k = 1 n 1 c k T k ( x )
is exact for all n roots of the Chebyshev polynomial T n ( x ) , see Theorem 6.7 of [14].

2.2. Extending the Domain: The Case of [−a, a]

In Equation (1), if we extend the domain [ 1 , 1 ] to [ c , d ] , with c < d , then the definition of a Chebyshev polynomial of the first kind of degree n is modified in the following manner:
T ^ n ( x ) = T n ( s ) , where s = 2 x ( c + d ) d c , c d ,
as in [14]. In the special case of [ a , a ] , with a > 1 ,
T ^ n ( x ) = T n x a .
Equation (10), via Equation (1), may be rewritten as:
T ^ n ( x ) = cos n arccos x a , x [ a , a ] , a > 1 , n N .
Equation (11) yields
T ^ 0 ( x ) = 1 , and T ^ 1 ( x ) = x a ,
whereas the recursive Relations (2b) imply
T ^ n + 1 ( x ) = 2 x a T ^ n ( x ) T ^ n 1 ( x ) , n 1 .
The roots of T ^ n ( x ) are located in
x ^ k ( n ) = a cos 2 k 1 n π 2 , k = 1 , 2 , , n .
Defining the n coefficients { c ^ j } j = 0 n 1 as
c ^ j = 2 n k = 1 n f ( x ^ k ( n ) ) T ^ j ( x ^ k ( n ) ) = 2 n k = 1 n f ( x ^ k ( n ) ) cos ( 2 k 1 ) j n π 2
for the domain [ a , a ] , Equation (8) becomes
f ( x ) 1 2 c ^ 0 + k = 1 n 1 c ^ k T ^ k ( x ) .

3. Chebyshev Interpolation with Almost Equally Spaced Points

We aim to approximate an unknown function f ( x ) , defined in [ 1 , 1 ] , via a polynomial of degree of at most q 1 , given that the values of f ( x ) are known at the following q equally space points:
x ˜ m = 1 2 ( m 1 ) q 1 = q + 1 2 m q 1 , m = 1 , 2 , , q , q > 1 .
Hence,
x ˜ 1 = 1 , x ˜ 2 = 1 2 q 1 , x ˜ 3 = 1 4 q 1 , , x ˜ q = 1 .
In Equation (16), the equally spaced points are in descending order, following the trend of the roots of the Chebyshev polynomial, as in Equations (5) and (13). Clearly, in this case, we cannot apply Chebyshev interpolation since the values of f ( x ) in the roots of T q ( x ) are unknown.
To overcome this difficulty, we extend the domain from [ 1 , 1 ] to [ a , a ] , a > 1 , as in Section 2.2. Let n denote the number of roots of the corresponding Chebyshev polynomial T ^ n ( x ) in [ a , a ] (see Equation (13)). We require that exactly q of the n Chebyshev roots in [ a , a ] are included in [ 1 , 1 ] . It is worth mentioning that the roots of T n ( x ) , defined in (5), and of T ^ n ( x ) , defined in Equation (13), are not equally spaced in the standard sense; however, they are equally spaced along the unit semicircle and a semicircle with radius a, respectively, (see Figure 1). In the following proposition, we provide a criterion for the choice of a, so that the q Chebyshev roots included in [ 1 , 1 ] are almost equally spaced.
Proposition 1.
Let a > 1 be such that exactly q > 1 of the n Chebyshev roots in [ a , a ] , defined in Equation (13), are included in [ 1 , 1 ] , and let ℓ be an odd integer, such that n = q . Then, for n sufficiently large, the equidistant points defined in Equation (16) are almost equal to the q Chebyshev nodes.
Proof. 
Since a is such that q of the n Chebyshev roots in [ a , a ] , defined in Equation (13), are included in [ 1 , 1 ] , it follows that there exists an integer λ > 1 , such that the number of Chebyshev roots in [ a , 1 ) , and similarly for ( 1 , a ] , is exactly λ (see Figure 2). Hence, λ is such that
2 λ + q = n λ , q > 1 .
Taking into account Equation (18), the q Chebyshev roots in [ 1 , 1 ] via Equation (13) may be rewritten as
x ^ m ( q ) = a cos 2 ( λ + m ) 1 n π 2 , m = 1 , 2 , , q , q > 1 .
Given that n = q , where is an odd integer, Equation (18) yields
λ = ( 1 ) 2 q , odd , > 1 ,
i.e., λ is a multiple of q. Inserting Equation (20) in Equation (19) implies
x ^ m ( q ) = a cos ( 1 ) q + 2 m 1 q π 2 = a sin q + 1 2 m q π 2 , m = 1 , 2 , , q .
Since q of the n roots x ^ k ( n ) are included in [ 1 , 1 ] , x ˜ 1 = x ˜ q = 1 and x ^ 1 ( q ) = x ^ q ( q ) , it is reasonable to require that x ^ 1 ( q ) x ˜ 1 = 1 . Thus, Equation (21) implies
a 1 sin q 1 q π 2 .
Furthermore, taking into account Equation (22), we may rewrite Equation (21) in the following form:
x ^ m ( q ) sin q + 1 2 m q π 2 sin q 1 q π 2 .
Taking the corresponding Taylor series yields
x ^ m ( q ) = sin q + 1 2 m q π 2 sin q 1 q π 2 = q + 1 2 m q 1 + O 1 2 .
Hence, for (and therefore n) sufficiently large, the Chebyshev nodes are almost equally spaced, since x ^ m ( q ) x ˜ m , for all m, as implied by Equations (16) and (24).    □
Figure 2. Chebyshev roots distribution in the extended domain [ a , a ] [ 1 , 1 ] , a > 1 .
Figure 2. Chebyshev roots distribution in the extended domain [ a , a ] [ 1 , 1 ] , a > 1 .
Mathematics 11 04757 g002
Since the approximation of the unknown function f ( x ) involves the interval [ 1 , 1 ] , and not [ a , a ] , we may assume that f ( x ) has compact support, namely f ( x ) = 0 , for | x | > 1 . Thus, we may employ Chebyshev interpolation in [ 1 , 1 ] , as in Equation (15), utilizing Equation (14), assuming f ( x ^ m ( q ) ) = f ( x ˜ m ) , for m = 1 , , q , and f ( x ) = 0 , for any other root of T ^ n ( x ) . A step-by-step explanation of our modified Chebyshev interpolation method, including the procedure for extending the domain, is presented in Algorithm 1.
Algorithm 1 Computational steps of the proposed Chebyshev-based interpolation.
Input: q, , f ( x ) given at x ˜ m with m = 1 , , q
Compute: a via Equation (22)
Assume: f ( x ^ m ( q ) ) = f ( x ˜ m ) , for m = 1 , , q , and f ( x ) = 0 , for any other root of T ^ n ( x )
Compute: c ^ j via Equation (14)
Output: approximation of f ( x ) via Equation (15)

4. Examples of Chebyshev Interpolation with Almost Equally Spaced Points

We implemented and evaluated the efficacy of our methodology by considering two interpolation examples, namely, an exponential function and a rational function. In particular, we considered the exponential function
f ( x ) = e 5 x 2 ,
and the rational function
f ( x ) = 1 1 + 16 x 2 ,
which is well known regarding the so-called Runge’s phenomenon.
In order to approximate f ( x ) defined on the interval [ 1 , 1 ] , we employed three different polynomial interpolation schemes, namely,
  • Lagrange interpolation with equally spaced points on [ 1 , 1 ] (e.g., [19]);
  • Chebyshev interpolation with Chebyshev nodes on [ 1 , 1 ] , according to Equations (7) and (8);
  • Chebyshev interpolation with almost equally spaced points, according to Proposition 1 and Equations (14) and (15).
For the first two interpolation methodologies, we considered q = 11 interpolation points on [ 1 , 1 ] . For the third, we further assumed an extended interval of the form [ a , a ] , a > 1 . This interval is such that exactly q = 11 almost equally spaced Chebyshev nodes are included in [ 1 , 1 ] . Furthermore, according to Proposition 1, we assumed n = q , odd. In the case of the exponential function, we chose = 15 , hence n = 165 Chebyshev nodes are lying on [ a , a ] , which, according to Equation (22), implies that a 10.52 . In the case of the rational function, choosing = 5 implies that n = 55 Chebyshev nodes are lying on [ a , a ] , which, in turn, implies that a 3.55 (see Equation (22)). We then performed Chebyshev interpolation on the interval [ a , a ] via Equations (14) and (15). It is worth mentioning that in Equation (14) we assumed f ( x ^ m ( 11 ) ) = f ( x ˜ m ) , m = 1 , , 11 , and f ( x ^ k ( n ) ) = 0 , for all other n q = 2 λ roots of T ^ n ( x ) on [ a , 1 ) ( 1 , a ] (see Equation (18) and Figure 2). All methods of interpolation investigated for both examples of f ( x ) are shown in Figure 3 and Figure 4, respectively, where blue represents the function itself, and red corresponds to its approximation.
We calculated the L norm of the error of the equally spaced points x ˜ , defined in Equation (16), and of the almost equally spaced Chebyshev points x ^ ( 11 ) , defined in (19), i.e.,
Exponential : x ˜ x ^ ( 11 ) = max | x ˜ 1 x ^ 1 ( 11 ) | , , | x ˜ 11 x ^ 11 ( 11 ) | 0.0006 , Rational : x ˜ x ^ ( 11 ) = max | x ˜ 1 x ^ 1 ( 11 ) | , , | x ˜ 11 x ^ 11 ( 11 ) | 0.0053 .
Upon designating the new Chebyshev nodes as “almost equally spaced”, it is evident, from the above L norms of their difference from uniformly spaced points, that their spacing is indeed nearly uniform. This is also apparent due to the fact that the L norm of the difference is 0.03% and 0.265% of the length of the integral [ 1 , 1 ] for the exponential and rational functions, respectively.
We calculated the L norm of the interpolation error, namely f p , where p denotes the interpolating polynomial, with p P q 1 for the Lagrange and standard Chebyshev cases, and p P n 1 for the almost equally spaced points Chebyshev case. This error norm was computed by taking 10,000 points along the interval [ 1 , 1 ] . The results of our investigation are presented in Table 1. For both examples investigated, the results demonstrate that the novel Chebyshev-based method consistently outperforms the traditional approach with respect to the L norm. In particular, for the case of the exponential function, the decrease in the L norm of the interpolation error from the traditional Chebyshev method to the modified Chebyshev approach was 50%. Similarly, for the case of the rational function, the corresponding decrease was 75%.
It must be emphasized that our analysis is, at this stage, preliminary. Furthermore, the limitations of global versus local or piecewise (splines) polynomial interpolation remain the same, since both the traditional and modified Chebyshev approaches are special cases of Lagrange interpolation. Specifically, the Chebyshev approach consists of a Lagrange polynomial interpolation employed at the Chebyshev nodes.

5. Application to Emission Tomography

It is well known [27] that the inversion of the Radon transform in two dimensions is given by
f ( x 1 , x 2 ) = 1 4 π 2 0 2 π H ( ρ , θ ) ρ d θ , ρ = x 2 cos θ x 1 sin θ ,
where H denotes the Hilbert transform of the sinogram, f ^ , i.e.,
H ( ρ , θ ) = f ^ ( ρ , θ ) ρ ρ d ρ = 1 1 f ^ ( ρ , θ ) ρ ρ d ρ ,
where, as commonly attributed to medical imaging (see for example [26]), we assume that f ^ has compact support on [ 1 , 1 ] , and that f ^ ( 1 , θ ) = f ^ ( 1 , θ ) = 0 [27,28].
In order to apply our method, we consider an extension of the usual sinogram ( f ^ ). Let f ˜ denote the extension of f ^ from [ 1 , 1 ] to [ a , a ] , a > 1 , defined as follows:
f ˜ ( ρ , θ ) = f ^ ( ρ , θ ) , ρ [ 1 , 1 ] , 0 , ρ [ a , 1 ) ( 1 , a ] , a > 1 .
Then, according to Equation (15), the Chebyshev-type expansion of the first-kind f ˜ is given by
f ˜ ( ρ , θ ) = 1 2 c ^ 0 ( θ ) + k = 1 n 1 c ^ k ( θ ) T ^ k ( ρ ) ,
where T ^ n ( ρ ) are defined in Equation (11). Hence, the Hilbert transform of f ^ via Equation (28) may be rewritten in terms of f ˜ as follows:
H ( ρ , θ ) = f ˜ ( ρ , θ ) ρ ρ d ρ = 1 1 f ˜ ( ρ , θ ) ρ ρ d ρ = 1 1 f ˜ ( ρ , θ ) f ˜ ( ρ , θ ) ρ ρ d ρ + 1 1 f ˜ ( ρ , θ ) ρ ρ d ρ = f ˜ ( ρ , θ ) ln 1 ρ 1 + ρ + k = 1 n 1 c ^ k ( θ ) 1 1 T ^ k ( ρ ) T ^ k ( ρ ) ρ ρ d ρ = 1 2 c ^ 0 ( θ ) + k = 1 n 1 c ^ k ( θ ) T ^ k ( ρ ) ln 1 ρ 1 + ρ + k = 1 n 1 c ^ k ( θ ) 1 1 T ^ k ( ρ ) T ^ k ( ρ ) ρ ρ d ρ .
We denote the integral of the last line of Equation (31) by I k ( ρ ) , namely,
I k ( ρ ) = 1 1 T ^ k ( ρ ) T ^ k ( ρ ) ρ ρ d ρ .
Taking into account Equation (12a), it is easy to establish that
I 0 ( ρ ) = 0 , and I 1 ( ρ ) = 2 a .
Inserting the recursive Relation (12b) in Equation (32) yields
I k + 1 ( ρ ) = 2 a J k ( ρ ) I k 1 ( ρ ) ,
where
J k ( ρ ) = 1 1 ρ T ^ k ( ρ ) ρ T ^ k ( ρ ) ρ ρ d ρ .
We modify J k ( ρ ) by adding and subtracting the term ρ T ^ k ( ρ ) on the numerator inside the integral, as follows:
J k ( ρ ) = A k + ρ I k ( ρ ) ,
where A k denotes the following definite integral
A k = 1 1 T ^ k ( ρ ) d ρ .
Inserting Equation (36) in Equation (34) yields the following recursive relation:
I k + 1 ( ρ ) = 2 ρ a I k ( ρ ) I k 1 ( ρ ) + 2 a A k .
We note that it can be shown that
A k = 2 1 k 2 cos k sec 1 a + k a 2 1 sin k sec 1 a , k even , 0 , k odd .
By differentiating Equation (38) with respect to ρ , we obtain
I k + 1 ( ρ ) = 2 a I k ( ρ ) + 2 ρ a I k ( ρ ) I k 1 ( ρ ) ,
where prime denotes differentiation. Equation (33) implies that I k ( ρ ) vanishes for k = 0 and k = 1 , i.e.,
I 0 ( ρ ) = I 1 ( ρ ) = 0 .
Taking into account all of the above, the derivative of the Hilbert transform of the sinogram may be written as
H ( ρ , θ ) ρ = c ^ 0 ( θ ) ρ 2 1 + k = 1 n 1 c ^ k ( θ ) ln 1 ρ 1 + ρ T ^ k ( ρ ) + 2 T ^ k ( ρ ) ρ 2 1 + I k ( ρ ) ,
where T ^ k ( ρ ) denotes the derivative of T ^ k ( ρ ) with respect to ρ , and I k ( ρ ) is defined in Equations (40) and (41).

6. Numerical Implementation for PET Image Reconstruction

For the numerical implementation of the new, Chebyshev-based image reconstruction technique, we considered the function
f ( x 1 , x 2 ) = 1 , x 1 2 + x 2 2 < 1 4 , 0 , otherwise ,
i.e., a uniform circular phantom. This phantom and its line profile for x 2 = 0 are shown on the right and left side of the upper row of Figure 5, respectively.
The corresponding sinogram of the uniform phantom f, defined in Equation (43), involves its Radon transform R , and is denoted by f ^ , namely,
f ^ ( ρ , θ ) = R f ( x 1 , x 2 ) = L f x 1 ( s , ρ , θ ) , x 2 ( s , ρ , θ ) d s ,
where ( x 1 , x 2 ) are the usual Cartesian coordinates, ρ represents the signed distance from the origin, i.e.,
ρ : = ρ ( x 1 , x 2 ; θ ) = x 2 cos θ x 1 sin θ ,
d s denotes an arc length differential along parallel lines L for all angles θ , and s represents a parameter along the line L. We note that the Radon transform of f can be analytically calculated, i.e.,
f ^ ( ρ , θ ) = 2 1 4 ρ 2 , | ρ | < 1 2 , 0 , otherwise .
By appropriately approximating the sinogram in the variable ρ for all θ , we applied two different reconstruction techniques, namely:
  • The Chebyshev-based reconstruction technique with Chebyshev nodes for ρ [ 1 , 1 ] , following the inversion formula described in Algorithm 1 of [22];
  • Our novel Chebyshev-based reconstruction technique with almost equally spaced points, according to Proposition 1 and the inversion approach described in Section 5.
For both reconstruction techniques, we evaluated the sinogram for 90 equally spaced angles θ in [ 0 , π ) . Furthermore, for the case of standard Chebyshev approximation of the sinogram where ρ [ 1 , 1 ] , we considered q=119 values of ρ , at the roots of the corresponding Chebyshev polynomials of degree q, T q ( ρ ) . For the case of the Chebyshev approximation with q almost equally spaced points on [ 1 , 1 ] , we assumed an interval of the form [ a , a ] for the sinogram, thus extending its domain of definition as in Equation (29). Then, based on the relation n = q and the fact that q = 119 points are required to lie on the interval [ 1 , 1 ] , we chose = 27 , thus n = 3213 , which, via Equation (22), implies a 17.34 . The reconstruction results for the standard Chebyshev approximation with Chebyshev nodes in [ 1 , 1 ] , and for the Chebyshev with almost equally spaced points in [ 1 , 1 ] are shown in the middle and lower rows of Figure 5, respectively. It is worth mentioning that in Figure 5, instead of [ 1 , 1 ] , we employed the set of pixel indices from 1 to q in both x 1 and x 2 directions, i.e., a mapping of the q Chebyshev roots in [ 1 , 1 ] onto their corresponding pixel index, i q = 1 , 2 , , q .
Regarding the error analysis of our reconstructions, we calculated the L -norm of the difference (error) of the equally spaced points ρ ˜ = ( ρ ˜ 1 , , ρ ˜ 119 ) T , defined in Equation (16), and of the almost equally spaced Chebyshev points ρ ^ ( 119 ) = ρ ^ 1 ( 119 ) , , ρ ^ 119 ( 119 ) T , defined in Equation (19), as follows:
ρ ˜ ρ ^ ( 119 ) = max ρ ˜ 1 ρ ^ 1 ( 119 ) , , ρ ˜ 119 ρ ^ 119 ( 119 ) 0.0002 , for all θ .
We also calculated the L norm of the pixel-by-pixel interpolation error of the reconstruction line profile at x 2 = x 2 * , i.e.,
f x 2 * r x 2 * : = max f ( x 1 , 1 , x 2 * ) r ( x 1 , 1 , x 2 * ) , , f ( x 1 , q , x 2 * ) r ( x 1 , q , x 2 * ) ,
where r ( x 1 , x 2 ) denotes the reconstruction, f x 2 * and r x 2 * denote the tabulated line profile of the phantom and its corresponding reconstruction at the fixed x 2 = x 2 * , respectively, in the sense that
f x 2 * = f x 1 , 1 , x 2 * , , f x 1 , q , x 2 * T , and r x 2 * = r x 1 , 1 , x 2 * , , r x 1 , q , x 2 * T .
For the purposes of our analysis, given that q = 119 , we set x 2 * in the middle of the x 2 -direction, i.e., x 2 * = 60 . Furthermore, we calculated the entrywise L 1 norm of the reconstruction error matrix, as in [30],
f r 1 : = i = 1 q j = 1 q f ( x 1 , i , x 2 , j ) r ( x 1 , i , x 2 , j ) = i = 1 q j = 1 q f i j r i j ,
as well as its corresponding entrywise L 2 norm, also known as the Frobenius norm [31], usually denoted by · F ,
f r 2 = f r F : = i = 1 q j = 1 q f ( x 1 , i , x 2 , j ) r i j ( x 1 , i , x 2 , j ) 2 = i = 1 q j = 1 q ( f i j r i j ) 2 ,
where we employed square q × q reconstructions, e.g., the tabulated forms of both the phantom (f) and its corresponding reconstructions (r) were considered as square matrices f i j i , j = 1 q and r i j i , j = 1 q , respectively.
The results of our reconstruction error investigation are presented in Table 2. For the uniform phantom investigated, the results demonstrate that the reconstruction provided by the novel Chebyshev-based method outperforms the traditional approach with respect to both the L 1 and L 2 norms. Furthermore, it outperforms the standard Chebyshev method with respect to the L norm for the line profile at the middle row ( x 2 * ) of the phantom. In particular, the decrease in the L , L 1 , and L 2 norms of the interpolation error from the traditional Chebyshev method to the modified Chebyshev approach was 18%, 49%, and 19%, respectively.

7. Conclusions

In the present work we have introduced a novel method for the polynomial interpolation of the Chebyshev type in almost equally spaced Chebyshev nodes. Although standard polynomial interpolation methods usually employ equally spaced points on an interval, in Chebyshev interpolation, this is not the case. Instead of equally spaced points along a line, Chebyshev interpolation, by definition, involves the Chebyshev nodes, i.e., the roots of Chebyshev polynomials, corresponding to equally spaced points along the unit semicircle. Given the non-uniform nature of the Chebyshev nodes and the fact that in most medical imaging modalities, the domain is uniformly partitioned, we extended the usual domain from [ 1 , 1 ] , where the function to be interpolated is known at q points, to [ a , a ] , a > 1 . The selection of a is such that not only exactly q Chebyshev nodes are included in 1 , 1 but also these Chebyshev nodes are almost equally spaced. Our preliminary results indicate that our novel Chebyshev-based method consistently outperforms the traditional approach with respect to the L norm. Our approach opens the way for using Chebyshev polynomials in the solution of the inverse problems arising in PET and SPECT image reconstruction.
The results of our phantom study demonstrate that the reconstruction provided by the novel Chebyshev-based method are superior to the traditional approach, with respect to the L 1 , L 2 , and L norms. In future studies, we aim to investigate the performance of our interpolation technique in clinical PET and SPECT studies, where specific potential limitations and challenges will be analyzed and discussed.

Author Contributions

Conceptualization, V.M. and N.E.P.; methodology, V.M., A.S.F., G.A.K. and N.E.P.; software and validation, V.M. and N.E.P.; writing—original draft preparation, V.M. and N.E.P.; writing—review and editing, V.M., A.S.F., G.A.K., and N.E.P.; visualization, V.M. and N.E.P.; supervision, N.E.P. All authors have read and agreed to the published version of the manuscript.

Funding

This work has been co-financed by the European Union and Greek national funds through the Operational Program Competitiveness, Entrepreneurship and Innovation, under the call RESEARCH–CREATE–INNOVATE (project code: TAEDK-06189).

Data Availability Statement

The data presented in this study are available upon request from the corresponding author.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Chebyshev, P.L. Théorie des mécanismes connus sous le nom de parallélogrammes. Mém. Acad. Imp. Sci. St.-Pétersbg. 1854, 7, 539–568. [Google Scholar]
  2. Szegő, G. Orthogonal Polynomials, 4th ed.; American Mathematical Society Colloquium Publications; American Mathematical Society: Providence, RI, USA, 1975; Volume 23. [Google Scholar]
  3. Rivlin, T.J. Chebyshev Polynomials; Dover Publications: Mineola, NY, USA, 2020. [Google Scholar]
  4. Hubert, E.; Singer, M.F. Sparse interpolation in terms of multivariate Chebyshev polynomials. Found. Comput. Math. 2022, 22, 1801–1862. [Google Scholar] [CrossRef]
  5. Pucanović, Z.; Pešović, M. Chebyshev polynomials and r-circulant matrices. Appl. Math. Comput. 2023, 437, 127521. [Google Scholar] [CrossRef]
  6. Piessens, R. Computing integral transforms and solving integral equations using Chebyshev polynomial approximations. J. Comput. Appl. Math. 2000, 121, 113–124. [Google Scholar] [CrossRef]
  7. Occorsio, D.; Ramella, G.; Themistoclakis, W. Lagrange-Chebyshev Interpolation for image resizing. Math. Comput. Simul. 2022, 197, 105–126. [Google Scholar] [CrossRef]
  8. Wang, M.; Au, F. Precise integration methods based on the Chebyshev polynomial of the first kind. Earthq. Eng. Eng. Vib. 2008, 7, 207–216. [Google Scholar] [CrossRef]
  9. Rattez, H.; Veveakis, M. Weak phases production and heat generation control fault friction during seismic slip. Nat. Commun. 2020, 11, 350. [Google Scholar] [CrossRef] [PubMed]
  10. Defferrard, M.; Bresson, X.; Vandergheynst, P. Convolutional Neural Networks on Graphs with Fast Localized Spectral Filtering. In Proceedings of the Advances in Neural Information Processing Systems, Barcelona, Spain, 5–10 December 2016; Lee, D., Sugiyama, M., Luxburg, U., Guyon, I., Garnett, R., Eds.; Curran Associates, Inc.: New York, NY, USA, 2016; Volume 29. [Google Scholar]
  11. Sobczyk, J.E.; Roggero, A. Spectral density reconstruction with Chebyshev polynomials. Phys. Rev. E 2022, 105, 055310. [Google Scholar] [CrossRef] [PubMed]
  12. Spier, N.; Nekolla, S.; Rupprecht, C.; Mustafa, M.; Navab, N.; Baust, M. Classification of polar maps from cardiac perfusion imaging with graph-convolutional neural networks. Sci. Rep. 2019, 9, 7569. [Google Scholar] [CrossRef] [PubMed]
  13. Huang, S.G.; Chung, M.K.; Qiu, A.; Alzheimer’s Disease Neuroimaging Initiative. Fast mesh data augmentation via Chebyshev polynomial of spectral filtering. Neural Netw. 2021, 143, 198–208. [Google Scholar] [CrossRef] [PubMed]
  14. Mason, J.C.; Handscomb, D.C. Chebyshev Polynomials; Chapman and Hall/CRC: Boca Raton, FL, USA, 2002. [Google Scholar]
  15. Boyd, J.P. Chebyshev and Fourier Spectral Methods; Dover Publications: New York, NY, USA, 2001. [Google Scholar]
  16. Trefethen, L.N. Approximation Theory and Approximation Practice; SIAM: Bangkok, Thailand, 2019. [Google Scholar]
  17. Runge, C. Über empirische Funktionen und die Interpolation zwischen äquidistanten Ordinaten. Z. Angew. Math. Phys. 1901, 46, 224–243. [Google Scholar]
  18. Atkinson, K.E. An Introduction to Numerical Analysis; John Wiley and Sons: Hoboken, NJ, USA, 1989. [Google Scholar]
  19. Hämmerlin, G.; Hoffman, K.H. Numerical Mathematics; Springer: Berlin/Heidelberg, Germany, 1991. [Google Scholar] [CrossRef]
  20. Bortfeld, T.; Oelfke, U. Fast and exact 2D image reconstruction by means of Chebyshev decomposition and backprojection. Phys. Med. Biol. 1999, 44, 1105. [Google Scholar] [CrossRef] [PubMed]
  21. Fokas, A.S.; Marinakis, V. Reconstruction algorithm for the brain imaging techniques of PET and SPECT. HERMIS 2003, 4, 45–61. [Google Scholar]
  22. Protonotarios, N.; Fokas, A.; Vrachliotis, A.; Marinakis, V.; Dikaios, N.; Kastis, G. Reconstruction of Preclinical PET Images via Chebyshev Polynomial Approximation of the Sinogram. Appl. Sci. 2022, 12, 3335. [Google Scholar] [CrossRef]
  23. Protonotarios, N.E.; Marinakis, V.; Dikaios, N.; Kastis, G.A. Chebyshev polynomials of the first kind and applications in tomography. In Mathematical Analysis, Differential Equations and Applications; Rassias, T.M., Pardalos, P.M., Eds.; World Sientific: Singapore, 2024. [Google Scholar] [CrossRef]
  24. Kaethner, C.; Erb, W.; Ahlborg, M.; Szwargulski, P.; Knopp, T.; Buzug, T.M. Non-Equispaced System Matrix Acquisition for Magnetic Particle Imaging Based on Lissajous Node Points. IEEE Trans. Med. Imaging 2016, 35, 2476–2485. [Google Scholar] [CrossRef] [PubMed]
  25. Droigk, C.; Maass, M.; Mertins, A. Direct multi-dimensional Chebyshev polynomial based reconstruction for magnetic particle imaging. Phys. Med. Biol. 2022, 67, 045014. [Google Scholar] [CrossRef]
  26. Clackdoyle, R.; Noo, F. A large class of inversion formulae for the 2D Radon transform of functions of compact support. Inverse Probl. 2004, 20, 1281–1291. [Google Scholar] [CrossRef]
  27. Fokas, A.; Iserles, A.; Marinakis, V. Reconstruction algorithm for single photon emission computed tomography and its numerical implementation. J. R. Soc. Interface 2006, 3, 45–54. [Google Scholar] [CrossRef] [PubMed]
  28. Kastis, G.A.; Kyriakopoulou, D.; Gaitanis, A.; Fernández, Y.; Hutton, B.F.; Fokas, A.S. Evaluation of the spline reconstruction technique for PET. Med. Phys. 2014, 41, 042501. [Google Scholar] [CrossRef] [PubMed]
  29. Gil, A.; Segura, J.; Temme, N.M. Numerical Methods for Special Functions; SIAM: Bangkok, Thailand, 2007. [Google Scholar]
  30. Song, Z.; Woodruff, D.P.; Zhong, P. Low rank approximation with entrywise 1-norm error. In Proceedings of the 49th Annual ACM SIGACT Symposium on Theory of Computing, Montreal, ON, Canada, 19–23 June 2017; pp. 688–701. [Google Scholar] [CrossRef]
  31. Liu, Z.; Vandenberghe, L. Interior-point method for nuclear norm approximation with application to system identification. SIAM J. Matrix Anal. Appl 2010, 31, 1235–1256. [Google Scholar] [CrossRef]
Figure 1. Chebyshev nodes (red) depicted as equally spaced points along the unit semicircle (blue).
Figure 1. Chebyshev nodes (red) depicted as equally spaced points along the unit semicircle (blue).
Mathematics 11 04757 g001
Figure 3. Interpolation of the exponential function f ( x ) = exp ( 5 x 2 ) on the interval [ 1 , 1 ] , with q=11, =15 and a = 10.52 .
Figure 3. Interpolation of the exponential function f ( x ) = exp ( 5 x 2 ) on the interval [ 1 , 1 ] , with q=11, =15 and a = 10.52 .
Mathematics 11 04757 g003
Figure 4. Interpolation of the rational function f ( x ) = 1 / ( 1 + 16 x 2 ) on the interval [ 1 , 1 ] , highlighting the occurrence of the Runge phenomenon, with q=11, =5, and a = 3.55 .
Figure 4. Interpolation of the rational function f ( x ) = 1 / ( 1 + 16 x 2 ) on the interval [ 1 , 1 ] , highlighting the occurrence of the Runge phenomenon, with q=11, =5, and a = 3.55 .
Mathematics 11 04757 g004
Figure 5. Uniform phantom line profiles and reconstructions. (Upper row) uniform phantom (ground truth). (Middle row) Chebyshev reconstruction with Chebyshev nodes, as in [22]. (Lower row) Chebyshev reconstruction with almost equally spaced points.
Figure 5. Uniform phantom line profiles and reconstructions. (Upper row) uniform phantom (ground truth). (Middle row) Chebyshev reconstruction with Chebyshev nodes, as in [22]. (Lower row) Chebyshev reconstruction with almost equally spaced points.
Mathematics 11 04757 g005
Table 1. Interpolation error results.
Table 1. Interpolation error results.
f p Lagrange [ 1 , 1 ] Chebyshev [ 1 , 1 ] Chebyshev [ a , a ]
Exponential0.02370.00180.0009
Rational1.17690.05880.0206
Table 2. Reconstruction error results.
Table 2. Reconstruction error results.
Chebyshev [ 1 , 1 ] Chebyshev [ a , a ]
f x 2 * r x 2 * 0.1840.151
f r 1 302.121153.160
f r 2 5.6514.597
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Marinakis, V.; Fokas, A.S.; Kastis, G.A.; Protonotarios, N.E. Chebyshev Interpolation Using Almost Equally Spaced Points and Applications in Emission Tomography. Mathematics 2023, 11, 4757. https://doi.org/10.3390/math11234757

AMA Style

Marinakis V, Fokas AS, Kastis GA, Protonotarios NE. Chebyshev Interpolation Using Almost Equally Spaced Points and Applications in Emission Tomography. Mathematics. 2023; 11(23):4757. https://doi.org/10.3390/math11234757

Chicago/Turabian Style

Marinakis, Vangelis, Athanassios S. Fokas, George A. Kastis, and Nicholas E. Protonotarios. 2023. "Chebyshev Interpolation Using Almost Equally Spaced Points and Applications in Emission Tomography" Mathematics 11, no. 23: 4757. https://doi.org/10.3390/math11234757

APA Style

Marinakis, V., Fokas, A. S., Kastis, G. A., & Protonotarios, N. E. (2023). Chebyshev Interpolation Using Almost Equally Spaced Points and Applications in Emission Tomography. Mathematics, 11(23), 4757. https://doi.org/10.3390/math11234757

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop