1. Introduction

Velocimetric techniques range from calculation of fluid speed via temperature-dependant resistance in a wire anemometer, to measurement of solar wind velocities in a passing comet via spectral-shape analysis [1]. Most comprehensive flow measurements however, are achieved using optical techniques based on imaging particles deliberately seeded within the flow. Often, a plane within the flow vessel is illuminated by a visible-light source (commonly a laser), focused into a light “sheet” with detectors arranged to capture images of light scattered from the particles. Seed particles are chosen for their neutral buoyancy and reflectivity. If the seeding density is low enough for individual particles to be identified from one image frame to the next, the technique is known as particle tracking velocimetry (PTV); for higher seeding densities one may calculate displacements via statistical correlation analysis of image frames, known as particle image velocimetry (PIV) [2]. Typically PIV yields velocity information with two vector components within the light sheet and over the two dimensional region of the sheet. Three components of velocity can be measured – but only within the two dimensional light sheet – by imaging from two or more locations with a technique known as stereo PIV [3]. Velocity information can be measured over a three dimensional volume using a number of PIV and PTV techniques [4]. These techniques vary in their ability to measure biological flows, with many limited by particle seeding density and high levels of complexity. A new technique known as volumetric correlation PIV [4] offers significant advantages over other volume techniques, but as with all visible-light PIV methods, shares the significant limitation of being unable to work with optically opaque flows.

In recent years X-rays have been shown to be a viable alternative to visible-light velocimetry, with the potential for true three-dimensional measurements, used with the sample and detectors in transmission geometry to simultaneously obtain flow information over the entire thickness of the vessel [5,6]. These X-ray-PIV/PTV methods have the advantage of being able to probe optically opaque samples. This has particular promise for biomedical applications which require high-resolution in-vivo measurements of movement and function [7], for example in the application to be investigated in this paper, namely that of blood flow measurement.

To motivate such an investigation, we note that diseases of the vascular system continue to constitute the world’s greatest cause of mortality and morbidity. A sometimes critical parameter in their diagnosis, monitoring and treatment evaluation is blood flow rate and velocity information. It is known that the formation of atherosclerosis is linked to the shear stress (proportional to the gradient in the flow velocity at the vessel wall). Extensive studies of the most effective methods for reliable, instantaneous blood-flow measurement at minimal levels of invasiveness have been carried out, concluding that none of the techniques currently in use meet all desired clinical requirements [8].

The potential of X-ray PIV for blood-flow measurements was discussed by Lee and Kim in 2005 [9]. Whilst in general X-ray PIV measurements may be achieved purely through absorption contrast given sufficient mass density difference between the tracer particles and surrounding medium, frequently a form of phase contrast may be used, for example utilizing the refractive index difference between the flow medium and bubbles of air within. In the case of blood, its composition may be exploited, with regards to the very high numbers of red blood cells within which allow velocimetric measurement without the addition of any tracer particles. The typical hematocrit, or red cell volume density, for an adult human male is approximately 45%; their small size (approx 2μm in one dimension, 6μm in the other) means that to X-rays they act in concert like a weak aberrated compound lens system [10]. Upon propagation of the phase-varying exit X-ray wavefield, self-interference results in a characteristic dynamic speckle pattern. The speckle pattern evolves via Fresnel diffraction with propagation distance and depends on the sample thickness and detector point spread function (PSF). Generally speaking, the dominant spatial frequency of the intensity speckles corresponds to a few µm for X-ray energies on the order of 20 keV, and optimum visibility occurs for propagation distances on the order of tens of centimeters. Figure 1 shows an example, with video illustrating the change of speckle pattern of flowing blood between two close points in time (14ms) at 45cm propagation distance, where the average displacement of blood cells between the frames is roughly equivalent to one cell length (6μm). More details of the imaging and flow conditions are described in Section 2.

 

Fig. 1 Example of dynamic X-ray phase-contrast blood speckle, obtained when a blood-filled vessel approximately 2mm thick (inner diameter) is illuminated by 25 keV X-rays, with the resulting exit wavefield propagating 45cm to the detector plane. Field of view is 128×128 pixels = 115×115 μm (0.9µm pixel size), with 14ms between exposures. 
(Media 1).

Download Full Size | PDF

