1. Introduction

A goal of many research efforts in the current x-ray diagnostic imaging is to greatly enhance the lesion-detection sensitivity. While the digital imaging has currently been a major effort toward to this goal in x-ray diagnostic imaging, but its potential for improving lesion-detection sensitivity is limited by the tissue’s radiological contrast per se. In the over100-years history of x-ray diagnostic imaging, the biological tissue contrast has always been based on the biological tissue’s differences in x-ray attenuation until quite recently. Recently a new class of tissue contrast, that is, the tissue phase-contrast brings new momentum into x-ray imaging research [1–13]. In fact when x-ray traverses body parts the tissues cause x-ray phase changes as well. The amount of the phase change is determined by biological tissue dielectric susceptibility, or equivalently, by the refractive index of the tissue. The refractive index n($\overrightarrow{r}$ _o ) of a tissue voxel at $\overrightarrow{r}$ _o is complex and equal to n($\overrightarrow{r}$ _o ) = 1-Δ($\overrightarrow{r}$ _o ) + iB($\overrightarrow{r}$ _o ), where Δ($\overrightarrow{r}$ _o ) is the refractive index decrement at $\overrightarrow{k}$ _o , and it is responsible for x-ray phase shift. In contrast, the imaginary part B($\overrightarrow{r}$ _o )of n($\overrightarrow{r}$ _o ) is responsible for x-ray attenuation. When x-ray traverses tissues, the x-ray phase change by tissue is given by [1–3]

where λ is the x-ray wavelength, and the integral is over the ray path. In other words, the phase change is proportional to the projection of refractive index decrement Δ. It is interesting to note that tissue’s Δ is much larger than B. We have estimated Δ and B for biological tissues, and found that the issue’s Δ (10^-6-10^-8) is about 1000 times greater thanB (10^-9-10^-11) for x-rays of 10 keV to 150keV[13]. Therefore, the differences in x-ray phase shifts between different tissues are about 1000 times greater than the differences in the projected linear attenuation coefficients. Hence phase-contrast imaging may greatly increase the lesion-detection sensitivity for x-ray imaging.

The x-ray phase tomography is a technique for three-dimensional imaging of tissue refractive index decrement Δ($\overrightarrow{r}$ _o ). With the high phase-sensitivity x-ray phase tomography has a great potential for three-dimensional imaging of soft tissues. In order to develop x-ray phase tomography, it is essential to develop the basic image reconstruction formula, and then reconstruction algorithms can be further developed from the basic reconstruction formula. Bronnikov first developed a basic reconstruction formula for phase tomography of objects with weak attenuations [14–15]. His reconstruction formula can be applied to phase tomography of very thin biological samples imaged by almost parallel x-ray beams from synchrotron-based sources. Assuming a parallel x-ray beam and pure phase objects, his reconstruction formulas require only a single phase-contrast image per tomographic view (projection). For objects with substantial attenuation, he assumed that the attenuation for a given tomographic view is a constant across the field of view, but the attenuation may vary for different views. To apply his reconstruction formula to these cases, one needs to acquire two images per view: one is the “contact print” of the object and another is the phase-contrast image for the same view [15].

In general for phase retrieval one need acquire at least two images that are separated by a distance in x-ray propagation direction [3, 5, 7–9, 11]. According to Eq. (1) the phase-image for a given tomographic view provides the ray-sum map of Δ($\overrightarrow{r}$_o ) needed for 3-D image reconstruction. While a phase-contrast image for a tomographic view contains contributions from both the ray sums of Δ($\overrightarrow{r}$_o ) and B($\overrightarrow{r}$_o ), but acquisition of the single phase-contrast image is not enough to disentangle the ray-sums of Δ($\overrightarrow{r}$_o ) and B($\overrightarrow{r}$_o ) for phase retrieval. For this reason in general at least two images separated by a distance in x-ray propagation direction are needed for the phase retrieval for a given view. This requirement of at least two acquired images for each tomographic view makes the implementation of phase tomography much more difficult than for the projectional phase-contrast imaging, because of the limitations imposed by body motion and radiation dose associated with the large number of tomographic views [9]. In addition this multiple-image requirement impedes the dynamical imaging such as cardiac imaging.

