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
P A. Snigirev, I. Snigireva, and V. Kohn, et al, “On the possibilities of x-ray phase contrast micro-imaging by coherent high-energy synchrotron radiation,” Rev. Sci. Instrum.
66, 5486–5492, (1995). [CrossRef]
]. 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(
o
) of a tissue voxel at
o
is complex and equal to
n(
o
) = 1-Δ(
o
) +
iB(
o
), where Δ(
o
) is the refractive index decrement at
o
, and it is responsible for x-ray phase shift. In contrast, the imaginary part B(
o
)of
n(
o
) is responsible for x-ray attenuation. When x-ray traverses tissues, the x-ray phase change by tissue is given by [
1–3
P A. Snigirev, I. Snigireva, and V. Kohn, et al, “On the possibilities of x-ray phase contrast micro-imaging by coherent high-energy synchrotron radiation,” Rev. Sci. Instrum.
66, 5486–5492, (1995). [CrossRef]
]
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
X. Wu, A. Dean, and H Liu, “X-ray diagnostic techniques,” in Biomedical photonics handbook, T. VoDinh ed., Chapter 26, p.26
26–34, (CRC Press, Tampa, Florida, 2003).
]. 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 Δ(
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
A. Bronnikov, “Reconstruction formulas in phase-contrast tomography,” Optics Commun.
171, 29–244 (1999) [CrossRef]
]. 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
A. Bronnikov, “Theory of quantitative phase-contrast computed tomography,” J. Opt. Soc. Am. A
19, 472–480 (2002). [CrossRef]
].
In general for phase retrieval one need acquire at least two images that are separated by a distance in x-ray propagation direction [
3
K. Nugent, T. Gureyev, D. Cookson, D. Paganin, and Z. Barnea, “Quantitative phase imaging using hard x rays,” Phys. Rev. Lett.
77,2961–2965, (1996). [CrossRef] [PubMed]
,
5
D. Paganin and K. Nugent, “Noninterferometric phase imaging with partialcoherent light,” Phys. Rev. Lett.
80, 2586–2589 (1998). [CrossRef]
,
7–9
X, Wu and H. Liu, “A general formalism for x-ray phase contrast imaging,” J. X-Ray Sci. Technol.
11, 33–42 (2003).
,
11
X. Wu and H. Liu, “A new theory of phase-contrast x-ray imaging based on Wigner distributions,” Med. Phys.
31, 2380–2384 (2004). [CrossRef]
]. According to Eq. (
1) the phase-image for a given tomographic view provides the ray-sum map of Δ(
o
) needed for 3-D image reconstruction. While a phase-contrast image for a tomographic view contains contributions from both the ray sums of Δ(
o
) and B(
o
), but acquisition of the single phase-contrast image is not enough to disentangle the ray-sums of Δ(
o
) and B(
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
X. Wu and H. Liu, “A dual detector approach for X-ray attenuation and phase imaging,” J. X-Ray Sci. Technol.
12, 35–42, 2004.
]. 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
μ(
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
μ(
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
D. Paganin, S. Mayo, T. Gureyev, P. Miller, and S. Wilkins, “Simultaneous phase and amplitude extraction from a single defocused image of a homogeneous object,” J. Microsc.
206, 33–40 (2002). [CrossRef] [PubMed]
]. 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
X. Wu, H. Liu, and A. Yan, “X-ray phase-attenuation duality and phase retrieval,” Optics Lett.
30, 379–381 (2005). [CrossRef]
], 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
X. Wu, H. Liu, and A. Yan, “X-ray phase-attenuation duality and phase retrieval,” Optics Lett.
30, 379–381 (2005). [CrossRef]
]. 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 Δ(
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
C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging , (IEEE Press, New York, 1987).
] 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
A. Katsevich, “Analysis of an exact inversion algorithm for spiral cone-beam CT,” Phys. Med. Biol.
47, 2583–2597 (2003). [CrossRef]
] 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
X. Wu and H. Liu, “A reconstruction formula for soft tissue X-ray phase tomography,” J. X-Ray Sci. Technol.
12, 273–279 (2004).
].
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
1,
x
2,
x
3) is attached with the object. As is the convention in the cone-beam CT literature [
18–19
C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging , (IEEE Press, New York, 1987).
], we denote a virtual detector plane passing the origin as the (
y
2,
y
3)-plane shown in
Fig. 1. This virtual detector plane represents the real 2-D imaging detector plane a distance of R
2 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
1, X
2) -plane. The fixed frame is attached with object, while the detector-plane denotes the 2-D imaging detector a distance of R2 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.
A cone-beam transform g(,) of Δ(o
) along an x-ray is given by:
where {
,
⊥,
3} 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(
,
) is equal to a measurement of the line integral of refractive index decrement Δ(
o
)along the ray. Obviously all the line integrals acquired for a view
represent the object’s phase-map for this view. According to Eq. (
1) the phase-map
ϕ(
;
) is given by:
Therefore the cone-beam transforms g(,) can be found by the phase retrieval of ϕ(,
) and a three-dimensional map of Δ(o
) can be reconstructed if ϕ(,
) 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
X. Wu and H. Liu, “Clinical implementation of phase contrast x-ray imaging: theoretical foundation and design considerations,” Med. Phys.
30, 2169–2179 (2003). [CrossRef] [PubMed]
,
11
X. Wu and H. Liu, “A new theory of phase-contrast x-ray imaging based on Wigner distributions,” Med. Phys.
31, 2380–2384 (2004). [CrossRef]
]. 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
X. Wu and H. Liu, “A new theory of phase-contrast x-ray imaging based on Wigner distributions,” Med. Phys.
31, 2380–2384 (2004). [CrossRef]
]. We found that
Ĩ(
;
θ), the Fourier transform of the phase-contrast
, satisfies the following equation [
11
X. Wu and H. Liu, “A new theory of phase-contrast x-ray imaging based on Wigner distributions,” Med. Phys.
31, 2380–2384 (2004). [CrossRef]
]:
Fig. 2. The moving frame attached to the virtual detector plane.
where the symbol
F̂ denotes the 2-D Fourier transform and
Iin
is the intensity of the incident x-ray upon the object, and
is the spatial frequency vector. Here we denote the distance from a x-ray source to the plane of the imaged object as R
1, and the distance from object plane to the detector plane as R
2, hence the geometric magnification factor M is equal to (R
1 + R
2)/ R
1. Interesting readers are referred to ref. [
11
X. Wu and H. Liu, “A new theory of phase-contrast x-ray imaging based on Wigner distributions,” Med. Phys.
31, 2380–2384 (2004). [CrossRef]
] for details of derivation of Eq. (
6). In Eq. (
6)
represents tissue’s attenuation:
where
μ(
R
1
+
t(
-
R
1
)/|
-
R
1
|), is the linear attenuation coefficient along the cone beam
R
1
+
t(
-
R
1
)/|
-
R
1
|) is the imaginary part of the refractive index. In Eq. (
6)
OTF
det(
) is the detector’s spatial frequency response, which is determined by detector pixel
11
X. Wu and H. Liu, “A new theory of phase-contrast x-ray imaging based on Wigner distributions,” Med. Phys.
31, 2380–2384 (2004). [CrossRef]
].
Fig. 3. Schematic of an undulator source.
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
in Eq. (
6). For the perfectly coherent incident x-ray one has
for all
. For partially coherent incident x-ray
and decreases with increasing
. X-ray tubes are incoherent sources, hence we can apply the van Cittert-Zernike theorem to them [
23
M. Born and E. Wolf, Principle of Optics , 6th ed. (Pergamon, Oxford
1980).
] and we found [
11
X. Wu and H. Liu, “A new theory of phase-contrast x-ray imaging based on Wigner distributions,” Med. Phys.
31, 2380–2384 (2004). [CrossRef]
]:
where
the optical transfer function (OTF) for the geometrical unsharpness associated with a finite focal spot [
23
M. Born and E. Wolf, Principle of Optics , 6th ed. (Pergamon, Oxford
1980).
]. For example, for an anode source with a round focus spot we found [
11
X. Wu and H. Liu, “A new theory of phase-contrast x-ray imaging based on Wigner distributions,” Med. Phys.
31, 2380–2384 (2004). [CrossRef]
]
where f is the diameter of the focus spot, and J
1(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
H. Wiedemann, Synchrotron Radiation , (Springer-verlag, Berlin Heidelberg
2003). [CrossRef]
]. 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
H. Wiedemann, Synchrotron Radiation , (Springer-verlag, Berlin Heidelberg
2003). [CrossRef]
]. 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
H. Wiedemann, Synchrotron Radiation , (Springer-verlag, Berlin Heidelberg
2003). [CrossRef]
]. Under the Gaussian-beam assumption we found that
Here
σ
1 denotes total effective horizontal-axis electron beam size, which accounts for the electron-beam’s divergence as well [
24
H. Wiedemann, Synchrotron Radiation , (Springer-verlag, Berlin Heidelberg
2003). [CrossRef]
]. The total effective vertical-axis electron beam size is denoted by
σ
2.
In many applications including clinical imaging one has
For these cases Eq. (
6) is reduced to
Since both the phase-map ϕ(
;) and attenuation-map (
;) are unknown, one cannot solve for ϕ(
;) by means of computing from the image intensity measurement for the view . That is why one needs in general at least two acquired images for retrieving a phase-map ϕ(
;). 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
X. Wu, A. Dean, and H Liu, “X-ray diagnostic techniques,” in Biomedical photonics handbook, T. VoDinh ed., Chapter 26, p.26
26–34, (CRC Press, Tampa, Florida, 2003).
,
25
N. A. Dyson, X-rays in atomic and Nuclear physics , (Longman, 1973).
] 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
Tissue Substitutes in Radiation Dosimetry and Measurement, Report 44 of the International Commission on Radiation Units and Measurements (Bethesda, MD, 1989).
]. 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
X. Wu, H. Liu, and A. Yan, “X-ray phase-attenuation duality and phase retrieval,” Optics Lett.
30, 379–381 (2005). [CrossRef]
], 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
X. Wu, H. Liu, and A. Yan, “X-ray phase-attenuation duality and phase retrieval,” Optics Lett.
30, 379–381 (2005). [CrossRef]
],
ϕ(
;
) and attenuation-map
(
;
) of soft tissues are really determined by the same map of the projected electron density
ρe,p
(
;
).
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
X. Wu, H. Liu, and A. Yan, “X-ray phase-attenuation duality and phase retrieval,” Optics Lett.
30, 379–381 (2005). [CrossRef]
]. 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
(
;
) is defined as
where
ρe
(
o
) is tissue electron density at position
o
. The phase-attenuation duality manifests as [
17
X. Wu, H. Liu, and A. Yan, “X-ray phase-attenuation duality and phase retrieval,” Optics Lett.
30, 379–381 (2005). [CrossRef]
]:
In Eq. (
13
X. Wu, A. Dean, and H Liu, “X-ray diagnostic techniques,” in Biomedical photonics handbook, T. VoDinh ed., Chapter 26, p.26
26–34, (CRC Press, Tampa, Florida, 2003).
)
re
is the classical electron’s radius, and
σKN
is the total cross section for x-ray Compton scattering from a single free electron:
where
, here we denote the photon energy of the primary x-ray beam by
Ephoton
, and
mec
2 is the rest electron energy and equal to 511 keV [
25
N. A. Dyson, X-rays in atomic and Nuclear physics , (Longman, 1973).
].
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;
) is the measured phase-contrast image intensity. Note that Eq. (
17) generalize the phase retrieval formula in [
17
X. Wu, H. Liu, and A. Yan, “X-ray phase-attenuation duality and phase retrieval,” Optics Lett.
30, 379–381 (2005). [CrossRef]
] 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 Δ(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
A. Bronnikov, “Reconstruction formulas in phase-contrast tomography,” Optics Commun.
171, 29–244 (1999) [CrossRef]
,
27
X. Wu and H. Liu, “A reconstruction formula for soft tissue X-ray phase tomography,” J. X-Ray Sci. Technol.
12, 273–279 (2004).
] 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
C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging , (IEEE Press, New York, 1987).
]. The idea of applying FDK reconstruction algorithm to phase tomography is to first calculate the contribution to Δ(
o
) from a given view
by means of the conventional fan-beam filtered backprojection. Apparently the contributions to Δ(
o
) from different views come from different planes, and these planes forms a sheaf with its vertex at
o
. Disregarding this fact the FDK algorithm sums all these contributions over all different views to reconstruct a good approximation of Δ(
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
1, 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
and
⊥ are defined in Eq. (
3). Applying the FDK cone-beam reconstruction formula to our case we found by means of Eq. (
6) [
18–19
C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging , (IEEE Press, New York, 1987).
]:
where (y
2,y
3) are related to object point o
as
Let us briefly explain the meaning of terms in Eq. (
19
F. Natterer, The Mathematics of Computerized Tomography , (SIAM, Philadelphia, 2001). [CrossRef]
) as follows. In phase tomography the ray-sums for each given tomographic view are given by the retrieved phase
ϕ(
;
) of Eq. (
18
C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging , (IEEE Press, New York, 1987).
) for that view. The retrieved phase-maps are first weighted by the factor
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´
2 from -
ρ to
ρ(the data are assumed be zero for |
y´
2|>
ρ). The convolution kernel
v(
y
2) in Eq. (
19) is the conventional ramp filter in 2D tomography [
18–19
C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging , (IEEE Press, New York, 1987).
]. Note that the filtering is applied view by view for all the views, though different views are really not co-planar. The front factor
in Eq. (
19) is the
x
3 -axis (
Fig. 2) away from the mid-plane the reconstruction is approximate, though the average along
x
3 -
direction is exact
18–19
C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging , (IEEE Press, New York, 1987).
]. 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.
FDK algorithm. In order to understand this note that a three-dimensional Radon transform of Δ(o
) is a plane-integral of Δ(o
):
where
n is the plane’s normal vector, and
δ(
o
∙
-
s) is a 1-D Dirac-function. As is shown in [
18–19
C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging , (IEEE Press, New York, 1987).
], Eq. (
21) represents a 3-D Radon transform of Δ(
o
) along the plane
o
∙
=
s. Hence a 3-D Radon transform
(
s,
) is represented by a point (
s,
)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
C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging , (IEEE Press, New York, 1987).
]. 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
A. Katsevich, “Analysis of an exact inversion algorithm for spiral cone-beam CT,” Phys. Med. Biol.
47, 2583–2597 (2003). [CrossRef]
], 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
A. Katsevich, “Analysis of an exact inversion algorithm for spiral cone-beam CT,” Phys. Med. Biol.
47, 2583–2597 (2003). [CrossRef]
] to derive a phase tomography formula valid for helical scanning of large objects. Consider now that the source traces a helix of radius
R
1 around the
x
3 -axis of
Fig. 1, and the source’s position is specified by
(
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(o
)can reconstructed from its cone beam projections acquired as the point-source traces a helix:
Fig. 5. Helical scanning and the PI-segment.
In Eq. (
23) we assume that the point-source traces a helical trajectory {
(
s)}as is shown in
Fig. 5. Let’s first briefly explain what Eq. (
23) does for reconstruction. We start with
g(
,
(
s)) which denotes a cone beam projection of
f(
o
) with a ray emitted from the source at
(
s) and the ray-direction unit-vector Θ→. Eq. (
23) says that one takes the derivative of
g(
,
(
s)) with respect to source angle s for the same
; and then convolve the derivatives
[
20–22
A. Katsevich, “Analysis of an exact inversion algorithm for spiral cone-beam CT,” Phys. Med. Biol.
47, 2583–2597 (2003). [CrossRef]
]. 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
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
IPI
(
o
), which is an interval [
sb
(
o
),
st
(
o
)]of source angle
s such that
sb
(
o
) and
st
(
o
) are the endpoints of the PI–line segment that passing
o
[
20–22
A. Katsevich, “Analysis of an exact inversion algorithm for spiral cone-beam CT,” Phys. Med. Biol.
47, 2583–2597 (2003). [CrossRef]
]. By definition a PI-line segment means the ‘bottom” endpoint
sb
(
o
) and the “top” endpoint.
st
(
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
(
s), and its orthogonal projection onto the(
x
1,
x
2)-plane is denoted by
R
1
(
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
3(
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
* (s,o
) is the projection of an object point o
onto the virtual detector plane when the source is at
(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
(
*(
s,
o
),
s) is the filtered phase defined as
where
ϕ´ (
y
2,
y
3,
s)is the derivative of the retrieved phase with respective to
s at constant ray direction, which is given by [
22
F. Noo, J. Pack, and D. Heuscher, “Exact helical reconstruction using native cone-beam geometry,” Phys. Med. Biol.
48, 3787–3818 (2003). [CrossRef]
]
Hence
ϕF
(
*(
s,
o
),
s) of Eq. (
27) is essentially the phase derivative convolved with an one-dimensional Hilbert filter in
y
2. Taking advantage of phase-attenuation duality, we calculate
ϕ´(
y
2,
y
3,
s) in terms of
as well:
In Eq. (
29)
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
2,
) 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
2,
y
3κ
)is a point on a
Κ–line [
20–22
A. Katsevich, “Analysis of an exact inversion algorithm for spiral cone-beam CT,” Phys. Med. Biol.
47, 2583–2597 (2003). [CrossRef]
]:
and angle
ψ is an offset to source angle
s. Therefore
κ–lines are labeled by offset angle
ψ. On the other hand, the angle
in
y
3κ
(
y
2,
) is of the smallest |
ψ| that labels the
κ-line passing through the object-point’s cone beam projection onto virtual detector plane. In other words,
is of the smallest |
ψ| satisfies the following equation [
22
F. Noo, J. Pack, and D. Heuscher, “Exact helical reconstruction using native cone-beam geometry,” Phys. Med. Biol.
48, 3787–3818 (2003). [CrossRef]
]:
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
C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging , (IEEE Press, New York, 1987).
]. In contrast to conventional computed tomography, x-ray phase tomography provides a three-dimensional map of the tissue refractive index decrement Δ(
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
A. Bronnikov, “Theory of quantitative phase-contrast computed tomography,” J. Opt. Soc. Am. A
19, 472–480 (2002). [CrossRef]
]. 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
A. Katsevich, “Analysis of an exact inversion algorithm for spiral cone-beam CT,” Phys. Med. Biol.
47, 2583–2597 (2003). [CrossRef]
] 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
A. Bronnikov, “Reconstruction formulas in phase-contrast tomography,” Optics Commun.
171, 29–244 (1999) [CrossRef]
,
27
X. Wu and H. Liu, “A reconstruction formula for soft tissue X-ray phase tomography,” J. X-Ray Sci. Technol.
12, 273–279 (2004).
]. 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
H. Wiedemann, Synchrotron Radiation , (Springer-verlag, Berlin Heidelberg
2003). [CrossRef]
], 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
H. Wiedemann, Synchrotron Radiation , (Springer-verlag, Berlin Heidelberg
2003). [CrossRef]
]. 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.
and the reduced complex degree of coherence
for any values of |
λR
2
u⊥|
14–15
A. Bronnikov, “Reconstruction formulas in phase-contrast tomography,” Optics Commun.
171, 29–244 (1999) [CrossRef]
,
27
X. Wu and H. Liu, “A reconstruction formula for soft tissue X-ray phase tomography,” J. X-Ray Sci. Technol.
12, 273–279 (2004).
]. 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
X. Wu and H. Liu, “A reconstruction formula for soft tissue X-ray phase tomography,” J. X-Ray Sci. Technol.
12, 273–279 (2004).
]. 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
X. Wu and H. Liu, “A reconstruction formula for soft tissue X-ray phase tomography,” J. X-Ray Sci. Technol.
12, 273–279 (2004).
]:
Note that Eq. (
33) is essentially the 3-D inverse Radon transform formula [
19
F. Natterer, The Mathematics of Computerized Tomography , (SIAM, Philadelphia, 2001). [CrossRef]
]. The phase map
ϕ(
;
) 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 =
o
∙
, the unit vectors
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
A. Bronnikov, “Reconstruction formulas in phase-contrast tomography,” Optics Commun.
171, 29–244 (1999) [CrossRef]
] 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.