Low radiation tomographic reconstruction with and without template information

2019·Arxiv

Abstract

Abstract

duced 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.

Keywords: low-dose tomographic reconstruction, compressed sensing, priors,

longitudinal studies.

1. Introduction

subjects [1] and biological specimens [2]. 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:

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 x is the underlying image representing the density values. The noise model for y is primarily Poisson in nature as this is a photon counting process [3], and the added Gaussian noise is due to the thermal effects [4]. 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 y is given as: Poisson(exp) + , 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

in areas outside of CT. [5] recovered images from Poisson-noise corrupted and blurred images using alternating direction method of multipliers(ADMM). Lowdose 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 [6, 7, 8]. For long, almost all of the commercial CT machines used FBP as the standard reconstruction technique [9]. Only recently are the iterative techniques being deployed for commercial use [10]. The power of iterative routines was reinforced by [11], where it was proved that iterative reconstructions from ultra-low dose 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.

imizes log-likelihood of the Poisson distribution and a patch-based spatially encoded non-local penalty. [13] 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 [14], 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 [15] 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. [16] 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) [17], Adaptive Dictionary based Statistical Iterative Reconstruction (ADSIR) [18], Gradient Projection Barzilai Borwein (GBPP) [19] and their method (L1-DL), in the same order. [20] used edge-based priors to reconstruct normal-dose CT along with Compressed Sensing (CS) sparsity prior. An iterative method [21] in a related area (electrical impedance tomography) reconstructs using Split Breg-

man 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.

construction. [22] proposed one such neural network to learn features of the image that is later imposed along with data-fidelity during iterative reconstruction. [23] 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 chal-

lenging 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.

literature tune the parameters omnisciently. A recent work [24] used the Lcurve 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

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.

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

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 y given and x (i.e., by log )). Other alternatives could include a simple least squares term , or a weighted version of the same. The second term R(x) is a regularizer (weighted by the regularization parameter ) representing prior knowledge about x. This could be in the form of the well-known total variation prior TV (x) = + 1+ (x(i, j + 1) 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 [25]). Hence, faster methods such as the Fast Iterative Soft Thresholding Algorithm (FISTA) [25] 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)

traditional CS reconstruction after linearizing the measurements [26]. The latter process is performed by computing the logarithm of the acquired measurements.

The linearized measurements are given by log

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 min(y) is zero or negative, is set to min(y) + 0.001. The cost function is

given by

is minimized by solver [27]. 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 X, we have log((log(X)) as

Figure 1: Histogram of statistical bias in post-log methods. The bias is computed as (refers to linearized post-log measurements. Here, the added Gaussian noise had a mean value of 0 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.

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 [28].

4.2. Non-linear Least Squares with CS

fidelity cost to mimic the non-linearity inherent in the acquisition process. The

cost function is then given by

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

earized measurements: log = . The slice or volume x is then re- constructed from the linearized measurements by filtered backprojection (FBP) in case of parallel beam projections or FeldKamp David Kress (FDK) algorithm [29] 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 [30],[31],[32] and the latter has been found to be well suited for low-dose reconstructions [33].

4.4. Negative Log Likelihood-Poisson with CS

part) and searches for a solution that minimizes the negative log-likelihood of the observed measurements. Given m measurements, the likelihood of is defined

as

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

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

) =

4.5. Negative Log Likelihood-Poisson-Gaussian with CS

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 V denote the Poisson random variable, i.e. . As seen earlier, the

Poisson likelihood of is given by

where . Poisson negative log-likelihood of is given by

Next, if the assumed Gaussian noise has a variance of , then Gaussian likeli-

hood of is given by

The Gaussian negative log-likelihood of is given by

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

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

i.e., it is convex in if v is kept fixed and vice versa. Such a cost-function was used in [34] 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

in Sec.4.2. Since, the variance of a Poisson random variable is proportional to its mean, the variance of y is directly proportional to exp(). Hence the data-fidelity cost must be rescaled as shown below:

Again, the cost is minimized using FISTA solver. This technique is in some sense similar to the Penalized Weighted Least Squares (PWLS) technique from [35]

which seeks to minimize

where W is a diagonal matrix of weights which are explicitly set (prior to running the optimization) based on the values in y. 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

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

We noticed that in [36], tomographic reconstruction was performed by minimizing the following cost function:

which is inspired by the approximation of Poisson(z) by N(z, z) 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 [37]. 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

on the fact that if a random variable Q is the sum of two random variables R and S, then the density function of Q is given by the convolution of the density functions of R and S. This scheme has been used earlier [38] 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:

where . The measurement is given as: Poisson() + , where . The probability density of the measurement is given by the following convolution:

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 m 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

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

Since l! is computationally intractable for large l, it has been approximated using Stirling’s approximation: . Further, in order to make the optimization numerically feasible, the value that l takes for a particular measurement is limited to the range max(0) to where K 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 [39].

4.8. Results on comparison of different methods

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 [39]. 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 156 156, and the size of image from Colon CT dataset was 154 154. The sum of the intensity values for the Walnut and Colon dataset images were 75 and 60 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: = 20, 40, 80, 160, 320 and 620. 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 1) of 25% for = 20, and 4.5% for = 620, for both the datasets. In addition, Gaussian noise of 0 mean and variance equal to 2% of average Poisson-corrupted measurement was added to measurements. The regularization parameter was chosen omnisciently.

Figure 2: Ground truth test slices used for comparison of low dose reconstruction techniques.

Figure 3: 2D Low-dose reconstructions of Walnut dataset for 320 and 620. Gaussian noise of 0 mean and variance equal to 2% of average Poisson-corrupted measurement was added to simulate the low-dose acquisition. The SSIM values are shown in Fig. 5.

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 Poissononly likelihood have very similar performance. However, at a theoretical level, the former is a more principled method, and can deal with negative-valued

Figure 4: 2D Low-dose reconstructions of Colon dataset for 320 and 620. Gaussian noise of 0 mean and variance equal to 2% 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.

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 = 620 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 [39]),

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 [39].

5. Reconstruction with prior

ing the reconstruction performance. However, when the x-ray dose is less, the performance can be further improved by incorporation of useful priors [42, 43]. 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 [44] 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.,

Poisson(Iexp) is not applicable to the spatial reconstructed image domain.

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:

bin in y is not part of the ‘new changes’, then the following is true: + 3/8 + + 3/8 + 4).

6. Based on the aforementioned fact, hypothesis testing is performed on + 3/8 + + 3/8 + to detect bins corresponding to new

as follows:

where the eigenvectors V and mean of the templates form the eigenspace which is built from the available high-dose reconstructions of the templates. denotes the coefficients of V , 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 W . Eq. 23 is solved by alternating minimization on and until convergence is reached.

5.1. Reconstruction results

low dose measurements. Fig. 6 shows a slice from each of the template and test volumes of the potato dataset. This dataset 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 150 150. The

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 [150

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); SSIM = 0.22 (c) Slice from unweighted prior reconstruction; SSIM = 0.42 . The new change is missing. (d) Slice from weighted prior reconstruction; SSIM = 0.69. 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 [39].

corresponding size of the reconstructed volume is 150 150 150. 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 900 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 1.5%. Mean of the added Gaussian noise was 0 and was set to 0.1% of the mean of Poisson-corrupted measurements. Fig 7 shows the same slice from each of the reconstructed volumes. A patch size of [5, 5] 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. 7e. The reconstructed volumes can be found in [39].

5.2. Re-irradiation to improve reconstruction

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:

The new measurement model is: Poisson(exp) + . 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 reimaged 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.

Figure 8: Dataset for illustrating re-irradiation: Templates (first four from the left) and test (extreme right). Size of each slice is (310 The RoI shows the region of difference between the test and the templates.

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 20% of the bins enable better reconstruction of new

Figure 10: Sprouts Dataset for illustrating re-irradiation: Templates (first row) and test (second row). Size of each slice is (156156). The RoI shows the region of difference between the test and the templates.

6. Tuning of parameters

weight for CS term and : weight for object-prior. Below are few of the ways

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

8 times the earlier low-dose X-ray intensity at 25% of the bins enable better reconstruction

to select these parameters.

6.1. Selection of weightage for CS term

ularization 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 fi-delity term. For the Poisson-Gaussian noise model, the cost function is given by ) = .