It should be pointed out there is some exception of the multiple-image requirement for phase retrieval. For a pure phase objects with negligible attenuation a single phase-contrast image is enough for the phase-retrieval. However in many practical cases, and especially in clinical application, object attenuation is significant and cannot be ignored at all for phase-retrieval. The average ray-sum of the linear attenuation coefficients μ($\overrightarrow{r}$ _o ) of an abdomen along an anterior-posterior projection is about 4 to 6 for x-rays of 60 keV, and even for small animal the average ray-sum of μ($\overrightarrow{r}$_o ) can be as high as 1 at this energy. Another special case is the case of homogeneous objects of single-material. In this case the ratio of Δ and B is a constant over the object. In this case hence the phase map and attenuation-map of the object are determined by the thickness-map, hence one single phase-contrast image is enough to retrieve the thickness-map [16]. However interesting objects for tomography are all inhomogeneous and this approach cannot be applied.

Quite recently in a brief communication we proposed a new concept of phase-attenuation duality and a new phase-retrieval approach for soft tissues [17], for which only a single phase-contrast image is needed for phase retrieval for objects with any strong attenuation such as thick body parts. The duality-based phase retrieval formula was presented without derivation in that brief communication due to space limitation in [17]. Besides, only the projectional imaging of soft-tissue was discussed in that work and it is natural to extend this new phase-retrieval approach to three-dimensional imaging of refraction index-decrement Δ($\overrightarrow{r}$_o ). In this work, we first present detailed derivation of the phase-retrieval formula based on the phase-attenuation duality. In addition we have included the effects of x-ray source coherence and detector resolution into the phase-retrieval formula. We then combine this formula with the FDK cone-beam reconstruction algorithm [18–19] to derive a three-dimensional phase tomography formula for soft tissue objects of relatively small sizes, such as small animals or small human body parts such as breast. For large objects the FDK reconstruction formula is inaccurate and one may have to use helical scanning for phase tomography. We combined our duality-based phase retrieval approach with the exact cone-beam reconstruction algorithm developed recently by Katsevich [20–22] to derive a phase tomography formula valid for helical scanning of large objects. In the final section we compare these new reconstruction formulas of cone-beam phase tomography with that for the parallel beams in a previous work [27].

2. Phase-attenuation duality and cone beam tomography for soft tissues

Consider a soft-tissue object irradiated by a monochromatic x-ray of wavelength λ emitted from a point source S that is moving on a circle, as is shown in Fig. 1. Note a fixed coordinate frame (x ₁,x ₂,x ₃) is attached with the object. As is the convention in the cone-beam CT literature [18–19], we denote a virtual detector plane passing the origin as the (y ₂,y ₃)-plane shown in Fig. 1. This virtual detector plane represents the real 2-D imaging detector plane a distance of R₂ downstream from the origin. In the cone beam scan geometry, the virtual detector plane is kept perpendicular to the line joining the origin and the source moving along the circle.

 

Fig. 1. Cone-beam tomography with a source orbiting on a circle in (X ₁, X ₂) -plane. The fixed frame is attached with object, while the detector-plane denotes the 2-D imaging detector a distance of R₂ downstream. The imaging detector is always kept perpendicular to the line joining the origin and the source. For definitions of the angles, see the text for details.

Download Full Size | PDF

A cone-beam transform g($\overrightarrow{\theta }$,$\overrightarrow{y}$) of Δ($\overrightarrow{r}$_o ) along an x-ray is given by:

where {$\overrightarrow{\theta }$,$\overrightarrow{\theta }$ _⊥,$\overrightarrow{e}$ ₃} are the unit basis vectors of the moving frame attached to the virtual detector plane D, which passes through the origin. The basis vectors are related to the source-angle β by (Figs. 1–2)

In addition, the position vector in virtual detector plane D is denoted by:

In fact the cone-beam transform g($\overrightarrow{\theta }$,$\overrightarrow{y}$) is equal to a measurement of the line integral of refractive index decrement Δ($\overrightarrow{r}$_o )along the ray. Obviously all the line integrals acquired for a view $\overrightarrow{\theta }$ represent the object’s phase-map for this view. According to Eq. (1) the phase-map ϕ($\overrightarrow{y}$;$\overrightarrow{\theta }$) is given by:

