RAPToR: A Resampling Algorithm for Pseudo-Polar based Tomographic Reconstruction

08/17/2017
by   Shahar Tsiper, et al.
Technion
0

We propose a stable and fast reconstruction technique for parallel-beam (PB) tomographic X-ray imaging, relying on the discrete pseudo-polar (PP) Radon transform. Our main contribution is a resampling method, based on modern sampling theory, that transforms the acquired PB measurements to a PP grid. The resampling process is both fast and accurate, and in addition, simultaneously denoises the measurements, by exploiting geometrical properties of the tomographic scan. The transformed measurements are then reconstructed using an iterative solver with total variation (TV) regularization. We show that reconstructing from measurements on the PP grid, leads to improved recovery, due to the inherent stability and accuracy of the PP Radon transform, compared with the PB Radon transform. We also demonstrate recovery from a reduced number of PB acquisition angles, and high noise levels. Our approach is shown to achieve superior results over other state-of-the-art solutions, that operate directly on the given PB measurements.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 2

page 3

page 5

page 9

09/17/2020

Improved recovery guarantees and sampling strategies for TV minimization in compressive imaging

In this paper, we consider the use of Total Variation (TV) minimization ...
06/15/2019

An Exact and Fast CBCT Reconstruction via Pseudo-Polar Fourier Transform based Discrete Grangeat's Formula

The recent application of Fourier Based Iterative Reconstruction Method ...
08/04/2017

Accelerated Image Reconstruction for Nonlinear Diffractive Imaging

The problem of reconstructing an object from the measurements of the lig...
10/09/2012

Level Set Estimation from Compressive Measurements using Box Constrained Total Variation Regularization

Estimating the level set of a signal from measurements is a task that ar...
12/24/2015

Fast Acquisition for Quantitative MRI Maps: Sparse Recovery from Non-linear Measurements

This work addresses the problem of estimating proton density and T1 maps...
05/29/2017

PCM-TV-TFV: A Novel Two Stage Framework for Image Reconstruction from Fourier Data

We propose in this paper a novel two-stage Projection Correction Modelin...
05/05/2018

Polar Wavelets in Space

Recent work introduced a unified framework for steerable and directional...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

I Introduction

Projectional radiography imaging is the most widely used imaging modality in medicine today. The computed tomography (CT) technique for 3D imaging was invented in the early 1970’s by Sir Godfrey Newbold Hounsefield, and since then CT scanners have come a long way. Today’s scanners can perform full body scans in under 30 seconds, incorporating detector arrays with hundreds of elements that acquire more than a thousand different readings in every revolution of the scanner. Despite significant technological progress in the manufacturing of CT scanners, the reconstruction algorithms in use have only slightly evolved from the very first methods developed more than 30 years ago. The vast majority of commercial algorithms are based on the filtered back projection (FBP) algorithm [1, 2] for parallel-beam (PB) and fan-beam scans, or the Feldkamp-Davis-Kress (FDK) algorithm [3, 1] for cone-beam/helical scans. Modern reconstruction approaches, which rely on iterative optimization based solvers, are seldom used in commercial machines due to the prohibitive computational complexity involved, leading to long reconstruction durations.

Recently, a paradigm for processing discrete tomographic measurements has been proposed, which relies on the pseudo-polar (PP) Radon transform (PPRT) [4, 5, 6]. PPRT is an algebraic mapping that relates an object to its discrete tomographic projections. This transform has many advantages, further described below, that make it an appealing framework for reconstruction of tomographic images. However, the PP framework assumes the projections are taken over a non-uniform set of angles and detectors, and therefore cannot be used directly on today’s tomographic scans.

There are three main advantages to using the PP framework. First, the PP grid points are closer to the polar grid points than the Cartesian grid, thus interpolating polar measurements to the PP grid is easier and more accurate than transforming them to a Cartesian grid, which is the standard practice today. Second, the PP Fourier transform (PPFT) and the PPRT can be computed with a fast accurate algorithm 

[5, 6]. The adjoint PPRT operation is computable at the same speed and accuracy as the forward PPFT, thus allowing a fast iterative inversion scheme. The third advantage, is that the algebraic system that describes the PPRT has a significantly lower condition number than an equivalent PB system. This allows to reconstruct from the PP measurements significantly faster and more accurately than from PB samples, with lower noise amplification. Although these are significant advantages, the PP framework is not used today for commercial tomographic reconstruction, since CT machines do not acquire measurements on this grid.

Fig. 1: Diagram of the different transforms. (a) The Shepp-Logan phantom . (b) The PB sinogram defined in creftype 1 and creftype 6. (c) The 1D-DFT of the PB sinogram creftype 8, which samples the 2D-CSFT of in polar coordinates. (d) The PP sinogram creftype 11. (e) The 1D-DFT of the PP sinogram creftype 10, which evaluates the 2D-CSFT of on the PP grid.

In this work, we propose a way to bridge the gap between PP theory and real-world CT scans, using a Resampling Algorithm for Pseudo-polar based Tomographic Reconstruction (RAPToR). Our main contribution is a resampling algorithm, which is based on modern sampling tools, that transforms actual CT measurements to comply with the PP model. This approach relies on a mathematical derivation of a subspace, in which the tomographic measurements lie. The subspace and its respective kernel take into account the geometrical structure of the tomographic scan, expanding the work in [7], and utilize tools from generalized sampling theory [8]. The resampling technique presented is shown to be numerically accurate, and to simultaneously denoise the measurements, while being computationally efficient. We demonstrate that a fast iterative solver (FISTA-TV [9, 10]) coupled with the PP framework, produces excellent reconstructions from the transformed PP measurements, surpassing other state-of-the-art solutions. Better performance is maintained for both low signal to noise ratios (SNR) and when given a reduced number of measurements. These translate to radiation dosage reduction for CT scanning.