. Let . Clearly, we have ) = . Hence we can state that )] = m. Furthermore, our

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

simulations (Fig. 12) have shown that

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 R = is very small compared to its mean. The expected value of R 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

and observe how D and relative MSE of reconstructions vary for different values

Figure 13: A method to choose the parameter in low-dose reconstruction: We expect D to

be minimum at the same for which relative MSE is minimum. Here, the and relative MSE are minimum are very close. Refer to Fig. 14 to observe the reconstruction results for different values of

of . At the optimum must be minimum. The test image (154 154) and the reconstructions are shown in Figure 14. For these reconstructions, 410 projection views were chosen and Gaussian noise = 0.3% 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 D 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 D.

6.2. Selection of weightage for object-prior term

for every dataset. However, for a variation of 300, 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.

Figure 14: Colon test data and its reconstructions for different values of is minimum for 2, shown in green, with a relative MSE of 0.1691. The reconstruction for shown in red, gives the minimum relative MSE of 0.1501.

7. Conclusions

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 spe-cific 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

[1] D. J. Brenner, E. J. Hall, Computed tomography an increasing source of radiation exposure, New England Journal of Medicine 357 (22) (2007) 2277–2284.

[2] M. Howells, T. Beetz, H. Chapman, C. Cui, J. Holton, C. Jacobsen, J. Kirz, E. Lima, S. Marchesini, H. Miao, D. Sayre, D. Shapiro, J. Spence, D. Starodub, An assessment of the resolution limitation due to radiation-damage in X-ray diffraction microscopy, Journal of Electron Spectroscopy and Related Phenomena 170 (1) (2009) 4 – 12.

[3] J. C. Dainty, R. Shaw, Image Science: principles, analysis and evaluation of photographic-type imaging processes, Academic Press, 1974.

[4] J. Xu, B. M. W. Tsui, Electronic noise modeling in statistical iterative recon- struction, IEEE Transactions on Image Processing 18 (6) (2009) 1228–1238.

[5] H. Zhang, Y. Dong, Q. Fan, Wavelet frame based poisson noise removal and image deblurring, Signal Processing 137 (2017) 363 – 372.

[6] O. Barkan, J. Weill, S. Dekel, A. Averbuch, A mathematical model for adaptive computed tomography sensing, IEEE Transactions on Computational Imaging 3 (4) (2017) 551–565.

[7] A. Fischer, T. Lasser, M. Schrapp, J. Stephan, P. B. Nol, Object specific trajec- tory optimization for industrial X-ray computed tomography, Scientific Reports 6 (19135).

[8] A. Dabravolski, K. J. Batenburg, J. Sijbers, Dynamic angle selection in x-ray computed tomography, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 324 (2014) 17 – 24, 1st International Conference on Tomography of Materials and Structures.

[9] X. Pan, E. Y. Sidky, M. Vannier, Why do commercial CT scanners still employ traditional, filtered back-projection for image reconstruction?, Inverse problems 25 (12).

[10] Philips, Philips IMR offers new capabilities to simultaneously reduce CT radiation and enhance image quality, https://www. philips.com/a-w/about/news/archive/standard/news/press/2013/ 20130617-Philips-IMR-offers-new-capabilities.html (news-article).

[11] M. Fujita, T. Higaki, Y. Awaya, T. Nakanishi, Y. Nakamura, F. Tatsugami, Y. Baba, M. Iida, K. Awai, Lung cancer screening with ultra-low dose ct using full iterative reconstruction, Japanese Journal of Radiology 35 (4) (2017) 179–189.

[12] K. Kim, G. El Fakhri, Q. Li, Low-dose CT reconstruction using spatially encoded nonlocal penalty, Medical physics 44 (10) (2017) e376–e390.

