Low radiation tomographic reconstruction with and without template information

12/23/2019 ∙ by Preeti Gopal, et al. ∙ Monash University IIT Bombay 33

Low-dose tomography is highly preferred in medical procedures for its reduced radiation risk when compared to standard-dose Computed Tomography (CT). However, the lower the intensity of X-rays, the higher the acquisition noise and hence the reconstructions suffer from artefacts. A large body of work has focussed on improving the algorithms to minimize these artefacts. In this work, we propose two new techniques, rescaled non-linear least squares and Poisson-Gaussian convolution, that reconstruct the underlying image making use of an accurate or near-accurate statistical model of the noise in the projections. We also propose a reconstruction method when prior knowledge of the underlying object is available in the form of templates. This is applicable to longitudinal studies wherein the same object is scanned multiple times to observe the changes that evolve in it over time. Our results on 3D data show that prior information can be used to compensate for the low-dose artefacts, and we demonstrate that it is possible to simultaneously prevent the prior from adversely biasing the reconstructions of new changes in the test object, via a method called “re-irradiation”. Additionally, we also present two techniques for automated tuning of the regularization parameters for tomographic inversion.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 13

page 14

page 21

page 22

page 25

This week in AI

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

1 Introduction

Reduction in radiation exposure is a critical goal, especially in CT of medical subjects Brenner and Hall (2007) and biological specimens Howells et al. (2009). One of the ways to reduce this radiation is to acquire projections from fewer views. An alternate way, which is the focus of this work, is to lower the strength (‘dose’) of X-ray beam. The CT imaging model that incorporates the strength of X-rays, , is non-linear and non-deterministic and is given by:

(1)

where

represents the zero mean additive Gaussian noise vector with a fixed signal-independent standard deviation

, where is the sensing matrix which represents the forward model for the tomographic projections, and is the underlying image representing the density values. The noise model for is primarily Poisson in nature as this is a photon counting process Dainty and Shaw (1974), and the added Gaussian noise is due to the thermal effects Xu and Tsui (2009). This Poisson-Gaussian noise model is quite common in optical or X-ray based imaging systems, but we consider it here explicitly for tomography, where it induces a non-linear inversion problem. Specifically, the index (for bin number and projection angle) in the measurement vector is given as: , where is the row of the sensing matrix

. The major effect of low-dose acquisition is the large magnitude (relative to the signal) of Poisson noise due to the low strength of X-ray beam. This is because the Signal-to-Noise-Ratio (SNR) of Poisson noise with mean

and variance

is given by . Due to the inherently low SNR, traditional low dose reconstructions are noisy.

2 Previous Work

Modelling of Poisson noise and recovery of images also finds applications in areas outside of CT. Zhang et al. (2017)

recovered images from Poisson-noise corrupted and blurred images using alternating direction method of multipliers(ADMM). Low-dose imaging and reconstruction (with dense projection view sampling) has been more widely studied than the few-views imaging. This is probably because the former does not involve a strategy for selection of the set of view angles, which in itself is an active field of research 

Barkan et al. (2017); Fischer et al. (2016); Dabravolski et al. (2014). For long, almost all of the commercial CT machines used FBP 111Filtered Backprojection as the standard reconstruction technique Pan et al. (2009). Only recently are the iterative techniques being deployed for commercial use Philips (news-article). The power of iterative routines was reinforced by Fujita et al. (2017), where it was proved that iterative reconstructions from ultra-low dose 222Typically, low-dose imaging is performed at 120 kVp and 30 mAs beam current, and ultra-low dose imaging is performed at 80-100 kVp and 20-30 mAs beam current settings. CT are of similar quality to those of FBP reconstructions from low-dose CT. Here, a commercial forward projected model-based algorithm was deployed and compared with FBP.

Among the other iterative methods, Kim et al. (2017)

presented a technique that minimizes log-likelihood of the Poisson distribution and a patch-based spatially encoded non-local penalty. 

Lyu et al. (2019) used a smoothness prior along with data-fidelity constraint and solved using ADMM. In order to improve the reconstruction further, various prior-based and learning-based methods have also been explored in literature. In these techniques, properties of available standard-dose CT images influence low-dose reconstruction of the test (i.e., the object which needs to be reconstructed from the current set of new tomographic projections). One such technique was described by Zheng et al. (2016), wherein the iterative reconstruction was formulated as a penalized weighted least squares problem with a pre-learned sparsifying transform. While the weights were set manually, the sparsifying transform was learned from a database of regular-dose CT images. Another technique presented by Jia et al. (2016)

clustered overlapping patches of previously scanned standard-dose CT images using Gaussian Mixture Model (GMM). The texture of the prior was learned for each cluster. Following this, patches from a pilot reconstruction of the test were classified using the learned GMM and depending on the class, the corresponding texture priors were imposed on patches of the reconstructed test image. The limitation here is– patches that correspond to new changes between the test and the templates will also be influenced by some inappropriate texture of patches from prior.

Zhang et al. (2016) solved a cost function with L1 norm for imposing similarity to a learned dictionary. They concluded that the number of measurements needed is progressively less for each of the four methods: Simultaneous Algebraic Reconstruction Technique (SART) Andersen and Kak (1984), Adaptive Dictionary based Statistical Iterative Reconstruction (ADSIR) Xu et al. (2012), Gradient Projection Barzilai Borwein (GBPP) Park et al. (2012) and their method (L1-DL), in the same order. Wu et al. (2013) used edge-based priors to reconstruct normal-dose CT along with Compressed Sensing (CS) sparsity prior. An iterative method Wang et al. (2012) in a related area (electrical impedance tomography) reconstructs using Split Bregman algorithm for L1 minimization. None of these methods explore optimizing a log-likelihood based cost-function that accurately reflects the Poisson-Gaussian noise statistics. In addition, they do not address the issue of the prior playing a role in the reconstruction of parts of the test that are dissimilar to the parts of the prior, which is undesirable. In contrast, this work focuses on applying a computationally fast global prior on only those regions of the test that are similar to the prior.