Prior works have already suggested ways to benefit from the PP framework. The authors in [11, 12, 13] proposed iterative reconstruction schemes based on the PP transform, yet assumed that the measurements were acquired directly on the PP grid. Unfortunately, this assumption is not valid for CT scanners in use today. In [4] the uniform PB measurements are interpolated to the non-uniform PP grid. The suggested interpolation process is performed in the frequency domain, and involves extrapolations, which suffer from numerical inaccuracy, usually leading to unwanted reconstruction artifacts. Furthermore, no results or comparisons were provided to assess the quality of the interpolation. In [14] the measurements are acquired over uniformly spaced detectors, and non-uniform angles, so that the angles comply with the PP paradigm, which might be challenging for commercial CTs. The measurements are then transformed to the PP grid, by performing frequency based interpolation solely over the detector axis.

Our approach offers a fundamentally different methodology that generates PP measurements from a PB projection, by exploiting the inherent structural properties of the tomographic scan. All the interpolations in our solution are performed in the spatial domain of the measurements, and operate with a low computational burden. In addition, we do not assume that the acquisition angles comply with the PP angles. The computational complexity of both resampling and reconstruction is on the order of , where is the number of pixels in the reconstructed CT axial scan. We show that our method outperforms state-of-the-art algorithms for tomographic reconstruction, both in of computational complexity and in quality.

While the method presented in this manuscript focuses on PB measurements, we have found it applicable as well to both fan-beam and cone-beam reconstructions, as acquired in modern CT devices. Our method can be used on these acquisition schemes by first applying a rebinning step [15, 16, 17, 18] on the fan-beam/cone-beam measurements that transforms them into PB measurements. Rebinning to PB measurements is used today in the reconstruction flow of commercial devices for enabling fast filtration and back-projecting using the PB FBP algorithm. Our approach can then replace the traditional FBP step and enable fast reconstruction via the PP grid, with improved noise resistance and accuracy, even when compared to direct reconstruction from the fan/cone-beam measurements. This application-oriented scheme for reconstructing from fan and cone-beam acquisitions is further discussed in the Appendix A.

The paper is organized as follows. Section II contains the mathematical prerequisites and defines the problem. In Section III we review the PP transform and discuss its properties, followed by Section IV where we present our subspace approach for resampling the sinogram. Section V suggests a reconstruction method, composed of resampling to the PP grid and then reconstructing using a modern and a fast iterative solver, regularized by the total variation (TV) [19] norm. An experiment is detailed in Section VI, followed by a conclusion in Section VII.

Ii Problem Formulation

We begin by describing the acquisition process of a PB tomographic scan, and then formulate the problem we aim to solve, first using continuous notation and then defining the discrete sampling process. The 2D axial object to be scanned is denoted by , where are spatial coordinates. The value of the scanned object in each point, represents the X-ray attenuation coefficient. We assume that the support of the object is confined to a circle with radius , so that , for . In addition, we define as the discrete image on a Cartesian grid, such that , where is the spacing between two pixels and . Requiring that ensures that the discrete image spans the entire support of the object.

The CT measurements are composed from projections, denoted by , where is the projection angle and is the detector location. These projections are commonly referred to as the sinogram. The relation between the sinogram and the scanned object is given by the PB Radon transform

(1)

The 1D and 2D continuous space Fourier transforms (CSFT) of the sinogram, and the 2D-CSFT of the object are denoted by , and respectively, and given by

(2)
(3)
(4)

Note that the definition of in creftype 3 is not intuitive, since a CSFT performed on an angular variable is not a common procedure. In Section IV we use creftype 3 for constructing the subspace of the measurements.

The function and the 2D-CSFT of the object are related by

(5)

This property is commonly referred to as the Fourier slice theorem [1]. Its intuitive meaning is that a PB projection of an object taken at an angle , is equivalent to a line in the frequency domain, sampling the 2D-CSFT of the object, and intersecting the horizontal frequency axis with an angle of .

Fig. 2: Presenting the object and its PB projections. (a) The ground-truth digitized sheep phantom , with full dynamic range (b) the polar sinogram creftype 1, (c) the 2D-CSFT of the sinogram creftype 3. The spectral support structure is shown to be limited to the region creftype 15.

The sampled PB sinogram is given by

(6)

Here, and are the spacings between two adjacent projection angles and detectors respectively. We use the short-hand notation

(7)

where the PB sampling operator operates on a continuous sinogram , producing its discrete samples . The full PB scan takes place over an angular range of radians and a detector axis of length . Note that the PB sampling process is uniform with respect to the variables of , namely .

The 1D discrete Fourier transform (DFT) of is defined as

(8)

and denoted . Here is the 1D-DFT operator operating on the columns of its argument. Using creftypeplural 8, 6 and 5, the values of can be related to by

(9)

for and . This establishes that a PB scan is equivalent to sampling the 2D-CSFT of the scanned object on a polar grid. An example of a polar grid in the frequency domain can be seen in Fig. 1(c).

One approach for reconstructing a PB scan is to recover from the polar measurements of the 2D-CSFT given in creftype 9, by first transforming the measurements to using creftype 8. This approach was recently demonstrated by [20, 21]. However, inverting the 2D-CSFT from a set of points known on a polar grid is a difficult task, both in terms of computational complexity and in terms of inversion accuracy. We will compare our method to the SParse Uniform ReSampling (SPURS) [21] algorithm, a state-of-the-art algorithm for performing non-uniform to uniform resampling in the frequency domain, that surpassed the non-uniform fast Fourier transform algorithm proposed in [20].