We have shown previously [11] that preprocessing of such speckle images, in particular the use of single-image phase retrieval algorithms [12] which “decode” the propagation-induced X-ray speckle to yield the projected blood-cell thickness maps which resulted in such speckle, leads to improved correlation measurements in subsequent PIV analyses. Three-dimensional velocimetric analysis begins with the understanding that the resulting optimized correlation peaks, for a given “line of sight” through the vessel, is a projection representing the probability density function (PDF) of particle displacement for all contributing particles across the full vessel thickness [13]. For the simplest possible axi-symmetric geometry, i.e. a cylindrical vessel with rotationally-symmetric flow, the resulting flow vectors are single-component (with flow in the axial direction only); in this case, one may use a method of scalar tomography to reconstruct the radial velocity profile via the inverse Abel transform [11].

Here, we expand the theory of X-ray phase-retrieval PIV flow reconstruction to include any axially-symmetric flow geometry, whose constituent vectors will be multi-component, resulting in the graduating step from scalar to vector velocimetric tomography. To achieve this, it has been necessary to develop a variant of the Abel transform which takes into account the inability to detect motion that occurs along the axis of propagation (see Sec. 3).

Experimentally, we consider a sample modeling a simple symmetric blood-vessel constriction (as seen in the presence of stenoses), and pump a periodic pulsatile flow so as to approximately model a beating heart. We reconstruct the full velocity profile across our in vitro model in all four dimensions. Our simple reconstruction technique allows for close to instantaneous measurements (at best needing only 2 consecutive frames, for a single orientation of the flow vessel, for tomographic velocimetry analysis). Our technique complements new computed tomographic X-ray velocimetry techniques [7,14], which are not limited by any assumptions of rotational symmetry, but require rotation of either sample or detector and thus place greater restrictions on time-resolved measurement.

2. Experiment/acquisition of blood speckle images

Experimental images of in-vitro blood flow were acquired at the undulator medical imaging beamline BL20XU at the SPring-8 synchrotron in Japan. We utilized the hutch 80m downstream of the source, for maximum flux at 25keV X-ray energy (using a Si-111 double crystal monochromator). The imaging setup is shown in Fig. 2 . The mechanism of propagation-based phase contrast [15] renders the transverse phase shifts, imparted upon the nominally planar X-ray waves as they traverse the blood-filled tube, visible as a flowing speckle pattern over the surface of the detector. In this context, the distance between the exit surface of the sample and the detector, which in the present experiment was 45 cm, must be sufficiently large for phase contrast to be observed.

 

Fig. 2 Schematic of X-ray phase contrast blood flow velocimetry setup. Monochromated shutter-synchronized X-ray synchrotron radiation illuminates a rotationally-symmetric vessel (in this case a 3D-printed acrylic model with constriction), through which blood is flowing in a pulsatile manner. The blood cells act as weak X-ray lenses, forming a flowing speckle pattern over the detector plane. X-ray velocimetric phase-retrieval tomography seeks to reconstruct the fully-three-dimensional time-dependent rotationally-symmetric flow velocity field within the vessel, given a time series of dynamic speckle patterns.

Download Full Size | PDF

For the detector, a PCO 4000 CCD was coupled to a beam monitor consisting of a phosphor and a 20× magnification lens (2×2 binning gave effective pixel size = 0.9μm, with power-spectrum analysis yielding a cut-off spatial frequency corresponding to 3.4 μm). The camera shuttering system was synchronised to a high-speed shutter via a Uniblitz VMM-T1 timer with exposures of 3.4ms and an inter-frame rate of 500ms. Whole blood was pumped through the sample using a micro-syringe pump (WPI Inc.) at an average rate of 5 μLs⁻¹.

The sample was a simple rotationally-symmetric stenotic model, manufactured by 3D printing (Objet Geometries Ltd). This may be represented as a cylindrical vessel of diameter 2.7mm, with axially-symmetric circular constriction giving a minimum diameter of 1.9mm (a cross-sectional area ratio of 2:1). As a result of this symmetry, measurements made over more than half of the vessel diameter are redundant. In order to maximize the spatial resolution, the lens magnification was chosen such that the field of view occupied just over half of the maximum vessel diameter.

The pump cycle is represented by a full-wave rectified sine wave, with sine-function period of 2.1s, as shown in Fig. 3 . Fifty image pairs were recorded in a total of 25 seconds (with our camera in “PIV” mode two consecutive images are recorded in quick succession, then after a longer period another pair is acquired, etc.). Whilst it is possible to achieve measurements from a single image pair at each desired time point, for the purposes of optimal signal-to-noise we employed phase-averaging of the pump cycle. To determine the phase of each measurement time-point, conventional PIV analysis of the images was performed using software described in reference [3], followed by spatial averaging over the central region to give a simple time trace [see Fig. 3(a)]. This time trace was then decomposed to give 8 representative points on the pulse cycle, as shown in Fig. 3b.

 

