Fast non-uniform Fourier transform based regularization for sparse-view large-size CT reconstruction

1.   Introduction

With rapidly increasing interest in investigating micro-scale internal structures of diverse objects in a nondestructive way, i.e., cells and macromolecules and their basic arrays, synchrotron-based micro-computed tomography (micro-CT) has been extensively applied to biomedical sciences and materials engineering [1–4]. The distinct synchrotron radiation (SR) offers a greater range of capabilities including higher spatial resolution, better phase contrast, elementally sensitive and high-speed data collection for observing dynamic processes in three-dimensions [5]. However, high-resolution micro-CT imaging usually demands a huge number of projections (more than 1000 views) for robust reconstruction under the condition of the Shannon-Nyquist sampling theorem [5].

On the other hand, sparse-view micro-CT imaging has been receiving lots of requirements from practical applications due to the following significant factors. First, a lower number of projections naturally minimize the accumulated exposure radiation to the living specimens, which is strongly desired in medical imaging [1]. Second, a short acquisition time owing to a small number of projections reduces the potential motion artefacts, as well as favors time-sensitive studies [5–7]. Besides, the reduced projection data decreases the computational cost and reconstruction time, enabling the reconstruction feasible for desktop computers or small workstations with limited capacity. However, insufficient projections result in severe streak artefacts in the reconstructed images when conventional reconstruction methods such as filtered back projection (FBP) are adopted [8].

Therefore, efficient reconstruction algorithms become crucial for sparse-view micro-CT imaging [9]. In earlier stage, several iterative algorithms [10–12] were proposed and outperformed FBP when the projections are incomplete. However, they suffer from expensive computational demand even for moderate-size CT imaging [7]. In the past decade, along with the advent of a few seminal works of [13–15], compressed sensing (CS) has received lots of attentions in signal processing including CT imaging reconstruction. CS can recover a satisfactory reconstruction from much fewer projections than required by the Shannon-Nyquist sampling theorem, by making use of the sparsity priors of an image. Specifically, the numerical results in [8] demonstrated that in a noiseless environment CS-based reconstruction methods might recover the exact image from approximately one tenth of the number of projections needed in FBP. Inspired by the CS theorem, a few CS-based algorithms have been proposed for sparse-view CT reconstructions. Examples include [7, 16–19]. Some of them have been applied to real moderate-size CT imaging. But none of them have been used to Synchrotron-based micro-CT imaging. All these algorithms need to involve a system matrix which represents the projection transform. The system matrix is linearly proportional to the product of the image size and the number of projections. The Synchrotron-based Micro-CT imaging commonly acquires hundreds of projections and produces a large image with millions of pixels. Hence, the resulting system matrix is extremely huge and impractical for workstations with limited capacity.

Yet, the classical Fourier slice theorem puts forward another straightforward way for CT reconstruction. If the projections are examined in the Fourier domain, the reconstruction can be solved directly in a 2-dimensioanl manner [20]. The huge system matrix is then avoided. Unfortunately, the projections are sampled radially while the Fourier transform needs to perform on a Cartesian grid [21]. So far, there are two outstanding ways to tackle such a problem: the pseudo-polar Fourier transform (PPFT) [22] and the non-uniform fast Fourier transform (NUFFT) [23–27]. PPFT has already been applied to medium-scale sparse-view CT reconstruction in [8, 28]. NUFFT has been proposed for the reconstruction of fully sampled CT imaging as well [29, 30].

