About this Journal Submit a Manuscript Table of Contents
International Journal of Biomedical Imaging
Volume 2011 (2011), Article ID 285130, 16 pages
http://dx.doi.org/10.1155/2011/285130
Research Article

Inverse Fourier Transform in the Gamma Coordinate System

Yuchuan Wei,1,2 Hengyong Yu,2,3,4 and Ge Wang2,4

1Department of Radiation Oncology, Wake Forest University School of Medicine, Winston-Salem, NC 27157, USA
2Biomedical Imaging Division, VT-WFU School of Biomedical Engineering and Sciences, Wake Forest University Health Sciences, Winston-Salem, NC 27157, USA
3Department of Radiology, Division of Radiologic Sciences, Wake Forest University Health Sciences, Winston-Salem, NC 27157, USA
4Biomedical Imaging Division, VT-WFU School of Biomedical Engineering and Sciences, Virginia Tech., Blacksburg, VA 24061, USA

Received 28 May 2010; Accepted 29 July 2010

Academic Editor: Yangbo Ye

Copyright © 2011 Yuchuan Wei et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

This paper provides auxiliary results for our general scheme of computed tomography. In 3D parallel-beam geometry, we first demonstrate that the inverse Fourier transform in different coordinate systems leads to different reconstruction formulas and explain why the Radon formula cannot directly work with truncated projection data. Also, we introduce a gamma coordinate system, analyze its properties, compute the Jacobian of the coordinate transform, and define weight functions for the inverse Fourier transform assuming a simple scanning model. Then, we generate Orlov's theorem and a weighted Radon formula from the inverse Fourier transform in the new system. Furthermore, we present the motion equation of the frequency plane and the conditions for sharp points of the instantaneous rotation axis. Our analysis on the motion of the frequency plane is related to the Frenet-Serret theorem in the differential geometry.

1. Introduction

In [1], we developed a general scheme for 2D and 3D, parallel- and divergent-beam computed tomography (CT). Different from traditional Radon's or Tuy's formulas starting with the inverse Fourier transform in the spherical coordinate system, our new framework is based on an instantaneous cylindrical coordinate system. However, in [1], the coordinate system was not formally defined, and the Jacobian of the transform was not explicitly calculated. Given the basic role of our general scheme in CT theory and some confusion in the literature [25], here, we treat the Γ coordinate system in a mathematically strict way.

The organization of this paper is as follows. In Section 2, we demonstrate how the inverse Fourier transform formula generates different reconstruction formulas, with a special attention on different properties of projection truncation. In Section 3, we define the Γ coordinate system and deduce its Jacobian factor. In Section 4, we study the motion of the frequency plane and present the condition for sharp points of the instantaneous rotation axis. In Section 5, a simplified model for helical CT reconstruction is given. In the appendixes, we present a weighted Λ reconstruction formula, a weighted Radon's formula, a deduction of Orlov's formula, and an alternative way to study the motion of the frequency plane. For convenience, we use the same notations as in [1].

2. Inverse Fourier Transform in Commonly Used Coordinate Systems

In this section, we write the inverse Fourier transform in commonly used coordinate systems, derive corresponding reconstruction formulas, and demonstrate their difference on truncation of projection data.

As shown in Figure 1, in the 3D space 𝑅3,Ψ(𝑟) is an object function to be reconstructed, whose Fourier transform is Ψ(𝑘). According to Fourier analysis, we have Ψ=𝑟Ψ𝑘exp2𝜋𝑖𝑘𝑟𝑑𝑘1𝑑𝑘2𝑑𝑘3,(1) where 𝑟=(𝑥,𝑦,𝑧) and 𝑘=(𝑘1,𝑘2,𝑘3) are 3D vectors in the real space and Fourier domain, respectively. The unit vectors along the three axes are i,j, and k. Note that the frequency vector 𝑘 and the unit vector k are different. This formula states that to reconstruct an 3D object Ψ(𝑟), we need its Fourier transform Ψ(𝑘) in the whole 3D space. In X-ray CT, the Fourier transform Ψ(𝑘) can be calculated when sufficiently many parallel projections are measured.

fig1
Figure 1: 3D reconstruction problem in parallel-beam geometry and the cylindrical coordinate system. (a) A parallel-beam scans an object along a semicircle locus on the unit sphere; (b) In the frequency domain, when 𝜙0 increases from 0 to 𝜋, the frequency plane 𝜙=𝜙0 scans every point exactly once except for the points on the 𝑘3 axis, and the end point of the normal vector 𝑂𝑆 moves along the equator from (0,1,0) to (0,1,0); (c) For simplicity, we suppose that the detector panel passes through the origin.

Now, let us specify the simplest 3D reconstruction problem. In Figure 1(a), 𝑆(𝜙)=𝑆(sin𝜙,cos𝜙,0) with 𝜙[0,𝜋] is a point on a great semicircle of the unit sphere. In the 3D space, we define the unit vectors 𝑒=𝑂𝑆=(sin𝜙,cos𝜙,0) and 𝑒=(cos𝜙,sin𝜙,0). Then, 𝑒,𝑒 and k are another set of orthogonal unit vectors. As shown in Figure 1(c), the new components of the vector 𝑟 can be denoted as 𝑢=𝑟𝑒,𝜌=𝑟𝑒,𝑧=𝑟k.(2) The projection at point 𝑆 is defined by 𝑃𝑆(𝜌,𝑧)=𝑃𝜙(𝜌,𝑧)=Ψ𝑢𝑒+𝜌𝑒k+𝑧𝑑𝑢.(3) Our purpose is to reconstruct Ψ(𝑟) when 𝑃𝜙(𝜌,𝑧) is known for all 𝜙[0,𝜋]. Let us see how the selection of a coordinate system will determine the reconstruction formula, paying an attention to longitudinal data truncation.

2.1. Signed Cylindrical Coordinate System

In a cylindrical coordinate system, the inverse Fourier transform can be rewritten as Ψ=𝑟𝜋0Ψ𝑘||𝑘exp2𝜋𝑖𝑘𝑟𝜌||𝑑𝑘𝜌𝑑𝑘3𝑑𝜙,(4) where (𝑘𝜌,𝑘3,𝜙) with <𝑘𝜌<,<𝑘3<,0𝜙<𝜋 is the signed cylindrical coordinates of 𝑘. As shown in Figure 1(b), we have 𝑘𝑘𝑘=𝜌,𝑘3,𝜙=𝑘𝜌𝑒(𝜙)+𝑘3k.(5) When 𝜙0 increases from 0 to π, the plane 𝜙=𝜙0 passes through every frequency point exactly once (except for the points on the axis 𝑂𝑘3 which have a zero measure), and the normal vector 𝑂𝑆 of the plane 𝜙=𝜙0 moves along a great semicircle on the unit sphere with the starting point A(0,1,0) and the end point B(0,1,0). Note that the coordinate transform from (𝑘1,𝑘2,𝑘3) to (𝑘𝜌,𝑘3,𝜙) is one to one (except for the points on the axis 𝑂𝑘3 with zero measure), and the absolute value of Jacobian is |𝑘𝜌|.

From the parallel projection 𝑃𝜙(𝜌,𝑧), we can calculate the Fourier transform of the object Ψ𝑘𝜌,𝑘3=𝑃,𝜙𝜙𝑘𝜌,𝑘3=𝑃𝜙𝑘(𝜌,𝑧)exp2𝜋𝑖𝜌𝜌+𝑘3𝑧𝑑𝜌𝑑𝑧.(6) Substituting (6) into (4), we have Ψ=𝑟𝜋0𝑃𝜙𝑘𝜌,𝑘3||𝑘exp2𝜋𝑖𝑘𝑟𝜌||𝑑𝑘𝜌𝑑𝑘3𝑑𝜙.(7) The reconstruction scheme behind this formula can be denoted as 𝑃𝜙Ψ𝑘(𝜌,𝑧)Ψ𝑟,(8) which indicates that we first calculate all the Fourier components from projections and then reconstruct the object using the inverse Fourier transform. In this process, the object is reconstructed as a whole body.