The main algorithms in use today, FBP [2] and FDK [3], approximate the solution to this inverse problem, yet lead to substantial reconstruction artifacts. They offer a direct one-shot reconstruction approach, without applying any interpolations, by performing one-dimensional filtering on the sinogram and then back-projecting it. However, these direct methods require a multitude of measurements to produce an image of clinical quality and produce more artifacts than recent iterative approaches as further discussed in [22]. In addition, their computational complexity is on the order of , which is higher than that of a single iterative step of our algorithm.

We propose a new reconstruction scheme composed of two steps: we first transform the given PB measurements creftype 6 to a modified PP sinogram , defined below and illustrated in Fig. 1(d), by the “spatial resampling” arrow. Next, we solve the linear PP system given by , obtaining the reconstructed object , where is the PPRT linear operator that takes an image , and returns the PP projection samples . By resampling our measurements from the polar grid to the PP grid, shown in Fig. 1(e), we gain all the advantages of the PP framework, including improved numerical accuracy, algebraic stability and fast iterative inversion. The fast inversion is due to the low condition number associated with the PPRT operator, compared to a PB Radon operator. These advantages eventually lead to more accurate reconstructions while reducing the overall computational complexity.

Iii The Pseudo-Polar Transform

In this section we define the PP transform and mathematically relate it to PB scanning. The relations established here, are used in the next sections for formulating a robust resampling scheme that transforms the measurements between their acquisition grid and the PP grid, alleviating the main drawback of the PP framework.

The PPFT evaluates the 2D discrete space Fourier transform (DSFT) of an image on a PP grid, seen in Fig. 1(e), over the frequency domain . This grid is also known as a “concentric squares” grid [4], and occasionally referred to as “equally sloped tomography”. The PPFT operator computes the 2D-DTFT values of on the PP grid, and is given by

(10)

where , and . We define the PPRT using its connection to the PPFT, by applying an inverse 1D-DFT on the columns of :

(11)
(12)

where and , such that is the PP sampled sinogram. The relations between the PPFT, PPRT and the PB tomographic scan are illustrated in Fig. 1.

Next, our goal is to obtain a relation between the PP sampled sinogram and the continuous sinogram , in order to define the PP sampling operator, which is used for resampling the samples onto the PP grid. By combining the discrete Fourier slice theorem (Theorem 1 in [4]), together with creftypeplural 11, 10, 9 and 8, we get

(13)

where is the PP sampling operator, and

(14)

define the non-uniform PP grid points in the domain. The above expressions creftypeplural 14 and 13 define the direct sampling method for the PP sinogram, as can be seen in Fig. 1(d).

To summarize, by relying on creftypeplural 14, 13 and 6 we formulate in the next section our first step; a subspace-based resampling from the given PB sinogram to a corresponding PP counterpart . Once we have the PP measurements, we reconstruct the object using the PP operators given in creftypeplural 11 and 10.

Fig. 3: The different acquisition grids in Fourier (frequency) domain. The Cartesian grid (left) is the grid used by the 2D-DFT. After parallel-beam scanning we obtain the values of the Fourier transform of the scanned object on the polar grid (right). The pseudo-polar grid (middle), acts as an intermediate step between the two grids.

Iv The Sinogram Subspace

Iv-a Mathematical Definition

We now formulate a subspace prior that will aid in the resampling of the PB sinogram to the PP grid, and can also be used for denoising tomographic measurements. A PB tomographic scan of any bounded object yields a sinogram that has an approximately compact support over its 2D-CSFT  creftype 3 [7], in the region

(15)

Here, is the maximal spatial frequency of the detector axis, determined by the distance by adjacent detectors , such that . The parameter determines the intersection point with the axis, as can be seen in Fig. 4. The transform of a sinogram to the 2D-CSFT domain is given in creftype 3. An illustration of can be seen in Fig. 4.

The authors in [7] showed that more than of the sinogram’s energy is confined to , when the intersection constant is set to . We extend their work by introducing as a variable parameter in the region , in order to empirically improve results. By choosing , we get a slightly bigger region , allowing for greater flexibility in our algorithm. A performance analysis for the selection of different values, is given in the next subsection.

Fig. 4: The support region of the sinogram in Fourier domain. More than 98% of the energy is located within this shape when . is the diameter of the scanned object and is the maximal spatial frequency on the detector axis , determined by the spacing between each adjacent detector .
(16)

 

The sinogram approximately lies in an SI subspace , defined by all the functions whose 2D-CSFT is limited to . Every continuous sinogram is spanned by , so that

(17)

for some (possibly infinite) matrix . The kernel is computed by performing an inverse 2D-CSFT over an indicator function of the region . The analytic expression for is shown in creftype 16 and its full derivation is found in Appendix B.

Since the region has finite support in the frequency domain, we deduce that its respective spanning kernel has an infinite support in the sinogram domain. Multiplying the kernel with a window function limits its support, which leads to a low computational burden when performing direct interpolation in the sinogram domain. In our implementation, we use a Hamming window given by:

(18)

with the continuous normalized radius function calculated by

Here, is a parameter which sets the effective radius of the window in natural units, and is the total number of X-ray detectors.

Once the kernel is multiplied by the window , the sinogram is spanned by a limited number of coefficients at every given point, such that

(19)

where represents the coefficients that span over the modified subspace , with its associated kernel function . An example of the different kernels and (with and without multiplying with a window function, respectively) and their 2D-CSFT is shown in Fig. 5 for .