Lately, artificial neural networks have also been designed for low-dose reconstruction. 

Wu et al. (2017) proposed one such neural network to learn features of the image that is later imposed along with data-fidelity during iterative reconstruction. Shan et al. (2019) showed that deep neural network based reconstructions are faster than iterative reconstructions for comparable reconstruction quality. All of these neural-network based techniques need large amount of data. This can be challenging in longitudinal studies where usually only a few of the previous scans of the same object are available. Hence, this paper focuses on analytical iterative techniques.

We also present a technique for parameter selection. Most techniques in literature tune the parameters omnisciently. A recent work Gong and Zeng (2019) used the L-curve method in which data-fidelity residue is plotted against regularization norm. The parameter can then be selected based on the performance required for the application at hand. However, this method does not utilize the available information about noise statistics in low-dose imaging. In this work, we use the noise-model for the purpose of automated parameter selection.

3 Contributions

This paper discusses the following:

  1. How the quality of reconstruction is affected by the manner in which Poisson and Gaussian noise are modelled within iterative routines.

  2. How a prior of the object being scanned can be effectively used to compensate for the noisy low dose measurements, while simultaneously identifying genuine structural changes between the currently scanned objects and its priors, and preventing undesirable influence of the prior in those regions. We also propose a technique called re-irradiation which improves the reconstruction quality in these regions of change at the cost of a small amount of added radiation.

  3. In addition, most of the iterative schemes involve a cost function with a data-fidelity term, a sparsity term and a data-prior term. We discuss a systematic way to tune the parameters involved with these terms.

Specifically, this work presents a few new reconstruction methods and their comparison with existing methods, each of which model noise in a slightly different way. In addition to this, a technique for detecting new changes (i.e., differences between the test and templates) directly in the measurement space is presented. This is applied for prior-based reconstruction in longitudinal studies.

Sec. 4 describes two new techniques and its comparison with a few methods in literature. Sec. 5 describes the new prior-based low-dose reconstruction and its validation on 3D tomographic data. In Sec. 6, we illustrate a systematic technique to parameter-tuning. Finally, key results are summarized in Sec. 7

4 Reconstruction without prior

A good low-dose reconstruction technique should make optimal use of noise statistics as well as appropriate signal priors. Most techniques will involve minimizing a cost function of the following form: . Here the first term involves a data-fidelity cost, and may possibly (though not necessarily) be expressed by the negative log-likelihood of given and (i.e., by . Other alternatives could include a simple least squares term , or a weighted version of the same. The second term is a regularizer (weighted by the regularization parameter ) representing prior knowledge about . This could be in the form of the well-known total variation prior or penalty on the norm of the coefficients in a sparsifying basis where Such cost functions are minimized by iterative shrinkage and thresholding algorithms such as ISTA. However, ISTA by itself is known to have slow convergence (as discussed in Sec.3 of (Beck and Teboulle, 2009)). Hence, faster methods such as the Fast Iterative Soft Thresholding Algorithm (FISTA) (Beck and Teboulle, 2009) may be used, which is the method adopted in this paper. Below are some of the existing reconstruction methods, or intuitive variants thereof, and two new proposed techniques.

4.1 Post-log Compressed Sensing (CS)

A preliminary approach is to ignore the presence of Poisson noise and apply traditional CS reconstruction after linearizing the measurements (Hou and Zhang, 2014). The latter process is performed by computing the logarithm of the acquired measurements. The linearized measurements are given by , where is a small positive constant added to the measurements to make them all positive and thus suitable for linearizing by applying a logarithm. For practical purposes, if is zero or negative, is set to . The cost function is given by

(2)

is minimized by solver (Koh et al., last viewed–July, 2016). This method is however not true to the Poisson-Gaussian statistics and suffers from an inherent statistical bias (as seen in Fig. 1

) as it is a so-called ‘post-log’ method. The bias arises because for any non-negative random variable

, we have as per Jensen’s inequality. Another way of viewing this is that the noise in (i.e. post-log) is being treated as if it were Gaussian with a constant variance). This is not true except at very high intensity () value. The adverse effects of computing post-log measurements is also discussed in (Fu and others, 2016).

Figure 1: Histogram of statistical bias in post-log methods. The bias is computed as , where refers to linearized post-log measurements. Here, the added Gaussian noise had a mean value of and average Poisson-corrupted projection value. The fact that every bin has a different bias, but is shifted by a constant is problematic. This results in poor reconstructions, as shown in a later Sec. 4.8.

4.2 Non-linear Least Squares with CS

An intuitive way to modify the previous cost is by allowing the data fidelity cost to mimic the non-linearity inherent in the acquisition process. The cost function is then given by

(3)

The FISTA routine is used for this minimization. Since the attenuation constant of an object is never negative, a non-negativity constraint is imposed on . It can be seen that this cost function is non-convex in . Moreover, it treats all measurements as though they had the same noise variance, which is not true of Poisson settings.

4.3 Filtered Backprojection

In this technique, the classic filtered backprojection is applied on the linearized measurements: . The slice or volume is then reconstructed from the linearized measurements by filtered backprojection (FBP) in case of parallel beam projections or FeldKamp David Kress (FDK) algorithm (Feldkamp et al., 1984) in case of cone beam projections. This method is called the post-log FBP. While it is computationally efficient, it suffers from a statistical bias for the same reasons as post-log CS, as described in 4.1. The performance of post-log FBP has been extensively compared with iterative schemes in (Pontana et al., 2011),(Wang et al., 2013),(Koyama et al., 2014) and the latter has been found to be well suited for low-dose reconstructions (Willemink and Noël, 2019).

4.4 Negative Log Likelihood-Poisson with CS

