X-ray Talbot moiré interferometers can now simultaneously generate two differential phase images of a specimen. The conventional approach to integrating differential phase is unstable and often leads to images with loss of visible detail. We propose a new reconstruction method based on the inverse Riesz transform. The Riesz approach is stable and the final image retains visibility of high resolution detail without directional bias. The outline Riesz theory is developed and an experimentally acquired X-ray differential phase data set is presented for qualitative visual appraisal. The inverse Riesz phase image is compared with two alternatives: the integrated (quantitative) phase and the modulus of the gradient of the phase. The inverse Riesz transform has the computational advantages of a unitary linear operator, and is implemented directly as a complex multiplication in the Fourier domain also known as the spiral phase transform.
© 2014 Optical Society of America
Recent advances in X-ray imaging allow the capture of differential phase images of the test specimen. Differential phase provides information about X-ray deflection rather that the X-ray absorption of conventional radiography. Clinical differential phase provides fine detail about the soft-tissue structure, detail simply not present in normal X-ray images. However, there is a catch: differential phase gives a highly asymmetric and anisotropic depiction of the optical path (or integrated refractive index) through a specimen. Many systems (using linear gratings) provide just one oriented differential phase image - the x derivative for example. More recent systems have been developed with crossed-gratings (or 2-D arrays) that simultaneously provide two differential phase images; the x and the y derivatives for example.
With two differential phase images there is a question of how best to view the information. The conventional approach is to integrate (literally and mathematically) the two images into one. Simple and efficient computational methods have already been devised for integrating optical differential phase images (such as Nomarski Differential Interference Contrast or DIC)  utilizing fast Fourier correspondences. Although integration has the important advantage of producing an intensity image that is quantitatively related to the original optical path length there are a number of disadvantages:
- • integration may encounter inconsistent gradient data between two components, resulting in damaging image artifacts,
- • integration typically increases dynamic range beyond the range of LCD displays,
- • integration introduces blur due to the intrinsic low-pass nature of 2-D integrator.
With just one differential phase image there is the related problem of how best to view the image. Again (one dimensional) integration is the conventional approach. However, it has been proposed that Hilbert transformation is a better approach . Hilbert transformation has two advantages. Firstly the streaking problem due to the unknown constants of integration is circumvented. Secondly, the fine detail and dynamic range are both maintained whilst the odd Fourier symmetry of the derivative is subverted by the odd symmetry of the Hilbert Fourier multiplier.
In this paper we propose extending the 1-D Hilbert symmetrizing method formally for 2-D images, using the Riesz transformation. Only the essential underlying mathematical theory will be presented in this preliminary exposition. The paper is structured as follows:
- Section 2 outlines the Fourier theory of differential imaging.
- Section 3 introduces the Riesz transform and the related complex embedding (known as the spiral phase transform)
- Section 4 outlines the main algorithmic definition of the process.
- Section 5 presents inverse Riesz visualization in comparison with conventional images.
- Section 6 concludes with a discussion of the main results.
2. Differential phase imaging prerequisites
For convenience we consider the differential phase images generated by a 2-D X-ray grating-based phase contrast imaging system . The essential principal behind the 1-D Talbot system was first proposed by Momose et al , then refined by Pfeiffer et al . The 2-D system multiplexes the different directional phase derivatives onto different Fourier spectral side-lobes. Alternatively, we could just assume that the two requisite differential phase images are provided from a black-box system (such as two separate, orthogonal, Nomarski DIC images).
2.1 Differential phase from 2-D X-ray grating Fourier side-lobes
First we define the spatial and Fourier coordinate systems, both Cartesian and polar:Figure 1(a) shows a simulated moiré pattern for an exaggerated, but simple object. The nine separate lobes are clear.
We denote the 2-D continuous Fourier transform (FT) as a double headed arrow and define it via the forward and inverse transforms:Eq. (2) is real, the FT of the moiré image in Eq. (2) is Hermitian, namely the sum of nine triple convolutions which correspond to four Hermitian pairs and a DC lobe:
The individual Fourier lobes are defined:Figure 1(b) shows the magnitude (in a log intensity scale) of the FT in Eq. (3). In practice, two orthogonal side-lobes P are selected and isolated, inverse Fourier transformed to obtain the X-ray phase and magnitude (absorption) information. The easiest example being the (1,0) and the (0,1) side-lobes:6].) can be applied separately or jointly to and , before proceeding to the next step.
2.2 Phase gradients and Fourier correspondence
A simple Fourier relation exists for the gradient of a function is defined as follows:7] to be the original Fourier transform multiplied by linear factors in the coordinates u and v:8] applied the Fourier method to wave-front reconstruction. Notably, Frankot and Chellapa  solve the problem in the context of shape from shading. Ghiglia and Romero  proposed a related (fast cosine transform) solution to phase unwrapping via the integration of the local phase differences. More recently Arnison et al  exploited a compact, complex scalar solution to this vector problem. An application of quantitative wavefront reconstruction to a multilateral shearing interferogram with 9 side-lobes, similar to our Fig. 1, was described by Velghe et al  in 2005. Their objective is to use all the side-lobes in the Fourier-based reconstruction and thereby reduce noise.
3. Riesz transform
The Riesz transform is widely accepted in the mathematics community (harmonic theory in particular) as the proper, isotropic generalization of the 1-D Hilbert transform to higher dimensions (or Euclidean spaces). In engineering, especially signal and image processing, anisotropic approximations like the orthant Hilbert transform remain popular, and the Riesz transform is little used.
3.1 Origin of the Riesz transform
The concept Riesz transform arises naturally in the analysis of harmonic functions, the Laplace and Poisson equations, and the solution of partial differential equations. The Riesz transform can be interpreted as a partial factorization of certain differential operators, namely the Laplacian and the gradient. Readers interested in some background on the Riesz transform should consult Stein  for the mathematics, or Larkin  and Felsberg  for related image processing perspectives.15]. For example, Smith & Keinert  define it as their (lambda) operator:17] then use the lambda operator defined above to generate a locally soluble tomographic inversion based on the known spatial localization of the Laplacian operator. In contrast our proposed application maintains the spatial non-locality of the Riesz operators. We factorize the gradient operator, applied to the phase in this instance, as follows:
3.2 Embedding of the vector Riesz transform in the complex scalar domain
It often happens that the Riesz transform needs to be implemented in a pre-existing image processing software application. Many applications work with gray-scale images and allow complex representation for compatibility with discrete Fourier transform (DFT) processing. Surprisingly, Riesz transform processing can be implemented in 2-D without the further extension to complex vector representation. A little more care has to be taken to account for cross-talk between somewhat overloaded channels, but in most cases everything just works out automatically. The embedded formalism is known as the complexified Riesz transform or spiral phase transform .
We define a complex gradient operator  D, related to the Cauchy or Wirtinger complex derivatives. We also define the spiral phase operator R as follows:19] for steerable wavelet research.
3.3 Discrepancy term of the inverse Riesz transform
Considering Eq. (19) it would seem that application of the properly scaled inverse Riesz transform results in a real image with lambda enhancement. However we have not considered random noise or systematic effects. In general each of the differential images will be a pure differential image with some unknown additive component. From our work  in general decompositions of 2-D phase fields and 2-D phase unwrapping we recognized immediately that there is a part of the phase that is not representable as the gradient of a scalar. This so-called phase “discrepancy” can be represented as the curl of a vector field, and leads inevitably from the Helmholtz decomposition (viz the Fundamental Theorem of Vector Calculus). Ghiglia and Pritt (chapter 2 ,) describe in detail the conditions necessary for a phase function in 2-D to generate a curl discrepancy. Even for non-phase differentials, noise can still induce apparent discrepancies.
One important systematic error can easily occur if the direction of the derivative operation is rotated. For example, if the x differential is really along a direction rotated an angle counterclockwise from the x-axis, and the y differential an angle counterclockwise from the y-axis, then the output of the inverse Riesz defined in Eq. (19) is simply multiplied by a complex constant . The end result is that if the x and y differentials are interchanged, for example, then the desired output switches from the real part to the imaginary part. So it is important to make sure that the derivative direction is known beforehand. In practice it is possible to automate algorithms to compensate this effect under normal circumstances (when the noise level is not excessive).
Note that the rotation of derivative directions gives the same effect as introducing a curl (discrepancy) component, and this fact ties in with the observation that the curl effect can be interpreted as 90º rotated gradient components.
More generally, the discrepancy term is a measure of how reliable the gradient data are. The conventional approach to integration by solving the Poisson equation is completely blind to the discrepancy data. Discrepancy is discussed further in section 5.4.
4. Algorithmic definition inverse Riesz transform applied to differential images
As outlined above, the complexified Riesz transform, and inverse Riesz transform are simple to implement in an image processing system that allows complex variables. Figure 2 shows the basic flow chart.
Two orthogonal differential images are input into the processor. The x-differential image is placed in the real part, and the y-differential in the imaginary part. The fast Fourier transform of the complexified gradient image is computed and the result is then multiplied by the inverse Riesz spiral phase factor. Complex constants are used to re-normalize the gradient Fourier multiplier. However this step may be omitted as it only multiplies the final image by a constant factor and swaps the real and imaginary parts. The inverse FFT (IFFT) produces an output with the desired image in the real part and a “discrepancy” image in the imaginary channel.
5. Application to differential images
Although our proposed method is applicable to any differential image pairs, whether phase or otherwise, we have chosen an example from X-ray Talbot moiré differential phase imaging.
5.1 Example from X-ray Talbot moiré interferometry
Our example is a molded plastic chess piece shown in Fig. 3. The chess-piece was imaged in an experimental X-Ray Talbot moiré interferometer . Some plastics have low X-ray absorption yet significant refractive index not unlike some animal soft tissues. A complete sequence of 16 (4x4) phase shifted interferograms was collected. Each phase shift is 90º in a predominantly horizontal or vertical direction. One frame of the sequence is shown in Fig. 4.
Two Fourier side-lobes were isolated using 16 step 2-D phase shifting algorithms (PSAs). The side-lobes could easily have been isolated using the windowed Fourier transform method , or the traditional Fourier transform method [22,23], but we prefer to illustrate our method with the highest isolation, resolution and noise reduction obtained from a PSA.
Figure 5(a) shows the corresponding Fourier transform magnitude, the extracted amplitude modulation (or absorption) image, and the extracted x and y differential phase images. The absorption image corresponds to the conventional X-ray image. The variation in background intensity of the absorption in Fig. 5(b) is due to the interferometer source geometry , but has no effect on the phase estimation.
The x and y differentials in Fig. 6(a) and 6(b) have the well-known bas-relief asymmetry of optical differential phase images. On their own the differential images are highly anisotropic, losing considerable detail at certain orientations. Together the images are difficult to interpret; requiring continual visual scanning and comparison.
5.2 Visual assessment
In Fig. 7(a) we see that integration produces low frequency background variation of intensity (mainly due to our imposing overly simplistic periodic boundary conditions). Furthermore the integrated image has a large dynamic range (hard to view properly on a monitor or in print). There is also a low-pass filtering effect apparent when compared to the original differentials.
The gradient magnitude images of Fig. 7(b) and 7(c) seem to maintain the required detail and isotropy with the negated magnitude of Fig. 7(b) suitable for print out and Fig. 7(c) perhaps more suitable for monitor display. What is not so obvious in this example is the potential aliasing that the modulus operation generates. It is well known in signal processing that modulus-squaring doubles the bandwidth required by an image without generating any more useful information (quite the contrary in fact). Furthermore the modulus operation generates even higher order aliasing than mod-squared. The absorption image and the integrated image both have the advantage that the (demodulation) processing is linear (in terms of the 16 input moiré images); and so no spurious (aliased) frequencies are generated.
The inverse Riesz transformed image in Fig. 7(d) we see an image with a wealth of fine detail. The image also retains some of the global intensity ordering of the idealized integral image in Fig. 7(a). Dark regions in the inverse Riesz image broadly correspond with dark regions in the integrated image – so intensity ordering is approximately maintained. Although the gradient image of Fig. 7(b) has fine detail it has completely lost the intensity ordering. For example all the air bubbles are rendered as brightly as their surrounding regions. The inverse Riesz image of Fig. 7(d) shows the classic edge enhancement common to optical phase contrast imaging systems. In this instance the edge effect results from two processes, the first being the differential phases from Talbot distance propagation in the interferometer, the second from the isotropic melding of the inverse Riesz transformer. However, the overshoot/undershoot is purely a mathematical characteristic of the lambda operator. Although we do not currently have a quantitative criterion for comparing the relevant information content of Fig. 7 (a), Fig. 7(b) and Fig. 7(d) we think that our initial visual comparison will be compelling for many readers. The main properties can be summarized in a Table 1:
5.3 Optical (Craik-O'Brien-Cornsweet) illusion
Careful measurement of local patch intensities shows that the intensity ordering in the Riesz image is to some extent an optical illusion, albeit serendipitous. For example the internal air bubble near the middle of the bottom edge appears dark in both Fig. 7(a) and 7(d). One of the reasons the Riesz processed image resembles the integrated image is that the sharp overshooting edges produce apparently darker or lighter regions via a human perceptual effect, known as the Craik-O'Brien-Cornsweet illusion . The quantitatively darker regions in the integral image coincide with the apparently darker regions in the Riesz image, likewise for the lighter regions. Adjacent regions with the same intensity will be perceived to be different if the connecting edges have an overshooting and undershooting intensity profile.
5.4 Discrepancy image
As mentioned in section 3.3 the gradient of a phase image has two components related to the Helmholtz decomposition of a vector field. In this initial exposition of inverse Riesz visualization we had not planned to dwell on what is essentially a side-product of the main method. Figure 8 is the inverse Riesz discrepancy image related to Fig. 7(d). The image is shown in linear gray-scale with zero represented by mid-gray. Indeed most of the image appears mid-gray, with significant positive and negative deviations mostly occurring at the specimen edges. The mid-gray actually contains wide-band noise. In terms of signal energy the discrepancy image has a variance about one third of primary inverse Riesz image, in this instance. Closer scrutiny reveals that the input images (Fig. 6(a) and Fig. 6(b)) contain phase wrapping errors on the more extreme edges. These wrap errors fail to satisfy the underlying gradient assumption and inevitably propagate through to the (curl) discrepancy image. Consequently the discrepancy image can be used to check whether or not the input images are consistent with some underlying assumptions. It may even be possible to go further and use the discrepancy image to improve the reliability or accuracy of the main inverse Riesz image (e.g. noise reduction), but we have not pursued this line of research further.
The critical reader might suggest that our proposed method is merely equivalent to proceeding with conventional integration, then applying a high-pass (i.e. lambda or ramp) filter. Whilst this is almost true, our method completely avoids the most difficult and unstable step of integration. The inverse Riesz transform is a stable, unitary linear operation. Conventional (Poisson equation) integration also omits the informative discrepancy image. The Riesz transform can be viewed, then, as a simple short-cut to a rather desirable result.
We would like to thank our Canon Inc. colleagues for providing access to an X-ray Talbot moiré interferometer , and experimental data summarized by Fig. 4. In particular, we thank Kentaro Nagai, Soichiro Handa, Takashi Nakamura, Takeshi Kondoh, Yutaka Setomoto, Takayuki Teshima, and Hidenosuke Itoh. Don Bone provided the chess piece appearing in Figs. 3-8 and collected the corresponding interferogram data.
References and links
1. M. R. Arnison, K. G. Larkin, C. J. R. Sheppard, N. I. Smith, and C. J. Cogswell, “Linear phase imaging using differential interference contrast microscopy,” J. Microsc. 214(1), 7–12 (2004). [CrossRef] [PubMed]
2. M. R. Arnison, C. J. Cogswell, N. I. Smith, P. W. Fekete, and K. G. Larkin, “Using the Hilbert transform for 3D visualization of differential interference contrast microscope images,” J. Microsc. 199(1), 79–84 (2000). [CrossRef] [PubMed]
3. H. Itoh, K. Nagai, G. Sato, K. Yamaguchi, T. Nakamura, T. Kondoh, C. Ouchi, T. Teshima, Y. Setomoto, and T. Den, “Two-dimensional grating-based X-ray phase-contrast imaging using Fourier transform phase retrieval,” Opt. Express 19(4), 3339–3346 (2011). [CrossRef] [PubMed]
4. A. Momose, S. Kawamoto, I. Koyama, Y. Hamaishi, K. Takai, and Y. Suzuki, “Demonstration of X-ray Talbot interferometry,” Jpn. J. Appl. Phys., Part 2 , 42, 866–868 (2003).
5. F. Pfeiffer, T. Weitkamp, O. Bunk, and C. David, “Phase retrieval and differential phase-contrast imaging with low-brilliance X-ray sources,” Nat. Phys. 2(4), 258–261 (2006). [CrossRef]
6. D. C. Ghiglia and M. D. Pritt, Two-dimensional phase unwrapping (John Wiley and Sons, New York, 1998).
7. R. N. Bracewell, The Fourier Transform and its Applications, McGraw-Hill Electrical and Electronic Engineering Series (McGraw Hill, New York, 1978).
8. K. R. Freischlad and C. L. Koliopoulos, “Modal estimation of a wave front from difference measurements using the discrete Fourier transform,” J. Opt. Soc. Am. A 3(11), 1852–1861 (1986). [CrossRef]
9. R. T. Frankot and R. Chellappa, “A method for enforcing integrability in shape from shading algorithms,” IEEE Trans. Pattern Anal. Mach. Intell. 10(4), 439–451 (1988). [CrossRef]
10. D. C. Ghiglia and L. A. Romero, “Robust Two-Dimensional Weighted and Unweighted Phase Unwrapping That Uses Fast Transforms and Iterative Methods,” J. Opt. Soc. Am. A 11(1), 107–117 (1994). [CrossRef]
11. S. Velghe, J. Primot, N. Guérineau, M. Cohen, and B. Wattellier, “Wave-front reconstruction from multidirectional phase derivatives generated by multilateral shearing interferometers,” Opt. Lett. 30(3), 245–247 (2005). [CrossRef] [PubMed]
12. E. M. Stein, Singular Integrals and Differentiability Properties of Functions, Princeton mathematical series (Princeton University Press, Princeton, N.J., 1970).
13. K. G. Larkin, D. J. Bone, and M. A. Oldfield, “Natural demodulation of two-dimensional fringe patterns. I. General background of the spiral phase quadrature transform,” J. Opt. Soc. Am. A 18(8), 1862–1870 (2001). [CrossRef] [PubMed]
14. M. Felsberg and G. Sommer, “The Monogenic Signal,” IEEE Trans. Signal Process. 49(12), 3136–3144 (2001). [CrossRef]
15. R. S. Strichartz, “Radon inversion - variations on a theme,” Am. Math. Mon. 89(6), 377–384 (1982). [CrossRef]
17. A. Faridani, E. L. Ritman, and K. T. Smith, “Local tomography,” SIAM J. Appl. Math. 52(2), 459–484 (1992). [CrossRef]
19. M. Unser and N. Chenouard, “A Unifying Parametric Framework for 2D Steerable Wavelet Transforms,” SIAM J. Imaging Sciences 6(1), 102–135 (2013). [CrossRef]
21. K. Nagai, H. Itoh, G. Sato, T. Nakamura, K. Yamaguchi, T. Kondoh, and S. H. T. Den, “New phase retrieval method for single-shot x-ray Talbot imaging using windowed Fourier transform,” Optical Modeling and Performance Predictions V, Proc. SPIE 8127, San Diego, California, August 21, 2011.
23. M. Takeda, H. Ina, and S. Kobayashi, “Fourier-transform method of fringe-pattern analysis for computer-based topography and interferometry,” J. Opt. Soc. Am. 72(1), 156–160 (1982). [CrossRef]
24. M. W. Levine and J. M. Shefner, Fundamentals of Sensation and Perception (Pacific Grove, CA: Brooks/Cole, 1991) p. 675.