The number of coefficients used for a direct computation of a single point using creftype 19 is proportional to . Modifying the choice of the window size offers the user a trade-off between improving the resampling accuracy and reducing the overall computational complexity. In our experiments, any value of achieves satisfactory results. Throughout the rest of the paper we fix to get a good balance between accuracy and speed when implementing the direct convolution computation as described by creftype 19. Increasing the size of further leads to diminishing results, as shown in a simulation study conducted for optimizing the window size parameter, which is presented in Fig. 6.

Fig. 5: (a-b) The infinite kernel that spans and its 2D-CSFT. (c-d) The kernel multiplied by the window creftype 18, with size , and its 2D-CSFT. This kernel spans the subspace and is denoted by .
Fig. 6: Comparison of different window sizes. We compare the interpolation results for various window sizes vs. the input sinogram’s noise ratio. A window size of 6 is sufficiently big for achieving satisfactory results, and the performance of using bigger windows diminishes significantly. The complexity of the resampling procedure is proportional to .

Iv-B Denoising with the Subspace

We now wish to demonstrate the sinogram denoising capabilities of the derived subspace . The uniformly sampled PB sinogram contains noise induced by various physical phenomena related to the X-ray projection process and properties of the detectors. By filtering the noisy sinograms with the kernel function we cancel the components of that lie outside of the region , therefore removing noise components, while preserving the informative parts of the sinogram required for reconstruction.

To test the subspace denoising performance we first add Poisson noise to a ground truth sinogram obtained by simulating a PB projection of a numerical sheep phantom. The simulation is performed with AIR Tools v1.3 [23]. We denote the contaminated measurements by . Next, we denoise the sinograms by convolving with the kernel function to obtain

(20)

Here are the denoised measurements that were filtered by the kernel .

The denoised sinogram is compared to the ground truth sinogram for multiple values of , under different noise levels. The comparison metric is chosen as SNR measured in dB, and calculated by the following function

(21)

Figure 7 shows that filtering with the subspace kernel successfully denoises the sinogram when compared to the input measured SNR, given by . The black diagonal line denotes equality between the input noise level and the output noise level . Under moderate noise conditions of SNR, we observe considerable gains of more than , when choosing the optimal value for the intersection parameter . According to the simulation in Fig. 7, the optimal value of is empirically found and set to be .

Fig. 7: Denoising performance and optimal choice for . In the left graph, the line plots show the output SNR calculated between and the vs. the input SNR, calculated between and for different values of . The diagonal black line is where the input SNR equals to the output SNR. In the right graph, the line plots show the output SNR vs. values of B, for different levels of input SNR. The denoising is effective for moderate levels of SNR, and less effective for low SNR. In addition, for a relatively clean input there is little to no gain in the denoising procedure, due to the small information loss it incurs.

Note that when trying to denoise images with high levels of SNR (very low noise energy) we observe the performance diminishes. That is because the region is an approximation that contains most, yet not all, of the necessary spectral components of the PB radon transform of an object. In Section V we show that in spite of this fact, reconstructing from a denoised sinogram fully maintains fine details of the scanned object.

We conclude that by simply restricting the 2D-CSFT support of a PB sinogram to , we can faithfully denoise it from additive noise.

Iv-C Resampling the Sinogram

Let be the SI subspace spanned by the windowed kernel . The discrete PB tomographic samples, , can then be modeled as

(22)

where is the finite dimensional underlying coefficient matrix that represents the continuous sinogram over the subspace , is the PB sampling operator creftype 6 that samples on the uniform PB grid, and is the subspace operator corresponding to convolution with the kernel function .

Our primary goal in this section is to resample the sinogram from PB to PP coordinates. We first recover the underlying coefficients, denoted by , from the PB measurements, , by inverting creftype 22 in the frequency domain. Once we recover we can compute the sinogram over any grid constellation, and specifically over the PP grid.

According to modern sampling theory [8], if a (continuous) signal lies within a subspace and sampled by , we can perfectly reconstruct it by applying the operator to its samples , as long as the operator is invertible. Otherwise, we can obtain a consistent reconstruction by applying the pseudo-inverse .

In practice, computing the exact inverse is a numerically difficult task. Therefore, we propose to solve the following regularized least squares (LS) convex deconvolution problem,

(23)

The operator can be written as a linear convolution such that

(24)

By using the known DFT identity for linear convolution, the above expression is evaluated in complexity by

(25)

where denotes an element-wise multiplication between matrices.

The LS objective is regularized with a quadratic norm in order to avoid an ill-posed formulation, which is generally the case when reconstructing from a small number of tomographic measurements. The regularization constant can be arbitrarily small ( in our implementation). By taking the LS formulation approximates the exact pseudo-inverse .

The optimal solution to creftype 23 has a closed form expression given by

(26)

where is the identity operator. We note that due to the spatial symmetry of the kernel function , the operator is Hermitian, and the calculation of its adjoint , is synonymous to applying . We can therefore transform creftype 26 to the frequency domain and calculate it using creftype 25. The LS solution is then computed by element-wise multiplications and divisions on the 2D-DFT of the discrete kernel and PB measurements , given by

Thus, the optimal solution to creftype 23 has a closed form solution, computed in a single step of complexity.

Once we recover the coefficients from the given PB measurements , we can resample the continuous sinogram onto the PP grid, resulting in

(27)

Writing explicitly we have

where creftype 13 is a sampling operator that samples on the PP grid points creftype 14.

We note that the above formulation is completely general and is not limited to point-wise sampling, or a specific subspace. The sampling operator can be generalized so it further takes into account X-ray detector characteristics.

