Isotropic scalar image visualization of vector differential image data using the inverse Riesz transform

: 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 .


Introduction
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 Xray 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) [1] 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 [2]. 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.

Differential phase imaging prerequisites
For convenience we consider the differential phase images generated by a 2-D X-ray gratingbased phase contrast imaging system [3]. The essential principal behind the 1-D Talbot system was first proposed by Momose et al [4], then refined by Pfeiffer et al [5]. 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).

Differential phase from 2-D X-ray grating Fourier side-lobes
First we define the spatial and Fourier coordinate systems, both Cartesian and polar: , , , , cos , sin , , , , cos , sin x y r x r y r The simplest system generates a moiré intensity image ( ) , f x y on the sensor with nine Fourier spectral lobes, as follows: The intensity image is real and positive. Figure 1(a) shows a simulated moiré pattern for an exaggerated, but simple object.
In practice we utilize discrete Fourier transforms (DFTs) and in particular fast Fourier transforms (FFTs), but for simplicity descriptions are in terms of continuous transforms. As the intensity in 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: Only Hermitian pairs of Fourier lobes represent real spatial functions. Single isolated lobes represent complex functions with both amplitude and phase. 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: Physically, there is the constraint that the coefficients b(x,y) are real and positive, so the phase can be extracted by subtracting the carrier phases ( ) 0 0 , u v and taking the complex logarithm or Arg function. Alternatively the side-lobes can be shifted back to the Fourier origin then inverse Fourier transformed and the Arg found: The result is a wrapped estimate of the desired differential phase where k and l are the wrapping integers. In this initial analysis we consider only phase gradients in the interval [ ] , π π − + so that phase unwrapping is unnecessary. If large gradients are expected, then a variety of phase unwrapping methods (e.g [6].) can be applied separately or jointly to x ψ and y ψ , before proceeding to the next step.

Phase gradients and Fourier correspondence
A simple Fourier relation exists for the gradient of a function ψ is defined as follows: The Fourier transform of the gradient is easily shown [7] to be the original Fourier transform multiplied by linear factors in the coordinates u and v: Given the gradient of a 2-D function it is then possible to invert this Fourier multiplier to find the original function. However, care must be taken with the singularity (division by zero) at the origin. Various methods and regularization schemes have been proposed over the years.
Notably, Frankot and Chellapa [9] solve the problem in the context of shape from shading. Ghiglia and Romero [10] proposed a related (fast cosine transform) solution to phase unwrapping via the integration of the local phase differences. More recently Arnison et al [1] exploited a compact, complex scalar solution to this vector problem. An application of quantitative wavefront reconstruction to a multilateral shearing interferogram with 9 sidelobes, similar to our Fig. 1, was described by Velghe et al [11] in 2005. Their objective is to use all the side-lobes in the Fourier-based reconstruction and thereby reduce noise.

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.

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 [12] for the mathematics, or Larkin [13] and Felsberg [14] for related image processing perspectives. We have already stated the Fourier multiplier property of the gradient operator in Eq. (9). It is not difficult to prove by differentiation of Eq. (2): Similarly the Laplacian operator . div grad Δ = has the following form: Written in more expressive shorthand we have: It looks like the Laplacian operator has a possible "square root" in the Fourier domain. The square root of the Laplacian occurs in the mathematical foundations of tomography [15]. For example, Smith & Keinert [16] define it as their Λ (lambda) operator: Faridani & Smith [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: The new vector operator R is the Riesz transform, and its Fourier multiplier has unit magnitude and odd symmetry, very like the Hilbert multiplier, and can be compactly expressed in terms of the Fourier polar angle φ : The Riesz transform has numerous useful and special properties that make it invaluable in proving important theorems in the foundations of harmonic analysis, but brevity prohibits further elaboration.

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 [13].
We define a complex gradient operator [18] D, related to the Cauchy or Wirtinger complex derivatives. We also define the spiral phase operator R as follows: The combination of Riesz transform and gradient operator gives the lambda (or frequency ramp) operator, and the combination can be achieved conveniently in complex variables, without vectors. The computational advantages of the complexified Riesz transform have been recognized recently [19] for steerable wavelet research.

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 [20] 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 socalled 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 [6],) 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 [ ] exp iβ . 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.

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 2 i π 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.

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.

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 [3]. 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. Fig. 4. One of sixteen phase-shifted moiré interferograms, f(x,y), from an experimental X-ray Talbot grating interferometer . The specimen is a plastic chess piece (a knight). The image is 1024 x 1024 pixels, pixel size 48μm. 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 [3], 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. Figures 5, 6, and 7 are available at the full 1024x1024 resolution as Media hyperlinks to allow readers detailed inspection of image quality.

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:

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 [24]. 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.

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.

Discussion
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.