Fig. 3 The pulsatile flow cycle of the blood pumped through the vessel is described by (a) a full-wave rectified sine wave function, recorded over 50 image pairs, where each point represents a single pair of images. The phase-averaging (8 sections per period, i.e. A-H) of these pairs is shown in (b). Measurements here are of the axial flow component, performed using conventional PIV analysis of the images with spatial averaging over the central region of the vessel (thus indicating the modal values for that region).

Download Full Size | PDF

3. Theory of vector tomographic reconstruction

For a vessel exhibiting rotational symmetry about its main axis, one may decompose the resulting flow vector field in terms of its radial and axial velocity components (see 
Fig. 4 ):

 

Fig. 4 Geometry for velocimetric vector tomography. (a) Indicative cross-correlation map, for the z-projection of the optical thickness through the vessel, taken at two closely-spaced times. The characteristic cross-correlation peak has both horizontal (radial) and vertical (axial) components as shown. These components may be considered separately, in a tomographic velocimetry analysis. (b, c) The vertical velocity component (shown in blue) of any particle displacement is fully discernable by the X-rays which have travelled in the out-of-page (z) direction. However, a sub-component of the radial component (shown in yellow) is in the z direction and thus only the perpendicular sub-component (x-component, shown in red) contributes to the measurements.

Download Full Size | PDF

Here, a caret denotes a unit vector, with the geometry as shown in Fig. 4b. The direction of propagation of the X-rays in this case is perpendicular to the vessel axis (i.e., in the z direction). These two vector components of the flow’s velocity field may then be considered separately in terms of both their measurement and reconstruction.

Note that in this section we will speak of vector components and their “projection” over the direction z of X-ray propagation, meaning the integral of a given vector component over z. Note that the vector components (together with their z projections) may be negative as well as positive. In such projections, any particle displacements in a positive z direction will be cancelled out in measurements by a displacement of equal magnitude in the negative z direction.

3.1 Axial velocity-field component

The axial component of flow vectors (taken to be along the y axis, in Cartesian xyz coordinates, cf. Figure 4b) means that the projection of said vectors in the z direction is equivalent to the standard forward scalar Abel transform:

If the projection V(x,y) is known, then v(r,y) may be calculated using Eq. (4). In addition there are a number of alternatives to this equation, for example:

Since this alternative expression features the derivative outside of the integral it has the advantage of being more computationally stable (integrating first has the effect of smoothing noise whereas differentiating will amplify it).

3.2 Radial velocity-field component

The net radial component integrated over the volume of the vessel is, by the object’s symmetry, zero. However, a projection of the radial component, along the direction of the X-ray’s propagation, is measurable – see Fig. 4. Since any particle displacement which occurs along the propagation axis is undetectable in projection, only the perpendicular component contributes:

The new forward transform is therefore very closely related to Eq. (3):

The corresponding inverse transform may be computed by manipulating the transform pairs in sections 8.10.2 to 8.10.3 of the book by Deans et al. [16], resulting in:

4. Phase-retrieval pre-processing of X-ray dynamic speckle images

As described by Irvine et al. [11], spatial filtering of the blood speckle images is desirable prior to performing a quantitative velocimetric tomography analysis. After average subtraction of the images in the sequence (which removes all stationary remnants such as some of the beam structure, dust, vessel walls etc.) the remaining image information is ideally “pure” speckle, with a characteristic Gaussian grayscale histogram and toroidally-shaped power spectrum of frequencies. One then applies a simple single-image phase retrieval algorithm [12,15] to decode the propagation-based phase contrast information in each frame of the dynamic X-ray speckle, taking each such frame and computing the projected thickness of the particles (in this case, red blood cells) in the vessel at each time step. The algorithm [11] is

