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
they are orthogonal, i.e., pairwise perpendicular, in the
-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
-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
. We extend the domain to
,
, and select a sufficiently large value of
a, such that
exactly q Chebyshev nodes are included in
, 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.
3. Chebyshev Interpolation with Almost Equally Spaced Points
We aim to approximate an unknown function
, defined in
, via a polynomial of degree of at most
, given that the values of
are known at the following
q equally space points:
Hence,
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
in the roots of
are unknown.
To overcome this difficulty, we extend the domain from
to
,
, as in
Section 2.2. Let
n denote the number of roots of the corresponding Chebyshev polynomial
in
(see Equation (
13)). We require that exactly
q of the
n Chebyshev roots in
are included in
. It is worth mentioning that the roots of
, defined in (
5), and of
, 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
are almost equally spaced.
Proposition 1. Let be such that exactly of the n Chebyshev roots in , defined in Equation (13), are included in , and let ℓ be an odd integer, such that . 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
, defined in Equation (
13), are included in
, it follows that there exists an integer
, such that the number of Chebyshev roots in
, and similarly for
, is exactly
(see
Figure 2). Hence,
is such that
Taking into account Equation (
18), the
q Chebyshev roots in
via Equation (
13) may be rewritten as
Given that
, where
ℓ is an odd integer, Equation (
18) yields
i.e.,
is a multiple of
q. Inserting Equation (
20) in Equation (
19) implies
Since
q of the
n roots
are included in
,
and
, it is reasonable to require that
. Thus, Equation (
21) implies
Furthermore, taking into account Equation (
22), we may rewrite Equation (
21) in the following form:
Taking the corresponding Taylor series yields
Hence, for
ℓ (and therefore
n) sufficiently large, the Chebyshev nodes are almost equally spaced, since
, for all
m, as implied by Equations (
16) and (
24). □
Figure 2.
Chebyshev roots distribution in the extended domain , .
Figure 2.
Chebyshev roots distribution in the extended domain , .
Since the approximation of the unknown function
involves the interval
, and not
, we may assume that
has compact support, namely
, for
. Thus, we may employ Chebyshev interpolation in
, as in Equation (
15), utilizing Equation (
14), assuming
, for
, and
, for any other root of
. 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, ℓ, given at with |
Compute: a via Equation (22) |
Assume: , for , and , for any other root of |
Compute: via Equation (14) |
Output: approximation of 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
and the rational function
which is well known regarding the so-called
Runge’s phenomenon.In order to approximate defined on the interval , we employed three different polynomial interpolation schemes, namely,
Lagrange interpolation with equally spaced points on
(e.g., [
19]);
Chebyshev interpolation with Chebyshev nodes on
, 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
interpolation points on
. For the third, we further assumed an extended interval of the form
,
. This interval is such that exactly
almost equally spaced Chebyshev nodes are included in
. Furthermore, according to Proposition 1, we assumed
,
ℓ odd. In the case of the exponential function, we chose
, hence
Chebyshev nodes are lying on
, which, according to Equation (
22), implies that
. In the case of the rational function, choosing
implies that
Chebyshev nodes are lying on
, which, in turn, implies that
(see Equation (
22)). We then performed Chebyshev interpolation on the interval
via Equations (
14) and (
15). It is worth mentioning that in Equation (
14) we assumed
,
, and
, for all other
roots of
on
(see Equation (
18) and
Figure 2). All methods of interpolation investigated for both examples of
are shown in
Figure 3 and
Figure 4, respectively, where blue represents the function itself, and red corresponds to its approximation.
We calculated the
norm of the error of the equally spaced points
, defined in Equation (
16), and of the almost equally spaced Chebyshev points
, defined in (
19), i.e.,
Upon designating the new Chebyshev nodes as “almost equally spaced”, it is evident, from the above
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
norm of the difference is 0.03% and 0.265% of the length of the integral
for the exponential and rational functions, respectively.
We calculated the
norm of the interpolation error, namely
, where
p denotes the interpolating polynomial, with
for the Lagrange and standard Chebyshev cases, and
for the almost equally spaced points Chebyshev case. This error norm was computed by taking 10,000 points along the interval
. 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
norm. In particular, for the case of the exponential function, the decrease in the
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
where
H denotes the Hilbert transform of the sinogram,
, i.e.,
where, as commonly attributed to medical imaging (see for example [
26]), we assume that
has compact support on
, and that
[
27,
28].
In order to apply our method, we consider an extension of the usual sinogram (
). Let
denote the extension of
from
to
,
, defined as follows:
Then, according to Equation (
15), the Chebyshev-type expansion of the first-kind
is given by
where
are defined in Equation (
11). Hence, the Hilbert transform of
via Equation (
28) may be rewritten in terms of
as follows:
We denote the integral of the last line of Equation (
31) by
, namely,
Taking into account Equation (
12a), it is easy to establish that
Inserting the recursive Relation (
12b) in Equation (
32) yields
where
We modify
by adding and subtracting the term
on the numerator inside the integral, as follows:
where
denotes the following definite integral
Inserting Equation (
36) in Equation (
34) yields the following recursive relation:
We note that it can be shown that
By differentiating Equation (
38) with respect to
, we obtain
where prime denotes differentiation. Equation (
33) implies that
vanishes for
and
, i.e.,
Taking into account all of the above, the derivative of the Hilbert transform of the sinogram may be written as
where
denotes the derivative of
with respect to
, and
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
i.e., a uniform circular phantom. This phantom and its line profile for
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
, and is denoted by
, namely,
where
,
are the usual Cartesian coordinates,
represents the signed distance from the origin, i.e.,
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.,
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
, 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
. Furthermore, for the case of standard Chebyshev approximation of the sinogram where
, we considered
q=119 values of
, at the roots of the corresponding Chebyshev polynomials of degree
q,
. For the case of the Chebyshev approximation with
q almost equally spaced points on
, we assumed an interval of the form
for the sinogram, thus extending its domain of definition as in Equation (
29). Then, based on the relation
and the fact that
points are required to lie on the interval
, we chose
, thus
, which, via Equation (
22), implies
. The reconstruction results for the standard Chebyshev approximation with Chebyshev nodes in
, and for the Chebyshev with almost equally spaced points in
are shown in the middle and lower rows of
Figure 5, respectively. It is worth mentioning that in
Figure 5, instead of
, we employed the set of pixel indices from 1 to
q in both
and
directions, i.e., a mapping of the
q Chebyshev roots in
onto their corresponding pixel index,
.
Regarding the error analysis of our reconstructions, we calculated the
-norm of the difference (error) of the equally spaced points
, defined in Equation (
16), and of the almost equally spaced Chebyshev points
, defined in Equation (
19), as follows:
We also calculated the
norm of the pixel-by-pixel interpolation error of the reconstruction line profile at
, i.e.,
where
denotes the reconstruction,
and
denote the tabulated line profile of the phantom and its corresponding reconstruction at the fixed
, respectively, in the sense that
For the purposes of our analysis, given that
, we set
in the middle of the
-direction, i.e.,
. Furthermore, we calculated the entrywise
norm of the reconstruction error matrix, as in [
30],
as well as its corresponding entrywise
norm, also known as the Frobenius norm [
31], usually denoted by
,
where we employed square
reconstructions, e.g., the tabulated forms of both the phantom (
f) and its corresponding reconstructions (
r) were considered as square matrices
and
, 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
and
norms. Furthermore, it outperforms the standard Chebyshev method with respect to the
norm for the line profile at the middle row (
) of the phantom. In particular, the decrease in the
,
, and
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 , where the function to be interpolated is known at q points, to , . The selection of a is such that not only exactly q Chebyshev nodes are included in 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 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 , , and 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.