Therefore the cone-beam transforms g($\overrightarrow{\theta }$,$\overrightarrow{y}$) can be found by the phase retrieval of ϕ($\overrightarrow{\theta }$,$\overrightarrow{y}$ ) and a three-dimensional map of Δ($\overrightarrow{r}$_o ) can be reconstructed if ϕ($\overrightarrow{\theta }$,$\overrightarrow{y}$ ) is retrieved for all the

In order to understand how to retrieve the phase of a soft-tissue object, let’s recall how tissue’s phase distribution determines the phase-contrast image of an object. As is pointed out in our previous work, in addition to tissue phase distribution, the limited spatial coherence of the incident x-ray, tissue attenuation and the detector resolutions are all the important factors determining the image intensity [8, 11]. In order to quantitatively account for the effects of spatial coherence of the incident x-ray, we developed a theory based on the Wigner distribution formalism [11]. We found that Ĩ($\overrightarrow{u}$ ; θ), the Fourier transform of the phase-contrast $\overrightarrow{\theta }$, satisfies the following equation [11]:

 

Fig. 2. The moving frame attached to the virtual detector plane.

Download Full Size | PDF

where the symbol F̂ denotes the 2-D Fourier transform and I_in is the intensity of the incident x-ray upon the object, and $\overrightarrow{u}$ is the spatial frequency vector. Here we denote the distance from a x-ray source to the plane of the imaged object as R₁, and the distance from object plane to the detector plane as R₂, hence the geometric magnification factor M is equal to (R₁ + R₂)/ R₁. Interesting readers are referred to ref. [11] for details of derivation of Eq. (6). In Eq. (6) $A_o^2$ represents tissue’s attenuation:

where μ(R ₁ $\overrightarrow{\theta }$ + t($\overrightarrow{y}$ - R ₁ $\overrightarrow{\theta }$)/|$\overrightarrow{y}$ - R ₁ $\overrightarrow{\theta }$|), is the linear attenuation coefficient along the cone beam R ₁ $\overrightarrow{\theta }$ + t($\overrightarrow{y}$ - R ₁ $\overrightarrow{\theta }$)/|$\overrightarrow{y}$ - R ₁ $\overrightarrow{\theta }$|) is the imaginary part of the refractive index. In Eq. (6) OTF ^det($\overrightarrow{u}$ ) is the detector’s spatial frequency response, which is determined by detector pixel 11].

 

Fig. 3. Schematic of an undulator source.

Download Full Size | PDF

It should be noted that the effects on image intensity of spatial coherence of the incident x-ray is accounted for by a multiplication factor ${\tilde{\mu }}_{\mathit{in}}\left(\frac{\lambda R_2\overrightarrow{u}}{M}\right)$ in Eq. (6). For the perfectly coherent incident x-ray one has $\mid {\tilde{\mu }}_{\mathit{in}}\left(\frac{\lambda R_2\overrightarrow{u}}{M}\right)\mid =1$ for all $\left|\frac{\lambda R_2\overrightarrow{u}}{M}\right|$. For partially coherent incident x-ray $\mid {\tilde{\mu }}_{\mathit{in}}\left(\frac{\lambda R_2\overrightarrow{u}}{M}\right)\mid \le 1$ and decreases with increasing $\left|\frac{\lambda R_2\overrightarrow{u}}{M}\right|$. X-ray tubes are incoherent sources, hence we can apply the van Cittert-Zernike theorem to them [23] and we found [11]:

where ${\mathit{OTF}}_{G.U.}\left(\frac{\overrightarrow{u}}{M}\right)$ the optical transfer function (OTF) for the geometrical unsharpness associated with a finite focal spot [23]. For example, for an anode source with a round focus spot we found [11]