Iv-D Resampling Step Performance

To assess the performance of our resampling step, we perform a series of simulated tests and comparisons. A sheep phantom, with resolution of , is selected as our scanned object . The digital phantom was computing by scanning a sedated sheep with a Siemens SOMATOM Definition Flash CT scanner. The projection was performed with the high dosage of 480 mAs over 1250 projection angles, and is therefore chosen to be used as a phantom, due to its relatively high quality. The PB projection is simulated by the AIR-Tools v1.3 [23] package, generating either or projections with equispaced angles, each one composed from detectors. The synthetic sinogram

is then corrupted by Poisson noise with a standard deviation of

generating the noisy sinogram . The interpolation constants are chosen as , and , and the window constant is .

We compare our method to bilinear, bicubic and a modern spline based interpolation methods, where for reference we use the PPRT creftype 11 of the ground-truth phantom. Each of the regular interpolations are tested for resampling directly from both the noisy PB samples , and the denoised sinogram creftype 20. In this way, we can isolate whether our performance is gained by the subspace prior, or due to the generalized resampling algorithm, which solves (27). The interpolation results are measured in SNR and structural similarity index (SSIM), displayed in Fig. 8, confirming that our method successfully denoises the sinogram. We achieve better results than a two-step approach that first denoises the data by convolving with the kernel (20), and then use a conventional resampling technique. This confirms that our gains are achieved by both the defined subspace prior, and the generalized resampling approach.

Fig. 8: Comparison of different interpolation methods. The horizontal axis depicts the SNR of the inputted polar sinogram. The left graph demonstrates the performance of our resampling method compared to other interpolation methods. In the left graph we test the performance of our method compared to other approaches, where we use the denoised PB sinogram creftype 20 as input to the other methods. By examining the right graph it is apparent that the performance of our method is not gained only due to the denoising process, but rather by performing the entire resampling algorithm as depicted in creftype 27. Our approach achieves superior results for the entire range of input SNR.

V Reconstruction

For reconstructing our scanned object from the PP measurements we first consider the following TV optimization problem:

(28)

where is the PPRT operator, and is the Total-Variation [19] isotropic norm. This is a standard optimization approach for reconstruction, shown to produce superior results when coupled with the proposed PP resampling technique, than other state of the art tomographic reconstruction algorithms.

Before approaching creftype 28, we would like to introduce several modifications to it, that allow for faster and more accurate convergence, with better robustness to noise. First, the PP operator can be preconditioned by introducing the following operator

(29)

where is an element-wise operator described by

(30)

where . The operator is equivalent to applying a high-pass filter to each of the columns in the sinogram it operates on. More information on how we calculate an element-wise preconditioner to the PP system is found in Appendix C. When examining the combined operator , its condition number is significantly lower than that of , which leads to faster convergence when using first-order optimization techniques. A complete comparison between the condition number of the PB and the PP systems, with and without a preconditioner, is shown in Fig. 9. Since the computation of both and its adjoint , has complexity on the order of , the overall algorithm complexity remains the same.

Fig. 9: Comparison between system condition numbers. The algebraic condition number of various CT measurement systems is plotted as a function of the number of detectors (logarithmic scale). Preconditioning either the PP or the PB system by a diagonal matrix improves the condition number by more than an order of magnitude in each case. In either case, the PP system still has a significantly lower condition number. A condition number of is marked by a black line.

Next, we would like to incorporate a statistical model for the measurement noise. By considering a monochromatic radiation source, the measured X-ray data follows a Poisson distribution 

[24]. The raw projections measured by the CT detectors are distributed according to

(31)

where is the number of photons emitted towards the th detector, and is the average electrical bias in each detector. The line measurements are extracted by calculating the mean of the Poisson distribution in creftype 31

(32)

By assuming that each X-ray path is i.i.d with the others, the authors in [25]

have showed that the first order approximation to the maximum likelihood (ML) estimator that recovers

from is now given by a weighted LS [25]

(33)

where is an element-wise weighting operator dependent on the raw detector measurements , and given by

(34)

Modifying the quadratic norm term in creftype 28 together with the operator , helps to cope better with the additive Poisson noise, and achieve improved reconstructions.

In addition, we would like to reconstruct from a reduced number of measurements. To that end, we define a decimating operator as

(35)

where is the index set that contains the acquired measurements.

Combining creftype 30, creftype 34 and creftype 35, we can define a set of modified operators and a modified system matrix, such that

(36)

where contains a limited number of measurements, as described by . Since the operators and are all multiplying the sinogram (or its 1D-DFT in the case of ) element-wise, we simply take the square root of each of the multipliers to obtain and . Both operators, and , exhibit a computational complexity of , since they are composed from a composition of element-wise operators with complexity and the PPRT, with complexity .

By writing creftype 28 with the new operators creftype 36, the modified optimization problem is given by

(37)

For solving creftype 37, we use the FISTA-TV algorithm [10, 9] shown to produce excellent reconstructions from a reduced number of measurements under high noise levels.

The primary gradient step of the solver is given by

(38)

where is a constant, bigger or equal than the Lipschitz constant of the modified PP system

. This constant is computed by first approximating the largest singular value of

using the power method.

A summary of the entire algorithm is brought in Algorithm 1, and the FISTA-TV algorithm is presented in Algorithm 2.