In fact, Formula (7) can be written in the filtered backprojection form as Ψ=𝑟𝜋0𝑃𝜙(𝜌,𝑧)𝑑𝜙,(9) where the filtered projection is given by 𝑃𝜙=(𝜌,𝑧)𝑃𝜙𝑘𝜌,𝑘3||𝑘exp2𝜋𝑖𝑘𝑟𝜌||𝑑𝑘𝜌𝑑𝑘3=𝑃𝜙𝑘𝜌,𝑘3exp2𝜋𝑖𝑘𝜌𝜌exp2𝜋𝑖𝑘3𝑧||𝑘𝜌||𝑑𝑘𝜌𝑑𝑘3=𝑃𝜙1(𝜌,𝑧)2𝜋2𝜌2.(10) The symbol denotes the 1D convolution operation about variable 𝜌. Since the filtration is along the horizontal direction only, the projection can be longitudinally truncated, that is, a plane with a specified 𝑧 value can be reconstructed independently. The reconstruction scheme behind (9) with (10) can be denoted by Ψ𝑘𝑃𝜙(𝜌,𝑧)Ψ𝑟,(11) which states that certain frequency components form the projection, and the object can be reconstructed directly from projections.

When projections are longitudinally truncated, the Fourier transform cannot be calculated. How and why can we reconstruct a part of the object?

The answer is as follows. What we want to reconstruct is the object Ψ(𝑟) instead of the Fourier transform Ψ(𝑘). Though we cannot calculate every single frequency component, we can calculate the contribution of the frequency components on the whole frequency plane (as well as its neighborhood) to a point to be reconstructed, see (10). In the aforementioned case, projections contain as all frequency components as the object does, but the ratio of the frequency components is not suitable. Then, the only thing we need to do is to perform a 1D filtration to adjust this ratio. An analog is that the Fourier components can be considered as raw material, projections a semifinished product, and the reconstructed object the final product. If we use the semifinished product to make a product, we do not need change it back to the raw material. The role of Fourier analysis is to tell us how far from the semifinished product to the final product. Simply speaking, the essence of image reconstruction from projections is filtration. Therefore, a suitable coordinate system should show the character of the filtration in an easy way.

2.2. Signed Spherical Coordinate System

In the signed spherical coordinate system shown in Figure 2, the inverse Fourier transform can be rewritten as Ψ=𝑟𝜋0𝜋/2𝜋/2Ψ𝑘𝑘exp2𝜋𝑖𝑘𝑟2cos𝜗𝑑𝑘𝑑𝜗𝑑𝜙,(12) where (𝑘,𝜗,𝜙) with <𝑘<,𝜋/2<𝜗<𝜋/2,0𝜙<𝜋 is the signed spherical coordinates of 𝑘.

fig2
Figure 2: The sphere coordinate system and the detector plane. (a) The sphere coordinate system in the Fourier domain, and (b) the relationship between projections and the Radon transform.

For convenience, we denote 𝑘=𝑘(𝑘,𝜗,𝜙)=𝑘𝑛, where 𝑛=𝑛(𝜗,𝜙)=𝑒(𝜙)cos𝜗+ksin𝜗 is a unit vector.

For a unit vector 𝑛=𝑛(𝜗,𝜙), and a real number 𝑙(,), the Radon transform of the object is defined as an integral 𝑅Ψ𝑙,𝑛=𝑅Ψ(𝑙,𝜗,𝜙)=𝑟𝑛=𝑙Ψ𝑑𝑟2𝑟,(13) on the plane described by 𝑟𝑛=𝑙.

The Fourier transform can be calculated from the Radon transform by Ψ𝑘=Ψ=𝑘𝑛𝑅Ψ𝑙,𝑛exp(2𝜋𝑖𝑘𝑙)𝑑𝑙.(14) Thus, the inverse Fourier transform (12) becomes Radon's formulaΨ1𝑟=4𝜋2𝜋0𝜋/2𝜋/2𝜕2𝜕𝑙2||||𝑙=𝑟𝑛1𝑅Ψ𝑙,𝑛cos𝜗𝑑𝜗𝑑𝜙=4𝜋22𝜋0𝜋/2𝜋/2||𝑅Ψ𝑙,𝑛𝑙=𝑟𝑛1cos𝜗𝑑𝜗𝑑𝜙=8𝜋22Ω||𝑅Ψ𝑙,𝑛𝑙=𝑟𝑛𝑑𝑛,(15) where 2 is the Laplace operator, and Ω is the unit sphere.

Radon's formula tells us that to reconstruct the object at a point, what we need is the Radon transform of the planes through the point and its neighborhood. However, what we measure by the detector is line integrals along X-rays instead of the Radon transform [6] but it can be calculated from projections by =𝑅Ψ𝑙,𝑛=𝑅Ψ(𝑙,𝜗,𝜙)𝑃𝜙(𝑙cos𝜗𝑣sin𝜗,𝑙sin𝜗+𝑣cos𝜗)𝑑𝑣.(16) We call (16) the first relation between the Radon transform and parallel projection.

Formulas (15) and (16) form a scheme to reconstruct the object from its parallel projections. Here, the projection is not allowed to be longitudinally truncated. From Section 2.1, however, we have seen that the object can be reconstructed slice by slice from parallel projections. Why is the truncation not allowed when using the Radon formula?

For comparison, we rewrite Radon’s formula in the filtered backprojection format Ψ=𝑟𝜋0𝑃𝜙(𝜌,𝑧)𝑑𝜙.(17) with 𝑃𝜙(𝜌,𝑧)=𝜋/2𝜋/2Ψ𝑘𝑘exp2𝜋𝑖𝑘𝑟2=1cos𝜗𝑑𝑘𝑑𝜗4𝜋2𝜋/2𝜋/2𝜕2𝜕𝑙2||||𝑙=𝑟𝑛=𝜌cos𝜗+𝑧sin𝜗𝑅Ψ𝑙,𝑛cos𝜗𝑑𝜗.(18) Comparing (18) with (10), we have the following relation between the Radon transform and parallel-beam projections 14𝜋2𝜋/2𝜋/2𝜕2𝜕𝑙2||||𝑙=𝑟𝑛=𝜌cos𝜗+𝑧sin𝜗𝑅Ψ𝑙,𝑛(𝜗,𝜙)cos𝜗𝑑𝜗=𝑃𝜙1(𝜌,𝑧)2𝜋2𝜌2.(19) We call this the second relation between the Radon transform and parallel-beam projections. Based on the relationship between fan-beam and parallel-beam projections [1], the right hand side of (19) can be expressed in terms of cone-beam projections as well.

Substituting (16) into (19), we get an identity of parallel projections 14𝜋2𝜋/2𝜋/2𝜕2𝜕𝑙2||||𝑙=𝜌cos𝜗+𝑧sin𝜗×𝑃𝜙(𝑙cos𝜗𝑣sin𝜗,𝑙sin𝜗+𝑣cos𝜗)𝑑𝑣cos𝜗𝑑𝜗=𝑃𝜙(𝜌,𝑧)12𝜋2𝜌2.(20) Now, we see that reconstruction of an object from parallel-beam projections does allow truncation of projection data but to verify the identity (20) needs the information on the whole projection domain. Reconstruction from parallel-beam projections using the Radon formula means numerically verifying this identity. We point out that if the projection 𝑃𝜙(𝜌,𝑧) is replaced by a general function with two variables, this identity holds true as well; see Appendix A.

Similarly, in cone-beam CT, reconstruction from cone-beam projections using the Radon formula, that is, the Grangeat's framework, means numerically verifying the cone-beam version of this identity. This is the reason why the Grangeat framework cannot deal with the long object problem.

We see that the Radon formula comes from the inverse Fourier transform in the signed spherical coordinate system. In this simplest reconstruction problem, there are two shortcomings with the Radon formula: (1) what we measure with X-ray is line integrals rather than planar integrals, (2) the projection cannot be longitudinally truncated.

2.3. Standard Spherical Coordinate System

In a standard spherical coordinate system, the inverse Fourier transform is written as Ψ=𝑟02𝜋𝜋/2𝜋/20Ψ𝑘𝑘exp2𝜋𝑖𝑘𝑟2cos𝜗𝑑𝑘𝑑𝜗𝑑𝜙,(21) where (𝑘,𝜗,𝜙) with  0𝑘<,𝜋/2<𝜗<𝜋/2,0𝜙<2𝜋 is the standard spherical coordinates of 𝑘, as shown in Figure 2.