where f is the diameter of the focus spot, and J₁(x) is a Bessel function of the 1^st kind. Currently synchrotron radiation such as the undulator radiation is the most common x-ray source for phase contrast imaging experiments. An undulator source is an insertion device installed along the electron beam path in synchrotron, in which a spatially-periodic magnetic field makes electrons to deflect periodically along the beam path as is shown schematically in Fig. 3 [24]. In undulators the electron’s deflection is in order of microns. Due to the interferences of x-ray waves emitted by wiggling electrons under each pole of the undulator, a bright and narrow forward-cone of x-ray beam are generated with photon-energies in resonance-harmonics [24]. Different from the incoherent anode sources, an undulator source is partially coherent at the source level, that is why undulator radiation is narrowly collimated. Hence the van Cittert-Zernike theorem cannot be applied to undulator radiation in general.

For a synchrotron a fully damped and stored electron beam can be well approximated by a Gaussian distribution in transverse position and velocity-divergence [24]. Under the Gaussian-beam assumption we found that

Here σ ₁ denotes total effective horizontal-axis electron beam size, which accounts for the electron-beam’s divergence as well [24]. The total effective vertical-axis electron beam size is denoted by σ ₂.

In many applications including clinical imaging one has $\frac{\mathit{\pi \lambda }{R}_2{\overrightarrow{u}}^2}{M}\ll 1.$ For these cases Eq. (6) is reduced to

Since both the phase-map ϕ($\overrightarrow{y}$ ;$\overrightarrow{\theta }$) and attenuation-map $A_o^2$($\overrightarrow{y}$ ;$\overrightarrow{\theta }$) are unknown, one cannot solve for ϕ($\overrightarrow{y}$ ;$\overrightarrow{\theta }$) by means of computing $\tilde{I}(\frac{\overrightarrow{u}}{M};\overrightarrow{\theta })$ from the image intensity measurement for the view $\overrightarrow{\theta }$. That is why one needs in general at least two acquired images for retrieving a phase-map ϕ($\overrightarrow{y}$ ;$\overrightarrow{\theta }$). As we pointed out earlier this requirement of multiple-image for each

As the tissue phase arises from the x-ray coherent scattering from tissue, the tissue attenuation arises from three x-ray-tissue interactions: the photoelectric absorption, coherent scattering and incoherent scattering [13, 25] in the x-ray energy range employed for medical imaging. Soft tissues are oftentimes imaged in clinical exams and experiments at synchrotron centers. Soft tissues are composed by atoms of predominantly the light elements with Z<10, except trace heavier elements in soft tissue [26]. Quite recently in a brief communication we proposed a new concept of phase-attenuation duality and a new phase-retrieval approach for soft tissues [17], by means of which only a single phase-contrast image is acquired for phase retrieval for objects with any strong attenuation such as thick body parts. As is shown in [17], ϕ($\overrightarrow{y}$ ;$\overrightarrow{\theta }$) and attenuation-map $A_o^2$($\overrightarrow{y}$ ;$\overrightarrow{\theta }$) of soft tissues are really determined by the same map of the projected electron density ρ_e,p ($\overrightarrow{y}$ ;$\overrightarrow{\theta }$).

Considering the different origins in x-ray-tissue interactions for soft-tissue phase (from coherent x-ray scattering) and attenuation (from incoherent x-ray scattering), we called this complementary relationship between phase and attenuation as the phase-attenuation duality [17]. The phase-attenuation duality roots deeply in the complementary quantum-mechanical relationship between the atomic form factor for x-ray coherent scattering and the atomic incoherent scattering function. More specifically, the projected electron density ρ_e,p ($\overrightarrow{y}$ ;$\overrightarrow{\theta }$) is defined as

where ρ_e ($\overrightarrow{r}$_o ) is tissue electron density at position $\overrightarrow{r}$_o . The phase-attenuation duality manifests as [17]:

In Eq. (13)r_e is the classical electron’s radius, and σ_KN is the total cross section for x-ray Compton scattering from a single free electron:

where $\eta =\frac{E_{\mathit{photon}}}{m_e{c}^2},$, here we denote the photon energy of the primary x-ray beam by E_photon , and m_ec ² is the rest electron energy and equal to 511 keV [25].

Taking advantage of phase-attenuation duality, we substitute the duality condition Eq. (13) back into Eq. (11) for simplification. Note that the duality condition Eq. (13) means that