For our experimental data, the full beam width was not much bigger than the vessel, and the associated illuminating-beam intensity variations in time manifested as streaks and striations which were not all removable with flat-field correction or in the average-subtraction. The phosphor and monochromator also contributed to these streak effects. The result of such artifacts is strong signals in the power spectrum, predominantly at low spatial frequencies, some of which share the same Fourier-space region as the speckle. Application of the phase retrieval in Eq. (8) to these images produced amplification of the artifact, visible as cloudiness which unacceptably occluded the image correlation peaks. To obviate this, we employed Fourier mask filtering, to pre-process the speckle data. As shown in Fig. 5 , Fourier-space streaks were removed with smooth directional masks with angle and position determined through inspection of the power spectrum. A smooth low-frequency mask, in the vicinity of the Fourier-space origin, was also employed. For samples with negligible absorption, the transport-of-intensity equation implies that propagation of the wavefunction for low spatial frequencies is approximated in Fourier space via multiplication by a factor proportional to $k_{\perp }^2$ [12,15], therefore a smooth parabolic profile was used for this filtering of low Fourier frequencies. The maximum size of the mask was determined through simulations of blood speckle patterns [18]. These simulations model the red blood cells as spherical (the “average” shape of a rotating cell), packed in representative volume concentrations from which the projected thickness is calculated and the projection approximation [15] applied. An angular-spectrum formalism [15,19] is used to simulate free-space propagation of the exit wavefield to the detector plane. See Kitchen et al. [20] for an in-depth discussion of a similar model.

 

Fig. 5 The magnitude of the Fourier transform of the speckle images showing the mask (in dark) used to suppress directional (horizontal, vertical and diagonal) artefacts, as well as low-frequency “noise”. The toroidal shape is the typical spread of frequencies comprising image speckle.

Download Full Size | PDF

Given the simultaneous occupation of certain regions of the spatial frequency spectrum by both some of the speckle and artifact, after the latter’s removal via masking it was not possible to completely remove the Fresnel-diffraction-induced dark rings surrounding the correlation peaks. In other words, a small missing portion of the lower frequencies resulted in some slight minima around the correlation peaks.

In Fig. 6 a typical example of cross-correlation peaks, both before and after the application of phase retrieval to the speckled images, is shown. The optimal phase-retrieval parameters, tuned in order to achieve minimal dark rings in the correlation peaks, used ε = 100 μ. Note that the optimization of the two balancing filters (one is low-pass, one is high-pass) was found to be fairly robust, in that small changes resulted in insignificant variations to the final velocimetric analysis.

 

Fig. 6 Example cross correlation peaks, before (top left) and after phase retrieval (top right). Cross-correlation windows have been cropped from 256 pixels to 100 for clarity. The black spots close to the base of each peak denote the cross-correlation center or origin corresponding to zero displacement. Cross-sectional slices of both peaks are shown on the line plot (bottom center) for direct comparison. This analysis strongly suppresses the diffraction induced dark rings surrounding the correlation peaks.

Download Full Size | PDF

5. Analysis for obtaining the forward transform

In most conventional visible-light PIV analyses, consecutive images of flow particles are subdivided into “interrogation windows”; a pair of images results in many pairs of such windows. Statistical correlation analysis of the interrogation windows results in correlation peaks, one for each window, whose coordinates represent the 2-D vector of the flow (cf. Figure 4a). However when the sample is illuminated with X-rays the corresponding windows cannot be said to contain a single representative flow vector as the information contained within is really the projection of variable flow, along a “line of sight” spanning the entire sample thickness. Cross-correlation analysis then yields a stretched peak representing the probability density function of particle displacement (like a histogram) convolved with the autocorrelation function whose width is a function of the particle size. The position of the cross-correlation maximum represents merely the coordinates of the modal flow vector [13]. The mean displacement for each window may be found via normalized integration of the cross-correlation function (as the centroid). For the sake of simplicity, this is calculated from the cross-correlation peak profile. One may then multiply by the sample thickness T(x) to give

This gives the projection V(x) of velocity components, i.e. the Radon transform of detectable (normal to the X-rays) velocities. For rotationally-symmetric cases, this reduces to the (forward) Abel transform, as discussed previously.

Note that for flows with both radial and axial components, the profile of the cross-correlation peaks will lie along a line whose angle depends on the position of the window relative to the vessel, according to the relative strength of the axial and radial components of the flow at that point. For the axial component in our representative geometry, the displacement parameter s represents the y coordinates of the points along the profile, whilst the x coordinates of the profile relate to the radial component.

6. 4D Velocimetry Results