From the well-known generalized function relations [7] 𝛿(𝑙)=𝑖exp(2𝜋𝑖𝑘𝑙)𝑑𝑘,=𝜋𝑙sgn(𝑘)exp(2𝜋𝑖𝑘𝑙)𝑑𝑘,(22) we have 12𝑖𝛿(𝑙)+=2𝜋𝑙01exp(2𝜋𝑖𝑘𝑙)𝑑𝑘,2𝛿𝑖(𝑙)+𝜋𝑙3=4𝜋20𝑘2exp(2𝜋𝑖𝑘𝑙)𝑑𝑘.(23) Based on the convolution theorem, we have 0Ψ𝑘𝑘exp2𝜋𝑖𝑘𝑟2=𝑑𝑘0Ψ𝑘𝑛exp(2𝜋𝑖𝑘𝑙)𝑘21𝑑𝑘=8𝜋2𝜕2𝜕2𝑙𝑖𝑅Ψ𝑙,𝑛4𝜋3𝑙3.𝑅Ψ𝑙,𝑛(24) Therefore, the inverse Fourier transform becomes Ψ1𝑟=8𝜋202𝜋𝜋/2𝜋/2𝜕2𝜕𝑙2||||𝑅Ψ𝑙,𝑛𝑙=𝑟𝑛𝑖cos𝜗𝑑𝜗𝑑𝜙4𝜋302𝜋𝜋/2𝜋/21𝑅Ψ𝑙,𝑛𝑙3||||𝑙=𝑟𝑛cos𝜗𝑑𝜗𝑑𝜙.(25) Based on the odd and even symmetry, we have the Radon formula again [6] Ψ1𝑟=4𝜋2𝜋0𝜋/2𝜋/2𝜕2𝜕𝑙2||||𝑅Ψ𝑙,𝑛𝑙=𝑟𝑛cos𝜃𝑑𝜃𝑑𝜙.(26) The difference between the signed and standard spherical coordinate systems is little, and the signed spherical coordinate system leads to Radon's formula more easily.

It is important to underline that the inverse Fourier transform may yield various reconstruct formulas in different coordinate systems and allow different degrees of data truncation. Among the commonly used coordinate systems, the signed cylindrical coordinate system is the most convenient one to solve the simplest 3D reconstruction problem. To handle the general 3D reconstruction problem, in the following we will introduce a variant of the signed cylindrical coordinate system, which is referred to as the Γ coordinate system.

3. Γ Coordinate System and Its Jacobian Factor

What plays an important role in the 3D reconstruction field is a variant of the cylindrical coordinate system, which can be named the Γ coordinate system.

As shown in Figure 3(b),ΓAB is a three-times differentiable curve connecting points A(0,1,0) and B(0,1,0) on the unit sphere, whose length is 𝜃0𝜋.𝑆(𝜃) is a point on ΓAB parameterized by 𝜃[0,𝜃0] the length of the segment AS. Let us introduce three orthogonal unit vectors 𝑒3=𝑂𝑆,𝑒1=𝑑𝑑𝜃𝑒3,𝑒2=𝑒3×𝑒1,(27) where 𝑒1 represents the tangential direction, 𝑒2 the instantaneous rotation axis of the frequency plane. Note that they are functions of 𝜃=𝑆𝐴|𝑑𝑒3| (see Figures 3(a) and 3(b)).

fig3
Figure 3: Nonuniqueness of the representation of a point in the Γ coordinates. (a) The Γ coordinate system, (b) a general curve ΓAB on the unit sphere, and (c) the intersections of ΓAB and the great circle orthogonal to a given vector 𝑘. When 𝑆 moves from A to B ΓAB, some region may be scanned more than once.

In the Fourier domain, we call the plane through the origin 𝑂 and orthogonal to the vector 𝑂𝑆 the frequency plane, denoted as (𝑂𝑆). Every point (𝑘𝑂𝑆) can be expressed as 𝑘=𝜔1𝑒1+𝜔2𝑒2,(28) with 𝜔1=𝑘𝑒1,𝜔2=𝑘𝑒2.

More generally, for every given triplet (𝜔1,𝜔2,𝜃) with 𝜔1(,),𝜔2(,),𝜃[0,𝜃0], we can define a 3D vector 𝑘𝜔𝑘=1,𝜔2,𝜃=𝜔1𝑒1(𝜃)+𝜔2𝑒2(𝜃),(29) and the triplet (𝜔1,𝜔2,𝜃) is a set of Γ coordinates of the point 𝑘.

However, a given point 𝑘𝑅3 may be represented by several sets of Γ coordinates, because when 𝑆(𝜃) moves along ΓAB, the rotation axis 𝑒2 keeps changing, and hence the frequency plane (𝑂𝑆) may scan the given point 𝑘 more than once. The number of times of the point being scanned, denoted as 𝐽(𝑘), equals the number of intersections between the curve ΓAB and the great circle orthogonal to the vector 𝑘 (Figure 3(c)). To find the all intersections, one can solve the equation 𝑘𝑒3(𝜃)=0,(30) to obtain the solutions 1𝜃,2𝜃,,𝐽(𝑘)𝜃. For every solution, such as 1𝜃, one can further have 1𝜔1=𝑘𝑒1(1𝜃),1𝜔2=𝑘𝑒2(1𝜃). Thus, one gets the 𝐽(𝑘) groups of coordinates 𝑘𝑘=1𝜔1,1𝜔2,1𝜃=𝑘2𝜔1,2𝜔2,2𝜃𝑘==𝐽𝑘𝜔1,𝐽𝑘𝜔2,𝐽𝑘𝜃.(31) We call this new system the Γ coordinate system, since the system is based on the curve ΓAB. This coordinate system is a variant of the cylindrical system. For a differential arc 𝑆𝑆, the new system is similar to a cylindrical coordinate system if the scanned volume is considered. When Γ𝐴𝐵 is a great semicircle, the Γ coordinate system becomes the cylindrical coordinate system.

In the Γ coordinate system, the inverse Fourier transform can be expressed as Ψ=𝑟𝜃00𝑊𝜔1,𝜔2×Ψ𝑘||𝜔,𝜃exp2𝜋𝑖𝑘𝑟1||𝑑𝜔1𝑑𝜔2𝑑𝜃,(32) with 𝑘=𝜔1𝑒1(𝜃)+𝜔2𝑒2(𝜃). In [1], we obtained the absolute value of the Jacobian, |𝜔1| for the coordinate transform from (𝑘1,𝑘2,𝑘3) to (𝜔1,𝜔2,𝜃), by similarity between the Γ coordinate system and the cylindrical system.

Based on (29) in the Γ coordinate system, the absolute value of the Jacobian factor can be calculated as follows: ||||𝜕𝑘𝜕𝜔1×𝜕𝑘𝜕𝜔2𝜕𝑘||||=|||𝜕𝜃𝑒1×𝑒2𝜔1𝑑𝑑𝜃𝑒1+𝜔2𝑑𝑑𝜃𝑒2|||=|||𝑒3𝜔1𝑑𝑑𝜃𝑒1+𝜔2𝑑𝑑𝜃𝑒2|||=||𝜔1||,(33) where we have used the relations discussed in the next section𝑒3𝑑𝑑𝜃𝑒1=1,𝑒3𝑑𝑑𝜃𝑒2=0.(34) Then, a weight function 𝑊(𝜔1,𝜔2,𝜃) is defined for all the coordinates (𝜔1,𝜔2,𝜃), such that the summation of the weight function on the 𝐽(𝑘) groups of coordinates of the same 𝑘𝑅3 remains one, that is, 𝐽𝑘𝑗=1𝑊𝑗𝜔1,𝑗𝜔2,𝑗𝜃=1.(35) An example of the weight functions is [1, 8] 𝑊𝜔1,𝜔2𝑘,𝜃=sgn2𝜔sgn1,(36) with 𝑘2=𝑘𝑗=[𝜔1𝑒1(𝜃)+𝜔2𝑒2(𝜃)]𝑗.

With this weight function, (32) becomes a parallel-beam reconstruction formula from the Fourier slice theorem, and further Ye-Wang’s cone-beam formula [8] via the relationship between parallel-beam and cone-beam projections [1]. With the weight function, the filtration is no longer along the tangential direction generally.

Clearly, the weight function is not unique. For instance, for frequency points with 𝐽(𝑘)=3, in addition to the weight function defined by (36), (11+1), we can also choose (1/3+ 1/3+1/3), (0+1+0), or a more complicated one, (1+11) for some points and (1+1+1) for the rest points, leading to different parallel-beam reconstruction formulas. Furthermore, using the approach described in [1], we can arrive at Katsevich’s helical cone-beam formulas [911], and so forth.