In this paper, we proposed a non-uniform Fourier transform based method for sparse-view synchrotron-based micro-CT reconstruction. There are two significances. First, since we target at high spatial resolution micro-CT reconstruction, computational efficiency becomes critical as well as image quality. Second, if the projections are inadequate, an optimization solution will be preferable. Direct non-uniform Fourier transform is computationally costly. So, we choose the optimized NUFFT algorithm [24], which performs better than PPFT [22] in both computational cost and accuracy based on our study. The total variation (TV) demonstrates excellent performance in representing medical images in a sparse manner [18] and is chosen as the penalty term. Thus, the proposed approach comes to minimize the TV of the image subject to the constraint that the estimated image satisfies the Fourier slice theorem within a specified noise tolerance of the available projection data and that the pixel values are real numbers. We then formulated the proposed approach into an unconstrained minimization problem. It faces two major challenges: 1) there are probably millions of variables, so the algorithm is desired to be memory-efficient and fast in convergence; 2) the Hessian matrix for the objective function is difficult to solve. Consequently, the Limited-memory Broyden–Fletcher–Goldfarb–Shanno algorithm (L-BFGS) [31] optimization method is adopted. However, we utilized the back-tracking linear search instead of the recommended Wolfe conditions to find the optimal step size to accelerate the algorithm. Accordingly, a large step size is initialized empirically. Both numerical simulations and real data are tested to demonstrate the proposed algorithm. The real data set was obtained from the Australian Synchrotron facility. All the reconstruction algorithms were run on a common laptop computer. In simulations, a Shepp-Logan phantom with 256 × 256 pixels is used as the target image, and three other representative algorithms are performed to compare with the proposed approach. In real data test in which the target images are larger than 2000 × 2000 pixels, since two of the algorithms fails due to memory issues, only FBP is compared to the proposed approach. The Australia synchrotron facility only provides parallel x-ray beams, so only parallel projection is taken into consideration in this paper. But, if the projection data is rebinned from fan beam to parallel beam geometry [32], the proposed algorithm will be applied as well.

2.   Theory and method

The synchrotron-based micro-CT imaging system provides parallel projection, and both the parallel x-ray lines and the detectors are equally spaced. In principle, each projection of an object at a given angle is seen as a set of line integrals. Each line integral denotes the total attenuation of the beam of x-rays along a straight-line path through the object and then measured by the detector (see Figure 1). Different issues (or materials) at different positions each have a different attenuation coefficient. As a result, a CT image comprises a gray scale representation of the attenuation coefficient in Hounsfield units for each pixel. Since each detector records an integral of attenuation along a projection line, reconstruction is necessary to obtain the underlying CT image. Many algorithms have been developed, but in this paper we focus on those techniques which can handle the insufficiency of projections.

2.1.   Reconstruction algorithms in space domain

In the Cartesian system ($x, y$-axes) of Figure 1, arbitrary line at angle $\theta$ and position $r$ can be expressed as

Let $f(x, y)$ denote the attenuation coefficient at point $(x, y)$, and then the $\theta$-view projection along such a line can be expressed by the line integral

where $\delta$ denote the 2-dimensional (2D) Dirac delta function, which is also known as Radon transform.

In real applications, the attenuation coefficients are digitized into the pixel-based representations. Each pixel represents a tiny grid with side length $\Delta$ in the Cartesian system, and $\Delta$ represents the spatial resolution of the micro-CT system. Suppose a weight coefficient ${a}_{ij}$ denote the intersection length between the projection line and the pixel ($i, j$), and then (2) is rewritten into discrete form as

where ${f}_{i, j}$ denotes the pixel value and the target image is made up of $N\times N$ pixels. The weight coefficient gauges how much scale each pixel contributes to each projection in terms of attenuation. It is noted that ${a}_{i, j}$ could be zero when the corresponding x-ray line and pixel do not intersect. Furthermore, the Radon transform can be expressed in a matrix-vector product:

where $A$ is conventionally named system matrix derived from the above weight coefficients, $f$ denote the vectorized target image and $y$ is the vectorized projection data.

Using (4) as a constraint, many CS-based algorithms [7, 16–18] have been developed to solve the following optimization problem:

where $\left\| \cdot \right\|_{n}$ denotes the $n$-norm, $\mathrm{\Psi }$ denote the sparse transform and $\varepsilon$ denotes the noise tolerance.

These algorithms have achieved superior performance over BFP in reconstructing the underlying CT image from sparse-view projections. However, the size of $A$ is of order $O\left(LM{N}^{2}\right)$ where $L$, $M$, $N$ respectively denote the number of projections, detectors and the image length in each dimension, and usually $M\ge N$. Thus, when $N$ is larger than 1000, the size of $A$ is likely to reach more than a terabyte. Such a huge matrix certainly imposes an extremely high demand on memory and computational cost.

2.2.   Reconstruction algorithm in Fourier domain