Once the axial and radial components have been reconstructed separately as scalar functions of radius (using the inverse Abel transforms described by Eqs. (4) or (4a), and (7) respectively), they may be recombined in vector form. The process is repeated at each row of windows across the vessel image. Figure 7 shows the results of the reconstruction both in scalar and vector form for the row 1mm above the throat of the stenosis. The axial component shows aspects of Womersley flow, which results from pulsatility of a sufficiently high frequency in a viscous material [21]. The typical Womersley flow is known to have a much flatter profile than the standard parabolic profile of steady flows. A slight increase in axial velocity in regions away from the central axis can even be seen. The radial velocity profile also agrees well with the expected results. There is a positive (outward) flow caused by the enlarging geometry found downstream of the stenosis throat. This increases in a nearly linear fashion away from the central axis, where continuity dictates zero axial flow, and again decreases towards zero near to the wall where both components of the flow must be zero.

 

Fig. 7 Reconstructed velocity profiles, for the radial (top-left), axial (top-right) components at 1.0 mm above the throat of the stenosis which has a radius of 0.95mm. A reduced susbset of the same data is shown in a combined vector plot (bottom). Measurements are from point F on the pump cycle (cf. Figure 3b).

Download Full Size | PDF

Built up in this way, a 2D vector map may be compiled which resembles a standard PIV map except that with r instead of x as the independent variable it is not simply flow across a plane, instead it contains all the 3-dimensional flow information of the rotationally symmetric flow-field. A representative example of this map at point F on the pump cycle (cf. Figure 3b) is shown below in Fig. 8 .

 

Fig. 8 Reconstructed 2D vector velocity map, obtained using X-ray phase retrieval tomographic particle image velocimetry, for point F on the pump cycle (cf. Fig. 3b). Although 2-dimensional, it contains all the 3-dimensional flow information for a rotationally-symmetric object. The axial position is relative to the throat of the stenosis, as illustrated by the vessel schematic on the left hand side (note: measurements begin at 0.84mm above the vessel throat).

Download Full Size | PDF

Finally, this process is repeated at each desired time point in the pump cycle. Figure 9 (Media 2) displays the full reconstructed 4-dimensional flow information for our rotationally symmetric stenotic model, displayed in Cartesian coordinates at a reduced spatial sampling rate for clarity. It is encouraging to note, while both axial and radial flow components are reconstructed entirely separately and that similarly each row was independently reconstructed, the resulting streamlines from their recombination follow quite closely the contours of the vessel. Additionally at the region of greatest constriction, the flow magnitude is the greatest, and predominantly in the axial direction.

 

Fig. 9 Movie of reconstructed reconstructed 3+1-dimensional velocity-field vector map for pulsatile stenotic blood flow, obtained using phase-retrieval X-ray tomographic particle image velocimetry. First frame for display is point E on the pump cycle 
(Media 2).

Download Full Size | PDF

As a test for the accuracy of our 3+1-dimensional velocity-field tomography, we applied the continuity equation to the data, which requires that at any given point in time, the flow through any transverse cross-section of the vessel should be constant, so that

Expression 10 was applied for all of the 15 resolved vertical (axial) coordinates, at each of the 8 different flow rates represented by the phase-averaged points on the pump cycle. The standard deviations were found to be consistent to within 5.9% of the mean flow.

7. Summary and conclusions

We have demonstrated the time-resolved tomographic reconstruction of the velocity field associated with pulsatile blood flow through a rotationally-symmetric stenotic vessel model. Raw flow data was acquired in the form of a dynamic speckled interference pattern which results from propagation-based X-ray phase contrast of the blood cells moving within the flow. A Fourier mask filter was used to enhance the signal of the speckle, after which regularized single-image phase retrieval was applied to decode the phase contrast from the speckle patterns and yield the projected-density maps of the flowing cells. We then used statistical correlation methods of particle image velocimetry (PIV) to calculate the probability density functions of cell displacement. The axial velocity-field component of the rotationally-symmetric flow was reconstructed with an inverse-Abel transform, with a modified inverse-Abel transform needed for reconstruction of the radial component. After performing this vector tomographic phase-retrieval velocimetry over the full pumping cycle, we have achieved complete characterization of the velocity field of the pulsatile blood flow in both space and time. This new tomographic technique for velocimetry is powerful in its simplicity and potential for very high temporal resolution, since with rotational symmetry it requires a minimum of only 2 consecutive frames to derive all the volumetric information of a single flow state. Our technique complements other computed tomographic X-ray velocimetry techniques which are not limited by assumptions of flow symmetry, but require rotation of either sample or detector and so suffer from greater restrictions on time-resolved measurement. The theory of reconstruction presented here would also be applicable to a number of other scenarios such as could be found in visible light PIV experiments.