Combining Eq. (11), Eq. (13) and Eq. (15), after a straight-forward calculation we found:

Rewriting Eq. (16) explicitly we found eventually the phase retrieval formula when phase-attenuation duality holds:

where I(M$\overrightarrow{y}$;$\overrightarrow{\theta }$) is the measured phase-contrast image intensity. Note that Eq. (17) generalize the phase retrieval formula in [17] by including the effects of partial coherence and detector resolution. Equivalently, we can rewrite Eq. (17) in terms of the retrieved phase:

Eq. (18) is a key result of this work, and it forms the basis of phase imaging based on the phase-attenuation duality.

3. Cone beam reconstruction formula of Δ($\overrightarrow{r}$_o )

Currently phase tomography experiments mostly involve small size objects. For small objects such as small animals or small human body parts such as breast, one can apply the Feldcamp-Davis-Kress (FDK) reconstruction algorithm for 3-D tomography. In this section we combine the new phase-retrieval formula of Eq. (18) with the FDK reconstruction algorithm to derive a reconstruction formula for 3-D cone-beam phase tomography.

Different from the parallel beam assumption made in previous works [14–15, 27] the object is illuminated by a cone beam of x-ray from source orbiting in a circle enclosing the object. In principle a planar circular orbit is an incomplete orbit for cone beam reconstruction. But for small cone angles the approximate FDK cone-beam reconstruction formula give satisfactory tomographic reconstruction [18–19]. The idea of applying FDK reconstruction algorithm to phase tomography is to first calculate the contribution to Δ($\overrightarrow{r}$_o ) from a given view $\overrightarrow{\theta }$ by means of the conventional fan-beam filtered backprojection. Apparently the contributions to Δ($\overrightarrow{r}$_o ) from different views come from different planes, and these planes forms a sheaf with its vertex at $\overrightarrow{r}$_o . Disregarding this fact the FDK algorithm sums all these contributions over all different views to reconstruct a good approximation of Δ($\overrightarrow{r}$_o ). Consider a x-ray source orbits a circle enclosing the object, as is shown in Fig. 1. The radius of the circle is R ₁, and the coordinates-frame and source angle β are defined in the same way as shown in Fig. 2. Note that the moving frame basis vectors $\overrightarrow{\theta }$ and $\overrightarrow{\theta }$ _⊥ are defined in Eq. (3). Applying the FDK cone-beam reconstruction formula to our case we found by means of Eq. (6) [18–19]:

where (y ₂,y ₃) are related to object point $\overrightarrow{r}$_o as

