2.1. Phase-contrast CT

In conventional X-ray imaging, image contrast arises from differences in attenuation of the X-ray beam by structures with differing linear attenuation coefficients. When imaging low-Z materials with similar linear attenuation coefficients, such as biological tissue, with X-rays in the 10–100 keV energy range, phase shift probabilities between different tissues are approximately 1000 times larger than for linear attenuation, thereby allowing visualization of structures which are not visible with conventional attenuation-based techniques (Zhou & Brahme, 2008; Wu et al., 2008). In propagation-based phase-contrast imaging the sample is irradiated with highly spatially coherent X-rays. The detector is set a certain distance downstream from the sample such that, when the X-ray wavefront is diffracted passing through the sample, the perturbed wave interferes upon free-space propagation after the sample according to Fresnel diffraction. This gives rise to interference fringes at the edges of the different sample structures consequently enhancing the edges (Diemoz et al., 2012). PC-CT uses these shifts in phase of the X-rays to provide high contrast. Both qualitative images, with strong edge enhancement but no phase retrieval, and quantitative (with phase retrieval) images can be obtained (Chen et al., 2010). Wu et al. (2008) terms these as direct `phase-contrast tomography' when the edge-enhancement without phase retrieval is utilized, and `X-ray phase tomography' where phase retrieval was performed. The direct phase-contrast tomogram is a mixture of three components [(i) the three-dimensional map of the sample's linear attenuation coefficients, (ii) the map of the three-dimensional Laplacians of its refraction indices, and (iii) artifacts related to the global distribution of tissue attenuation coefficients and refraction indices], thus there is need for phase retrieval for accurate and quantitative tomography for medical applications (Wu et al., 2008).

Burvall et al. (2011) identified seven methods of in-line phase retrieval which satisfy two conditions making them applicable to tomography. First, that they require only one image at each projection angle, and second, that they utilize analytical methods rather than iterative methods to minimize the reconstruction time. Burvall et al. also provide a scheme for choosing which reconstruction method to use based on the characteristics of the data set. Phase retrieval is recognized as an important aspect of phase-contrast imaging, but was not considered a main focus of this study for two reasons: (i) of the four programs assessed, only two (PITRE, H-PITRE) are capable of performing phase retrieval, and (ii) Mohammadi et al. (2013) provide a comparison of the free software programs ANKAPhase (Weitkamp et al., 2011, 2014) and PITRE (Chen et al., 2012a, 2013a) for implementing phase-retrieval algorithms with propagation-based X-ray tomography. As such this work only reviews the performance of the three-dimensional reconstruction software. We use phase-contrast images as the input data sets, but the results are equally valid for absorption contrast images.

Synchrotrons provide a high-intensity and highly coherent X-ray beam with parallel geometry and monochromatic energy. This allows for the properties of diffraction and interference to be exploited (Meuli et al., 2004) and avoids problems that occur with images produced using conventional polychromatic X-ray cone-beams such as low contrast, high noise and low dimensional accuracy (Kastner et al., 2010). Therefore, when in-line PC-CT is used in conjunction with synchrotron radiation, it dramatically increases the ability to image soft tissues, such as the prostate. An important advantage of the in-line PC-CT is that this method imparts low absorbed radiation dose to tissues (due to the use of high-energy X-rays) as compared with using only linear beam attenuation for contrast (Meuli et al., 2004). In addition, in-line PC-CT is a useful imaging technique due to its relatively simple set up. It does not require any crystals or additional equipment besides a high-intensity high-energy monochromatic beam and a detector (Guo et al., 2012).

2.2. Tomography algorithm

