Truncated Hilbert Transform: Uniqueness and a Chebyshev series Expansion Approach

2020·Arxiv

Abstract

Abstract

I. Background and introduction

The investigation of the finite Hilbert transform (FHT) has been a research topic for a long time in history from both mathematical interests and practical needs in many applications. Some early works on the theoretical investigation and practical applications in fluid mechanic can be found in [1-7]. The relation between the line integral and the Hilbert transform was first derived by Gelfand and Graev in [8]. The introduction of using the FHT to reconstruct images from partial data in single photon computed tomography (SPECT) and computed tomography (CT), started from [9-11]. We published a short paper toward a simple inversion of the FHT in [12] for the use in the second step of [11]. Based on one explicit inversion formula of the FHT and the analytical continuation on an interval, one uniqueness result and the stability of the analytical continuation were obtained in [13]. With the breakthrough work on the cone-beam data reconstruction [14], the arguments of [13] can be extended to the 3D cone-beam (region of interest) ROI data reconstruction [15-19]. Readers should keep in mind that the image reconstruction in these works is reduced to the inversion of the FHT or the truncated Hilbert transform (THT). Therefore, the uniqueness of the inverse THT and methods to find solutions of the inverse THT become very important. Besides CT and SPECT, the FHT and THT can be useful in the context of exponential Radon transform for many other applications such as positron emission tomography (PET), nuclear magnetic resonance (NMR) imaging, ultrasonic tomography and Doppler tomography as discussed in [20]. The projection onto convex sets (POCS) was used in the numerical experiment in [13]. Recently the authors of [18, 21, 22] investigated the single value decomposition (SVD) methods to find solutions of the inverse THT. In this paper, we will focus on the theoretical study of the uniqueness of the THT and numerical procedures to find solutions of the inverse THT without further explaining the application background.