1:Input: A tomographic PB scan, creftype 1, the kernel function , given in creftypeplural 18 and 16, a TV constant , and a regularization term
2:Output: Reconstructed scanned object
3:Find the underlying coefficients creftype 26:
4:
5:Calculate the sinogram on the PP grid by
6: creftype 27 such that
7:
8:Precondition the input measurements using the operators defined in creftypeplural 30, 35 and 34:
9:
10:Solve the following problem using Algorithm 2:
11:
12:where is defined in creftype 36
Algorithm 1 PP Tomography Resampling and Reconstruction
1: , ,
2:Initialize , and
3:for  do
4:     
5:     
6:     
7:end for
8:return
Algorithm 2 Fast Proximal Gradient Descent for RAPToR

Vi Simulation Setup and Results

To assess the performance of our algorithm, we perform a series of simulated tests and comparisons. We divide the tests into two parts, the first compares the quality of the resampling step, and the second evaluates the quality of the scanned object reconstruction.

Fig. 10: Comparison of algorithms for 60 (left) and 180 (right) viewing angles. The PSNR and SSIM metrics are compared for different state-of-the-art algorithms. Our method, which operates on the transformed PP measurements achieves significantly higher results when the input SNR is greater than .

We select the high dosage sheep phantom (digitized, with resolution ) as our scanned object . The PB projection is simulated by the AIR-Tools v1.3 [23] package, generating only projections with equispaced angles, each one composed from detectors. The synthetic sinogram is then corrupted by Poisson noise with a standard deviation of . The interpolation constants are chosen as in accordance with Fig. 7, and .

Fig. 11: Comparison of reconstruction results for 180 viewing angles with the same iteration budget. The reconstructions are presented here for the different algorithms. The reconstruction is performed from 180 viewing angles, with measurement noise of 28dB. (a) ground truth image, (b) FBP followed by TV denoising (FBP+TV), (c) Polar+TV, (d) ART+TV, (e) SPURS+TV, (f) Our method. All iterative methods were given the same iteration budget (of 8 iterations), demonstrating similar computational complexity.
Fig. 12: Comparison of reconstruction results for 120 viewing angles with full convergence. The reconstructions are presented here for the different algorithms. The reconstruction is performed from 120 viewing angles, with measurement noise of 24dB. (a) ground truth image, (b) FBP followed by TV denoising (FBP+TV), (c) Polar+TV, (d) ART+TV, (e) SPURS+TV, (f) Our method. All methods were allowed to run until convergence, specifically Polar+TV and ART+TV used 500 iterations.

For the reconstruction procedure we perform a simulated PB projection again with only 60 acquisition angles contaminated with noise, and then resample it to the PP domain using our method. This scenario is considered challenging for most state-of-the-art solvers available today. Throughout these simulations the TV coefficient is chosen as . After interpolating to the PP grid, the FISTA-TV algorithm runs for only 20 iterations until converging to the displayed results. The entire algorithm is compared to FBP, FBP followed by TV denoising, ART+TV [26], Polar+TV [27], and SPURS+TV [21]. The parameters in every algorithm were fine-tuned to run optimally, using an exhaustive grid search for each of them.

For SPURS+TV, we first used SPURS to resample the polar grid samples, given by the 1D-DFT of the PB projections , to the 2D-DFT of the object , recovering . In other words, we resampled in the frequency domain between a polar grid to a Cartesian one. To make the comparison fair, we then denoised the result by solving the following optimization problem:

(39)

where an optimal value for was searched and used.

The final results are shown in Fig. 12.The other solvers were allowed to run for a significantly higher number of iterations (), and a longer time until converging to a stable solution. To demonstrate the relative improvement in time complexity, we conducted an experiment that restricts the iterative methods to the same number of iterations as in our method. The results of this experiment are shown in Fig. 11, where it is evident that competing approaches are still far from a satisfactory solution. Consequently, we conclude that the time complexity of our combined solution (resampling and reconstruction) is significantly lower than compared solutions. Due to differences in the hardware implementation of each algorithm, we cannot provide an exact comparison of time complexity, except stating that the asymptotic complexity of our algorithm is slightly smaller than other approaches that has complexity of , while it runs for a number of iterations smaller by more than an order of magnitude.

For testing whether simply denoising the PB sinogram using the subspace prior has a profound effect on the results, we inputted a denoised sinogram into all other reconstruction approaches.

Vii Discussion

To conclude, we presented a robust reconstruction scheme for PB tomographic measurements, that performs better than state-of-the-art algorithms with a low computational complexity of . The reconstruction first relies on a resampling process that uses a prior condition on the sinogram, in the form of a subspace, and accurately transforms the PB sinogram to the PP grid, while reducing noise. The resampling is performed directly on the acquired sinogram measurements, and is shown to reduce noise and achieve high accuracy. It is shown to perform significantly better than known resampling techniques, either operating in the spatial or frequency domain. In addition, we showed that a PP-based reconstruction approach, that uses a modern solver which is fed with the resampled PP measurements, produces superior results than state-of-the-art algorithms, with a lower computational complexity and in less time. The method can also be adapted to modern helical and fan-beam scan methodologies, used in commercial CT scanners, as further discussed in Appendix A. We hope this work will pave the way for a CT scanner that uses the PP paradigm for producing clinical images that are taken with significantly less radiation dosage.