On the other hand, a few algorithms based on Fourier transform were proposed to reduce the computation time [29] or improve the image quality for sparse-view projections [8]. To better understand those techniques, the Fourier slice theorem is briefly reviewed in this subsection, and a detailed explanation refers to [20]. The Fourier slice theorem states that the 1-dimensional (1D) Fourier transform of a parallel projection of an object $\mathrm{f}(\mathrm{x}, \mathrm{ }\mathrm{y})$ obtained at angle $\theta$ equals a line in a 2-dimensional Fourier transform of $\mathrm{f}(\mathrm{x}, \mathrm{ }\mathrm{y})$ taken at the same angle.

First, we consider a simple case that a projection of $\mathrm{f}(\mathrm{x}, \mathrm{ }\mathrm{y})$ is taken parallel to the $\mathrm{y}$ axis. Thus, we have $\theta =0$ and $r=x.$ By substituting them into (1) and (2), we obtain

Then, taking a 1D Fourier transform with respect to $x$ on both sides of (6), we obtain

On the other hand, a 2D Fourier transform is performed on $\mathrm{f}(\mathrm{x}, \mathrm{ }\mathrm{y})$ directly, and then only the line at $v=0$ (corresponding to $\theta =0$) is evaluated,

(7) equals (8) if we set $u=\omega$. Then, the Fourier slice theorem is proven for $\theta =0$.

Now, the case of $\theta \ne 0$ is taken into consideration. Let us rotate the coordinate system counterclockwise by angle $\theta$ such that one of the axes, ${x}^{'}$, is parallel to the line path of the projection at $\theta$-view. Let ${f}^{'}({x}^{'}, {y}^{'})$ represent the object $\mathrm{f}(\mathrm{x}, \mathrm{ }\mathrm{y})$ in the rotated system, where the two coordinate systems can be converted to each other by