CT is a valuable medical technique that is often used to obtain parallel cross-sectional images of the affected body part or organ. The technique involves acquiring numerous `projection' images at angular interval steps around a sample and then using a computer-based geometric algorithm to reconstruct a view of a cross slice through the tissue (Meuli et al., 2004). The FBP method is one of the most commonly used algorithms to construct CT cross-sectional images (Herman, 2010) and makes use of the Radon transform (Radon, 1986). A projection image at a certain angle represents the line integral of the X-ray attenuation through a sample and the Radon transform of this integral represents the cross-sectional image of the sample along that path. The inverse Radon transform is applied across all the line integrals at different angular positions to reconstruct a single cross-sectional view (Kak & Slaney, 2001). One way to represent the projection data for filtering is by using sinograms which are plots of the absorption data for the reconstruction of a single slice (Stock, 2009). Sinograms are produced by collecting the linear signal along a specific detector row in the imaging plane for all the angular views along the 180° rotation and arranging them side by side (Chen et al., 2012a). FBP approximates the inverse Radon transform by a single convolution followed by a back projection of the data thus producing the cross-section image of the sample (Herman, 2010). The convolution involves the Fourier transform of the projection followed by multiplying it by a weighting function to filter the frequencies of the projections to produce higher quality reconstructions. Then the inverse Fourier transforms are summed over the image plane to back project the data over the cross-sectional plane (Kak & Slaney, 2001). FBP is one of the most frequently used methods for it only requires one convolution integral calculation reducing computation time compared with iterative methods that approximate and improve image quality with multiple calculations. In a comparison of the adaptive statistical iterative reconstruction (ASIR) algorithm, the model-based iterative reconstruction (MBIR) algorithm and the FBP algorithm, the MBIR method provided the best image quality of the three techniques in reducing noise and improving overall image quality (Deák et al., 2013). However, most iterative reconstruction methods require high computational power and long reconstruction times. ASIR is an exception with its reduced reconstruction time; however, it requires an initial building block of an FBP reconstructed image to accomplish this (Singh et al., 2010). Currently FBP is the algorithm used by many CT software systems and is the main algorithm used in all four software systems described here.

3.1. Data collection

Projection data used to test these programs were collected at the Canadian Light Source (CLS) Biomedical Imaging and Therapy Bending Magnet (BMIT-BM) beamline (Wysokinski et al., 2007). This is a bend magnet beamline with a field strength of 1.354 T. The ring energy is 2.9 GeV and the storage ring current is a maximum of 250 mA operating in decay mode. The critical energy of the bend magnet source is 7.57 keV. The beamline uses a Si (2,2,0) Bragg double-crystal monochronometer at a distance of 13.2 m from the source. An ex vivo formalin-fixed canine prostate was suspended in a plastic specimen tube within Knox Gelatine (Associated Brands Lp, Toronto, ON, Canada) positioned on a mechanical rotating stage 25 m from the 30 keV monochromatic synchrotron X-ray source. A Photonic Science VHR-90 (Photonic Science, Millham, Mountfield, UK) array radiation camera (FOP coupled CCD detector, 18.67 µm × 18.67 µm pixel size) was placed 5 m behind the sample to allow for free space propagation. The sample was then rotated over 180.096° and 3751 projections were taken each at a rotational step size of 0.048°. Exposure time for each projection was 200 ms with total imaging time of 44 min with total exposure of 12.5 min. Each projection was digitized to a 3968 × 251 pixel 12-bit image (TIFF greyscale format). The reconstruction software requires 16-bit TIFFs, and this change in format was made using Convert (Canadian Light Source, Inc., Saskatoon, SK, Canada). The effective pixel resolution of 12.4 µm was measured by acquiring two projection images with a needle moved a known distance perpendicular to the detector on the sample stage, adding the images together, counting the number of pixels in the known distance and dividing the distance by the number of pixels. Calibration images were taken before and after the tomographic projections were completed. Further steps in processing the data depended on the software used and are described below.

3.2. Software parameters for image calibrations and artifact removal

In-line PC-CT reconstruction requires that certain parameters be adjusted to suit an individual dataset. The following parameters were adjusted (when available) by the four programs reviewed here: dark-field correction, flat-field correction, centre of rotation alignment, beam-hardening correction, ring-artifact correction, and intensity or beam power normalization. Calibrations or corrections to the projection files usually occur prior to conversion to sinograms or slices.

3.2.2. Projection alignment

Practically, it is difficult for the centre of rotation of the specimen to fall exactly at the centre of the sinograms or the projection images. Thus, the projections need to be aligned so that problems in the slice reconstruction during the FBP algorithm can be avoided. These include streaking and alias artifacts that appear, especially around very dense structures within the tissue and along edges of objects (Fig. 1). Among the methods this can be accomplished are to align the image taken at 0° and the image taken at 180° degrees by using their histograms (as is done by NRecon), or by adjusting the centre of the sinogram horizontally so it falls on the centre of rotation (Chen et al., 2012a).

Figure 1 Previews reconstructed by PITRE. (a) Completely unaligned image, (b) partially aligned image and (c) aligned image.

3.2.5. Ring artifact removal

Ring artifacts are concentric rings that appear around the axis of rotation in an image (Fig. 2a). They come from an invalid measurement made along the projection line that occurs in every view. This can be caused by a defective pixel on the detector, a slightly mis-filtered beam, or a defective detector line parallel to the axis of rotation (Chen et al., 2013b). Dark and flat calibrations suppress, but often do not completely remove, the ring artifact due to imperfections in the dark- and flat-field images. Thus, many programs improve image quality by using a ring-removal algorithm (Bourland, 2012).

Figure 2 Reconstructed images of a canine prostate by NRecon showing images (a) without ring-removal algorithm applied and (b) with ring-removal correction.

3.2.6. Efficiency

Reconstructions of high-quality CT images are extremely computationally demanding tasks (Flores et al., 2013). Recently, image processing using the GPU of a computer has become more popular as the computational performance of GPUs has surpassed that of the CPU with the advent of powerful GPUs for computer games (Eklund et al., 2013). One report showed that the use of the GPU over the CPU accelerated reconstruction with the FBP algorithm by up to 100 times (Schiwietz et al., 2006). The average amount of time spent pre-processing, reconstructing and post-processing images was recorded for each of the four programs (Table 1).

Table 1 Timing comparison for reconstruction of 3751 projections Comparative chart of average time spent on pre-processing and reconstructing of images with all four programs on a Gateway FX Desktop (Intel Core i7-2600, CPU 3.40 GHz, 16 GB RAM, 1.5 TB HDD) with NVIDIA Quadro 6000 Video Card (384-bit memory interface, memory bandwidth 144 GB s−1). NRecon was the only program requiring calibration and pre-processing of the data prior to importing into the software for image reconstruction. No pre-processing of the data was required for the other three programs. Program Pre-processing (h) Reconstruction (h) NRecon 0.333 <1 Athabasca Recon 0 6 PITRE 0 4 H-PITRE 0 1

Table 2 lists the freeware programs and their features which were investigated for their ability to deliver high-quality images and complete different facets of PC-CT reconstruction. For image quality assessment, a reference image was generated with extra post-processing with the program Ring Remover (Z. Wei, Canadian Light Source, Inc.) in MATLAB (MATLAB R2013A, MathWorks, Natick, MA, USA). The quality of each output image (Fig. 3) was assessed by calculating the peak signal to noise ratio (PSNR) and structural similarity index (SSIM). The PSNR expresses the ratio between the maximum signal value and the value of noise that can distort the quality of the image, while the SSIM provides a measure of the similarity of two images (Wang et al., 2004). Both values provide a quantitative measure of digital image quality.

Table 2 Comparative list of the four freeware programs outlining their available data-processing components Program Required platform Dark-field calibration Flat-field calibration Beam normalization Alignment Beam hardening correction Ring artifact correction Athabasca Recon Windows Powershell Yes Yes Yes No No No NRecon None No No No Yes Yes Yes PITRE IDL Virtual Machine Yes Yes Yes Yes No Yes H-PITRE IDL Virtual Machine Yes Yes Yes Yes No Yes

Figure 3 Output images for reconstructions performed with (a) NRecon, (b) Athabasca Recon, (c) PITRE and (d) H-PITRE.

4.1. NRecon

NRecon, version 1.6.8.0, is a freeware program developed by SkyScan [now Bruker MicroCT (Bruker-microCT, Kontich, Belgium)]. A commercial version of NRecon is available that distributes the reconstruction calculations across GPUs of many computers thereby decreasing the computation time. The no-cost version of NRecon uses the GPU of a single workstation. The program was developed to work with micro-CT scanners and thus lacks some features for manipulating raw data from synchrotron CT; e.g. the program does not apply the dark-field or flat-field correction to the projections, nor does it allow beam intensity normalization. The program requires 16-bit TIFF files that have been calibrated for dark-field and flat-field correction. Dark-field and flat-field corrections on our dataset were performed with a macro plugin for ImageJ (D. Cooper, College of Medicine, University of Saskatchewan, Canada). NRecon had a reconstruction time of less than 1 h which was comparable with H-PITRE's processing time, but was faster than both PITRE and Athabasca Recon. NRecon required some pre-processing, such as dark- and flat-field calibration that took 20 min, and there was no post-processing.

After uploading calibrated projection data, NRecon automatically identifies the optimized settings for projection alignment, beam-hardening correction, ring artifact correction, and smoothing filter strength. Previews of one reconstructed slice can be generated on the main page in the Start tab; the default preview is the central slice. The program approximates the post-alignment value by comparing and aligning the histograms of the first and last projections at angular positions of 0° and 180°, respectively. The strength of the ring artifact, smoothing, beam-hardening filters and the value of the post-alignment factor can also be set manually. Final values for these parameters can be fine-tuned as the software allows the user to select from a range of previews before processing the final dataset. Fig. 2(a) shows a slice reconstructed with NRecon without ring-removal correction and Fig. 2(b) shows the same slice with ring removal applied. The reconstruction of the dataset set can be completed individually or queued into batch jobs. Reconstructed images can be saved in various formats including TIFF, JPEG and BMP.

NRecon offers some advanced settings to improve image quality. The FBP algorithm can be adjusted for use with under-sampled datasets, when objects are larger than the field of view, and for defective pixel masking. Multiple smoothing algorithms are available which alter each pixel on the projection images with a M × N (horizontal × vertical dimension, respectively) neighbourhood using either a box symmetric, box asymmetric or Gaussian kernel (Skyscan, 2011). The dynamic range of the image histogram may be edited to allow image enhancement.

NRecon has a useful manual and extensive help section; however, there is a lack of manual control to alter algorithms to optimize reconstruction precision as can be accomplished when using more direct algorithm application methods. NRecon is a user-friendly program which edits images quickly and efficiently and has an intuitive interface, but lacks features needed for many synchrotron-based datasets such as dark-field or flat-field correction and beam intensity normalization.

4.2. PITRE

PITRE (Phase-sensitive X-ray Image processing and Tomography REconstruction) can process data from high-energy X-ray source modalities including PC-CT, diffraction-enhanced imaging processing, and phase retrieval for additional PC-CT data processing (Chen et al., 2012a). PITRE accepts a number of file types including 8-, 16- and 32-bit TIFF, 8-bit PNG and 8-bit BMP. PITRE uses the CPU of the workstation rather than the GPU and so has a slower reconstruction time than NRecon and H-PITRE which use the GPU. With no pre- or post-processing for our test dataset, reconstruction took approximately 4 h using PITRE (Table 1). PITRE_BM, a batch manager program that is provided with PITRE, allows maximization of CPU usage.

For in-line PC-CT reconstruction the `Projection to Slice' data processing function was used. The other processes available are for other synchrotron-based data such as diffraction enhanced imaging. Prior to loading data one can choose between four methods used to calibrate images for dark-field and flat-field correction depending on whether the dark and flat images were taken before, after or during data acquisition, and this calibration occurs as the data is loaded. PITRE then converts the projections into sinograms for filtering and image processing. A single sinogram can be produced to create preview slices and define the main parameters (alignment, normalization and ring artifact removal) prior to producing all of the sinograms and completing the reconstruction. There are two reconstruction algorithms available to use: FBP (Kak & Slaney, 2001) and Gridrec (Dowd et al., 1999). Gridrec improves computational efficiency with negligible artifacts through the use of the gridding method for resampling the Fourier space from polar to Cartesian coordinates (Marone & Stampanoni, 2012). Fig. 4 shows output images with the Gridrec algorithm (Fig. 4a) and the FBP algorithm (Fig. 4b). PITRE provides an alignment tool that defines the centre of rotation of the sinograms by their pixel coordinates. A dynamic aligning method can also be used if there were shifts of the vertical axis during data collection. The dynamic adjustment extrapolates the centre of rotation for each sinogram by using pixel coordinates of the centre of rotation for the first and last sinogram. Beam intensity normalization is available whereby a range of projections are used to calculate and normalize an intensity curve which is then applied to the dataset. PITRE also provides a ring-removal algorithm that acts on the sinograms with an adjustable filter width. The algorithm finds the value of an average row of the sinogram by summing down each column and dividing by the number of rows, subtracts a smoothed version of this row to identify detector anomalies, and then subtracts these anomalies from each row of the sinogram (Chen et al., 2012a).