This technique accounts for only the Poisson noise (ignoring the Gaussian part) and searches for a solution that minimizes the negative log-likelihood of the observed measurements. Given measurements, the likelihood of is defined as

(4)

where . Thus, the negative log likelihood of is given by

(5)

The cost function combines the likelihood and the CS term as shown below.

(6)

4.5 Negative Log Likelihood-Poisson-Gaussian with CS

A natural extension of the earlier method is one wherein both the Poisson and Gaussian noise processes are accounted for in the design of the cost function. Here, given the measurements, the solution that minimizes the sum of negative likelihood terms of both Poisson and Gaussian noise models, is selected. Let denote the Poisson random variable, i.e. . As seen earlier, the Poisson likelihood of is given by

(7)

where . Poisson negative log-likelihood of is given by

(8)

Next, if the assumed Gaussian noise has a variance of , then Gaussian likelihood of is given by

(9)

The Gaussian negative log-likelihood of is given by

(10)

We minimize the sum of the two negative log-likelihoods:

(11)

and are solved for alternately. Note that is integer-valued, but a typical gradient-based method will not restrict to remain in the domain of integers. For computational convenience, needs to be ‘softened’ to real values. Consequently must be replaced by the gamma function.

This cost function is non-convex. However it can be shown to be bi-convex, i.e., it is convex in if is kept fixed and vice versa. Such a cost-function was used in (Xie et al., 2017) as a method of pre-processing/denoising of projections prior to tomographic reconstruction. In contrast, we directly use it as a data-fidelity term for tomographic reconstruction. This appears more principled because denoising of a projection induces some ‘method noise’ which cannot be accurately modelled and which may affect subsequent reconstruction quality.

4.6 Proposed Rescaled non-linear Least Squares (RNLLS) with CS

This new method integrates Poisson noise model into the technique described in Sec.4.2. Since, the variance of a Poisson random variable is proportional to its mean, the variance of is directly proportional to . Hence the data-fidelity cost must be rescaled as shown below:

(12)

Again, the cost is minimized using FISTA solver. This technique is in some sense similar to the Penalized Weighted Least Squares (PWLS) technique from (Fessler, 1994) which seeks to minimize

(13)

where is a diagonal matrix of weights which are explicitly set (prior to running the optimization) based on the values in . In contrast to PWLS, in RNLLS, no weights are set as a prior. Rather the weights are equal to the underlying noiseless measurement values, and are explicitly inferred on the fly. In fact, a major motivation for our proposed technique is based on the fact that

(14)

This technique can be used for the case of Poisson-Gaussian noise as well, as in

(15)

We noticed that in (Ding et al., 2019), tomographic reconstruction was performed by minimizing the following cost function:

(16)

which is inspired by the approximation of by and treating it as a maximum quasi-likelihood problem. On the other hand, the proposed method (RNLLS) can be interpreted as a weighted form of the well-known LASSO problem (Hastie et al., 2015). We also note that the cost function for RNLLS is convex in the case of Poisson noise, as shown in the supplemental material. In the case of Poisson-Gaussian noise, our numerical simulations reveal that the cost function is not convex in the worst case. However, this non-convexity did not affect the numerical results significantly.

4.7 Proposed Poisson-Gaussian Convolution

This new technique models both the Poisson and Gaussian noise. It is based on the fact that if a random variable is the sum of two random variables and , then the density function of is given by the convolution of the density functions of and . This scheme has been used earlier (Chouzenoux et al., 2015) for image restoration from linear degradations such as blur, followed by Poisson-Gaussian corruption of the signal. In contrast, in CT, the measured signal is a non-linear function of the underlying image (i.e. its attenuation coefficients) as per Beer’s law. Eq. 17 refers to the Beer’s law along with the Poisson and Gaussian noise. The measurement is the sum of a Poisson random variable and a Gaussian random variable:

(17)

where = . The measurement is given as: , where . The probability density of the measurement is given by the following convolution:

(18)

The running variable does not take on negative values because the Poisson is a counting process and hence the corresponding random variable is always positive. Because all the measurements are independent (i.e., the noise in the sensor at any one pixel is independent of the noise at any other pixel on it), we have

(19)

The that maximizes the above probability needs to be computed. This is equivalent to minimizing the negative log-likelihood of the above probability. Hence, our cost function is given by

(20)

Since is computationally intractable for large , it has been approximated using Stirling’s approximation: . Further, in order to make the optimization numerically feasible, the value that takes for a particular measurement is limited to the range to where

is an integer that is usually set to 3. It is assumed here that some estimate of the variance

of the Gaussian noise is already known. This is usually feasible by recording the values sensed by the detector during an empty scan (without any object), usually before the actual scan is taken. Among the methods discussed here, the ones that model both Poisson and Gaussian noise are non-convex. A few of the methods that model Poisson noise alone are convex and their convexity is proved in Sec.1 of Supplementary material (video, proofs and comparison of methods).

4.8 Results on comparison of different methods

In order to compare the performance of various methods, 2D reconstructions of two datasets (Walnut and Colon CT) shown in Fig. 2 were computed for varying low-dose intensities. Reconstructions of two other datasets (Pelvis and Shoulder CT) are shown later in the supplemental material Supplementary material (video, proofs and comparison of methods). Following are the details of the datasets and the conditions used for simulating low-dose imaging: The size of the image from Walnut dataset was , and the size of image from Colon CT dataset was . The sum of the intensity values for the Walnut and Colon dataset images were and respectively. Measurements were simulated using parallel beam geometry. The Cosine filter was applied for filtered backprojection. While the number of projection views was large (200 views for all datasets) and kept constant, the beam strength was varied as follows: and . Based on the intensity (attenuation coefficients) of the images, the above values of correspond to a Poisson noise-to-signal ratio (i.e. average value of ) of for , and for , for both the datasets. In addition, Gaussian noise of mean and variance equal to of average Poisson-corrupted measurement was added to measurements. The regularization parameter was chosen omnisciently.