Let us briefly explain the meaning of terms in Eq. (19) as follows. In phase tomography the ray-sums for each given tomographic view are given by the retrieved phase ϕ($\overrightarrow{y}$ ;$\overrightarrow{\theta }$) of Eq. (18) for that view. The retrieved phase-maps are first weighted by the factor $\frac{R_1}{\sqrt{R_1^2+{\left({y\text{'}}_2\right)}^2+y_3^2}}$ to account for the fan-beam coordinates conversion from the parallel beam coordinates. The weighted data are then filtered row by row along horizontal lines in by means of one-dimensional convolution in y´₂ from -ρ to ρ(the data are assumed be zero for |y´₂|> ρ). The convolution kernel v(y ₂) in Eq. (19) is the conventional ramp filter in 2D tomography [18–19]. Note that the filtering is applied view by view for all the views, though different views are really not co-planar. The front factor $\frac{R_1^2}{{\left(R_1^2-{\overrightarrow{r}}_o\bullet \overrightarrow{\theta }\right)}^2}$ in Eq. (19) is the x ₃ -axis (Fig. 2) away from the mid-plane the reconstruction is approximate, though the average along x ₃ - direction is exact 18–19]. The reason of this inexactness roots in the incomplete 3-D Radon-transform data limited by the co-planar source orbit associated with

 

Fig. 4. Sampled Radon space for the FDK reconstruction algorithm.

Download Full Size | PDF

FDK algorithm. In order to understand this note that a three-dimensional Radon transform of Δ($\overrightarrow{r}$_o ) is a plane-integral of Δ($\overrightarrow{r}$_o ):

where n is the plane’s normal vector, and δ($\overrightarrow{r}$_o ∙ $\overrightarrow{n}$ - s) is a 1-D Dirac-function. As is shown in [18–19], Eq. (21) represents a 3-D Radon transform of Δ($\overrightarrow{r}$_o ) along the plane $\overrightarrow{r}$_o ∙ $\overrightarrow{n}$ = s. Hence a 3-D Radon transform $\widehat{\Delta }$ (s,$\overrightarrow{n}$ ) is represented by a point (s,$\overrightarrow{n}$ )in the space of Radon transforms, or the Radon space. Fig. 4 plots the all such data-points in Radon space that aare sampled by a source in a circular trajectory and the figure shows clearly where the Radon transform data are missing [18–19]. In spite of this the FDK algorithm is the most widely used algorithm for cone beam attenuation-based tomography of small-size objects, such as small animal and human breast. Therefore by means of our duality-based phase retrieval approach Eqs. (18–19) provide good reconstruction formulas of cone beam phase tomography for small animal and human breast.

Although current phase tomography experiments mostly involve small-size objects only, but for sake of completeness let us briefly comment on cone-beam phase-tomography for large objects. As is well known, currently clinical CT imaging with multislice scanners undergoes explosive development with advent of commercially available 64-slice CT scanners. Obviously the helical scanning approach for large objects can be applied to the phase tomography as well. Moreover, recently an exact filtered backprojection reconstruction formula for a helical cone beam CT has been discovered by Katsevich [20–22], and its implementation and generalization is currently an rapidly evolving research area of the attenuation-based tomography. In this work we combine our duality-based phase retrieval approach with Katsevich’s formulas [20–22] to derive a phase tomography formula valid for helical scanning of large objects. Consider now that the source traces a helix of radius R ₁ around the x ₃ -axis of Fig. 1, and the source’s position is specified by $\overrightarrow{a}$ (s)as is shown in Fig. 5:

where P is the pitch of this helical scanning. Obviously s denotes the source-rotation angle. According to Katsevich’s original formula, a three-dimensional object function f($\overrightarrow{r}$_o )can reconstructed from its cone beam projections acquired as the point-source traces a helix:

 

Fig. 5. Helical scanning and the PI-segment.

Download Full Size | PDF

In Eq. (23) we assume that the point-source traces a helical trajectory {$\overrightarrow{a}$ (s)}as is shown in Fig. 5. Let’s first briefly explain what Eq. (23) does for reconstruction. We start with g($\overrightarrow{\Theta }$ ,$\overrightarrow{a}$(s)) which denotes a cone beam projection of f($\overrightarrow{r}$_o ) with a ray emitted from the source at $\overrightarrow{a}$ (s) and the ray-direction unit-vector Θ→. Eq. (23) says that one takes the derivative of g($\overrightarrow{\Theta }$ ,$\overrightarrow{a}$ (s)) with respect to source angle s for the same $\overrightarrow{\Theta }$; and then convolve the derivatives $\frac{1}{\sin \gamma }$[20–22]. Therefore the expression inside the square bracket in Eq. (23) accounts for the differentiating and filtering of the cone-beam projection data. The filtered projection data are then multiplied by a factor of $\frac{1}{\mid {\overrightarrow{r}}_o-\overrightarrow{a}\left(s\right)\mid }$ for distance weighting to account for transverse magnification associated with the cone beam projection. Finally the integral over source angle s is to sum the backprojection of the filtered projections over the different views. Note that the integral interval is I_PI ($\overrightarrow{r}$_o ), which is an interval [s_b ($\overrightarrow{r}$_o ),s_t ($\overrightarrow{r}$_o )]of source angle s such that s_b ($\overrightarrow{r}$_o ) and s_t ($\overrightarrow{r}$_o ) are the endpoints of the PI–line segment that passing $\overrightarrow{r}$_o [20–22]. By definition a PI-line segment means the ‘bottom” endpoint s_b ($\overrightarrow{r}$_o ) and the “top” endpoint.s_t ($\overrightarrow{r}$_o ) are separated by less than a helix turn. Putting together, the Katsevich formula Eq. (23) provides an exact filtered backprojection reconstruction formula for cone beam reconstruction.

We can apply Eq. (23) for cone beam phase tomography as well. The source’s position in the rest frame is specified by $\overrightarrow{a}$ (s), and its orthogonal projection onto the(x ₁,x ₂)-plane is denoted by R ₁ $\overrightarrow{\theta }$(s) in according to similar notation adopted in Fig. 1. In this way the coordinates frame on the virtual detector is set up such that the orthogonal projection of the source point is the origin and {θ _⊥(s),e ₃(s)} is the basis in the same way as shown in Fig. 2. By means of Katsevich formula Eq. (23) we found the reconstruction formula for cone beam phase tomography valid for helical scanning:

where $\overrightarrow{y}$ ^* (s,$\overrightarrow{r}$_o ) is the projection of an object point $\overrightarrow{r}$_o onto the virtual detector plane when the source is at $\overrightarrow{a}$ (s) :

Generalizing Eq. (20) for the circular source-orbit to the case for the helical orbit we found that: for helical scanning

Note that in Eq.(24) ϕ^F ($\overrightarrow{y}$ ^*(s,$\overrightarrow{r}$_o ), s) is the filtered phase defined as

where ϕ´ (y ₂,y ₃, s)is the derivative of the retrieved phase with respective to s at constant ray direction, which is given by [22]

Hence ϕ^F ($\overrightarrow{y}$ ^*(s,$\overrightarrow{r}$_o ),s) of Eq. (27) is essentially the phase derivative convolved with an one-dimensional Hilbert filter in y ₂. Taking advantage of phase-attenuation duality, we calculate ϕ´(y ₂,y ₃,s) in terms of $A_o^2$ as well:

In Eq. (29) $A_o^2$ as can be determined directly from Eq. (17) we have:

It should be noted that Eq. (29) is more computationally efficient that Eq. (28), since unlike Eq. (18) there is no need to calculate logarithms in Eq. (30). Finally we should explain how to define y _3κ(y ₂, $\widehat{\psi }$) in Eq. (27). As is explained by Katsevich, the one-dimensional convolution is performed along a so-called κ–line in virtual detector plane and (y ₂,y _3κ)is a point on a Κ–line [20–22]:

and angle ψ is an offset to source angle s. Therefore κ–lines are labeled by offset angle ψ. On the other hand, the angle $\widehat{\psi }$ in y _3κ(y ₂,$\widehat{\psi }$) is of the smallest |ψ| that labels the κ-line passing through the object-point’s cone beam projection onto virtual detector plane. In other words, $\widehat{\psi }$ is of the smallest |ψ| satisfies the following equation [22]:

4. Discussion and conclusions

The conventional computed tomography a provides three-dimensional map of tissue’s linear attenuation coefficients, and it finds widely spread clinical applications [18–20]. In contrast to conventional computed tomography, x-ray phase tomography provides a three-dimensional map of the tissue refractive index decrement Δ($\overrightarrow{r}$_o ). With the notion of high sensitivity associated with tissue phase changes, as discussed at the beginning of this paper, x-ray phase tomography is expected to be much more sensitive than the conventional CT, hence it has a great potential for clinical applications.

To develop x-ray phase tomography it is important to establish the phase-tomography reconstruction formulas, since such formulas form the basis for various reconstruction algorithms that can be developed. As is shown earlier, for phase tomography at least two acquired images for each tomographic view are needed. This requirement increases the radiation dose to body parts, and makes the scanner design much more difficult as explained earlier. That is why in Bronnikov’s approach of phase tomography one has to acquire two images per view, since one need acquire a “contact print” image in addition to the phase-contrast image [15]. To alleviate this difficulty for phase tomography implementation, we applied the idea of phase-attenuation duality to phase tomography, and in this new approach one needs only one image per view for three-dimensional tomography reconstruction. In this work we presented for the first time the detailed derivation of the phase-retrieval formula based the duality. Furthermore in this derivation we have explicitly included the effects of partial coherence of the incident x-rays. We hope this inclusion of the partial coherence effects can be more useful for phase retrieval works involving the synchrotron and anode sources. Combining the derived phase-retrieval formula with the FDK cone-beam reconstruction algorithm we derived the reconstruction formula Eqs. (18–19) for the three-dimensional phase tomography of small objects. For phase tomography of large objects we combined our duality-based phase retrieval approach with Katsevich’s reconstruction formulas [20–22] to derive the phase tomography formulas, that is, Eq. (24), Eq. (27) and Eqs. (29–30) for helical scanning of large objects.

Another difference of this work from the previous ones is that we deal with the cone beam reconstruction rather that the parallel beams in [14–15, 27]. It should be noted that even though synchrotron radiation, though narrowly forward collimated, is not really parallel, and synchrotron radiation beams have smaller divergence in the vertical direction than that in horizontal direction [24], the cone beam reconstruction will be more accurate than the reconstruction based on parallel beam assumption. For tomography with x-ray tube sources, obviously the cone beam reconstruction is necessary and the parallel beam reconstruction cannot be used for tomography reconstruction. When an x-ray source is far away from the object, the beam can be treated parallel. For example, the Spring-8 synchrotron in Japan has a beam line of length of 1.3km and phase contrast imaging experiments has been performed in this beam line [24]. For parallel beams the reconstruction formulas Eq. (18–19) still holds but the magnification factor M becomes 1,

 

Fig. 6. Parallel beam tomography with a source orbiting along a circle.

Download Full Size | PDF

and the reduced complex degree of coherence $\mid {\tilde{\mu }}_{\mathit{in}}\left(\frac{\lambda R_2\overrightarrow{u}}{M}\right)\mid =1$for any values of |λR ₂ u⊥| 14–15, 27]. This is because with parallel beams one can find all the plane integrals of Δ(r⊥_o ) by the phase retrieval for all tomographic views. With this complete data set one can reconstruct the object by means of the 3-D inverse Radon transform [27]. In contrast with a cone-beam source orbiting in a circle, some plane integrals of Δ(r⊥_o ) are not measured and the data set is incomplete for the exact reconstruction, as is shown in Fig. 4. That is why the FDK reconstruction is approximate and good only for cases small cone angles. For sake of comparison we rewrite the phase tomography reconstruction formula in following, which was derived in previous work for cases with a parallel-beam source orbiting in a circle [27]:

Note that Eq. (33) is essentially the 3-D inverse Radon transform formula [19]. The phase map ϕ($\overrightarrow{y}$ ;$\overrightarrow{\theta }$) in Eq. (33) is retrieved according to Eq. (18) in the same way as in cone-beam phase tomography described earlier. In Eq. (33) n = (cosθ∙sinω,sinθ∙sinω,cosω), s = $\overrightarrow{r}$_o ∙ $\overrightarrow{n}$ , the unit vectors $\overrightarrow{n}$ _D = (sin ω,cos ω), as is depicted in Fig. 6. Note also that while the FDK reconstruction is a filtered back projection algorithm, Eq. (33) is not. But following the approach of Bronnikov [14–15] one can rewrite the 3-D inverse Radon transform in a filtered backprojection form by means of the integration over angle ω,

where the 2-D convolution kernel (filter) is defined as

Obviously Eq. (35) can be rewritten as

However note that in Eq. (36) the convolution symbol ⊗ denotes the two-dimensional convolution, rather the convolutions in Eq. (19) and Eq. (27) for cone-beam tomography are one-dimensional convolutions.

In summary in this work we first presented a detailed derivation of the phase-retrieval formula based on phase-attenuation duality. Moreover we have included the effects of x-ray source coherence and detector resolution into the phase-retrieval formula as well. As only a single image is needed for performing the phase retrieval in this approach, this new approach has great advantages for implementation of phase tomography. We then combine our phase-retrieval formula with the FDK cone beam reconstruction algorithm to provide a three-dimensional phase tomography formula for soft tissue objects of relatively small size, such as small animals or human breast. For large objects we briefly discussed how to apply Katsevich’s cone beam reconstruction formula to the helical phase tomography for large objects as well.