[13] Q. Lyu, D. Ruan, J. Hoffman, R. Neph, M. McNitt-Gray, K. Sheng, Iterative reconstruction for low dose CT using plug-and-play alternating direction method of multipliers (ADMM) framework, SPIE Medical Imaging 10949.

[14] X. Zheng, Z. Lu, S. Ravishankar, Y. Long, J. A. Fessler, Low dose CT image reconstruction with learned sparsifying transform, in: 2016 IEEE 12th Image, Video, and Multidimensional Signal Processing Workshop, 2016, pp. 1–5.

[15] X. Jia, Z. Bian, J. He, Y. Wang, J. Huang, D. Zeng, Z. Liang, J. Ma, Texture- preserved low-dose CT reconstruction using region recognizable patch-priors from previous normal-dose CT images, in: IEEE Nuclear Science Symposium, 2016, pp. 1–4.

[16] C. Zhang, T. Zhang, M. Li, C. Peng, Z. Liu, J. Zheng, Low-dose CT recon- struction via L1 dictionary learning regularization using iteratively reweighted least-squares, BioMedical Engineering OnLine 15 (1) (2016) 66.

[17] A. Andersen, A. Kak, Simultaneous algebraic reconstruction technique (SART): A superior implementation of the ART algorithm, Ultrasonic Imaging 6 (1) (1984) 81 – 94.

[18] Q. Xu, H. Yu, X. Mou, L. Zhang, J. Hsieh, G. Wang, Low-dose x-ray CT recon- struction via dictionary learning, IEEE Transactions on Medical Imaging 31 (9) (2012) 1682–1697.

[19] J. C. Park, B. Song, J. S. Kim, S. H. Park, H. K. Kim, Z. Liu, T. S. Suh, W. Y. Song, Fast compressed sensing-based CBCT reconstruction using barzilai-borwein formulation for application to on-line IGRT, Med. Phy. 39 (3) (2012) 1207–1217.

[20] J. Wu, F. Liu, L. Jiao, X. Wang, Multivariate pursuit image reconstruction using prior information beyond sparsity, Signal Processing 93 (6) (2013) 1662 – 1672, special issue on Machine Learning in Intelligent Image Processing.

[21] J. Wang, J. Ma, B. Han, Q. Li, Split Bregman iterative algorithm for sparse re- construction of electrical impedance tomography, Signal Processing 92 (12) (2012) 2952 – 2961.

[22] D. Wu, K. Kim, G. El Fakhri, Q. Li, Iterative low-dose CT reconstruction with priors trained by artificial neural network, IEEE Transactions on Medical Imaging 36 (12) (2017) 2479–2486.

[23] H. Shan, A. Padole, F. Homayounieh, U. Kruger, R. D. Khera, C. Nitiwarangkul, M. K. Kalra, G. Wang, Competitive performance of a modularized deep neural network compared to commercial algorithms for low-dose CT image reconstruction, Nature Machine Intelligence 1 (6) (2019) 269–276.

[24] C. Gong, L. Zeng, Adaptive iterative reconstruction based on relative total vari- ation for low-intensity computed tomography, Signal Processing 165 (2019) 149 – 162.

[25] A. Beck, M. Teboulle, A fast iterative shrinkage-thresholding algorithm for linear inverse problems, Imaging Sciences 2 (1) (2009) 183–202.

[26] W. Hou, C. Zhang, A compressed sensing approach to low-radiation CT recon- struction, 2014 9th International Symposium on Communication Systems, Networks Digital Sign (CSNDSP) (2014) 793–797.

[27] K. Koh, S.-J. Kim, S. Boyd, l1-ls: Simple matlab solver for l1-regularized least squares problems, https://stanford.edu/~boyd/l1_ls/ (solver).

[28] L. Fu, et al., Comparison between pre-log and post-log statistical models in ultra- low-dose CT reconstruction, IEEE transactions on medical imaging 36 (3) (2016) 707 – 720.