References

  • [1] J. Hsieh, Computed Tomography: Principles, Design, Artifacts, and Recent Advances.   Bellingham: WA: SPIE, 2003.
  • [2] A. C. Kak and M. Slaney, Principles of computerized tomographic imaging.   IEEE press, 1988.
  • [3] L. A. Feldkamp, L. C. Davis, and J. W. Kress, “Practical cone-beam algorithm,” Journal of the Optical Society of America A, vol. 1, no. 6, p. 612, 1984.
  • [4] A. Averbuch, R. R. Coifman, D. L. Donoho, M. Israeli, and J. Waldén, “Fast slant stack: A notion of Radon transform for data in a Cartesian grid which is rapidly computible, algebraically exact, geometrically faithful and invertible,” SIAM Journal on Scientific Computing, 2001.
  • [5] A. Averbuch, R. Coifman, and D. Donoho, “A framework for discrete integral transformations I - the pseudopolar Fourier transform,” SIAM Journal on Scientific Computing, vol. 30, no. 2, pp. 764–784, 2008.
  • [6] A. Averbuch, R. Coifman, D. Donoho, Y. Shkolnisky, and I. Sedelnikov, “A framework for discrete integral transformations II - the 2D discrete Radon transform,” SIAM Journal on Scientific Computing, vol. 30, no. 2, pp. 785–803, 2008.
  • [7] P. A. Rattey and A. G. Lindgren, “Sampling the 2-D Radon transform,” IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 29, no. 5, pp. 994–1002, 1981.
  • [8] Y. C. Eldar, Sampling Theory: Beyond Bandlimited Systems.   Cambridge University Press, 2015.
  • [9] A. Beck and M. Teboulle, “Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems,” IEEE Transactions on Image Processing, vol. 18, no. 11, pp. 2419–2434, 2009.
  • [10] D. P. Palomar and Y. C. Eldar, Convex Optimization in Signal Processing and Communications.   Cambridge university press,, 2010.
  • [11] E. Lee, B. P. Fahimian, C. V. Iancu, C. Suloway, G. E. Murphy, E. R. Wright, D. Castaño-Díez, G. J. Jensen, and J. Miao, “Radiation dose reduction and image enhancement in biological imaging through equally-sloped tomography,” Journal of Structural Biology, vol. 164, no. 2, pp. 221–227, 2008.
  • [12] B. P. Fahimian, Y. Mao, P. Cloetens, and J. Miao, “Low-dose x-ray phase-contrast and absorption CT using equally sloped tomography.” Physics in Medicine and Biology, vol. 55, no. 18, pp. 5383–5400, 2010.
  • [13] Y. Mao, B. P. Fahimian, S. J. Osher, and J. Miao, “Development and Optimization of Regularized Tomographic Reconstruction Algorithms,” IEEE Transactions on Image Processing, vol. 19, no. 5, pp. 1259–1268, 2010.
  • [14] S. Hashemi, S. Beheshti, and P. Gill, “Efficient Low Dose X-ray CT Reconstruction through Sparsity-Based MAP Modeling,” arXiv preprint, vol. arXiv:1402, pp. 1–10, 2014.
  • [15] H. Kudo, M. Defrise, M. Defrise, and I. Reconstruction, “Single-slice rebinning method for helical cone- beam CT,” Physics in Medicine & Biology, no. 44, pp. 561–570, 1999.
  • [16] K. Taguchi, G. L. Zeng, G. T. Gullberg, and M. Defrise, “3D cone-beam CT reconstruction for circular trajectories,” Physics in Medicine & Biology, no. 45, pp. 329–347, 2000.
  • [17] S. Schaller, F. Sauer, K. C. Tam, G. Lauritsch, and T. Flohr, “Exact radon rebinning algorithm for the long object problem in helical cone-beam CT.” IEEE transactions on medical imaging, vol. 19, no. 5, pp. 361–375, 2000.
  • [18] A. Louis, F. Natterer, and E. Todd Quinto, Mathematical Methods in Tomography, 2006.
  • [19] L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Physica D: Nonlinear Phenomena, vol. 60, no. 1-4, pp. 259–268, 1992.
  • [20] J. A. Fessler, S. Member, and B. P. Sutton, “Nonuniform Fast Fourier Transforms Using Min-Max Interpolation,” IEEE Transactions on signal processing, vol. 51, no. 2, pp. 560–574, 2003.
  • [21] A. Kiperwas, D. Rosenfeld, and Y. C. Eldar, “The SPURS Algorithm for Resampling an Irregularly Sampled Signal onto a Cartesian Grid,” IEEE transactions on medical imaging, vol. X, no. X, pp. 1–16, 2016.
  • [22] X. Pan, E. Y. Sidky, and M. Vannier, “Why do commercial CT scanners still employ traditional, filtered back-projection for image reconstruction?” Inverse Problems, vol. 25, no. 12, p. 123009, 2009.
  • [23] P. C. Hansen and M. Saxild-Hansen, “AIR tools - A MATLAB package of algebraic iterative reconstruction methods,” Journal of Computational and Applied Mathematics, vol. 236, no. 8, pp. 2167–2178, 2012.
  • [24] I. I. A. I. I. A. Elbakri and J. J. A. J. Fessler, “Statistical image reconstruction for polyenergetic X-ray computed tomography,” IEEE Transactions on Medical Imaging, vol. 21, no. 2, pp. 89–99, 2002.
  • [25] I. A. Elbakri and J. A. Fessler, “Efficient and accurate likelihood for iterative image reconstruction in {X}-ray computed tomography,” Proc. SPIE Med. Imag., vol. 5032, pp. 1839–1850, 2003.
  • [26] G. H. Chen, J. Tang, and S. Leng, “Prior image constrained compressed sensing (PICCS): a method to accurately reconstruct dynamic CT images from highly undersampled projection data sets.” Medical Physics, vol. 35, no. 2, pp. 660–663, 2008.
  • [27] E. Y. Sidky and X. Pan, “Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization.” Physics in Medicine and Biology, vol. 53, no. 17, pp. 4777–4807, 2008.
  • [28] C. C. Shaw, Cone Beam Computed Tomography.   CRC Press, 2014, vol. 1, no. 1.
  • [29] B. Alfano, M. Comerci, M. Larobina, A. Prinster, J. P. Hornak, S. E. Selvan, U. Amato, M. Quarantelli, G. Tedeschi, A. Brunetti, and M. Salvatore, “An MRI digital brain phantom for validation of segmentation methods,” Medical Image Analysis, vol. 15, no. 3, pp. 329–339, 2011.