If one takes the weight function 𝑊(𝜔1,𝜔2,𝜃)=1 on the right side of (32), what we will reconstruct is actually another function Φ=𝑟𝜃00Ψ𝑘||𝜔exp2𝜋𝑖𝑘𝑟1||𝑑𝜔1𝑑𝜔2=𝑑𝜃𝑅3Ψ𝑘𝐽𝑘𝑑exp2𝜋𝑖𝑘𝑟3𝑘.(37) From (37), we observe that the reconstructed image Φ(𝑟) equals the true object Ψ(𝑟) if and only if 𝐽(𝑘)=1, that is, the trajectory ΓAB is a great semicircle. The condition for Φ(𝑟) to be an approximation of the object function, that is, Φ(𝑟)Ψ(𝑟), is that 𝐽(𝑘)=1 in a major portion of the Fourier domain and 𝐽(𝑘) is close to 1 in the rest region. Based on the relation between 𝐽(𝑘) and the curve length of ΓAB, a necessary condition is 𝜃0𝜋, or equivalently 𝜀=𝜃0/𝜋10. Otherwise, Φ(𝑟) and Ψ(𝑟)will have a great difference. Unfortunately this formula was once thought as exact, which was essentially Theorem 3 in [2] and Theorem 4.3 in [3]. The details can be seen in [1].

Generally speaking, the criterion to identify exact formulas from approximate ones is to evaluate if every frequency component contributes exactly once to the final reconstruction, and the difference between exact formulas comes from different forms of weight functions.

4. Motion of the Frequency Plane and Condition for Sharp Points

In this section, we will study motion of the frequency plane, present the condition for sharp points of the locus of the rotation axis, and introduce the Euler angles. Let us simplify notations ̇𝑒3=(𝑑/𝑑𝜃)𝑒3,̈𝑒3=(𝑑2/𝑑𝜃2)𝑒3, and so on. In this section, we assume that Γ is a three times differentiable curve on the unit ball, which may be not complete.

Obviously we have the following.

Proposition 1. For any point 𝑆(𝜃)Γ, 𝑒1,𝑒2,𝑒3 defined by (27) are orthogonal unit vectors.

Definition 1. For every point 𝑆(𝜃)Γ, let us define 𝜎=𝑒3̇𝑒3×̈𝑒3,(38) which is the signed volume spanned by the three vectors 𝑒3,̇𝑒3,̈𝑒3. If we consider the curve length 𝜃 as time, 𝑒3,̇𝑒3,̈𝑒3 are the position, velocity, and the acceleration of point S.

Proposition 2. The motion of 𝑒1,𝑒2,𝑒3 can be described by the following differential equations 𝑑𝑑𝜃𝑒3=𝑒1,𝑑𝑑𝜃𝑒1=𝜎𝑒2𝑒3,𝑑𝑑𝜃𝑒2=𝜎𝑒1.(39)

Proof. The first equation is the definition. Now, let us prove the other two. Since 𝑒1,𝑒2,𝑒3 are orthogonal unit vectors, we can write 𝑒1𝑒1=1,𝑒2𝑒2=1,𝑒3𝑒3=1,𝑒1𝑒2=0,𝑒2𝑒3=0,𝑒3𝑒1=0.(40) Applying 𝑑/𝑑𝜃 on the above equations, one has ̇𝑒1𝑒1̇=0,𝑒2𝑒2̇=0,𝑒3𝑒3̇=0,𝑒1𝑒2+𝑒1̇𝑒2̇=0,𝑒2𝑒3=0,𝑒3̇e1=1.(41) Using the definition of 𝜎, we have ̇𝑒1𝑒2=̇𝑒1𝑒3×𝑒1=𝑒3̇𝑒3×̈𝑒3=𝜎,𝑒1̇𝑒2=𝜎.(42) Using the orthogonal unit vectors 𝑒1,𝑒2,𝑒3, ̇𝑒1 and ̇𝑒2 can be expressed as 𝑑𝑑𝜃𝑒1=𝜎𝑒2𝑒3,𝑑𝑑𝜃𝑒2=𝜎𝑒1,(43) which completes the proof.

From these differential equations, we recognize that the motion of 𝑒1,𝑒2,𝑒3 is a rotation in the 3D space.

Now, we consider the frequency plane Π(𝑂𝑆) fixed with 𝑒1,𝑒2, 𝑒3 as a rigid body. When 𝑆(𝜃) moves along the curve Γ, the rigid body rotates in the 3D space.

Let us define Ω=𝑒2+𝜎𝑒3.(44) We will see that Ω is the angular velocity of the rigid body by the following proposition.

Proposition 3. The motion of the vectors 𝑒1,𝑒2,𝑒3 can be expressed as 𝑑𝑑𝜃𝑒3=Ω×𝑒3,𝑑𝑑𝜃𝑒1=Ω×𝑒1,𝑑𝑑𝜃𝑒2=Ω×𝑒2.(45)

Proof. It is easy to verify the three equations one by one with Proposition 2 and the definition of Ω.

Proposition 4. The angular acceleration of the rigid body is ̇Ω=̇𝜎𝑒3(𝜃).(46)

Proof. In fact, we have ̇𝑑Ω=𝑑𝜃𝑒2+𝜎𝑒3=̇𝑒2+̇𝜎𝑒3̇+𝜎𝑒3=𝜎𝑒1+̇𝜎𝑒3+𝜎𝑒1=̇𝜎𝑒3.(47)

Proposition 5. For any point 𝑆(𝜃)Γ, we have ̈𝜎=0𝑒3=𝑒3,(48)̇𝜎=0𝑒3𝜎=2̇+1𝑒3.(49)

Proof. Let us first prove (48).
If 𝜎=0, by Proposition 2 we have ̈𝑒3=𝑒3. On the other hand, if ̈𝑒3=𝑒3, we have 𝜎=𝑒3̇[𝑒3×̈𝑒3]=0.
Now, let us prove (49).
By Proposition 2, we have ̇𝑒1=𝜎𝑒2𝑒3.(50) Applying 𝑑/𝑑𝜃 to both sides of (50) and using ̇𝜎=0, we have ̈𝑒1=̇𝜎𝑒2̇+𝜎𝑒2̇𝑒3=𝜎𝜎𝑒1𝑒1𝜎=2+1𝑒1.(51) On the other hand, if 𝑒3=(𝜎2̇+1)𝑒3, we have 𝑑̇𝜎=𝑑𝜃𝑒3̇𝑒3×̈𝑒3=̇𝑒3̇𝑒3×̈𝑒3+𝑒3̈𝑒3×̈𝑒3+𝑒3̇𝑒3×𝑒3=𝑒3̇𝑒3×𝑒3=0.(52) This completes the proof. Physically, 𝑒3 is the change rate of the acceleration ̈𝑒3.

Proposition 6. 𝜎=0 for every point 𝑆(𝜃)Γ, if and only if Γ is a great circle or a great arc.

Proof. Since 𝜎=0 for every point 𝑆(𝜃)Γ, the motion equation of the three unit vectors, (39), is simplified as 𝑑𝑑𝜃𝑒3=𝑒1,𝑑𝑑𝜃𝑒1=𝑒3,𝑑𝑑𝜃𝑒2=0.(53) Since (𝑑/𝑑𝜃)𝑒2=0 for every point, 𝑒2 is a fixed unit vector. Since (𝑑/𝑑𝜃)𝑒3=𝑒1=𝑒2×𝑒3, 𝑒3 rotates about the fixed axis 𝑒2. Since 𝑒3𝑒2, 𝑆, the endpoint of 𝑒3, draws a great circle or a great arc; see Figure 4(a).
On the other hand, if Γ is a great circle or a great arc, one can verify that ̈𝑒3=𝑒3. Therefore, we have 𝜎=𝑒3̇(𝑒3×̈𝑒3)=0. This completes the proof.

fig4
Figure 4: The circles on unit sphere and a helix in 3D space. (a) The illustration of 𝜎=0 anḋ𝜎=0, (b) a circular locus on the unit sphere and the motion of the frequency plane, and (c) the helix in 3D space and the motion of its TNB frame. The figure is not drawn to scale.

Proposition 7. ̇𝜎=0 for every point 𝑆(𝜃)Γ, if and only if Γ is a circle or a circular arc.

Proof. Since ̇𝜎=0, by Proposition 4 we have ̇Ω=̇𝜎𝑒3(𝜃)=0. We know that Ω is a fixed vector. By Proposition 3, we have 𝑑𝑑𝜃𝑒3=Ω×𝑒3.(54) Therefore, 𝑒3 rotates about the fixed axis Ω, and 𝑆, the endpoint of 𝑒3, draws a circular arc. Furthermore, the radius of the circle is 1/1+𝜎2; see Figure 4(b).
On the other hand, if Γ is a circular arc, the line connecting the circle center and the origin 𝑂 is the rotation axis. We know that ̈𝑒3 is pointing to the axis. By circular symmetry, 𝜎=𝑒3̇(𝑒3×̈𝑒3) is constant at every point of the circular arc. Therefore, ̇𝜎=0 on the circular arc. This completes the proof.