(a) walnut
(b) colon
Figure 2: Ground truth test slices used for comparison of low dose reconstruction techniques. A slice from (a) (Walnut, uni-muenster) dataset is of size , (b) (Colon, idash.ucsd.edu) dataset is of size
(a) Convolution
(b) Log-Likelihood Poisson-Gaussian
(c) Rescaled Non-Linear Least Squares
(d) Post-Log FBP
Figure 3: 2D Low-dose reconstructions of Walnut dataset for and . Gaussian noise of mean and variance equal to of average Poisson-corrupted measurement was added to simulate the low-dose acquisition. The SSIM values are shown in Fig. 5.
(a) Convolution
(b) Log-Likelihood Poisson-Gaussian
(c) Rescaled Non-Linear Least Squares
(d) Post-Log FBP
Figure 4: 2D Low-dose reconstructions of Colon dataset for and . Gaussian noise of mean and variance equal to of average Poisson-corrupted measurement was added to simulate the low-dose acquisition. The SSIM values are shown in Fig. 5.
Figure 5: SSIM of the reconstructions for Walnut and Colon datasets shown in Fig. 4 for varying values of X-ray doses. A higher SSIM implies better reconstruction. Here, the reconstructions by Poisson-likelihood and Poisson-Gaussian likelihood methods were very similar. Hence, their SSIM plots (blue and yellow respectively) overlap.