Appendix A Application to Fan-Beam Acquisition

In this section, we examine the applicability of our method to fan-beam acquistion. Known algorithms for direct cone-beam or fan-beam reconstruction, namely variants of filtered back-projection (FBP) [2, 1] and the Feldkamp-Davis-Kress (FDK) [3] algorithms are known to suffer from inaccuracies due to their involved weighting process and interpolations, that compensate for the scan geometry. These phenomena and others are thoroughly discussed in [22, 28], where the advantages of iterative reconstruction algorithms are proposed and highlighted. In practice, due to these limitations, many commercial CTs employ a rebinning step on the measured cone-beam data. The rebinning process aims to transform the fan-beam projections to an approximately equivalent PB projection, enabling reconstruction by a standard FBP operating on the rebinned measurements.

The problem of accurately rebinning has been already addressed in the past [15, 16, 17, 18], and scanners today are reconstructing by employing various methods for accurate rebinning. Many devices use carefully calibrated look-up tables to enable on-the-fly rebinning while the scan is performed to generate PB measurements. After rebinning the reconstruction process is performed by a FBP.

To further asses the applicability of our method to rebinned PB data, we have conducted a simple experiment. In the experiment we have generated noisy fan-beam projection of the digital brain phantom [29] phantom, using the AIR-Tools v1.3 [23] CT simulator. We then rebinned the projections to PB equivalent measurements using a simple linear interpolation, and on those measurements we applied our resampling algorithm to transform them to pseudo-polar grid. For comparison, we reconstructed the object by direct FBP from the fan-beam projections, an FBP used on the rebinned PB projections and a eventually a FBP performed on the pseudo-polar resampled data.

As can be seen in Fig. 13, our method achieves better reconstructions than both fan-beam and parallel-beam FBP by simply applying the resampling step to the rebinned measurements. This result is maintained throughout a range of scanning profiles, including different SNR values and the number of angular projections. After a thorough analysis, we can conclude that even when performing rebinning to PB, which is a lossy transformation, the gains achieved by reconstructing with our method are still superior to other reconstruction approaches.

Fig. 13: RAPToR Performance Comparison. The numerical phantom is first projected using a fan-beam simulator. Results are compared for direct reconstruction using an FDK [3] based algorithm, a parallel-beam FBP performed on rebinned data and an equivalent pseudo-polar FBP performed on data resampled with RaPTOR.

Appendix B Derivation of The Subspace Kernel

To find the analytic expression for the kernel we perform a 2D inverse continuous Fourier transform over an indicator function of the region creftype 15:

Since this is a continuously differentiable function with respect to and , we can find its limits for each of the irregular points. For the case where , the limit is given by

The limits for the lines where are given by

For the second transition above we employed L’Hôpital’s rule. The limit for the line is computed by

For computing the limit of the point we choose the trajectory when :

thus concluding the derivation of the subspace kernel. Its derived expression is given in creftype 16.

Appendix C The Pseudo-Polar Preconditioner

We would like to find an algebraic preconditioner to the PP system, that can be computed with low computational complexity. The ideal algebraic preconditioner to the PP system is given by the following expression

(A1)

This preconditioner (if exists), leads to a condition number of for the pseudo-polar tomographic system, since by plugging it in we get:

where is the identity operator.

Fig. 14: The optimally calculated preconditioner (after applying the operator on it) is shown on the left for a small dimensional PP transform (). It was calculated by numerically evaluating creftype A1. More than of the energy is found on the main diagonal. On the right, for clarity we plotted only the first 4 repetitions from the diagonal of , showing that the preconditioning is the same for every column in the PP sinogram. For bigger dimensions of , we apply an element-wise version of the operator , as depicted in creftype 30.

For finding a preconditioner that efficiently approximates , we first numerically evaluate the following matrix inversion problem:

(A2)

Here and are the matrix form of the PPRT operators and its adjoint,

is the identity matrix, and

is a regularization constant chosen arbitrarily small ( in our case), that makes sure the inverse to creftype A2

uniquely exists. These operate on vectorized scanned objects

and PP sinograms accordingly. In practice, due to large dimensions of the matrix versions of the PP operators, the ideal preconditioner can only be evaluated for small dimensions. In our experiment, it was computed for object dimensions of pixels ().

Even if we could compute and store the ideal preconditioner for conventional dimensions, applying it in every iteration is a prohibitive step with computational complexity of . In an attempt to alleviate this computational burden, we apply the 1D-DFT operator to , such that

Most of the energy in is concentrated on its main diagonal, as can be seen in Fig. 14. In addition, by plotting the main diagonal of , we observe it is periodic with a periodicity equal to the number of X-Ray detectors of the examined system. We can therefore approximately substitute the ideal preconditioner by taking a diagonal matrix, with its main diagonal taken from and applying an inverse 1D-DFT operator on it. This procedure corresponds to filtering each of the PP projections for every acquisition angle with the same filter, and has computational complexity of .

To summarize, we numerically show that in the case of the PP system, a good approximation to the optimal algebraic preconditioner is given by a 1D filter that operates on the columns (detector axis wise) of the sinogram. The exact filter can be deduced by examining and extracting the diagonal of for small dimensions, for which its exact computation is numerically feasible.