Clearly, Proposition 6 is a special case of Proposition 7; see Figure 4(a).

One important case is that 𝜎=𝑒3̇(𝑒3×̈𝑒3)=0, for a certain point 𝑆(𝜃)Γ. At such a moment, we have (𝑑/𝑑𝜃)𝑒2=𝜎𝑒1=0, the rotation axis of the frequency plane does not move, the point 𝑆 moves along a great arc, and the Γ coordinate system is more similar to a cylindrical coordinate system. We say that 𝜃 is a stationary point of 𝑒2 if 𝜎(𝜃)=0. If 𝜎 has different signs before and after 𝜎=0, the stationary point will become a sharp point, since in this case 𝑒2 goes forward, stops and goes back according to (𝑑/𝑑𝜃)𝑒2=𝜎𝑒1. Therefore, we have the following straightforward proposition.

Proposition 8. A necessary condition for the curve 𝑒2(𝜃) to have a sharp point at 𝜃1(0,𝜃0) is 𝜎(𝜃1)=0; a sufficient condition for the curve 𝑒2(𝜃) to have a sharp point at 𝜃1(0,𝜃0) is that 𝜎(𝜃) has different signs for 𝜃<𝜃1 and for 𝜃>𝜃1.

Based on the definition of 𝜎 and the equation of the osculating plane, 𝜎(𝜃1)=0 means that the osculating plane at 𝑆(𝜃1)Γ passes through the origin since (𝑒3̇0)(𝑒3×̈𝑒3)=0.𝜎(𝜃) changes its sign before and after 𝜃1 means that the osculating plane sweeps the origin, or relatively speaking, the origin goes from one side of the osculating plane to the other. An example can be seen in the next section. Sharp points serve as the landmarks for defining weight functions.

Proposition 9. For any given point on the frequency plane 𝑘=𝜔1𝑒1+𝜔2𝑒2, the motion equation is 𝑑𝑘𝑑𝜃𝑘=Ω×=𝜎𝜔2𝑒1(𝜃)+𝜎𝜔1𝑒2(𝜃)𝜔1𝑒3(𝜃).(55)

Proof. It is straightforward from Propositions 2 or 3.
We call (55) the motion equation of the frequency plane.
There are two applications of Proposition 9. First, it can be used to compute the Jacobian of the Γ coordinate system ||||𝜕𝑘𝜕𝜔1×𝜕𝑘𝜕𝜔2𝜕𝑘||||=||||𝜕𝜃𝑒3𝜕𝑘||||=||𝜕𝜃𝜔1||=||𝜔1||.(56) Second, since 𝑒1𝜕𝜕𝜃𝑘=𝜎𝜔2,𝑒2𝜕𝜕𝜃𝑘=𝜎𝜔1,(57) we know that among the three vectors 𝜕𝑘/𝜕𝜔1,𝜕𝑘/𝜕𝜔2,𝜕𝑘/𝜕𝜃, the last one is not orthogonal to the first two. Hence, a general Γ coordinate system is not orthogonal. Notice that 𝜕𝑘/𝜕𝜃=𝜕𝑘/𝜕𝜃|𝜔1,𝜔2=𝑑𝑘/𝑑𝜃
Note that the frequency plane rotates about the axis 𝑒2 with the unit angular velocity while the vectors 𝑒1 and 𝑒2 themselves rotate about the axis 𝑒3 with the angular velocity 𝜎; see Figure 5(a). We see that the Γ coordinate system is different from the cylindrical coordinate system in the motion process. However, since the spin about the axis 𝑒3 of the frequency plane does not have an effect on the scanned volume, their Jacobians have the same expression. Figure 5(b) shows the shape of a differential volume of the Γ coordinate system.

fig5
Figure 5: Differential motion of the frequency plane. (a) Two decomposed rotations and (b) one combined rotation.

Proposition 10. When S moves along a curve Γ on the unit sphere, the scanned volume in the unit ball by the frequency plane is 𝑉=4𝜃03=||𝑘||<1𝐽𝑘𝑑3𝑘,(58) where 𝜃0 is the length of Γ, and 𝐽(𝑘) is the number of times point 𝑘 is scanned.

Proof. In fact, the scanned volume can be calculated in both the Γ and the Cartesian coordinate systems as follows: 𝑉=𝜃00𝑘2=𝜔21+𝜔22<1||𝜔1||𝑑𝜔1𝑑𝜔2=𝑑𝜃𝜃001𝑘=02𝜋𝜙=0𝑘2||||=4cos𝜙𝑑𝑘𝑑𝜙𝑑𝜃3𝜃004𝑑𝜃=3𝜃0=𝑘21+𝑘22+𝑘23<1𝐽𝑘𝑑3𝑘.(59) This completes the proof.

This concise formula provides a basic relation between the curve length, the scanned volume by the frequency plane, and the J function.

When we study the rotation problem, Euler angles (𝜗,𝜙,𝛿) are helpful. As shown in Figure 6(a), (𝜗,𝜙) with 𝜗(𝜋/2,𝜋/2),𝜙(,) are the spherical coordinate of 𝑆=𝑆(𝜗,𝜙),𝛿 is the angle between 𝑒𝜗 and 𝑒2. 𝑒𝜙, and 𝑒𝜗 are the latitude and longitude unit vectors at point 𝑆 in the spherical coordinate system.

fig6
Figure 6: Euler angles. (a) The Euler angles of the frequency plane and (b) the rotation angle of frequency plane relative to the detector plane.

The Euler angles are functions of the curve length 𝜗=𝜗(𝜃),𝜙=𝜙(𝜃),𝛿=𝛿(𝜃).(60) The three unit vectors can be represented as follows: 𝑒3=𝑒3(𝜗,𝜙),𝑒1=sin𝛿𝑒𝜗+cos𝛿𝑒𝜙=̇𝜗𝑒𝜗+̇𝜙cos𝜗𝑒𝜙,𝑒2=cos𝛿𝑒𝜗+sin𝛿𝑒𝜙̇=𝜙cos𝜗𝑒𝜗+̇𝜗𝑒𝜙.(61) Furthermore, one can verify the motion equations of the three unit vectors in Proposition 2.

Based on the motion theory of rigid body, the angular velocity of a rigid body can be written as ̇Ω=𝜗𝑒𝜙+̇𝜙̇𝑘+𝛿𝑒3̇=𝜗𝑒𝜙+̇𝜙cos𝜗𝑒𝜗+̇̇𝛿𝜙sin𝜗+𝑒3=𝑒2+𝜎𝑒3.(62) Therefore, we have the following relationships:𝑒2̇=𝜗𝑒𝜙+̇𝜙cos𝜗𝑒𝜗=Ω0Ω0𝑒3𝑒3,̇̇Ω𝜎=𝜙sin𝜗+𝛿=0𝑒3̇Ω0𝑒1,̇̇Ω𝛿=0𝑒1=Ω0̇𝑒1,(63) where we have used a new notation Ω0̇=𝜗𝑒𝜙+̇𝜙𝑘,(64) which is the first two terms of the total angular velocity Ω of (62).

In fact, Ω0 is the angular velocity of the detector with the unit vectors 𝑒𝜙,𝑒𝜗 when 𝑆 moves along Γ. To be specific, we consider the detector for parallel-beam projection as a rectangle of finite width and length through the origin orthogonal to direction 𝑒3=𝑂𝑆. During any motion, one side of the detector always keeps horizontal, that is, parallel to the direction 𝑒𝜙. The pose of the detector is completely determined by the point 𝑆=𝑆(𝜗,𝜙), or two Euler's angles (𝜗,𝜙). The relative motion of the frequency plane to the detector is described by the angle 𝛿; see Figure 6(b).