Sample reconstructions are shown in Figs. 3 and 4. The corresponding SSIM values of the reconstructions are shown in Fig. 5. From these plots, the following can be inferred: the convolution method and the Poisson-Gaussian likelihood reconstructions were comparable and gave the best reconstructions for a majority of dose levels and datasets.The Poisson-Gaussian Likelihood and the Poisson-only likelihood have very similar performance. However, at a theoretical level, the former is a more principled method, and can deal with negative-valued measurements which have to be weeded out for the Poisson-only method. The non-linear least squares method (Sec. 4.2) performed poorly. This is because the data-fidelity term assumes constant variance for all signal values. In reality, the variance of Poisson noise increases as signal intensity increases. The post-log linear least squares (Sec. 4.1) failed because the linear model fails to approximate the highly non-linear low-dose acquisition. The post-log FBP yielded poor results, especially at slightly higher dose levels (for example at in Fig. 3. This could be due to the absence of iterative optimization when compared to the other methods and due to the post-log approximation. For all datasets except Walnut (Colon as discussed here, and Pelvis, Shoulder as discussed in Supplementary material (video, proofs and comparison of methods)), the performance of rescaled non-linear least squares (RNLLS) is inbetween the performance of likelihood-based methods and those of all other methods. For the Walnut dataset though, the RNLLS gives the best quality for many dosage levels. The performance of the above methods across multiple noise instances is discussed in Sec.2.1 of Supplementary material (video, proofs and comparison of methods).

5 Reconstruction with prior

As seen so far, principled data fidelity terms play a significant role in improving the reconstruction performance. However, when the x-ray dose is less, the performance can be further improved by incorporation of useful priors Chen et al. (2008); Rashed and Kudo (2016). These priors could be previous high-quality reconstructions of the same object in longitudinal studies, or high-quality reconstructions of similar objects. We refer to such prior data as templates. Here, our aim is to reconstruct an object from its low-dose measurements, using templates which are previous high-dose reconstructions of the same object in a longitudinal study. However, there is a danger of the templates overwhelming the current reconstruction and adversely affecting reconstruction of new regions in the test (i.e., the object which needs to be reconstructed from the current set of new tomographic projections) that are absent in any of the templates. In the case of reconstruction from few projection views, the above problem was tackled Gopal et al. (2018)

by generating a map (known as ‘weights-map’) that shows an estimate of the regions of new changes and their magnitude. This map was then used to modulate the influence of the prior on the reconstruction of the test. The weights-map was computed based on the difference between the pilot reconstruction from the test measurements and its projection onto an eigenspace spanned by representative templates. However, in the low-dose case, this is not a preferable method because

all information about the noise model is valid for the measurement space alone. The noise model (i.e., ) is not applicable to the spatial reconstructed image domain.

Hence, in this work, we propose a new algorithm to compute the weights-map (i.e to detect differences between the test and the templates) directly in the measurement space. The aim is to identify those measurement bins which correspond to the new changes in the test. Following are the steps followed in order to accomplish this:

  1. Let be high quality template volumes, i.e. template volumes reconstructed from their standard dose measurements.

  2. Simulate noiseless measurements from template volumes using the same used for imaging the test i.e. , where .

  3. Let be the tomographic projection of the template from the angle, where . Let represent the set of eigenspaces, where is the eigenspace built from the tomographic projections of each of the templates in the angle, i.e. built from

  4. Let be the noisy tomographic projection of the test volume from the angle. For each , project onto , i.e., compute the eigen-coefficients of the measurements

    , along the set of eigenvectors

    :

    (21)

    where denotes the mean tomographic projection of all templates in the angle. The in the suffix denotes that the eigenspace is computed in the measurement space (We will contrast this with another eigen-space computed in image domain, used later in Eq. 23). Next, compute the resultant projection , i.e.,

    (22)
  5. Note that if a random variable , where , then is approximately distributed as . The quality of the approximation is known to improve as increases. In the absence of Gaussian noise (equivalent to the case where ), this transform is called the Anscombe transform (Anscombe, 1948; Curtiss, 1943), and has been widely used in image processing. In the presence of Gaussian noise, it is referred to as the generalized Anscombe transform (Murtagh et al., 1995). Now consider the bin in the test measurement as well as in , which we shall denote as and respectively. If represents the same underlying structure as in , barring the effect of Poisson-Gaussian noise, i.e. if the bin in is not part of the ‘new changes’, then the following is true:

    .

    For bins falling in the regions of change in the test (compared to the template projections), the above hypothesis is false. The same argument can be extended for entire segments or 2D regions.

  6. Based on the aforementioned fact, hypothesis testing is performed on to detect bins corresponding to new changes in the measurement space. We use Z-test for hypothesis testing on 2D patches in the measurement space (note that since the volume is in 3D, the measurement space is in 2D for every imaging view). This

    test computes the probability that the given sample is likely to be drawn from a population as specified by the null hypothesis. In this case, the null hypothesis is that the intensity values of the patches are drawn from

    . The confidence level was set to , i.e. for null hypothesis to be false, the probability

    that the sample is drawn from Normal distribution must lie in the

    tail-end of the Normal distribution on either side. A lower -value denotes the presence of new changes i.e., presence of differences between the test and the templates in the measurement bins.

  7. Once the new changes are detected in the measurement space, filtered backprojection of the vectors (containing p-values) resulting from the hypothesis test gives the location of the new changes (which we denote ) in the original (3D) spatial domain. The Cosine filter was used in the filtered backprojection process.

  8. The final weights-map 333An alternate method to compute a weights-map (a simpler binary weights-map) is discussed in Sec.3 of Supplementary material (video, proofs and comparison of methods) is computed from by the following steps: (a) Inversion: . This step is just for inversion so that new regions get lower weight/intensity than prior-similar regions, (b) Linear stretching: Perform linear stretching on so that the weights lie between 0 and 1.

Finally, the computed weights-map is used in a reconstruction optimization as follows:

(23)

where the eigenvectors and mean of the templates form the eigenspace which is built from the available high-dose reconstructions of the templates. denotes the coefficients of , when the pilot reconstruction of the test is projected onto this eigenspace. Information about the location and magnitude of new changes in the test is present in the weights-map . Eq. 23 is solved by alternating minimization on and until convergence is reached.

5.1 Reconstruction results

(a)
(b)
(c)
(d)
Figure 6: Potato 3D dataset: One of the slices from template volumes (first four from the left) and test volume (extreme right). Size of each volume is [].

The above algorithm was validated by reconstructing a 3D volume from its low dose measurements. Fig. 6 shows a slice from each of the template and test volumes of the potato dataset. This dataset 444We are grateful to Dr. Andrew Kingston for facilitating data collection at the Australian National University. consisted of four scans of the humble potato, chosen for its simplicity. Measurements from each scan consisted of cone-beam projections from 900 views, each of size . The corresponding size of the reconstructed volume is . While the first scan was taken of the undistorted potato, subsequent scans were taken of the same specimen, each time after drilling a new hole halfway into the potato. The ground truth consists of FDK reconstructions from the full set of acquired measurements from projection views. Low dose cone-beam measurements were simulated from full-view FDK reconstructions of the test volume. was set to 4000, a value corresponding to Poisson noise of . Mean of the added Gaussian noise was and was set to of the mean of Poisson-corrupted measurements. Fig 7 shows the same slice from each of the reconstructed volumes. A patch size of was used for hypothesis testing and the location of new changes (marked in red RoI in test) was accurately detected in the weights-map as seen in Fig. 6(e). The reconstructed volumes can be found in (Supplementary material, video, proofs and comparison of methods).

(a) Test
(b) No prior
(c) Unweighted
(d) Our
reconstruction
(e) Weights
Figure 7: Prior-based low-dose reconstruction on 3D potato dataset. (a) Slice from test volume (b) Reconstruction using no prior (using RNLLS of Sec. 4.6); (c) Slice from unweighted prior reconstruction; . The new change is missing. (d) Slice from weighted prior reconstruction; . The new change is detected here and its reconstruction is guided by the low-dose measurements. (e) Weights map showing the location and intensity of the new changes. All SSIM values are averaged over 14 slices of the reconstructed volume in the red RoI region. The reconstructed volumes can be seen in (Supplementary material, video, proofs and comparison of methods).

5.2 Re-irradiation to improve reconstruction

Once the regions of new changes are detected by the weights map, this information can be used to re-irradiate them with standard-dose rays and further improve the quality of their reconstruction. Following are the steps of the re-irradiation process:

  1. Let the X-rays passing through the new regions have their source points denoted by , and the corresponding bins at the detector be denoted by . Let the X-rays passing through the other regions (i.e. regions where the test and the templates are not structurally different) have their source points denoted by , and the corresponding bins at the detector be denoted by .

  2. Block and re-irradiate the object by passing standard-dose rays from . This will generate measurements of high quality for regions of new changes. If the regions of new change are small in area, this process incurs only a small cost for the extra amount of radiation, since the latter is restricted to only specific regions.

  3. In the measurement matrix captured for pilot reconstruction, replace all the bins in set by their new measurements. Therefore, the final measurement matrix consists of standard-dose measurements corresponding to new regions of the object and low-dose measurements corresponding to the other regions of the object.

The new measurement model is: . Here denotes a diagonal matrix with denoting the strength of the X-ray incident on the bin of the sensor. Fig. 8 shows the templates and test images, and Fig. 9 shows the reconstructions illustrating the benefit of re-irradiation. The new changes within the RoI are reconstructed very well after they are re-imaged with standard-dose X-rays. This is also reinforced by results on the sprouts data (Fig. 10), shown in Fig. 11. The selection of bins for re-irradiation and the choice of new X-ray intensity can also be chosen in a supervised manner by the physician or scientist based on the particular clinical or non-clinical setting.

(a)
(b)
(c)
(d)
(e)
Figure 8: Dataset for illustrating re-irradiation: Templates (first four from the left) and test (extreme right). Size of each slice is . The RoI shows the region of difference between the test and the templates.
(a) Test
(b) Pilot
(c) weights
(d) Weighted
Prior
(e) After
re-irradiation
Figure 9: Improving reconstruction by re-irradiation in Okra 2D dataset. (a) test (b) pilot (c) weights-map; the lower the intensity, the higher the magnitude of new changes. (d) weighted prior reconstruction; the quality of reconstruction of new regions is poor because it is guided by the measurements alone. (e) re-irradiated reconstruction; new measurements with twice the earlier low-dose X-ray intensity at of the bins enable better reconstruction of new regions (as shown in RoI).
(a)
(b)
(c)
(d)
(e)
(f)
Figure 10: Sprouts Dataset for illustrating re-irradiation: Templates (first row) and test (second row). Size of each slice is . The RoI shows the region of difference between the test and the templates.
(a) Test
(b) Pilot
(c) weights
(d) Weighted
Prior
(e) After
re-irradiation
Figure 11: Improving reconstruction by re-irradiation in Sprouts 2D dataset. (a) test (b) pilot (c) weights-map; the lower the intensity, the higher the magnitude of new changes. (d) weighted prior reconstruction; the quality of reconstruction of new regions is poor because it is guided by the measurements alone. (e) re-irradiated reconstruction; new measurements with times the earlier low-dose X-ray intensity at of the bins enable better reconstruction of new regions (as shown in RoI).

6 Tuning of parameters

Two parameters were used in the techniques presented in this chapter: : weight for CS term and : weight for object-prior. Below are few of the ways to select these parameters.

6.1 Selection of weightage for CS term

In a large body of work on tomographic reconstruction (Zheng et al., 2016)(Liu et al., 2016), the regularization parameter is chosen in an “omniscient fashion”. That is, the optimization problem is solved separately for many different values of . The particular result which yields the least MSE with respect to a ground truth image is chosen to be the correct result. Such a method requires knowledge of the ground truth, and hence is infeasible in practice. Other alternatives include visual inspection or cross-validation. However none of these techniques are fully practical. Instead, we propose a method to choose based on sound statistical principles pertaining to the Poisson or the Poisson-Gaussian noise model. The method is shown here in conjunction with the rescaled non-linear least squares method, however in principle, it can be used with any data fidelity term. For the Poisson-Gaussian noise model, the cost function is given by .

(a)
(b)
(c)
(d)
Figure 12: Mean and variance of the data-fidelity term for different number of measurements (projection views) and beam strength . (a) Expected value of exactly coincides with , (b) Variance of is insignificant for any number of measurements, (c) mean of is independent of the beam-strength, and (d) Variance of is insignificant for all values.

Let denote the total number of bins, the reconstruction with optimal . Let . Clearly, we have . Hence we can state that . Furthermore, our simulations (Fig. 12) have shown that

(24)

where denotes element-wise division. We also observed that the variance of the above quantity is very small. This is illustrated in Fig. 12, which shows that the variance of is very small compared to its mean. The expected value of varies with the number of measurements (is equal to , and is independent of . Hence we conclude that the quantity should be as close to as possible. Therefore, we consider

(25)
Figure 13: A method to choose the parameter in low-dose reconstruction: We expect to be minimum at the same for which relative MSE is minimum. Here, the for which and relative MSE are minimum are very close. Refer to Fig. 14 to observe the reconstruction results for different values of .

and observe how and relative MSE of reconstructions vary for different values of . At the optimum , must be minimum. The test image () and the reconstructions are shown in Figure 14. For these reconstructions, 410 projection views were chosen and Gaussian noise = was added to the measurements. The dose of X-rays resulted in a Poisson NSR of 0.018. As shown in Fig. 13, the for which and relative MSE are minimum, are very close. In a real-life setting, when relative MSE cannot be computed because of absence of ground-truth, a brute force search needs to be done followed by selecting the value of that minimizes .

(a) Test
(b)
(c)
(d)
(e)
(f)
(g)

(h)
(i)
(j)

(k)
(l)
(m)
(n)
(o)
Figure 14: Colon test data and its reconstructions for different values of . is minimum for , shown in green, with a relative MSE of . The reconstruction for , shown in red, gives the minimum relative MSE of .

6.2 Selection of weightage for object-prior term

The weightage for the object prior, term needs to be chosen omnisciently for every dataset. However, for a variation of , there was no significant effect on the reconstructions. Lower values indicate that the reconstructions are primarily guided by the measurements, and higher values will strengthen the effect of the prior.

7 Conclusions

In the low-dose CT imaging regime, the noise in the measurements becomes significant and needs to be accounted for during the reconstruction. Two new techniques: Poisson-Gaussian convolution and rescaled non-linear least squares (RNLLS) were presented and extensively compared with many of the existing methods. RNLLS was further used in low-dose reconstruction for longitudinal studies to specifically detect new regions in the test and simultaneously reduce noise in the other reconstructed regions. The results were validated on both 2D and 3D biological data. We demonstrated that the reconstructions of the regions of new changes can be significantly improved by re-irradiating these specific regions by standard-dose X-rays. Further, different methods for choosing the parameters were also discussed, which has not been dealt with in literature. Our technique can possibly be extended to the case where templates of a similar class of objects are available, as against previous scans of the same object. This may further increase the utility of the technique in clinical settings.

References

  • A.H. Andersen and A.C. Kak (1984) Simultaneous algebraic reconstruction technique (SART): a superior implementation of the ART algorithm. Ultrasonic Imaging 6 (1), pp. 81 – 94. External Links: ISSN 0161-7346 Cited by: §2.
  • F. J. Anscombe (1948) The transformation of poisson, binomial and negative-binomial data. Biometrika 35 (3/4), pp. 246–254. External Links: ISSN 00063444 Cited by: item 5.
  • O. Barkan, J. Weill, S. Dekel, and A. Averbuch (2017) A mathematical model for adaptive computed tomography sensing. IEEE Transactions on Computational Imaging 3 (4), pp. 551–565. Cited by: §2.
  • A. Beck and M. Teboulle (2009) A fast iterative shrinkage-thresholding algorithm for linear inverse problems. Imaging Sciences 2 (1), pp. 183–202. Cited by: §4.
  • D. J. Brenner and E. J. Hall (2007) Computed tomography — an increasing source of radiation exposure. New England Journal of Medicine 357 (22), pp. 2277–2284. Cited by: §1.
  • G. Chen, J. Tang, and S. Leng (2008) Prior image constrained compressed sensing (PICCS): a method to accurately reconstruct dynamic CT images from highly undersampled projection data sets. Medical Physics 35 (2), pp. 660–663. Cited by: §5.
  • E. Chouzenoux, A. Jezierska, J. Pesquet, and H. Talbot (2015) A convex approach for image restoration with exact poisson–gaussian likelihood. SIAM Journal on Imaging Sciences 8 (4), pp. 2662–2682. Cited by: §4.7.
  • Colon (idash.ucsd.edu) CT Colonography. Note: https://idash.ucsd.edu/data-collections Cited by: Figure 2.
  • J. H. Curtiss (1943) On transformations used in the analysis of variance. Ann. Math. Statist. 14 (2), pp. 107–122. Cited by: item 5.
  • A. Dabravolski, K. J. Batenburg, and J. Sijbers (2014) Dynamic angle selection in x-ray computed tomography. Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 324, pp. 17 – 24. Note: 1st International Conference on Tomography of Materials and Structures Cited by: §2.
  • J. C. Dainty and R. Shaw (1974) Image science: principles, analysis and evaluation of photographic-type imaging processes. Academic Press. Cited by: §1.
  • Q. Ding, Y. Long, X. Zhang, and J. A. Fessler (2019) Statistical image reconstruction using mixed Poisson-Gaussian noise model for X-ray CT. submitted in Inverse Prob. and Imaging, axXiv. Note: Submitted Cited by: §4.6.
  • L. Feldkamp, L. C. Davis, and J. Kress (1984) Practical cone-beam algorithm. J. Opt. Soc. Am 1, pp. 612–619. Cited by: §4.3.
  • J. A. Fessler (1994) Penalized weighted least-squares image reconstruction for positron emission tomography. IEEE Transactions on Med. Imaging 13 (2), pp. 290–300. Cited by: §4.6.
  • A. Fischer, T. Lasser, M. Schrapp, J. Stephan, and P. B. Noël (2016) Object specific trajectory optimization for industrial X-ray computed tomography. Scientific Reports 6 (19135). Cited by: §2.
  • L. Fu et al. (2016) Comparison between pre-log and post-log statistical models in ultra-low-dose CT reconstruction. IEEE transactions on medical imaging 36 (3), pp. 707 – 720. Cited by: §4.1.
  • M. Fujita, T. Higaki, Y. Awaya, T. Nakanishi, Y. Nakamura, F. Tatsugami, Y. Baba, M. Iida, and K. Awai (2017) Lung cancer screening with ultra-low dose ct using full iterative reconstruction. Japanese Journal of Radiology 35 (4), pp. 179–189. External Links: ISSN 1867-108X Cited by: §2.
  • C. Gong and L. Zeng (2019) Adaptive iterative reconstruction based on relative total variation for low-intensity computed tomography. Signal Processing 165, pp. 149 – 162. Cited by: §2.
  • P. Gopal, S. Chandran, I. D. Svalbe, and A. Rajwade (2018) Learning from past scans: tomographic reconstruction to detect new structures. CoRR abs/1812.10998. Cited by: §5.
  • T. Hastie, R. Tibshirani, and M. Wainwright (2015) Statistical learning with sparsity. CRC Press, Taylor & Francis Group. Cited by: §4.6.
  • W. Hou and C. Zhang (2014) A compressed sensing approach to low-radiation CT reconstruction. 2014 9th International Symposium on Communication Systems, Networks Digital Sign (CSNDSP), pp. 793–797. Cited by: §4.1.
  • M.R. Howells, T. Beetz, H.N. Chapman, C. Cui, J.M. Holton, C.J. Jacobsen, J. Kirz, E. Lima, S. Marchesini, H. Miao, D. Sayre, D.A. Shapiro, J.C.H. Spence, and D. Starodub (2009) An assessment of the resolution limitation due to radiation-damage in X-ray diffraction microscopy. Journal of Electron Spectroscopy and Related Phenomena 170 (1), pp. 4 – 12. Cited by: §1.
  • X. Jia, Z. Bian, J. He, Y. Wang, J. Huang, D. Zeng, Z. Liang, and J. Ma (2016) Texture-preserved low-dose CT reconstruction using region recognizable patch-priors from previous normal-dose CT images. In IEEE Nuclear Science Symposium, Vol. , pp. 1–4. Cited by: §2.
  • K. Kim, G. El Fakhri, and Q. Li (2017) Low-dose CT reconstruction using spatially encoded nonlocal penalty. Medical physics 44 (10), pp. e376–e390. Cited by: §2.
  • K. Koh, S.-J. Kim, and S. Boyd (last viewed–July, 2016) L1-ls: simple matlab solver for l1-regularized least squares problems. Note: https://stanford.edu/~boyd/l1_ls/ Cited by: §4.1.
  • H. Koyama, Y. Ohno, M. Nishio, S. Matsumoto, N. Sugihara, T. Yoshikawa, S. Seki, and K. Sugimura (2014) Iterative reconstruction technique vs filter back projection: utility for quantitative bronchial assessment on low-dose thin-section MDCT in patients with/without chronic obstructive pulmonary disease. European Radiology 24 (8), pp. 1860–1867. External Links: ISSN 1432-1084 Cited by: §4.3.
  • J. Liu, Y. Hu, J. Yang, Y. Chen, H. Shu, L. Luo, Q. Feng, Z. Gui, and G. Coatrieux (2016) 3D feature constrained reconstruction for low dose CT imaging. IEEE Transactions on Circuits and Systems for Video Technology PP (99), pp. 1–1. Cited by: §6.1.
  • Q. Lyu, D. Ruan, J. Hoffman, R. Neph, M. McNitt-Gray, and K. Sheng (2019) Iterative reconstruction for low dose CT using plug-and-play alternating direction method of multipliers (ADMM) framework. SPIE Medical Imaging 10949. Cited by: §2.
  • F. Murtagh, J. Starck, and A. Bijaoui (1995) Image restoration with noise suppression using a multiresolution support. Astronomy and Astrophysics, Suppl. Ser 112, pp. 179–189. Cited by: item 5.
  • X. Pan, E. Y. Sidky, and M. Vannier (2009) Why do commercial CT scanners still employ traditional, filtered back-projection for image reconstruction?. Inverse problems 25 (12), pp. . Cited by: §2.
  • J. C. Park, B. Song, J. S. Kim, S. H. Park, H. K. Kim, Z. Liu, T. S. Suh, and W. Y. Song (2012) Fast compressed sensing-based CBCT reconstruction using barzilai-borwein formulation for application to on-line IGRT. Med. Phy. 39 (3), pp. 1207–1217. Cited by: §2.
  • Philips (news-article) Philips IMR offers new capabilities to simultaneously reduce CT radiation and enhance image quality. Note: https://www.philips.com/a-w/about/news/archive/standard/news/press/2013/20130617-Philips-IMR-offers-new-capabilities.html Cited by: §2.
  • F. Pontana, A. Duhamel, J. Pagniez, T. Flohr, J. Faivre, A. Hachulla, J. Remy, and M. Remy-Jardin (2011) Chest computed tomography using iterative reconstruction vs filtered back projection (part 2): image quality of low-dose ct examinations in 80 patients. European Radiology 21 (3), pp. 636–643. External Links: ISSN 1432-1084 Cited by: §4.3.
  • E. A. Rashed and H. Kudo (2016) Probabilistic atlas prior for CT image reconstruction. Computer Methods and Programs in Biomedicine 128, pp. 119 – 136. Cited by: §5.
  • H. Shan, A. Padole, F. Homayounieh, U. Kruger, R. D. Khera, C. Nitiwarangkul, M. K. Kalra, and G. Wang (2019) Competitive performance of a modularized deep neural network compared to commercial algorithms for low-dose CT image reconstruction. Nature Machine Intelligence 1 (6), pp. 269–276. Cited by: §2.
  • Supplementary material (video, proofs and comparison of methods) videos, more results. Note: https://www.dropbox.com/s/x4jwzjpcw8x7dz7/supplementary_material.pdf?dl=0 Cited by: §4.7, §4.8, §4.8, Figure 7, §5.1, footnote 3.
  • Walnut (uni-muenster) Walnut CT dataset. Note: https://www.uni-muenster.de/Voreen/download/workspaces_and_data_sets.html Cited by: Figure 2.
  • H. Wang, B. Tan, B. Zhao, C. Liang, and Z. Xu (2013) Raw-data-based iterative reconstruction versus filtered back projection: image quality of low-dose chest computed tomography examinations in 87 patients. Clin. Imaging 37 (6), pp. 1024 – 32. External Links: ISSN 0899-7071 Cited by: §4.3.
  • J. Wang, J. Ma, B. Han, and Q. Li (2012) Split Bregman iterative algorithm for sparse reconstruction of electrical impedance tomography. Signal Processing 92 (12), pp. 2952 – 2961. Cited by: §2.
  • M. J. Willemink and P. B. Noël (2019)

    The evolution of image reconstruction for CT—from filtered back projection to artificial intelligence

    .
    European Radiology 29 (5), pp. 2185–2195. External Links: ISSN 1432-1084 Cited by: §4.3.
  • D. Wu, K. Kim, G. El Fakhri, and Q. Li (2017) Iterative low-dose CT reconstruction with priors trained by artificial neural network. IEEE Transactions on Medical Imaging 36 (12), pp. 2479–2486. External Links: ISSN 0278-0062 Cited by: §2.
  • J. Wu, F. Liu, L. Jiao, and X. Wang (2013) Multivariate pursuit image reconstruction using prior information beyond sparsity. Signal Processing 93 (6), pp. 1662 – 1672. Note:

    Special issue on Machine Learning in Intelligent Image Processing

    Cited by: §2.
  • Q. Xie, D. Zeng, Q. Zhao, D. Meng, Z. Xu, Z. Liang, and J. Ma (2017) Robust low-dose CT sinogram preprocessing via exploiting noise-generating mechanism. IEEE Transactions on Medical Imaging 36 (12), pp. 2487–2498. Cited by: §4.5.
  • J. Xu and B. M. W. Tsui (2009) Electronic noise modeling in statistical iterative reconstruction. IEEE Transactions on Image Processing 18 (6), pp. 1228–1238. Cited by: §1.
  • Q. Xu, H. Yu, X. Mou, L. Zhang, J. Hsieh, and G. Wang (2012) Low-dose X-ray CT reconstruction via dictionary learning. IEEE Transactions on Medical Imaging 31 (9), pp. 1682–1697. Cited by: §2.
  • C. Zhang, T. Zhang, M. Li, C. Peng, Z. Liu, and J. Zheng (2016) Low-dose CT reconstruction via L1 dictionary learning regularization using iteratively reweighted least-squares. BioMedical Engineering OnLine 15 (1), pp. 66. Cited by: §2.
  • H. Zhang, Y. Dong, and Q. Fan (2017) Wavelet frame based poisson noise removal and image deblurring. Signal Processing 137, pp. 363 – 372. Cited by: §2.
  • X. Zheng, Z. Lu, S. Ravishankar, Y. Long, and J. A. Fessler (2016) Low dose CT image reconstruction with learned sparsifying transform. In 2016 IEEE 12th Image, Video, and Multidimensional Signal Processing Workshop, Vol. , pp. 1–5. Cited by: §2, §6.1.