Propagation-based Phase Contrast
From the complex refractive index to single-distance PAD phase retrieval
This chapter connects the complex refractive index, projected attenuation and phase, free-space propagation, the transport-of-intensity equation (TIE), and the single-distance phase–attenuation duality (PAD) retrieval formula. The derivation uses angular spatial frequencies unless stated otherwise.
Starting from the Complex Refractive Index
For X-rays, the complex refractive index of a material is generally written as
\[ n=1-\delta+i\beta, \]
where \(\delta\) describes the phase shift induced by the material and \(\beta\) describes X-ray absorption. The wavenumber in vacuum is
\[ k=\frac{2\pi}{\lambda}, \]
where \(\lambda\) is the X-ray wavelength.
Consider an incident monochromatic plane wave:
\[ \Psi_{\mathrm{in}}(z)=A_{\mathrm{in}}e^{ikz}. \]
If the X-ray propagates through a homogeneous material of thickness \(t\), the transmitted wave becomes
\[ \Psi=A_{\mathrm{in}}e^{iknt}. \]
Substituting the complex refractive index gives
\[ \begin{aligned} \Psi &=A_{\mathrm{in}}\exp\left[ik(1-\delta+i\beta)t\right] \\ &=A_{\mathrm{in}}e^{ikt}e^{-ik\delta t}e^{-k\beta t}. \end{aligned} \]
The three exponential factors represent:
- \(e^{ikt}\): ordinary propagation in vacuum;
- \(e^{-ik\delta t}\): the material-induced phase shift;
- \(e^{-k\beta t}\): attenuation of the wave amplitude.
After removing the common vacuum propagation phase, which does not affect the measured intensity, the exit wave immediately downstream of the sample is
\[ \Psi_0=A_{\mathrm{in}}e^{-k\beta t}e^{-ik\delta t}. \]
Obtaining Intensity and Phase from the Complex Wavefield
The complex exit wave can be written in amplitude–phase form:
\[ \Psi_0(x,y)=\sqrt{I_0(x,y)}e^{i\phi(x,y)}. \]
Comparing the amplitude–phase representation with the exit-wave expression gives
\[ \sqrt{I_0}=\sqrt{I_{\mathrm{in}}}e^{-k\beta t}, \]
and therefore
\[ I_0=I_{\mathrm{in}}e^{-2k\beta t}. \]
The linear attenuation coefficient is
\[ \boxed{\mu=2k\beta=\frac{4\pi\beta}{\lambda}}. \]
Consequently,
\[ \boxed{I_0=I_{\mathrm{in}}e^{-\mu t}}, \]
which is the Beer–Lambert law. At the same time, the material-induced phase shift is
\[ \boxed{\phi=-k\delta t}. \]
Thus, the same material modifies both the amplitude and phase of the incident X-ray wave.
Extension to an Inhomogeneous Three-Dimensional Object
For a general three-dimensional object, \(\delta\) and \(\beta\) vary with spatial position. Along an X-ray path \(L\), the projected phase shift is
\[ \phi(x,y)=-k\int_L\delta(x,y,\ell)\,\mathrm{d}\ell. \]
The normalized exit intensity satisfies
\[ \frac{I_0(x,y)}{I_{\mathrm{in}}(x,y)} = \exp\left[ -2k\int_L\beta(x,y,\ell)\,\mathrm{d}\ell \right]. \]
Taking the negative logarithm gives the absorption projection
\[ \begin{aligned} p(x,y) &=-\ln\frac{I_0(x,y)}{I_{\mathrm{in}}(x,y)} \\ &=2k\int_L\beta(x,y,\ell)\,\mathrm{d}\ell \\ &=\int_L\mu(x,y,\ell)\,\mathrm{d}\ell. \end{aligned} \]
The two projected quantities are therefore
\[ \phi=-k\int_L\delta\,\mathrm{d}\ell, \qquad p=2k\int_L\beta\,\mathrm{d}\ell. \]
The first is the projected phase shift. The second is the conventional attenuation line integral required for CT reconstruction.
Phase–Attenuation Duality
The key assumption of the phase–attenuation duality (PAD) method is that
\[ \gamma=\frac{\delta}{\beta} \]
is approximately constant throughout the object. Thus, \(\delta=\gamma\beta\), and the projected phase becomes
\[ \phi=-k\gamma\int_L\beta\,\mathrm{d}\ell. \]
Because the attenuation projection gives
\[ \int_L\beta\,\mathrm{d}\ell=\frac{p}{2k}, \]
the phase and attenuation projections are related by
\[ \boxed{ \phi=-\frac{\gamma}{2}p =-\frac{\delta}{2\beta}p }. \]
Using
\[ p=-\ln\frac{I_0}{I_{\mathrm{in}}}, \]
we obtain
\[ \boxed{ \phi(x,y) = \frac{\delta}{2\beta} \ln\frac{I_0(x,y)}{I_{\mathrm{in}}(x,y)} }. \]
Since \(I_0<I_{\mathrm{in}}\) for an attenuating object, the logarithm is generally negative. The recovered phase is therefore normally negative, consistent with \(\phi=-k\delta t\).
The PAD relation connects exit-surface intensity and phase, reducing two unknown functions to one. Without an additional constraint of this kind, a single measured intensity image is generally insufficient to recover both wavefield amplitude and phase.
How Phase Variations Become Intensity Variations
An X-ray detector measures intensity,
\[ I=|\Psi|^2, \]
but does not directly measure phase. Immediately after the sample, part of the object information is stored in \(\phi(x,y)\). During free-space propagation, transverse phase gradients redirect the local energy flow:
- some regions gain intensity and become brighter;
- other regions lose intensity and become darker.
This conversion is described by the transport-of-intensity equation:
\[ \boxed{ \frac{\partial I}{\partial z} = -\frac{1}{k} \nabla_\perp\cdot \left(I\nabla_\perp\phi\right) }, \]
where
\[ \nabla_\perp = \left( \frac{\partial}{\partial x}, \frac{\partial}{\partial y} \right) \]
is the transverse gradient operator.
Derivation of the Transport-of-Intensity Equation
In free space, the scalar X-ray wavefield satisfies the Helmholtz equation:
\[ \nabla^2\Psi+k^2\Psi=0. \]
Separate the rapidly oscillating carrier from the slowly varying envelope:
\[ \Psi(x,y,z)=\psi(x,y,z)e^{ikz}. \]
The longitudinal derivative is
\[ \frac{\partial^2\Psi}{\partial z^2} = e^{ikz} \left[ \frac{\partial^2\psi}{\partial z^2} +2ik\frac{\partial\psi}{\partial z} -k^2\psi \right], \]
and the transverse Laplacian is
\[ \nabla_\perp^2\Psi=e^{ikz}\nabla_\perp^2\psi. \]
Substitution into the Helmholtz equation gives
\[ \nabla_\perp^2\psi + \frac{\partial^2\psi}{\partial z^2} + 2ik\frac{\partial\psi}{\partial z} =0. \]
Under the paraxial approximation,
\[ \left| \frac{\partial^2\psi}{\partial z^2} \right| \ll \left| 2ik\frac{\partial\psi}{\partial z} \right|, \]
so the second-order longitudinal derivative is neglected. The paraxial wave equation is then
\[ 2ik\frac{\partial\psi}{\partial z} +\nabla_\perp^2\psi=0, \]
or
\[ \boxed{ \frac{\partial\psi}{\partial z} = \frac{i}{2k}\nabla_\perp^2\psi }. \]
Writing the Wavefield in Amplitude–Phase Form
Let
\[ \psi=\sqrt{I}e^{i\phi}, \qquad I=\psi^*\psi. \]
Differentiating the intensity with respect to \(z\) gives
\[ \frac{\partial I}{\partial z} = \psi^*\frac{\partial\psi}{\partial z} + \psi\frac{\partial\psi^*}{\partial z}. \]
Using the paraxial wave equation and its complex conjugate,
\[ \frac{\partial I}{\partial z} = \frac{i}{2k} \left( \psi^*\nabla_\perp^2\psi - \psi\nabla_\perp^2\psi^* \right). \]
The vector identity
\[ \psi^*\nabla^2\psi-\psi\nabla^2\psi^* = \nabla\cdot \left( \psi^*\nabla\psi-\psi\nabla\psi^* \right) \]
gives
\[ \frac{\partial I}{\partial z} = \frac{i}{2k} \nabla_\perp\cdot \left( \psi^*\nabla_\perp\psi - \psi\nabla_\perp\psi^* \right). \]
For \(\psi=\sqrt{I}e^{i\phi}\),
\[ \psi^*\nabla_\perp\psi - \psi\nabla_\perp\psi^* = 2iI\nabla_\perp\phi. \]
Therefore,
\[ \frac{\partial I}{\partial z} = \frac{i}{2k} \nabla_\perp\cdot \left(2iI\nabla_\perp\phi\right) = -\frac{1}{k} \nabla_\perp\cdot \left(I\nabla_\perp\phi\right), \]
which recovers the transport-of-intensity equation.
Physical Interpretation of the TIE
The TIE can be written like a continuity equation:
\[ \frac{\partial I}{\partial z} + \nabla_\perp\cdot\mathbf{J}_\perp=0, \]
where
\[ \mathbf{J}_\perp = \frac{I}{k}\nabla_\perp\phi \]
is the transverse energy-flow density.
When \(\nabla_\perp\phi\ne0\), the local energy flow is deflected:
- if \(\nabla_\perp\cdot\mathbf{J}_\perp<0\), energy converges and intensity increases;
- if \(\nabla_\perp\cdot\mathbf{J}_\perp>0\), energy diverges and intensity decreases.
Material boundaries generally produce large phase gradients. After propagation, they appear as bright–dark edge fringes in the measured image.
From the Differential Equation to a Single-Distance Model
Let \(I_0(x,y)\) be the exit-surface intensity and \(I_z(x,y)\) the intensity measured after propagation distance \(z\). A first-order Taylor expansion gives
\[ I_z(x,y) = I_0(x,y) + z\left. \frac{\partial I}{\partial z} \right|_{z=0} + O(z^2). \]
Neglecting higher-order terms and substituting the transport-of-intensity equation,
\[ \boxed{ I_z = I_0 - \frac{z}{k} \nabla_\perp\cdot \left(I_0\nabla_\perp\phi\right) }. \]
At this stage, the equation contains two unknown functions: \(I_0(x,y)\) and \(\phi(x,y)\). The PAD relation is used to eliminate \(\phi\).
Substituting the PAD Relation into the TIE
From the PAD intensity–phase relation,
\[ \phi = \frac{\delta}{2\beta} \ln\frac{I_0}{I_{\mathrm{in}}}. \]
After flat-field correction, assume that \(I_{\mathrm{in}}\) is approximately uniform over the local image region. Then
\[ \begin{aligned} \nabla_\perp\phi &= \frac{\delta}{2\beta}\nabla_\perp\ln I_0 \\ &= \frac{\delta}{2\beta} \frac{\nabla_\perp I_0}{I_0}. \end{aligned} \]
Thus,
\[ I_0\nabla_\perp\phi = \frac{\delta}{2\beta}\nabla_\perp I_0, \]
and
\[ \nabla_\perp\cdot \left(I_0\nabla_\perp\phi\right) = \frac{\delta}{2\beta}\nabla_\perp^2 I_0. \]
Substitution into the first-order propagation model gives
\[ \boxed{ I_z = I_0 - \frac{z\delta}{2k\beta} \nabla_\perp^2 I_0 }. \]
Define
\[ \boxed{ \alpha = \frac{z\delta}{2k\beta} = \frac{\lambda z}{4\pi} \frac{\delta}{\beta} }. \]
The propagation model becomes
\[ \boxed{ I_z = \left(1-\alpha\nabla_\perp^2\right)I_0 }. \]
This is the central result of the single-distance PAD derivation.
Why Free-Space Propagation Enhances Edges
The transverse Laplacian is
\[ \nabla_\perp^2 = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}. \]
It is small in slowly varying regions and large near sharp boundaries. The forward propagation model can therefore be interpreted as the sum of an exit-surface absorption image and a propagation-induced edge term:
\[ I_z=I_0-\alpha\nabla_\perp^2I_0. \]
For a step-like boundary, the second derivative changes sign across the edge. One side becomes brighter and the other darker. The observed bright–dark pair is not evidence of two physical boundaries; it is an intensity modulation generated by free-space propagation.
Normalized Form of the Propagation Model
Define the normalized exit-surface and propagated intensities:
\[ f(x,y)=\frac{I_0(x,y)}{I_{\mathrm{in}}(x,y)}, \qquad g(x,y)=\frac{I_z(x,y)}{I_{\mathrm{in}}(x,y)}. \]
Experimentally, \(g\) is obtained through dark-field and flat-field correction:
\[ \boxed{ g(x,y) = \frac{P(x,y)-D(x,y)} {F(x,y)-D(x,y)} }, \]
where \(P\) is the sample projection, \(D\) is the dark-field image, and \(F\) is the flat-field image.
If the incident intensity is locally uniform,
\[ \boxed{ g = \left(1-\alpha\nabla_\perp^2\right)f }. \]
The measured quantity is \(g\), while the desired exit-surface intensity is \(f\). Phase retrieval is therefore the inverse problem \(g\rightarrow f\).
Solving the Inverse Problem in the Fourier Domain
Let
\[ \widehat{f}(q_x,q_y)=\mathcal{F}\{f(x,y)\}, \]
where \(q_x\) and \(q_y\) are angular spatial frequencies in \(\mathrm{rad/m}\). The Fourier derivative properties are
\[ \mathcal{F} \left\{ \frac{\partial^2f}{\partial x^2} \right\} =-q_x^2\widehat{f}, \qquad \mathcal{F} \left\{ \frac{\partial^2f}{\partial y^2} \right\} =-q_y^2\widehat{f}. \]
Therefore,
\[ \mathcal{F}\{\nabla_\perp^2f\} = -\left(q_x^2+q_y^2\right)\widehat{f}. \]
Taking the Fourier transform of the normalized forward model gives
\[ \widehat{g} = \left[ 1+\alpha\left(q_x^2+q_y^2\right) \right]\widehat{f}. \]
Solving for \(\widehat{f}\),
\[ \boxed{ \widehat{f} = \frac{\widehat{g}} {1+\alpha(q_x^2+q_y^2)} }. \]
The inverse transform gives
\[ \boxed{ f(x,y) = \mathcal{F}^{-1} \left[ \frac{\mathcal{F}\{g(x,y)\}} {1+\alpha(q_x^2+q_y^2)} \right] }. \]
Substituting the definition of \(\alpha\) produces the single-distance PAD phase-retrieval formula:
\[ \boxed{ f(x,y) = \mathcal{F}^{-1} \left[ \frac{\mathcal{F}\{g(x,y)\}} { 1+ \displaystyle \frac{\lambda z}{4\pi} \frac{\delta}{\beta} (q_x^2+q_y^2) } \right] }. \]
Alternative Formula Using Cycles per Unit Length
Some papers use ordinary spatial frequencies \(u\) and \(v\) in \(\mathrm{cycles/m}\). The two conventions are related by
\[ q_x=2\pi u, \qquad q_y=2\pi v. \]
Hence,
\[ q_x^2+q_y^2=4\pi^2(u^2+v^2), \]
and the denominator in the angular-frequency formula becomes
\[ \boxed{ 1+ \pi\lambda z \frac{\delta}{\beta} (u^2+v^2) }. \]
The retrieval formula is therefore
\[ \boxed{ f(x,y) = \mathcal{F}^{-1} \left[ \frac{\mathcal{F}\{g(x,y)\}} { 1+ \pi\lambda z \displaystyle\frac{\delta}{\beta} (u^2+v^2) } \right] }. \]
The angular-frequency and ordinary-frequency formulas are physically identical. Their apparent difference arises only from the Fourier-transform convention and spatial-frequency definition.
Recovering the Absorption Projection
After phase retrieval,
\[ f(x,y)=\frac{I_0(x,y)}{I_{\mathrm{in}}(x,y)}. \]
From the Beer–Lambert law, \(f=e^{-p}\), where
\[ p(x,y)=\int_L\mu(x,y,\ell)\,\mathrm{d}\ell. \]
Therefore,
\[ \boxed{p(x,y)=-\ln f(x,y)}. \]
For a homogeneous material with projected thickness \(t(x,y)\),
\[ \boxed{ t(x,y)=-\frac{1}{\mu}\ln f(x,y) }, \qquad \mu=\frac{4\pi\beta}{\lambda}. \]
The phase projection can also be recovered:
\[ \boxed{ \phi(x,y) = \frac{\delta}{2\beta}\ln f(x,y) }. \]
A practical phase-retrieval package may output the retrieved exit-surface intensity \(f\), attenuation projection \(p\), projected thickness \(t\), or projected phase \(\phi\). Confirm the definition before applying an additional logarithmic transformation.
Complete Single-Distance Processing Formula
For each projection angle \(\theta\), first normalize the raw image:
\[ g_\theta(x,y) = \frac{P_\theta(x,y)-D(x,y)} {F(x,y)-D(x,y)}. \]
Retrieve the exit-surface intensity:
\[ f_\theta(x,y) = \mathcal{F}^{-1} \left[ \frac{\mathcal{F}\{g_\theta(x,y)\}} { 1+ \displaystyle \frac{\lambda z}{4\pi} \frac{\delta}{\beta} (q_x^2+q_y^2) } \right]. \]
Then calculate the attenuation projection:
\[ p_\theta(x,y)=-\ln f_\theta(x,y). \]
For a selected detector row, arrange the projections over all rotation angles into a sinogram \(p(\theta,s)\) and reconstruct the attenuation distribution using filtered backprojection or an iterative method:
\[ p(\theta,s)\longrightarrow\mu(x,y). \]
The complete processing chain is
\[ \boxed{ P_\theta \rightarrow g_\theta \rightarrow f_\theta \rightarrow p_\theta \rightarrow p(\theta,s) \rightarrow \mu(x,y) }. \]
Why the Retrieval Operator Is a Low-Pass Filter
The Fourier-domain transfer function is
\[ \boxed{ H(q_x,q_y) = \frac{1}{1+\alpha(q_x^2+q_y^2)} }. \]
At zero spatial frequency, \(H(0,0)=1\), so the direct-current component and slowly varying information are preserved. At high spatial frequency,
\[ q_x^2+q_y^2\rightarrow\infty \quad\Rightarrow\quad H(q_x,q_y)\rightarrow0, \]
so high-frequency components are suppressed.
It is important to distinguish the forward propagation operator from the phase-retrieval operator. Forward propagation is
\[ \widehat{g} = \left[ 1+\alpha(q_x^2+q_y^2) \right]\widehat{f}, \]
which enhances high spatial frequencies and generates sharp bright–dark fringes. Phase retrieval applies the inverse:
\[ \widehat{f} = \frac{\widehat{g}} {1+\alpha(q_x^2+q_y^2)}. \]
The smoothing is therefore not an arbitrary image-processing operation. It follows from inversion of the physical propagation model.
Influence of \(\lambda\), \(z\), and \(\delta/\beta\)
The retrieval strength is controlled by
\[ \alpha = \frac{\lambda z}{4\pi} \frac{\delta}{\beta}, \]
so
\[ \alpha\propto\lambda, \qquad \alpha\propto z, \qquad \alpha\propto\frac{\delta}{\beta}. \]
Propagation Distance \(z\)
Increasing \(z\) increases \(\alpha\) and strengthens the retrieval filter. Typical consequences are:
- stronger suppression of phase fringes;
- reduced image noise;
- smoother homogeneous regions;
- increased loss of high-frequency detail.
If the value used in the algorithm is smaller than the true distance, residual phase fringes may remain. If it is too large, thin walls, small pores, and fine cracks may be excessively blurred.
Ratio \(\delta/\beta\)
Increasing \(\delta/\beta\) makes the algorithm assume a stronger phase effect relative to absorption and applies a stronger correction.
If \(\delta/\beta\) is too small:
- the bright–dark edge pair may not be fully removed;
- residual halos may remain;
- reconstructed slices may contain false boundary structures.
If \(\delta/\beta\) is too large:
- noise is strongly reduced;
- boundaries become broader;
- small pores and cracks may disappear;
- spatial resolution is degraded.
The ratio therefore controls the trade-off between fringe suppression and spatial-resolution preservation.
X-Ray Wavelength \(\lambda\)
Since
\[ \lambda=\frac{hc}{E}, \]
a higher photon energy corresponds to a shorter wavelength. However,
\[ \delta=\delta(E), \qquad \beta=\beta(E), \]
so the energy dependence of phase retrieval cannot be determined from \(\lambda\) alone. The correct \(\delta(E)\) and \(\beta(E)\) values must be used for the selected X-ray energy.
Assumptions Required by the Derivation
Monochromatic or Narrow-Band X-Rays
The derivation uses a single set of parameters \(\lambda\), \(\delta\), and \(\beta\). A broad polychromatic spectrum contains energy-dependent wavelengths and refractive-index parameters, so a rigorous model requires integration over energy. Monochromatic synchrotron radiation is particularly suitable for this method.
Homogeneous Material or Constant \(\delta/\beta\)
The central assumption is
\[ \frac{\delta(\boldsymbol{r})}{\beta(\boldsymbol{r})} \approx\text{constant}. \]
This is approximately valid for an object dominated by one material. For a strongly multi-material object, one value of \(\delta/\beta\) may restore one interface correctly while over-smoothing or under-correcting another.
Paraxial Propagation
The wavefield must propagate mainly along the optical axis with small transverse angles. Large-angle scattering and strong multiple scattering are not represented.
Near-Field or First-Order Propagation
The single-distance model uses
\[ I_z \approx I_0+z\frac{\partial I}{\partial z}. \]
It is most appropriate in the near-field regime. When the propagation distance is large enough to produce complicated Fresnel fringes, the first-order TIE model may no longer be accurate.
Projection Approximation
The derivation assumes approximately straight X-ray paths and represents the object by projected \(\delta\) and \(\beta\). Strong diffraction, refraction-induced ray bending, and multiple scattering may violate this approximation.
Mathematical Meaning of PAD in the Processing Workflow
The essential PAD assumption is
\[ \boxed{ \frac{\delta}{\beta}\approx\text{constant} }. \]
It produces the phase–attenuation relation
\[ \boxed{ \phi=-\frac{\delta}{2\beta}p }. \]
Combining this relation with the TIE gives
\[ \boxed{ g=\left(1-\alpha\nabla_\perp^2\right)f }, \qquad \boxed{ \alpha = \frac{\lambda z}{4\pi} \frac{\delta}{\beta} }. \]
Fourier-domain inversion gives
\[ \boxed{ f = \mathcal{F}^{-1} \left[ \frac{\mathcal{F}\{g\}} {1+\alpha(q_x^2+q_y^2)} \right] }, \qquad \boxed{p=-\ln f}. \]
The derivation can be summarized as
\[ \boxed{ n=1-\delta+i\beta \;\Longrightarrow\; \phi=-\frac{\delta}{2\beta}p \;\xRightarrow{\mathrm{TIE}}\; g=(1-\alpha\nabla_\perp^2)f \;\xRightarrow{\mathcal{F}^{-1}}\; f,\;p }. \]
This is the mathematical basis for converting a propagation-based phase-contrast projection into an attenuation projection suitable for conventional CT reconstruction.