In summary, there are three sets of orthogonal vectors (beside ki,j,, the unit vectors of the Cartesian coordinate system) and three important planes. The orthogonal unit vectors 𝑒𝜙,𝑒𝜗, 𝑒3 are related to the movement of the detector and 𝑒3, ̇𝑒3=𝑒1, 𝑒2=𝑒3×𝑒1 are related to the movement of the frequency plane, see Figure 6(b). It is easy to verify that Ω=𝑒2+𝜎𝑒3=𝑒1×̇𝑒1 and hence the vectors 𝑒1, ̇𝑒1, Ω are the third set of orthogonal vectors. The direction of Ω is the normal of the osculating plane and its magnitude, 1+𝜎2, is the reciprocal of the radius of the osculating circle. Here the osculating plane is spanned by 𝑒1, ̇𝑒1 [15] and the osculating circle is referred to as the intersection between the osculating plane and the unit ball. When S moves along a curve Γ on the unit ball, it can be viewed that the point S moves along an osculating circle, whose orientation and radius are described by Ω. From Ω=𝑒2+𝜎𝑒3, we can see that 𝜎 is an angular velocity of the spin of the frequency plane. The locus lengths of the endpoints of 𝑒3,𝑒1, 𝑒2 are 𝜃0, 𝜃001+𝜎2𝑑𝜃, 𝜃00|𝜎|𝑑𝜃, respectively, according to Proposition 2.

5. Example of a Curve with 𝐽(𝑘)=1,3

Here, we discuss the reconstruction based on the simplified model of the 3D helical cone beam reconstruction, that is, 3D helical parallel-beam reconstruction.

The curve 𝐶 in 𝑅3 defined by 𝑆𝐶=𝑅3𝑂𝑆=𝑦(𝜑)=𝑟cos𝜑,𝑟sin𝜑,𝜑,𝜋2𝜋𝜑2,𝜋2(65) is a segment of a helix with radius 𝑟 and pitch , which is shown in Figure 7 for 𝑟=1,=3.

285130.fig.007
Figure 7: Motion process of the frequency plane. The red curve is a segment of a helical scanning trajectory. The blue curve is its projection on the unit sphere. The black curve is the great circle through 𝑆1,𝑆3, and 𝑆5. Five positions of the frequency plane are color coded. The points inside the triangle are scanned 3 times while the rest points are scanned only once.

In fact, it is formed by a combination of a uniform circular motion in the XOY plane and a uniform straight motion along the 𝑍 direction.

The projection of the curve 𝐶 on the unit sphere Ω is Γ=𝑆Ω𝑂𝑆=𝑒3(𝜑)=𝑦(𝜑)||||𝜋𝑦(𝜑),𝜑2,𝜋2,(66) or Γ=𝑆Ω𝑂𝑆=𝑒3(𝜗,𝜙),𝜗=tg1𝜙𝜋2𝜋,𝜙2,𝜋2.(67) The curve length 𝜃 can be given by ||𝑑𝜃=𝑑𝑒3||=||||𝑑(𝜑)𝑦(𝜑)||||||||,𝑦(𝜑)𝜃(𝜑)=𝜑𝜋/2𝑑𝜃𝑑𝜑𝑑𝜑=𝜑𝜋/2|||||𝑑𝜑𝑦||𝑦(𝜑)|||||||.(68) For example, we can calculate 𝜃0=𝜃(𝜋/2)3.177 and 𝜀=0.011, when 𝑟=1 and =3, indicating that the overlapping degree in the Fourier domain is small.

For every 𝜑[𝜋/2,𝜋/2], the three orthogonal unit vectors are𝑒3=𝑂𝑆=𝑦(𝜑)||||,𝑦(𝜑)𝑒1=(𝑑/𝑑𝜑)𝑒3||(𝑑/𝑑𝜑)𝑒3||,𝑒2=𝑒3×𝑒1.(69) The inverse Fourier transform can be expressed with either the Γ coordinates (𝜔1,𝜔2,𝜃) or the new coordinates (𝜔1,𝜔2,𝜑)Ψ=𝑟𝜃00𝑊𝜔1,𝜔2Ψ𝑘×||𝜔,𝜃exp2𝜋𝑖𝑘𝑟1||𝑑𝜔1𝑑𝜔2=𝑑𝜃𝜋/2𝜋/2𝑊𝜔1,𝜔2×Ψ𝑘||𝜔,𝜃(𝜑)exp2𝜋𝑖𝑘𝑟1||𝑑𝜃𝑑𝜑𝑑𝜔1𝑑𝜔2𝑑𝜑,(70) with 𝑘𝜔𝑘=1,𝜔2,𝜃(𝜑)=𝜔1𝑒1(𝜃(𝜑))+𝜔2𝑒2(𝜃(𝜑)).(71) To visualize the motion of the frequency plane (𝑂𝑆), its five positions are marked in Figure 7 at 𝜑1=𝜋/2, 𝜑2=𝜋/4, 𝜑3=0,𝜑4=𝜋/4 and 𝜑5=𝜋/2, respectively. Note that the initial position Π(𝑂𝑆1) and the final position Π(𝑂𝑆5) coincide with each other.

In Figure 7, the two curved sides of the red triangle on the top of the unit sphere is the locus of 𝑒3, the instantaneous rotation axis of the frequency plane. Similarly, 𝑒2 and Π(𝑂𝑆1) form the counterpart triangle on the bottom of the unit sphere. The locus of 𝑒2 is the division line of the regions with different 𝐽(𝑘) values. Based on (58), the total area of the two triangles is only 0.55% of the surface area of the unit sphere.

Note that a sharp point on the locus of 𝑒2 appears at 𝜑=𝜑3=0. The reason is as follows.

Let us define𝜎0𝑑(𝜑)=𝑦(𝜑)𝑑𝑑𝜑𝑦(𝜑)×2𝑑𝜑2𝑦(𝜑).(72) It can be verified that 𝜎0(0)=0, 𝜎0(𝜑)<0 for 𝜑<0, and 𝜎0(𝜑)>0 for 𝜑>0.

It is not difficult to prove 𝜎=𝑒3𝑒1×̇𝑒1𝑑=𝑦(𝜑)𝑑𝑑𝜑𝑦(𝜑)×2𝑑𝜑21𝑦(𝜑)||||×1𝑦(𝜑)||(||𝑑/𝑑𝜑)𝑦(𝜑)2𝑑𝜑𝑑𝜃3=𝜎0(𝜑)𝐾(𝜑),(73) where 𝐾(𝜑) is positive. Therefore, 𝜎 and 𝜎0 have the same zeros and sign.

Based on the motion process of Π(𝑂𝑆), we have 𝐽𝑘=1,𝑘outsidethetwotriangles,3,𝑘insidethetwotriangles.(74) The overlapped region in the frequency space is due to the motion of the instantaneous rotation axis 𝑒2 of the frequency plane, or equivalently the discrepancy between the curve Γ and half a great circle.

As shown in Figure 8, the weight factor can be (11+1), (1/3+ 1/3+1/3), or (0+1+0), and so forth, for the overlapped region. By dividing the overlapped region into 2 symmetric parts, one can define the weight function as (1+1+1) on the left side and (1+11) on the right side, or (11+1) on the left and (1+11) on the right, and so forth. Hence, one can design various parallel-beam reconstruction formulas and derive the associated cone-beam formulas using the approach described in [1].

fig8
Figure 8: Examples of possible weight factors. (a) 11+1, (b) 1/3+1/3+1/3, (c) 0+1+0, and (d) 1+1+1 for the left half, 1+11 for the right half, (e) 11+1 for the left half, and 1+11 for the right half. The straight line represents a great circle on the unit sphere. The weight factor in (a) for (𝑂𝑆2) is further illustrated in (f) as an example.

Appendices

A. Weighted Λ Reconstruction Formula

𝑓(𝑥,𝑦) is a function in the 2D space with the Fourier transform 𝑓𝑘1,𝑘2=𝑓(𝑥,𝑦)exp𝑖2𝜋𝑘1𝑥exp𝑖2𝜋𝑘2𝑦𝑑𝑥𝑑𝑦.(A.1) Its parallel-beam projection is denoted as 𝑃𝜗(𝑙)=𝑓(𝑙cos𝜗𝑣sin𝜗,𝑙sin𝜗+𝑣cos𝜗)𝑑𝑣(A.2) with 𝑙(,),𝜗[𝜋/2,𝜋/2].