Figure 4 Reconstruction of the PC-CT data by (a) PITRE with the Gridrec algorithm, (b) PITRE with the FBP algorithm.

PITRE also offers a number of additional reconstruction options with a choice of window filters to be used with either the FBP or Gridrec algorithm, including Shepp-Logan, Hanning, Hamm, Ram-bak and Logarithm, which can be used to make data more uniform, reduce outlier pixels and reduce image noise. When a dataset is loaded, PITRE automatically defaults to the optimal filter based on the type of dataset. An extended view option is also included in this program which uses a dataset acquired over 360° of rotation rather than 180° to allow for the field of view to be increased. The program also allows conversion of 12-bit TIFF images to 16-bit TIFF format using a GraphicsMagick plugin. This is a valuable plugin as many detectors produce images in compressed 12-bit TIFF format rather than 16-bit to reduce imaging acquisition time.

PITRE is a useful detailed program that provides the ability for acute adjustments to increase the quality of volume slices; however, the program can be difficult to work with for an inexperienced user. The manual and the associated journal article are helpful, but lack explanations on how certain inputs should be defined. For example, it is not explicitly outlined as to what the units are for certain parameters. User experience could be improved with the addition of an expanded user's manual, but, regardless, there is no need for additional programs or plugins and there are a large number of manual adjustment options.