Similarly, we get ${x}^{'}=r$. Then the projection $p\left(r, \theta \right)$ is simply the integral of ${f}^{'}({x}^{'}, {y}^{'})$ along the ${y}^{'}$-axis:

Then, a 1D Fourier transform is performed on (10) over variable ${x}^{'}$, and we get

By differentiating (9), we have

Then, substituting (9) and (12) into (11), we obtain

On the other hand, the direct 2-D Fourier transform of $\mathrm{f}(\mathrm{x}, \mathrm{ }\mathrm{y})$ is

If we let $u= \omega \cos \theta$ and $v=\omega \sin \theta$, we can obtain

The left hand of (15) denotes a straight line in the 2D Fourier space passing through the origin and simultaneously forming an angle $\theta$ with respect to the $u$ axis; and the right hand indicates a 1D Fourier representation of the projection at $\theta$-view. So, the proof of the Fourier slice theorem is completed.

The Fourier slice theorem provides a potential way to reconstruct the CT image directly in the Fourier domain as: the entire 2D Fourier space will be filled up if enough projections are collected over the rang from 0 to $\pi$; once its Fourier representation is obtained, the target image can be directly reconstructed by the inverse Fourier transform. However, (12) shows that the Fourier samples obtained from the 1D Fourier transform of the projections are distributed in the polar system. In other words, the sampling frequencies are non-uniform from the perspective of the Cartesian system. The popular 2D FFT requires to be performed on uniform frequencies in the Cartesian system and cannot be applied to it directly. Therefore, the implementations of non-uniform 2D Fourier transform and its inverse play a key role for Fourier-based reconstruction algorithms.

3.   Fast non-uniform Fourier transform based approach

3.1.   2D non-uniform DFT for CT imaging

In digital image processing, the pixels are uniformly sampled along a row or column of a raster image, which is typically in the Cartesian system. The conventional discrete Fourier transform (DFT) converts a finite list of equally spaced samples of a function into the list of coefficients of a finite combination of complex sinusoids with ordered frequencies. The sampling frequencies of the output sinusoids are integer multiples of a fundamental frequency, whose corresponding period is the spatial interval of the pixels. The 2D uniform DFT is denoted as

where $0\le u, v\le N-1$ are uniform frequencies. And $f$ can be exactly recovered by the inverse DFT (IDFT) as

Furthermore, the FFT algorithm has been developed to compute (16) and (17) in an efficient way. It can reduce the arithmetical operations from $O\left({N}^{2}\right)$ to $O\left(N\mathrm{log}N\right)$ and largely speed up the computational speed for large-scale functions.

However, the proof of the Fourier slice theorem shows that the Fourier samples from the CT projections lie on the polar system rather than the Cartesian system (Figure 2). In other words, the sampling frequencies are non-uniform in the Cartesian system. So, FFT cannot be applied directly. One must consider a case where the samples are regularly sampled in the space domain but irregularly taken in the frequency domain. Accordingly, the 2D non-uniform DFT (NDFT) is defined as follows:

where $u\left(\omega , \theta \right)={\omega \cos \theta }$ and $v\left(\omega , \theta \right)=\omega \sin\theta$ ($0\le \omega \le N-1$, $\theta \in [0, \pi )$) denote the sampling frequencies in the Cartesian coordinates, which are real numbers. (18) can be converted in terms of the form of matrix-vector product into

where the NDFT is characterized by $G$ and $F$ is the Fourier representation of $f$. It is noted that $G$ only depends on the choice of the sampling point, and could be singular [33]. In a similar way to (17), the inverse non-uniform DFT (INDFT) is defined as

Figure 2 reveals that the sampling frequencies for CT imaging over-covers at the low frequency part but lacks coverage at the high frequency part. Thus, some high-frequency information is inevitably missing during transform. Due to the same reason, the INDFT of (19) cannot recover the exact $f$. Besides, the computational complexity of (18) without optimization is $O\left({LN}^{2}\right)$ where $L$ denotes the number of views, which is undesirably slow for large-scale images.

In view of the above drawbacks of direct NDFT, a few algorithms have been developed to approximate it, such as PPFT and NUFFT. PPFT utilises the fractional Fourier transform to approximate a pseudo-polar system based Fourier transform, but suffers from a high complexity of $O\left( {{N^2}\log N} \right)$ [22]. NUFFT adheres to the Cartesian system such that the efficient fast Fourier transform (FFT) can be incorporated to obtain a set of Fourier samples with uniform frequencies, and then evaluates each target Fourier sample at non-uniform frequency by interpolation with a few neighbouring Fourier samples. By contrast, the optimized NUFFT algorithm has a complexity of $O\left( {K\log N} \right)$ where $K=2N$ is the number of fundamental sampling frequencies of FFT [24]. The recent advance in the NUFFT literature even achieved $K\approx N$ with slight sacrifice of image quality [26]. NUFFT outperforms PPFT in terms of computational cost and is adopted in the proposed method.

Because direct INDFT of (16) cannot recover the target image and simultaneously has the same high computational complexity as NDFT, an alternative way is to solve the following linear least squares problem [24],

When the measurements are sufficient, (20) can solve a $\tilde f$ faithful to $f$. Otherwise, one must add a penalty function to attain an optimal solution. Most optimization algorithms may involve an adjoint operator, e.g., ${\mathrm{G}}^{\mathrm{T}}$ where $T$ denotes Hermitian transpose, which is also provided by the NUFFT algorithm.

3.2.   Problem formulation

By concluding section 3.1, the Fourier slice theorem for CT imaging can be specified as

where $\mathrm{G}$ denotes the 2D NDFT, $\mathcal{F}$ denotes 1D FFT along each projection, and $R$ denotes the Radon transform. An example is demonstrated in Figure 3. The target image is a standard Shepp-Logan phantom of 256 × 256 pixels. The projections are taken by Radon transform with 128 views. Then, the result shows that the Fourier samples obtained by 1D FFT on the projection data are nearly the same as the ones from direct NDFT on the target image. In real applications, the projection data is taken by the Micro-CT system and denoted by $p$, then we have

Unfortunately, the reconstruction from direct INDFT is time-consuming and low-quality.

Provided that $\mathrm{y}=\mathcal{F}\left(f\right)$ is the Fourier representation of projections and then $\mathrm{y}$ is constant, we reformulate (20) into a constraint with a noise tolerance $\varepsilon$:

where $A$ represents a generalized linear transform function which could be NDFT, PPFT or a system matrix, and its adjoint operate is denoted by $At$. Because it better preserves the edge information, TV is one of the popular sparse transforms for image signals,

By adding a TV-based penalty term to the constraint, our proposed method becomes to solve the following unconstrained minimization problem,

where $\lambda$ is the penalty parameters. Given $f$ be the vector form and $D$ denotes the matrix form of the TV operator (see [34]), (22) can be rewritten as

It is noted that $D$ is a very sparse matrix, so it remains computationally efficient even for large scale problems.

3.3.   Framework of proposed algorithm

Although (26) is convex and differentiable, solving it faces a computational obstacle. The target images for micro-CT imaging usually consist of millions of pixels, so a few optimization algorithms might be computationally intractable. Besides, the calculation of the Hessian matrix of (26) is difficult when $A$ is a linear function rather than a matrix. So, we select the Limited-memory BFGS (L-BFGS) method [31] to solve (26). L-BFGS can avoid the pursuit of the Hessian matrix by calculating one approximation to it. Besides, L-BFGS has a good convergence speed and requires a limited amount of computer memory.

The objective function is denoted by

The gradient of the objective function is important to the L-BFGS algorithm. The gradient of (24) is

where $ϵ$ is a scalar to remove the singularity of the TV, but should be small enough to preserve the shape of the gradient function [35]. In practice, $ϵ$ was set to ${10}^{-12}$. The $real$ operator constrains the gradient to be consistent with the image signal which must be real. In order to avoid the involvement of large-scale matrix-vector product, we do the 2D NDFT (or $A\left(f\right)$) and the TV (or $Df$) both in a 2-D manner that the operations are performed along one dimension after another. $D$ is a quite sparse matrix, so it has little impact on memory and computational cost even for large-scale problems. Yet, it is still recommended to pre-compute $D'D$ to speed up the algorithm. $\lambda$ needs be tuned manually to achieve the best performance.

The proposed method is then assembled as in Algorithm 1 in Figure 4. The noise tolerance $\varepsilon$ is set to ${10}^{-5}$. maxIteration denotes the allowed maximum iterations, which is set to 100 for medium-size images and 20 for large-scale ones. In the algorithm, ${s}_{k}$ and ${y}_{k}$ separately record the position difference and the gradient difference at iteration $k$, but at most $m$ latest values are preserved for calculating an approximation to the Hessian matrix. The bigger $m$, the better the approximation obtained, but the more memory the algorithm consumes. So, we select $m= 100$ for medium-size images and $m= 15$ for large-size images. The loss of quality due to smaller maxIteration or $m$ is negligible. The remaining parameters are tuned empirically.

4.   Experimental results

In this section, both numerical simulations and real data are tested. All the programs are implemented with MATLAB and run on a Toshiba Laptop with Intel Core-i5-2520M CPU and 8GB memory. Numerical simulations are used to compare the proposed method with three other competing algorithms: the classical filter back projection (FBP), and another two recent approaches which are the space based algorithm adapted from Zhu et al. [7] and the PPFT based algorithm adapted from Hashemi et al. [8]. Because the latter two algorithms face computation memory challenges on our test bed, a medium-size Shepp-Logan phantom of 256 × 256 pixels are chosen as the target image, and the bin number is assumed to be 256. The simulated projections are obtained by standard Radon transform at specific angles uniformly distributed over 180°. The real data were taken from the Australian Synchrotron facility which provides parallel projections. The image size is larger than 2000 × 2000 pixels. It is too large for the algorithms of Zhu et al. and Hashemi et al. to deal with, so only FBP is performed to compare with the proposed method. The raw projection data totally contain 1800 views at 0.1° steps. When sparse-view test cases are executed, a few views with equal interval are selected.

In radiology, signal to noise ratio (SNR) is a common criterion that measures the level of a desired signal to the level of background noise in a particular image, which results in a grainy appearance. The SNR is calculated as

The Root Mean Square Error (RMSE) is another frequently used measure of the residual between the ground truth and the reconstruction, defined as the square root of the mean squared error:

Besides, the structural similarity (SSIM) index [36] is widely used to measure the similarity between two images, which has proven to be better consistent with human subjective evaluation. These three criteria are selected to evaluate the reconstruction quality from quantitative perspective. The bigger the SNR or SSIM is, the better quality the reconstruction has, but the RMSE evaluates in an opposite way. Computational complexity is significant for large-scale CT imaging, so the reconstruction time is compared as well.

The numerical simulation results are shown in Figure 5, 6 and Table 1–3. Figure 5 (a, b and c) illustrates the trend of reconstruction quality with different number of views. The quality performances for all the algorithms are improving with the increase of views. However, for a specific number of views, the proposed method has the best performance in terms of SNR and RMSE, and the second-best performance in terms of SSIM, although the gain decreases with increasing views. A related visual comparison is shown in Figure 6, and a corresponding quantitative assessment is shown in Table-1. When the projections are insufficient, lots of streak artifacts and grainy noise are observed in the reconstructions of FBP and Hashemi et al. while the methods of Zhu et al. and ours produce much clearer results owing to the regularization.

Regarding reconstruction time, FBP performs the best and the proposed method outperforms the other two methods. The reason is that FBP is a linear method whereas the Hashemi et al. and Zhu et al. methods and ours involve solving an optimization problem which needs a lot of computational cost. The computational time of the Zhu et al. method increases exponentially with the amount of projection data, which approximately follows its computational complexity of $O\left(LM{N}^{2}\right)$. The methods of Hashemi et al and ours process the data in a 2D manner, such that the computational complexity is much reduced. However, the Hashemi et al. method needs a larger-support interpolation, and thus has a higher complexity than ours. The experiment results in Figure 5 and Table-1 roughly testify the theoretical analysis.

To better demonstrate the practice of the proposed method, two different specimens were tested. One specimen is an adult mouse head, and the other is an adult mouse chest. The specimen is held in a falcon tube with thickness of 5 mm. The two target images are respectively 2020 × 2020 pixels and 2496 × 2496 pixels. The sparse-view case was examined. The sparse case adopts 180 views (see Figure 7), occupying 10% of the total views. One can estimate that the acquisition time and storage space would decrease accordingly at a similar proportion. In some studies, researchers may hope to examine more micro-scale details if sufficient projections are provided. Therefore, a third test case is executed as all the 1800 projections (see Figure 8) are completely utilized. The reconstructed images are shown in Figures 7 and 8.

First, the sparse-view test results are compared. For 180 views, the proposed method recovers nearly all the desired structure. A few interesting regions are enumerated in Figure 7. The lungs, lesser airway, furs, and most soft tissues are evident. Lots of high frequency textures of the sinuses are distinguishable. However, those fine details fail to recognize in the reconstructions via FBP. Besides, it is noted that the performance improvement of the proposed method over FBP is much better in real applications than in numerical simulations. In fact, the noise is inevitable in real applications. FBP is a linear method and not robust to noise while the proposed method can take advantage of the regularization to greatly reduce the noise effect. In conclusion, the proposed method is much superior to FBP for sparse-view reconstruction.

Then, the reconstructions with total projections are compared. Both FBP and the proposed method recover all the desired textures. But lots of tiny grainy noise still exists in the FBP reconstructions, making the fine structures in the tissues blurry. Owing to the sufficient measurements, the reconstructions of the proposed method present high-contrast contents so that a few fine-scale details in the tissues are visible, e.g., the soft tissues near the backbone in Figure 8.

Finally, it is worthwhile to mention the reconstruction time. The proposed method is able to reconstruct the target images with full projections in around 3 minutes on a laptop computer, though FBP costs less than 30 seconds. The proposed method is still potential to recover the CT image in real-time once a powerful workstation is utilized.

5.   Conclusions

A fast non-uniform Fourier transform based approach is studied in this paper, targeting at high-resolution Synchrotron-based micro-CT reconstruction. The non-uniform Fourier transform enables the proposed method to reconstruct the target image based on the Fourier slice theorem and to avoid the involvement of a huge system matrix. As a result, the proposed method can handle high-resolution CT images consisting of millions of pixels, which are difficult for many conventional and recent methods due to high memory demand. To remedy the insufficiency of projections and simultaneously improve the robustness to noise, the TV-based penalty is incorporated into the proposed optimization framework. Real data from the Australian Synchrotron facility are tested to present a practical demonstration of the proposed method. The results show that the proposed method can recover a satisfactory image with 180 views, but its quality is comparable to the reconstructions via the classical FBP with 1800 views. Therefore, the proposed method is significant to sparse-view high-resolution CT imaging.

Author′s biography

Dr. Zhenglin Wang is a research fellow of School of Engineering and Technology with Central Queensland University, Australia. His research interests include computer vision, machine learning, deep learning, and precision agriculture.