Let 𝑓(𝑥,𝑦) pass through the filter |𝑘1|, we get a new function =Φ(𝑥,𝑦)||𝑘1||𝑓𝑘1,𝑘2exp𝑖2𝜋𝑘1𝑥exp𝑖2𝜋𝑘2𝑦𝑑𝑘1𝑑𝑘21=𝑓(𝑥,𝑦)2𝜋2𝑥2.(A.3) Defining the signed polar coordinates (𝑘𝜌,𝜗) with 𝑘𝜌(,),𝜗(𝜋/2,𝜋/2), which satisfy 𝑘1=𝑘𝜌cos𝜗,𝑘2=𝑘𝜌sin𝜗, we haveΦ(𝑥,𝑦)=𝜋/2𝜋/2||𝑘1||𝑓𝑘𝜌,𝜗×exp𝑖2𝜋𝑘𝜌||𝑘(𝑥cos𝜗+𝑦sin𝜗)𝜌||𝑑𝑘𝜌=𝑑𝜗𝜋/2𝜋/2𝑓𝑘𝜌,𝜗exp𝑖2𝜋𝑘𝜌(𝑥cos𝜗+𝑦sin𝜗)×𝑘2𝜌cos𝜗𝑑𝑘𝜌=1𝑑𝜗4𝜋2𝜋/2𝜋/2𝜕2𝜕𝑙2||||𝑙=𝑥cos𝜗+𝑦sin𝜗𝑓𝑘𝜌,𝜗×exp𝑖2𝜋𝑘𝜌𝑙𝑑𝑘𝜌=1cos𝜗𝑑𝜗4𝜋2𝜋/2𝜋/2𝜕2𝜕𝑙2||||𝑙=𝑥cos𝜗+𝑦sin𝜗𝑃𝜗(𝑙)cos𝜗𝑑𝜗.(A.4) Comparing (A.3) and (A.4), we have1𝑓(𝑥,𝑦)2𝜋2𝑥2=14𝜋2𝜋/2𝜋/2𝜕2𝜕𝑙2||||𝑙=𝑥cos𝜗+𝑦sin𝜗𝑃𝜗(𝑙)cos𝜗𝑑𝜗.(A.5) That is, we have rediscovered the identity (20) in a general way. In contrast to the classical Λ reconstruction formula [1], a weight factor cos𝜗 is involved in (A.5). Hence, we can call (A.5) a weighted Λ reconstruction formula.

B. A Weighted Radon Formula and the Orlov Formula

Let us define 𝑢=𝑟𝑒3,𝑡=𝑟𝑒1,𝑠=𝑟𝑒2,(B.1) and the parallel-beam projection is 𝑃𝜃(𝑡,𝑠)=Ψ𝑢𝑒3+𝑡𝑒1+𝑠𝑒2𝑑𝑢,(B.2) as shown in Figure 6(b). In the frequency plane Π(𝑂𝑆), we define a signed polar coordinate system (𝜔,𝛼) with 𝜔(,),𝛼[0,𝜋] for a point (𝜔1,𝜔2) with 𝜔1=𝜔cos𝛼,𝜔2=𝜔cos𝛼. Similarly, let us define the signed polar coordinates (𝑋,𝛽) with 𝑋(,),𝛽[𝜋/2,𝜋/2] for a point (𝑡,𝑠) with 𝑡=𝑋cos𝛽,𝑠=𝑋sin𝛽.

Let the weight function in (32) be 𝑊(𝜔1,𝜔2,𝜃)=𝑊(𝛼,𝜃), we have Ψ=𝑟𝜃00𝑊𝜔1,𝜔2Ψ𝑘×||𝜔,𝛼exp2𝜋𝑖𝑘𝑟1||𝑑𝜔1𝑑𝜔2=𝑑𝜃𝜃00𝜋0𝑊(𝛼,𝜃)|cos𝛼|exp2𝜋𝑖𝜔1𝑡+2𝜋𝑖𝜔2𝑠×Ψ𝜔1,𝜔2𝜔,𝜃2=𝑑𝜔𝑑𝛼𝑑𝜃𝜃00𝜋0𝑊(𝛼,𝜃)|cos𝛼|𝜔2×=exp(2𝜋𝑖𝜔(𝑡cos𝛼+𝑠sin𝛼))Ψ(𝜔,𝛼,𝜃)𝑑𝜔𝑑𝛼𝑑𝜃14𝜋2𝜃00𝜋0𝜕𝑊(𝛼,𝜃)2𝜕𝑙2||||𝑙=𝑟𝑛=𝑡cos𝛼+𝑠sin𝛼=𝑅Ψ𝑙,𝑛×|cos𝛼|𝑑𝛼𝑑𝜃14𝜋22𝜃00𝜋0×𝑊(𝛼,𝜃)𝑅Ψ𝑡cos𝛼+𝑠sin𝛼,𝑛|cos𝛼|𝑑𝛼𝑑𝜃,(B.3) with 𝑛=𝑒1cos𝛼+𝑒2sin𝛼.

We call it the weighted Radon formula. When ΓAB is a great semicircle, we have 𝑊(𝛼,𝜃)=1 and 𝜃0=𝜋, and the above formula becomes the classic Radon formula [12, 13].

Moreover, it is easy to rederive Orlov's formula [14] from (32) Ψ=𝑟𝜃00𝑊𝜔1,𝜔2Ψ𝑘×||𝜔,𝜃exp2𝜋𝑖𝑘𝑟1||𝑑𝜔1𝑑𝜔2=𝑑𝜃𝜃00𝑊(𝛼,𝜃)exp2𝜋𝑖𝜔1𝑡+2𝜋𝑖𝜔2𝑠×Ψ𝜔1,𝜔2||𝜔,𝜃1||𝑑𝜔1𝑑𝜔2=𝑑𝜃𝜃0014𝜋2𝜕2𝜕𝑡2+𝜕2𝜕𝑠2𝑊(𝛼,𝜃)×exp2𝜋𝑖𝜔1𝑡+2𝜋𝑖𝜔2𝑠Ψ𝜔1,𝜔2||𝜔,𝜃1||𝜔2𝑑𝜔1𝑑𝜔21𝑑𝜃=4𝜋22𝜃000𝑥0200𝑑|cos𝛼||Ψ𝜔𝜔|𝑊(𝛼,𝜃)1,𝜔2,𝜃×exp2𝜋𝑖𝜔1𝑡+2𝜋𝑖𝜔2𝑠𝑑𝜔1𝑑𝜔21𝑑𝜃=4𝜋22𝜃00𝑝𝜃(𝑡,𝑠)𝑔𝜃(𝑡,𝑠)𝑑𝜃,(B.4) where 𝑔𝜃=(𝑡,𝑠)|cos𝛼||𝜔|𝑊(𝛼,𝜃)exp2𝜋𝑖𝜔1𝑡+2𝜋𝑖𝜔2𝑠𝑑𝜔1𝑑𝜔2=|cos𝛼|||𝑋||||𝑊(𝛼,𝜃)𝛼=𝜋/2+𝛽.(B.5) The sign denotes the 2D convolution operation.

We recall that a frequency point has 𝐽(k) sets of Γ coordinates 𝑘𝑘=1𝜔1,1𝜔2,1𝜃=𝑘2𝜔1,2𝜔2,2𝜃𝑘==𝐽𝑘𝜔1,𝐽𝑘𝜔2,𝐽𝑘𝜃.(B.6) If a given triplet (𝜔1,𝜔2,𝜃) is the 𝑗th set of Γ coordinates of 𝑘, 𝜔1,𝜔2=,𝜃𝑗𝜔1,𝑗𝜔2,𝑗𝜃,(B.7) a new normalized weight function [1] for the triplet can be defined by 𝑊𝜔1,𝜔2||,𝜃=𝑊(𝛼,𝜃)=1/𝑘𝑒1𝑗𝜃||𝐽(𝑘)𝑗=1||1/𝑘𝑒1𝑗𝜃||=1/|cos𝛼|𝑘𝐽(0)𝑗=1||𝑘1/0𝑒1𝑗𝜃||,(B.8) and in this case we have 𝑔𝜃(1𝑡,𝑠)=||𝑋||1𝑘𝐽(0)𝑗=1||𝑘1/0𝑒1𝑗𝜃||||||||𝑘0=𝑒3×𝑋0.(B.9) Here, 𝑘0=𝑘/|𝑘|, and 𝑋0=(𝑡𝑒1+𝑠𝑒2)/𝑡2+𝑠2=𝑒1cos𝛽+𝑒2sin𝛽 are unit vectors.

Equation (B.4) with (B.9) is the Orlov formula for a general curve on the unit sphere ((2a) with (2c) in [14]). Note that in [14] different letters were used such as 𝑘𝑚0,𝜏𝑒3,̇𝜏𝑒1,𝑥(𝑡𝑒1+𝑠𝑒2),|𝑥||𝑋|,𝐺𝑑𝐺𝜏𝜃00𝑑𝜃, and so forth, and ̇𝜏𝑘 should be printed as ̇𝜏𝑘.

C. Link between the Motion of the Frequency Plane and Differential Geometry