4.3. H-PITRE

H-PITRE (High performance Phase-sensitive X-ray Image processing and tomography REconstruction) provides a graphic user interface which is very similar to PITRE but is designed for use with in-line PC-CT datasets only (Chen et al., 2012b). It is almost the same program with a few extra features devoted to PC-CT and all other features removed to create a simpler more efficient program. The actual set-up and functionality of the program are almost identical to PITRE with a main `Projection to Sinogram' window and a window for viewing the reconstructed slices. The program accepts 8-, 16- and 32-bit TIFF, PNG and BMP files with a naming index and produces 32-bit TIFF images. The program makes use of the GPU of the computer to increase the efficiency and speed of the reconstruction. This does result in strict computer requirements including: a PC with at least Windows XP, an NVidia video card and a large RAM size to accommodate the sinogram production. Without any pre- or post-processing the reconstructions of the test dataset took approximately 1 h on the workstation used in this study (Table 1).

H-PITRE offers the same flat and dark calibration options, alignment method, normalization process, and ring artifact removal algorithm as PITRE and provides only the FBP algorithm to reconstruct data. Similar to PITRE, H-PITRE offers some advanced features such as an extended view option and window filters. An additional logarithmic filter is available in H-PITRE compared with PITRE which can be used alongside the other window filters. A region of interest can also be defined to crop the images, which can reduce reconstruction time and increase efficiency. PITRE_BM also works for H-PITRE to provide memory management for increased productivity. Overall the program is very similar to PITRE with a few extra features devoted to PC-CT and all other features removed to create a simpler more efficient program.

Overall H-PITRE is an excellent program for parallel-beam PC-CT. Its ability to use the GPU increases speed compared with the standard CPU version in PITRE. The program allows great freedom to alter parameters for optimal reconstructions. As the program is in its beta stage, the user's manual lacks details and explanations which can cause difficulty for in­experienced users.

4.4. Athabasca Recon

Athabasca Recon, version 1.3, was designed by the Bone Imaging Laboratory at the University of Calgary to work with parallel beam synchrotron tomographic data and provides all required features and data modifications needed to create good quality images. The program runs on a Windows Powershell command interface and uses ImageJ plug-ins to perform a number of tasks. This program does not have a graphic user interface, but is relatively simple to use. Athabasca Recon requires ITK MetaImage file (mhd) format to handle image sequences and uses an ImageJ plugin to create these files. The program allows access to some memory management settings, but the lack of GPU usage for computations makes reconstructions slower. Reconstruction of the test dataset took approximately 6 h.

The reconstruction is defined using configuration text files that outline the parameters chosen to run and build a reconstructed set of slices with Athabasca Recon. A first configuration file is used to complete dark- and flat-field correction, bad pixel removal, and variable beam power normalization while the second configuration file applies the FBP algorithm. The Beam Power Correction has four different methods for correction depending on the dataset used. Alignment to the centre of rotation is completed by an ImageJ plugin, the Align Projections Tool. This tool was used on the test data following the application of the first configuration file. Ring artifact removal is not provided by Athabasca Recon so post-processing with additional software is needed to accomplish this task.

Athabasca Recon includes some advanced parameters to improve reconstruction quality. Bad Pixel Correction is a filter which identifies bad pixels and replaces them with the average value of the four nearest neighbours. Also available are smoothing filters with definable filter radius (Gaussian or tapered cosine filters), and a low-pass filter (removes high frequency noise). Athabasca Recon allows a choice in the method of pixel interpolation including nearest neighbour (less accurate; uses the value of the nearest pixel), bi-linear interpolation (interpolation from the four nearest pixels' centres) and bi-linear interpolation with fall back (used for volume data near to the edges of the projection; bi-linear interpolation loses a half pixel width of information so the edges are extrapolated using nearest neighbour). Lastly, Athabasca Recon allows for the threshold level of the dark-field and flat-field correction to be defined which allows for any pixel below the base flat-field threshold and above the dark-field threshold to be ignored.

Athabasca Recon provides control over the parameters of reconstruction through variation in the configuration file. The manual and the program are both very explicit on the technical details allowing a better understanding of how the data are manipulated. However, the lack of a graphic user interface can be challenging to those with less exposure to command-line interface programs. The usability of the program is further reduced as various parameters must be chosen prior to reconstruction with no available method to preview the slices. Overall, Athabasca Recon may be appealing due to the availability of a variety of filters and settings for fine-tuning the reconstructions. Slow computational speed, lack of a GUI and no ring artifact removal feature are its major limitations.

Images created using H-PITRE (21.37 dB) and NRecon (20.23 dB) provided the highest PSNR when compared with the reference image, while PITRE (18.71 dB) and Athabasca Recon (19.05 dB) provided the lowest. When compared with the reference image the SSIM was the highest with PITRE (0.26) followed by NRecon (0.25), H-Pitre (0.23) and Athabasca Recon (0.18). PITRE, NRecon and H-PITRE thus have very close structural similarity indices, whereas Athabasca Recon has lower similarity with the reference image which is not surprising considering that Athabasca Recon does not have native ring removal.