5.1. Implementations of X-ray speckle-based phase-contrast tomography¶
From version 1.1, methods for speckle-based phase contrast imaging have been implemented into
Algotom. This is the result of the collaboration between the author, Nghia Vo, and his
collaborators: Hongchang Wang, Marie-Christine Zdora, Lingfei Hu, and Tunhe Zhou; who are
experienced with developing and using speckle-based imaging methods. This technical note
is to summarize the work done (for more detail, see [R22]).
When a sample interacts with a coherent X-ray beam, it will cause reduction in the intensity
and change in the direction of the beam. The latter effect comes from the phase shift of the
X-ray wave. Using the first effect to image samples in tomography, known as X-ray absorption-contrast
tomography, is a widely used technique. However, using the second effect for imaging samples
is much more challenging in practice. The resulting images, i.e. phase-shift images,
can be used for tomography to visualize internal features of samples having small differences in
densities at high contrast quality. This is the advantage of the technique, known as X-ray
phase-contrast tomography (X-PCT), over conventional X-ray tomography.
To retrieve a phase-shift image, there are two basic approaches: measuring change in the
direction of a beam then performing surface reconstruction; or matching a wave-propagation model
to measured intensities. Speckle-based techniques use the
first approach. They are very simple to use in terms of set-up, data acquisition, and data processing.
The idea of the techniques is to measure the shift of each speckle-image pixel, caused by a sample,
by comparing the images with and without the sample. Because the speckle-pattern size is larger
than the pixel-size of a detector, to resolve the shift for each pixel we use a stack of speckle-images
scanned at different positions, with and without sample, and analyze the images using a small window
around each pixel.
For producing a random speckle-pattern, a diffuser made of a sandpaper can be used. The purpose of
a diffuser is to provide a reference-image used to detect local displacements caused by a sample.
It is crucial to get a high-contrast reference-image with the average feature size of a few pixels, e.g. 5-9 pixels.
A reference-image with the visibility, calculated as the ratio between the standard deviation value and
the mean value of pixel intensities, above 10% is good enough for use. For high-energy X-ray sources,
users can stack sandpapers together or using a box of material powder to improve the contrast of a
speckle-image. Other practical considerations for setting up an experiment are as follows:
The diffuser can be positioned before or after a sample depending on experiment conditions. For example,
in a parallel-beam system with a highly spatial-coherent source the diffuser can be placed closer to
the source than the sample to make use of the edge-enhancement effect which helps to improve
the contrast of the speckle-image.
The sample-detector distance should be chosen as a compromise between increasing the displacement
effect and reducing the edge-enhancement effect caused by highly spatial-coherent sources.
There are different ways of shifting a diffuser and acquiring a tomogram at each position. However, using
a spiral path has been proven to be practical and efficient [R25].
The distance between two speckle positions should be larger than the analysis window, e.g. 5-11 pixels,
to ensure that each speckle-pattern in the analysis window is completely different. This improves the
robustness of methods measuring pixel shifts. Using this acquisition scheme, 20 positions of a diffuser
is enough to retrieve a high-quality phase-shift image. However, for tomography systems with fluctuating
sources, higher number of positions, e.g. 30-50, is needed.
Due to mechanical error, moving a sample in and out of the field-of-view repeatedly for each diffuser
position can cause small shifts between the same projections of different tomograms. This problem has
a significant impact to the quality of processed data and is difficult to correct. To avoid it, the best
scanning approach is to scan all positions of a diffuser first, then collect tomograms of a sample at each
diffuser position. This approach may result in small displacements between the same speckle positions due
to mechanical error. However, it is correctable by image alignment using a small area of empty space within
the sample image.
The core idea of the technique is to find the shift of each pixel of a speckle-image caused by a sample.
This is done by: selecting a small window (5-11 pixels) around each pixel of the sample-stack image; sliding
this window around a slightly larger window (margin ~10 pixels) taken from the reference-stack image and
calculating the cost function [R25] or the correlation coefficient [R2] between
two windows at each position. The resulting correlation-coefficient/cost-function map is used to locate the
maximum/minimum point where sub-pixel accuracy can be achieved by using a differential approach or a polynomial
fitting approach. The shift of a pixel is the distance from the maximum/minimum point to the center of the map.
The procedure of finding the shift of each pixel is depicted in Fig. 5.1.1.
Performing 2D searching for every pixel of a 2k×2k image is computationally very expensive which is why
using multicore-CPU and GPU for computing is crucially needed. An approximate approach to reduce the
computational cost is to perform 1D search [R24] using middle slices in vertical and horizontal
direction of image stacks, to find shifts in x and y-direction separately.
Fig. 5.1.1 Demonstration of how to find the shift of each speckle-pixel.¶
The result of the previous step is separated into an x-shift image and a y-shift image, i.e. gradient images.
A phase-shift image is then retrieved by applying a method of surface reconstruction, or normal integration
(Fig. 5.1.2). There are many available options for implementing this step. However, Fourier-transform-based
methods [R7, R15] are preferred over least-squares methods due to their low computational
cost which is critical for tomography. The disadvantage of these Fourier methods is that the DC-component (average value of an image)
is undefined resulting in the fluctuations in background between phase-retrieved images. This effect, however,
can be corrected (Fig. 5.1.3) by using the double-wedge filter as described in [C1]
Fig. 5.1.2 Phase-shift image (c) is retrieved by normal integration using two gradient
images: (a) x-direction; (b) y-direction.¶
Fig. 5.1.3 (a) Fluctuation of grayscale values in a sinogram caused by the FT-based surface-reconstruction
method. (b) Corrected image after using the double-wedge filter.¶
5.1.3.3. Extracting transmission and dark-field signals¶
Another interesting capability of the speckle-based technique is that transmission image
(absorption-contrast image) and dark-field image (small-angle scattering signal, not to be confused
with dark-noise of a camera) can be extracted from data together with the phase-shift image (Fig. 5.1.4).
There are several ways to determine dark-signal image for correlation-based methods.
For the cost-based approach [R25], dark-signal image is easily to be obtained as a part of the model
equation.
Fig. 5.1.4 All imaging signal retrieved by the speckle-based technique can be used for tomography. (a) Phase-shift image.
(b) Transmission image. (c) Dark-field image.¶
Practical design-principles have been followed in the implementation:
To ensure that the software can work across platforms and is easy-to-install; dependencies are minimized,
and only well-maintained Python libraries are used.
For high performance computing, making use of GPU, but ease of understanding and use; Numba library is used instead
of Cupy or PyCuda.
Methods are broken down into building blocks to be able to run on either small or large memory RAM/GPU.
More importantly, this design allows users to customize methods or build data processing pipeline.
Top layer methods, API reference, for the software are as follows:
Reading images from multiple datasets, in tif or hdf format, and stacking them.
Finding local shifts between two images or two stacks of images.
Performing surface reconstruction from gradient images.
Retrieving phase-shift image given two stacks of images.
Extracting transmission image and dark-field image.
A dedicated module in Algotom, named correlation,
is a collection of methods as the building blocks for the top layer methods described in the previous section.
The first block is a method to generate correlation-coefficient map between two 2D/3D images (Fig. 5.1.1).
This is the core method to find the shift between images. It works by sliding the second image over the reference
image and calculating the correlation coefficient at each position. There are many formulas to calculate this
coefficient. Here, we use Pearson’s coefficient
as it has been proven to be one of the most reliable metrics. The method includes low-level implementations for
specific cases: 2D or 3D input, using CPU or GPU.
The second block is a method to locate the maximum/minimum point of a correlation-coefficient/cost-function map
with sub-pixel accuracy where there are two approaches selected: either a differential approach [R6]
or a polynomial fitting approach [R3]. At low-level are implementations to handle different cases: 1D
or 2D input, using the differential method or fitting method.
The above blocks are for finding the shift of each pixel using a small window around it. This operation is
applied to ~ 2k × 2k pixel. In practice, input data for retrieving a phase-shift image is two stacks of images;
each stack is around 20 images (20 speckle-positions); each image has a size of 2k × 2k. Total shape of the input
is 2 × 20 × 2k × 2k with the size of ~300MB (16-bit image). As can be seen, many strategies can be used to
parallelize workload. Here we find that processing data using chunk-by-chunk of image-rows in one go and
calculating the shifts row-by-row in parallel is the most efficient way in term of memory management, performance,
and code readability.
As mentioned above, the next building block is a method for finding the shift of each pixel in a chunk of
image-rows in parallel. The method includes low-level implementations for different cases: 1D or 2D search, 2D or
3D input, CPU or GPU computing. For GPU, to reduce the overhead of transferring data and compiling functions,
the first two blocks are implemented at GPU-kernel level.
The top building block is a method for processing full-size images.
It includes many options for processing at the lower-level blocks. The chunk-size option enables the method to run on
either small memory or large memory of RAM or GPU. Other top-layer methods listed in the previous section are
straightforward to implement either directly or by making use of the methods in the correlation module.
Fig. 5.1.6 Building blocks of the correlation module.¶
Fig. 5.1.7 Speckle-based tomographic experiment at beamline I12.¶
Data collected at beamline I12
at Diamond Light Source are used for demonstration. Details of how data were acquired are as follows:
A box of fine sand was used as a speckle generator and can achieve a visibility of 13% at 53keV X-rays with the
detector-sample distance of 2.2m. A detector with the pixel size of 7.9 µm was used. Image-size is 2560 and 2160
in height and width. The speckle-size is around 8 pixels. The intensity of the beam profile is very stable which
is an important advantage of beamline I12. The sample is a picrite basaltic rock
from Iceland.
20 speckle positions following a spiral path
with the step of 30 times of the pixel size were used for scanning.
Speckle images without the sample were acquired at all positions first. Then for each speckle position a
tomographic scan of the sample, 1801 projections, was acquired. This strategy ensures that projections at
the same angle are not shifted between speckle positions. Due to mechanical error, the diffuser positions
were not the same between the first scan (without the sample) and the second scan (with the sample).
This problem can be solved by image alignment
using free-space areas in each image (Fig. 5.1.8).
Fig. 5.1.8 Demonstration of the impact of image alignment. (a) Small area of an image which is the result of dividing
between speckle-image with sample and without sample. (b) Same as (a) but after image alignment.¶
The following presents how the data were processed:
Reference-images for each position are loaded, averaged, normalized (flat-field corrected), aligned,
and stacked.
For each angle of tomographic datasets, projections at different speckle-positions are loaded, normalized, and stacked.
Phase-shift image is retrieved from two previous image-stacks (Fig. 5.1.9) using a single function.
Full options for choosing back-end methods, surface reconstruction methods, and searching parameters are at
the API reference page.
Fig. 5.1.9 Speckle-image stack (a). Sample-image stack (b). Phase-shift image (c) retrieved from (a) and (b).¶
Algotom implements three approaches: the correlation-based method using 1D [R24] and 2D [R2] search,
and the cost-based approach [R25], known as the UMPA (Unified Modulated Pattern Analysis) method.
A summary of computing time for retrieving a single phase-shift image using different options is shown
in Fig. 5.1.10 where the window size is 7 and the margin is 10. As can be seen, there is a huge
speed-up of computing time if using GPU.
Fig. 5.1.10 Comparison of computing time using different approaches.¶
For tomographic reconstruction, phase retrieval is applied to all projections then the sinograms are generated
for reconstructing slice-by-slice. This step can be manually parallelized for multiple-CPUs or multiple-GPUs to
reduce computing time. In practice, users may want to tweak parameters and check the results before running full
reconstruction. This can be done by performing phase retrieval on a small area of projection-images.
Fig. 5.1.11 shows reconstructed images in horizontal and vertical direction from the three approaches
where ring artifact removal methods [R10, R19] and the FBP reconstruction method were used. There
are several interesting findings from the results. Firstly, the 1D-search method gives less-sharp images than
other methods but with better contrast and clearer features. There is not much different between the 2D-search
method and the UMPA method out of the low-pass component. However, the main advantage of the UMPA approach over
the others is that three modes of image can be retrieved at the same time as shown in Fig. 5.1.12.
This figure is also a showcase for the speckle-based tomography technique where phase-shift images give
better contrast than transmission-signal images (red arrows). The technique reveals interesting features
of the sample which are mineral olivine. Because the olivine is a crystal it can enhance dark signal as
shown in Fig. 5.1.12 (c,f). Making use of dark-signal images to gain deeper understanding of materials
is a very promising application of the technique.
Fig. 5.1.11 Horizontal slice and vertical slice of reconstructed volumes from the 3 approaches: the 1D-search method (a,d);
the 2D-search method (b,e); and UMPA (c,f).¶
Fig. 5.1.12 Horizontal slice and vertical slice of reconstructed volumes from 3 imaging modes: phase-shift image (a,d);
transmission image (b,e); and dark-signal image (c,f).¶
Python codes used to process data for this report are at here.
Detailed references can be found in [R22].