[29] L. Feldkamp, L. C. Davis, J. Kress, Practical cone-beam algorithm, J. Opt. Soc. Am 1 (1984) 612–619.

[30] F. Pontana, A. Duhamel, J. Pagniez, T. Flohr, J.-B. Faivre, A.-L. Hachulla, J. Remy, M. Remy-Jardin, 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) (2011) 636–643.

[31] H. Wang, B. Tan, B. Zhao, C. Liang, Z. Xu, Raw-data-based iterative reconstruc- tion versus filtered back projection: image quality of low-dose chest computed tomography examinations in 87 patients, Clin. Imaging 37 (6) (2013) 1024 – 32.

[32] H. Koyama, Y. Ohno, M. Nishio, S. Matsumoto, N. Sugihara, T. Yoshikawa, S. Seki, K. Sugimura, 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) (2014) 1860–1867.

[33] M. J. Willemink, P. B. No¨el, The evolution of image reconstruction for CT— from filtered back projection to artificial intelligence, European Radiology 29 (5) (2019) 2185–2195.

[34] Q. Xie, D. Zeng, Q. Zhao, D. Meng, Z. Xu, Z. Liang, J. Ma, Robust low-dose CT sinogram preprocessing via exploiting noise-generating mechanism, IEEE Transactions on Medical Imaging 36 (12) (2017) 2487–2498.

[35] J. A. Fessler, Penalized weighted least-squares image reconstruction for positron emission tomography, IEEE Transactions on Med. Imaging 13 (2) (1994) 290–300.

[36] Q. Ding, Y. Long, X. Zhang, J. A. Fessler, Statistical image reconstruction using mixed Poisson-Gaussian noise model for X-ray CT, submitted in Inverse Prob. and Imaging, axXivSubmitted.

[37] T. Hastie, R. Tibshirani, M. Wainwright, Statistical Learning with Sparsity, CRC Press, Taylor & Francis Group, 2015.

[38] E. Chouzenoux, A. Jezierska, J. Pesquet, H. Talbot, A convex approach for image restoration with exact poisson–gaussian likelihood, SIAM Journal on Imaging Sciences 8 (4) (2015) 2662–2682.

[39] Supplementary material, videos, more results, https://www.dropbox.com/s/ x4jwzjpcw8x7dz7/supplementary_material.pdf?dl=0 (video, proofs and comparison of methods).

[40] Walnut, Walnut CT dataset, https://www.uni-muenster.de/Voreen/download/ workspaces_and_data_sets.html (uni-muenster).

[41] Colon, CT Colonography, https://idash.ucsd.edu/data-collections (idash.ucsd.edu).

[42] G.-H. Chen, J. Tang, S. Leng, Prior image constrained compressed sensing (PICCS): A method to accurately reconstruct dynamic CT images from highly undersampled projection data sets, Medical Physics 35 (2) (2008) 660–663.

[43] E. A. Rashed, H. Kudo, Probabilistic atlas prior for CT image reconstruction, Computer Methods and Programs in Biomedicine 128 (2016) 119 – 136.

[44] P. Gopal, S. Chandran, I. D. Svalbe, A. Rajwade, Learning from past scans: Tomographic reconstruction to detect new structures, CoRR abs/1812.10998.

[45] F. J. Anscombe, The transformation of poisson, binomial and negative-binomial data, Biometrika 35 (3/4) (1948) 246–254.

[46] J. H. Curtiss, On transformations used in the analysis of variance, Ann. Math. Statist. 14 (2) (1943) 107–122.

[47] F. Murtagh, J. luc Starck, A. Bijaoui, Image restoration with noise suppression using a multiresolution support, Astronomy and Astrophysics, Suppl. Ser 112 (1995) 179–189.

[48] J. Liu, Y. Hu, J. Yang, Y. Chen, H. Shu, L. Luo, Q. Feng, Z. Gui, G. Coatrieux, 3D feature constrained reconstruction for low dose CT imaging, IEEE Transactions on Circuits and Systems for Video Technology PP (99) (2016) 1–1.