1 )(1)( dt ts tfsF . (1.0)

Throughout this paper, the integral at the singular point should be understood in the sense of Cauchy principal value and we always use the pair of )(to stand for a function and its Hilbert transform. Certainly, for a function )(can be defined on ),( which is the conventional Hilbert transform. The FHT concerns the mapping between functions defined on an finite interval such as )1,1(. The FHT has many nice properties, for example,

hold in the weighted ]. In the case of the inverse THT, )(sF is only known on an interval ),( ba which is different from )1,1(, then the goal is to find )(tf using available )(. The existence of )(tf is readily available from the FHT, but the uniqueness of

the partial uniqueness of )(tf on certain subinterval of )1,1(, for example, the result of [13] is that if )(tf has support on ],1[ is known on )1,(uniquely determined in ]( . The main arguments used in [13] include the explicit inversion formula of the FHT and the analytical continuation on an interval. These arguments can be extended to study the interior problem of CT image reconstruction in [15-18]. The analytical continuation is by and large a mathematical procedure and lacks a numerical procedure to be realized. These uniqueness results are not applicable to the data settings assumed in [21, 22] to find )(tf in its entire support interval.

In this paper, using the Sokhotski–Plemelj formulas, we derive a stronger uniqueness result for )(tf on its support interval )1,1(weaker conditions than the ones found in previously mentioned works. The result in this paper confirms that the solution by the SVD schema of [21, 22] can be unique under the assumptions in [13, 15-18]. Moreover, our uniqueness results are valid to the exponential Radon transform. We also describe one example of using such stronger uniqueness result to resolve the typical data truncation problem in tomographic imaging applications when the object is out of the field of view. We continue to explore the Chebyshev polynomials to evaluate the Hilbert transform and its inversion mentioned in [5, 25] and further studied in our recent work [24]. Surprisingly an explicit inversion can be achieved through Lagrange interpolations if the functions only include finite Chebyshev polynomials. Using the explicit Chebyshev polynomials can avoid the estimate of the eigenvalues and eigenfunctions needed in the SVD schema of [21, 22]. Numerically, the evaluation of the Chebyshev polynomials can be implemented by fast sine and cosine transformations. To numerically find a solution to the THT problem, we express )(in the finite Chebyshev polynomial series, and then try to estimate the coefficients using the available data from both )(Two methods are proposed to estimate the coefficients, one is to minimize a cost function and the other is an extrapolation procedure. Finally, we present computer simulation results to show the feasibility of the proposed method for practical applications.

II. Uniqueness of THT and cosh-weighted Hilbert transform

For a function with compact support in 1R , through scaling and shift, we can always make it have the support on ]1,1[. Without loss of generality in this paper, we assume that )(tf belongs

to )1,1(which is defined as

follows:

t tftfLd (2.2)

with the inner product

1 )()(1, dt t tgtfgf d

Hereafter, we use the superscript star for the conjugate operation of a function. Define the closed subspace )1,1(2

We will use C for the complex plane and 1for the imaginary unit. The Sokhotski–Plemelj formula is a historical result and its detailed proof can be found in many graduate text books such as [4] on the integral operator theory. In Proposition 1, we reformulate the results on the Sokhotski-Plemelj formulas under the conditions for the use in this paper.

)( 2 1)( dt zt tf i z . (2.6)

Proof Formulas (2.7) and (2.8) are originally derived by Sokhotski and Plemelj for Hölder continuous functions. The existence of )(lim

was obtained in Privalov’s series of papers. Detailed expositions and derivations of these results for integrable functions can be found from [4] and [26].

determined in )1,1(under either one of the following conditions:

Let )(be the function defined by (2.6). Notice that in case 1, there is no jump in (2.7) and (2.8) because of 0)( is continuous in ),( ba in both cases. It follows that )(is analytical in ),(]}1,1[\{ . Recall the fact that an analytic function on an open connected set is completely determined by its values on any set of points containing a limit point

Remark 1. Theorem 1 reveals that )(can be obtained from their partial data while it is still lacking a numerical procedure to find them. Nonetheless, Theorem 1 yields the uniqueness of )(from very little data. We give one example in using Theorem 1 as follows. Let

branch cut for 1, respectively. On the branch cuts, we have

Remark 2. We point out that the condition that )(is known on a same interval cannot be reduced in Theorem 1. We provide one counterexample if the condition of Theorem 1 is not met. For a fixed positive number 0, we claim that there are infinite numbers of )(such that their Hilbert transforms are equal to

01 . Through shift and scaling, we change (2.9) to the following pair of functions

s sss s sF (2.13)

transform. From the definitions of (2.12) and (2.13), it is easy to see that there are many different pairs )(are known when

In SPECT, PET and other applications, the projection data in some cases can be expressed

as the exponential Radon transform of certain unknown function to be solved. To deal with the truncated data as investigated in [29], one approach is to invert the cosh-weighted Hilbert transform )(

)1,1(under either one of the following conditions:

Using the same arguments in the proof of Theorem 1, we prove the corollary.

Out-of-field-of-view in CT/SPECT scan. To conclude this section, we apply Corollary 1 to one of the typical truncation problems in CT and SPECT when the object to be scanned is out of field of view. For SPECT, the exponential Radon transform (ERT) ),,( can be obtained by multiplying the measurement with an exponential factor that is decided by the attenuation coefficient and the distance between the center of the rotation and the object’s edge of the SPECT detector. Assume that the detector size is 2 and the object is measured within the unit circle, as shown as a thick circle in Figure 1. A typical truncation problem in SPECT can be shown in Figure 1 when the support of ),( yxf is inside an ellipse whose major-axis (i.e., the longest diameter) is longer than 2 and minor-axis (i.e., the shortest diameter) is less than 2. Notice that the object is fully covered when the detector surface is parallel to x-axis while the outer part is not covered when the detector surface is parallel to y-axis.

Figure 1. A SPECT imaging configuration with an elliptical object. The object is measured within the unit circle.

Let

. We refer to [29] for detailed data formation of ),( in SPECT scanning. Let ),( as the weighted

backprojection of the derivative of the ERT:

Let us consider an arbitrary thick horizontal line between the two dashed lines in Fig. 1. On this thick line, ),( does not include truncation, so one is able to accurately reconstruct ),( yxf by inverting an FHT defined by (2.20). There are a lot of numerical results for the inverse FHT in our paper [29]. Let 0y be the distance between the dashed line and x-axis, then ),( known for 0|| is still unknown on the points away from two dashed lines

but inside the ellipse. We consider

1 is met, it follows that ),( yxf can be uniquely determined for 0|| . In summary, we first use the inverse FHT to reconstruct ),( along lines parallel to x-axis, then use the THT to uniquely determine the remaining part along lines parallel to y-axis. If 0example is reduced to a typical X-ray CT truncation data because the bordering part of the object is out of the field of view. We summarize the preceding scanning settings to Corollary 2.

Proof. First we apply Corollary 1 to (2.20), the uniqueness can be derived on the strip defined by , then apply Corollary 1 to (2.21) to derive the uniqueness for . It follows that the ellipsoid is fully covered.

Proof. First we apply Corollary 1 to (2.20), the uniqueness can be derived on the strip defined by ),min(|| 21 , then apply Corollary 1 to (2.21) to derive the uniqueness for ),min(|| 21

III. Chebyshev polynomial expansion of f(t) and F(s)

In this section, we only consider the THT without cosh weights. Based on the analyticity of )(sF , we give a formal statement on the THT problem with three conditions as follows:

We point out that conditions C2 and C3 can be converted into condition C1 through shift and scaling operations, but such conversion will lose the analytical properties of )(sF in different intervals. It is this reason that we use three conditions to represent the truncation problems studied in the literature. According to Theorem 1, mathematically there exists a unique solution to above truncation problem. However, to numerically find the solution remains a challenge since the analytical function )(exists but a numerical procedure is lacking to compute the solution. To the authors’ knowledge, no any sort of explicit formula has been discovered to find

methods in [21, 22] lack of uniqueness and are virtually impossible to lead to a numerically stable procedure since both eigenvalues and eigenfunctions do not have explicit expressions.

In this paper, we propose two methods to numerically find the solution. The key technique is to use the Chebyshev polynomial series expansion. We express )(in the series of Chebyshev polynomials, then convert finding )(tf to estimating the coefficients through the

minimization principle and an extrapolation procedure. First we introduce some results on the Chebyshev polynomials. Let 0)(1

))sin(acos( ))acos()1sin(()( t

Here )(are the Chebyshev polynomials of the first and second kind, respectively. Over ]1,1[, we have those orthogonal relations:

nm nm dtttUtU nm

From [5], the FHT of the Chebyshev polynomials on ]1,1[satisfies the following relations

seen the exact expression. Here we provide a proof based on the recurrence relations. Lemma 1 (Extension of Chebyshev polynomials).

Proof. The recurrence relations for the Chebyshev polynomials are

When 1, using the orthogonal relations in (3.4) we have

Combining equations (2.9, 3.10-3.12), we have proven (3.7). Similar arguments lead to (3.8).

A recent comprehensive review on the Chebyshev polynomials can be found in [25]. We collect several known results in Lemma 2 for the use in this paper. Lemma 2 (Properties of Chebyshev polynomials).

In the sense of approximation, on ]1,1[, we may have the following relations

1 2 )(1)( , (3.17)

For the polynomial interpolation, we introduce the following coordinate transformation

expressed as a finite series

1 2 )(1)( , (3.20)

Proof. With the assumption of (3.20), )(sF can be expressed as

expression

interpolation polynomial on these points

Rewriting (3.23) to (3.22), we obtain the coefficients }{ . It follows that }{ can be obtained by expressing in the Chebyshev polynomial series.

Here we used the following identity

Selecting N+1 distinct sampling points in )1,1( 22 , and then by the Lagrange interpolation polynomial, we can obtain }{ . The condition C3 can be converted into either C1 or C2. This completes the proof.

The Lagrange interpolations are notoriously ill-posed for high-order polynomials. In the finite dimension, the explicit procedure to find )(tf may be only meaningful in the perfect world. In practical applications such as CT/SPECT imaging, )(sF may contain measurement and computation errors and will be approximated by (3.17). In order to fully use the available data from both )(, we suggest a minimization criterion to estimate the coefficients for all three conditions of C1, C2 and C3.

, the Chebyshev polynomials have the orthogonal relations (3.3) and (3.4). Outside of ]1,1[, we only have a global Plancherel formula for )}(~),(~{

in

Here we point out that the last term of (3.27) outside of ]1,1[takes polynomials whose term is contracting due to (3.25). Also we use the Stieltjes type integral through )acos(td to avoid the singularity at the end points of ]1,1[norm is used while the 2L norm is

adopted outside of ]1,1[. Equation (3.27) is one way to define a cost function. There may exist

other optimized ways to define the cost function concerning a prior information.

be approximated by

Looking for }{ to minimize })({ leads to an estimate of )(tf in the expression of (3.17). Cost function (3.28) can be used to find the solution in the intervals described in [30]. Numerical algorithms in [31] can be used to find }{ to minimize (3.28). Let }ˆ{ be one set of

coefficients to reach the minimum of })({ , then the estimated solution is

t tftf . (3.30)

Extrapolation procedure. Condition C1 is one of the typical truncation problems arising in CT/SPECT as previously described in Section II for the out-of-field-of-view in CT/SPECT scanning. Motivated by the iterative schema for the band-limited signal extrapolation in [32, 33], we construct an iterative procedure to find )(tf under condition C1. We denote by H the finite

Hilbert transform (1.0), and let 1be the following inversion formula

Then we suggest the following iterative procedure

For the band-limited signal extrapolation, the operator H is the Fourier transform [32, 33]. In [34], the author found one minor issue of [33] in using that schema. This iteration also falls in the concept of the POCS. Using the Plancherel formula in (2.5), if )()()( , we have the

following strictly decreasing relations

However we are not able to prove the convergence in )1,1(2 . In order to use the fast cosine

and sine transformations of [31, 35] to calculate the Hilbert transform and its inversion, for

According to the trigonometric expression of the FHT and its inversion in [24], on the CGL sampling points, we have

be expressed as the combination of sine and cosine transforms. Then the discrete versions of (3.34-3.35) become

finite dimension of Euclid space is compact.

IV. Computer simulation results

As pointed out in the beginning of Section III, conditions C2 and C3 can be always converted to condition C1, and C1 is the common truncation problem arising in CT/SPECT scanning. Also the extrapolation procedure for condition C1 is very simple to numerically realize. In this paper, we will perform the computer simulations for the iterative procedure of (3.42, 3.43) for condition C1. Through scaling and shift on (2.9), we construct the following pair of functions

otherwise tttf (4.1)

Restricted on ]1,1[construct a pair of a function and the Hilbert transform. We will use (4.1) and (4.2) to verify the iterative procedure of (3.42, 3.43). The CGL sampling grid of }{ ) will be used in the computer simulations. The number of sampling points is N=256 in all numerical simulations. Since }{ is not evenly sampled on ]1,1[display purpose, we will resample )(resampling method is to first find the coefficients }{and then use the following recurrence

formulas

following recurrence formulas

In the numerical experiemnts, )( will be used as the initial data in (3.40, 3.41). Then we repeat 30 iterations of (3.42, 3.43) to extrapolate )( in the outer region. The extrapolated )(are shown in Figs. 2 and 3.

Fig 2. Original and extrapolated )(

Fig 3. Original and extrapolated )(

We mention that the above truncation data setting is the condition used in [22]. As shown in Figs 2 and 3, the iterative procedure (3.42, 3.43) does extrapolate both )(to the outer region in an accurate manner. Even under moderate level of noise, the extrapolation procedure still works well as shown in the Figs 4 and 5.

Fig 4. Original and extrapolated noisy )(

Fig 5. Original and extrapolated noisy )(

In conclusion, the iterative procedure (3.42, 3.43) is easy to carry out and seems to provide a reasonable solution for the truncation problem under condition C1.

V. Discussion and Conclusion

In summary, we have obtained a stronger uniqueness result under weaker conditions compared with the existing works [13, 15-18]. The arguments used in [13, 15-18] cannot handle the uniqueness in the context of exponential Radon transform. Because of the powerful Sokhotski– Plemelj formulas, our uniqueness result can be applicable to the exponential Radon transform. As an application of our uniqueness result, the solution by the SVD schema in [21, 22] can be unique under the conditions used in [13, 15-18]. Following the early idea in [5], we have found that the Chebyshev polynomials can be used to construct an SVD for the truncated Hilbert transform from functions of ]1,1[to functions of 1R . The Chebyshev polynomials have explicit expressions and can be computed using fast algorithms such as sine and cosine transforms [31, 35, 36]. For condition C1, the implementation of iterative procedure (3.42, 3.43) is very easy and the computer simulation results are very promising. From the numerical realization standpoint, using the Chebyshev polynomial series expansion has an advantage over other existing methods.

The SVD schema in [21, 22] was obtained from )1,1(2 to a subspace. As pointed out in [23], for 1)( , its inverse FHT 21/ does not belong to )1,1(2 . This implies

that the SVD scheme in [22] may not produce a convergent series. In this paper, we consider

Chebyshev polynomial series expansions (3.15) and (3.16), a prior knowledge of )(tf can be considered in both minimization criterion (3.27) and the iterative procedure (3.42, 3.43).

The minimization criterion is more general and can be applicable to the THT data settings considered in [13, 15-18, 21, 22] through selecting the location of ],[ ba , for example the case in [13] is the interval with )1,1(, the case in [15-18, 22] is the interval with )1,1(, , and the case in [21] is the interval with ),1(, . However, finding the minimum point of (3.28) may not be numerically easy as the scheme (3.42, 3.43). Theorem 2 may be not very useful to the numerical realization, but it can be used to obtain an initial estimate of }{ using the low-order polynomial interpolation toward finding the final solution of }{ the minimum of })({ and providing a reasonable initial guess of )(sF . Due to the evaluation of high-order polynomials, to numerically find the minimum point of (3.28) is not trivial. For example, the double type only has 15~16 exact digits, then for 15, the coefficients of may not be accurate because of round-off error. We used equally spaced sampling points in (3.28) to show a simple discrete version of (3.27), however this type of sampling causes the Runge Phenomenon near the boundary inside ]1,1[. We mention that [36] includes a graphic tutorial on the numerical behavior of the Chebyshev polynomials. It may be interesting to study the numerical behaviors of the minimization criterion compared with the extrapolative procedure (3.42, 3.43).

To conclude the paper, we like to mention two papers [37, 38] on the explicit inversion formulas of cosh-weighted Hilbert transform in [37] and numerical approximation in [38]. Currently we are investigating the numerical characteristics of the inversion formulas of [37].

Reference

[2] Koppleman W and Pincus JD, “Spectral representations for finite Hilbert transform,” Math Z. vol. 71, pp. 399-407, 1959.

[4] Muskhelishvili NI, Singular Integral Equations: boundary problems of function theory and their applications to mathematical physics, Wolters-Noordhoff, The Netherlands, 1953.

[8] Gelfand IM and Graev MI, “Crofton function and inversion formulas in real integral geometry,” Funct. Anal. Appl. vol. 25 pp. 1–5, 1991

[9] Rullgård H, “An explicit inversion formula for the exponential Radon transform using data from 180 degrees,” Ark. Math., vol. 42, pp. 353-362, 2004.

[10] Zou Y and Pan X, “Exact image reconstruction on PI-lines from minimum data in helical cone-beam CT,”

[11] Noo F, Clackdoyle R, and Pack JD, “A two-step Hilbert transform method for 2D image reconstruction”,

[13] Defrise M, Noo F, Clackdoyle R and Kudo H. “Truncated Hilbert transform and image reconstruction from

[14] Katsevich A, “Theoretically exact FBP-type inversion algorithm for spiral CT,” SIAM J. Appl. Math., vol. 62,

[15] Ye Y, Yu H, Wei Y and Wang G, “A general local reconstruction approach based on a truncated Hilbert

[16] Kudo H, Courdurier M, Noo F and Defrise M, “Tiny a priori knowledge solves the interior problem in

[17] Courdurier M, Noo F, Defrise M and Kudo H, “Solving the interior problem of computed tomography using

[18] Ye Y B, Yu H Y and Wang G, “Interior reconstruction using the truncated Hilbert transform via singular

[19] Yang J, Yu H, Jiang M and Wang G, “High-order total variation minimization for interior SPECT,” Inverse

[20] You J, “The attenuated Radon transform with complex coefficients,” Inverse Problems, vol. 23, pp. 1963-

[21] Katsevich A, “Singular value decomposition for the truncated Hilbert transform,” Inverse Problems, vol. 26,

[22] Katsevich A, “Singular value decomposition for the truncated Hilbert transform: part II,” Inverse Problems,

[24] You J, “Finite Hilbert transform in weighted L2 spaces,” arXiv: 2002-02071 (Cubic-Imaging LLC, Preprint-

[26] Bagemihl F, “The Plemelj formulas with unrestricted approach and the continuity of Cauchy-type integrals,”

[29] Huang Q, You J, Zeng GL, and Gullberg GT, “Reconstruction from uniformly attenuated SPECT projection

[30] Alaifari R, Pierce L and Steinerberger S, "Lower bounds for the truncated Hilbert transform", arXiv:

[31] Press WH, Teukolsy SA, Vetterling WT, and Flannery BP, Numerical Recipes in C, Cambridge University

[32] Papoulis A, A new algorithm in spectral analysis and band-limited extrapolation, IEEE Trans. on Circuits

[33] Cadzow JA, “An extrapolation procedure for band-limited signals,” IEEE Trans. Circuits Sys., vol. 27, pp. 4-

[34] You J, “A note on an extrapolation procedure for band-limited signals by Cadzow,” IEEE Trans. Signal

[35] Yip P and Rao KR, “Fast decimation-in-time algorithms for a family of discrete Sine and Cosine transforms,”

[36] Sarra SA, “Chebyshev Interpolation: An Interactive Tour,” Journal of Online Mathematics and its

[37] Bertola M, Katsevich A, and Tovbis A, "Inversion formulae for the cosh-weighted Hilbert transform",

[38] Boche H and Pohl V, "Limits of calculating the finite Hilbert transform from discrete samples," Applied and