To study the property of a curve Γ on the unit sphere and the associated motion of the frequency plane, one can introduce another curve, Γ, defined as 𝑟(𝜃)=𝜃0𝑒3𝜃𝑑𝜃+𝑟(0),(C.1) where 𝑟(0) is an arbitrary vector in the 3D space, and 𝜃 (as well as 𝜃) is an arc-length parameter of the curve Γ. Clearly, ̇𝑑𝑟=𝑑𝜃𝑟(𝜃)=𝑒3,̈𝑟=𝑒1,̇𝑟=𝑒1.(C.2) Suppose that the arc-length parameter of the curve Γ is 𝑠. The tangential direction of the curve Γ is 𝑑T=𝑑𝑠𝑟=𝑑𝜃𝑑𝑑𝑠𝑑𝜃𝑟(𝜃)=𝑑𝜃𝑑𝑠𝑒3(𝜃).(C.3) Furthermore, we have |||𝑑|||=|||𝑑𝑠𝑟𝑑𝜃|||||𝑑𝑠𝑒3||,|||(𝜃)1=𝑑𝜃|||,𝑑𝑠(C.4) where |(𝑑/𝑑𝑠)𝑟|=|𝑑𝑟/𝑑𝑠|=|𝑑𝑟|/𝑑𝑠=𝑑𝑠/𝑑𝑠=1 has been used.

When the positive directions are the same, we have 𝑑𝜃𝑑𝑠=1,𝜃=𝑠+𝑠0,T=𝑒3(𝜃).(C.5) Picking up the constant 𝑠0=0, we have 𝜃=𝑠. Therefore, the symbol 𝜃 represents the arc length parameters of both Γ and Γ. Based on the tangential direction T, the norm and binormal unit vectors [15] can be defined as ̇TN=|||̇T|||=𝑒1||𝑒1||=𝑒1,B=T×N=𝑒3×𝑒1=𝑒2.(C.6) The curvature 𝜅 and the torsion 𝜏 of the curve Γ are ||̈||=||𝜅=𝑟𝑒1||̇̈=1,𝜏=𝑟𝑟×𝑟=𝑒3𝑒1×̇𝑒1=𝜎,(C.7) and the Darboux vector is Ω=𝜏T+𝜅B=𝑒2+𝜎𝑒3.(C.8) The Frenet-Serret Equations of the curve Γ̇̇̇T=N,N=T+𝜎B,B=𝜎N,(C.9) and ̇̇̇T=Ω×T,N=Ω×N,B=Ω×B,(C.10)

are (39) and (45). In this frame, 𝑟(𝜃) is the original curve (function), and 𝑒3(𝜃) becomes its derived curve (function). In practice, it is difficult to find the physical meaning of the 3D curve 𝑟(𝜃).

Example C. A circle on the unit sphere and its corresponding curve in the 3D space.

As shown in Figure 4(b), a circle on the unit sphere, marked as 𝐶, can be described as 𝑒3=sin𝜗0k+cos𝜗0cos𝜙i+cos𝜗0jsin𝜙(C.11) with a fixed parameter 𝜗0(𝜋/2,+𝜋/2) and a variable 𝜙[0,2𝜋].

Using the arc-length parameter 𝜃[0,2𝜋cos𝜗0], we have 𝑒3(𝜃)=sin𝜗0k+cos𝜗0𝜃coscos𝜗0i+cos𝜗0𝜃sincos𝜗0j.(C.12) Furthermore, we have 𝑒1𝜃(𝜃)=sincos𝜗0𝜃i+coscos𝜗0̇j,𝑒11(𝜃)=cos𝜗0𝜃coscos𝜗01icos𝜗0𝜃sincos𝜗0j,𝜎=𝑒3𝑒1×̇𝑒1=tg𝜗0,Ω=𝑒2+𝜎𝑒3=1cos𝜗0k.(C.13) Based on the definition (C.1), the corresponding 3D curve 𝐶 is 𝑟(𝜃)=𝜃sin𝜗0k+cos2𝜗0𝜃sincos𝜗0icos2𝜗0𝜃coscos𝜗0j,(C.14) with the choice 𝑟(0)=𝑗cos2𝜗0.

Now, we see that 𝐶 is a helix, whose standard form is 𝜋𝑟(𝜙)=𝑏𝜙k+𝑎cos𝜙2𝜋i+𝑎sin𝜙2j,(C.15) with 𝑎=cos2𝜗0,𝑏=sin𝜗0cos𝜗0,𝜙[0,2𝜋] (see Figure 4(c)). Its curvature and torsion are 𝜅=1,𝜏=tg𝜗0.(C.16) The length of the helix 𝐶 for 𝜙[0,2𝜋] is 2𝜋𝑎2+𝑏2=2𝜋cos𝜗0, which is equal to the circumference of circle 𝐶 on the unit sphere. The radius of the helix 𝐶(𝑎=cos2𝜗0) is smaller than that of the circle 𝐶(cos𝜗0). The figure is not drawn to scale. The TNB frame rotates around the 𝑧 direction at a constant velocity Ω=2𝜋2𝜋cos𝜗01k=cos𝜗0k.(C.17)

Acknowledgment

The works of H. Yu and G. Wang are partially supported by the NSF/MRI Program (no. CMMI-0923297) and NIH/NIBIB Grant (no. EB011785). The authors also thank a reviewer for a valuable suggestion.

References

  1. Y. Wei, H. Yu, J. Hsieh, W. Cong, and G. Wang, “Scheme of computed tomography,” Journal of X-Ray Science and Technology, vol. 15, no. 4, pp. 235–270, 2007. View at Google Scholar · View at Scopus
  2. V. P. Palamodov, “Reconstruction from ray integrals with sources on a curve,” Inverse Problems, vol. 20, no. 1, pp. 239–242, 2004. View at Google Scholar · View at Scopus
  3. V. P. Palamodov, Reconstructive Integral Geometry, Birhauser, Boston, Mass, USA, 2004.
  4. H. Yu, Y. Ye, S. Zhao, and G. Wang, “Studies on Palamodov's algorithm for cone-beam CT along a general curve,” Inverse Problems, vol. 22, no. 2, pp. 447–460, 2006. View at Publisher · View at Google Scholar · View at Scopus
  5. H. Yu, Y. Ye, S. Zhao, and G. Wang, “Reply to the comment on 'studies on Palamodov's algorithm for cone-beam CT along a general curve',” Inverse Problems, vol. 22, no. 4, pp. 1505–1506, 2006. View at Publisher · View at Google Scholar · View at Scopus
  6. A. G. Ramm and A. I. Katsevich, The Radon Transform and Local Tomography, CRC Press, New York, NY, USA, 1996.
  7. I. M. Gel'fand and G. E. Shilov, Generalized functions: Volume 1, vol. 1, Academic Press, New York, NY, USA, 1977.
  8. Y. Yangbo and G. Wang, “Filtered backprojection formula for exact image reconstruction from cone-beam data along a general scanning curve,” Medical Physics, vol. 32, no. 1, pp. 42–48, 2005. View at Publisher · View at Google Scholar · View at Scopus
  9. A. Katsevich, “A general scheme for constructing inversion algorithms for cone beam CT,” International Journal of Mathematics and Mathematical Sciences, vol. 2003, no. 21, pp. 1305–1321, 2003. View at Publisher · View at Google Scholar · View at Scopus
  10. A. Katsevich, “Theoretically exact filtered backprojection-type inversion algorithm for spiral CT,” SIAM Journal on Applied Mathematics, vol. 62, no. 6, pp. 2012–2026, 2002. View at Publisher · View at Google Scholar · View at Scopus
  11. A. Katsevich, “Improved exact FBP algorithm for spiral CT,” Advances in Applied Mathematics, vol. 32, no. 4, pp. 681–697, 2004. View at Google Scholar
  12. J. Radon, “Über die Bestimmung von Funktionen durch ihre Integralwerte längs gewisser Mannigfaltigkeiten,” Berichte über die Verhandlungen der Sächsische Akademie der Wissenschaften, vol. 69, pp. 267–277, 1917. View at Google Scholar · View at Scopus
  13. P. C. Parks, “On the determination of functions from their integral values along certain manifolds,” IEEE Transactions on Medical Imaging, vol. 5, no. 4, pp. 170–176, 1986. View at Google Scholar · View at Scopus
  14. S. Orlov, “Theory of three dimensional reconstruction: II. The recovery operator,” Soviet Physics—Crystallography, vol. 20, no. 4, pp. 429–433, 1975. View at Google Scholar · View at Scopus
  15. A. N. Pressley, Elementary Differential Geometry, Springer Undergraduate Mathematics Series, Springer, New York, NY, USA, 2001.