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 Xray beam. The CT imaging model that incorporates the strength of Xrays, , is nonlinear and nondeterministic and is given by:
(1) 
where
represents the zero mean additive Gaussian noise vector with a fixed signalindependent 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 PoissonGaussian noise model is quite common in optical or Xray based imaging systems, but we consider it here explicitly for tomography, where it induces a nonlinear 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 lowdose acquisition is the large magnitude (relative to the signal) of Poisson noise due to the low strength of Xray beam. This is because the SignaltoNoiseRatio (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 Poissonnoise 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 fewviews 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 ^{1}^{1}1Filtered Backprojection as the standard reconstruction technique Pan et al. (2009). Only recently are the iterative techniques being deployed for commercial use Philips (newsarticle). The power of iterative routines was reinforced by Fujita et al. (2017), where it was proved that iterative reconstructions from ultralow dose ^{2}^{2}2Typically, lowdose imaging is performed at 120 kVp and 30 mAs beam current, and ultralow dose imaging is performed at 80100 kVp and 2030 mAs beam current settings. CT are of similar quality to those of FBP reconstructions from lowdose CT. Here, a commercial forward projected modelbased algorithm was deployed and compared with FBP.Among the other iterative methods, Kim et al. (2017)
presented a technique that minimizes loglikelihood of the Poisson distribution and a patchbased spatially encoded nonlocal penalty.
Lyu et al. (2019) used a smoothness prior along with datafidelity constraint and solved using ADMM. In order to improve the reconstruction further, various priorbased and learningbased methods have also been explored in literature. In these techniques, properties of available standarddose CT images influence lowdose 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 prelearned sparsifying transform. While the weights were set manually, the sparsifying transform was learned from a database of regulardose CT images. Another technique presented by Jia et al. (2016)clustered overlapping patches of previously scanned standarddose 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 (L1DL), in the same order. Wu et al. (2013) used edgebased priors to reconstruct normaldose 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 loglikelihood based costfunction that accurately reflects the PoissonGaussian 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 lowdose reconstruction.
Wu et al. (2017) proposed one such neural network to learn features of the image that is later imposed along with datafidelity 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 neuralnetwork 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 Lcurve method in which datafidelity 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 lowdose imaging. In this work, we use the noisemodel for the purpose of automated parameter selection.
3 Contributions
This paper discusses the following:

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

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 reirradiation which improves the reconstruction quality in these regions of change at the cost of a small amount of added radiation.

In addition, most of the iterative schemes involve a cost function with a datafidelity term, a sparsity term and a dataprior 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 priorbased reconstruction in longitudinal studies.
Sec. 4 describes two new techniques and its comparison with a few methods in literature. Sec. 5 describes the new priorbased lowdose reconstruction and its validation on 3D tomographic data. In Sec. 6, we illustrate a systematic technique to parametertuning. Finally, key results are summarized in Sec. 7
4 Reconstruction without prior
A good lowdose 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 datafidelity cost, and may possibly (though not necessarily) be expressed by the negative loglikelihood 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 wellknown 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 Postlog 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 PoissonGaussian statistics and suffers from an inherent statistical bias (as seen in Fig. 1
) as it is a socalled ‘postlog’ method. The bias arises because for any nonnegative random variable
, we have as per Jensen’s inequality. Another way of viewing this is that the noise in (i.e. postlog) 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 postlog measurements is also discussed in (Fu and others, 2016).4.2 Nonlinear Least Squares with CS
An intuitive way to modify the previous cost is by allowing the data fidelity cost to mimic the nonlinearity 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 nonnegativity constraint is imposed on . It can be seen that this cost function is nonconvex 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 postlog FBP. While it is computationally efficient, it suffers from a statistical bias for the same reasons as postlog CS, as described in 4.1. The performance of postlog 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 lowdose reconstructions (Willemink and Noël, 2019).
4.4 Negative Log LikelihoodPoisson with CS
This technique accounts for only the Poisson noise (ignoring the Gaussian part) and searches for a solution that minimizes the negative loglikelihood 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 LikelihoodPoissonGaussian 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 loglikelihood 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 loglikelihood of is given by
(10) 
We minimize the sum of the two negative loglikelihoods:
(11) 
and are solved for alternately. Note that is integervalued, but a typical gradientbased 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 nonconvex. However it can be shown to be biconvex, i.e., it is convex in if is kept fixed and vice versa. Such a costfunction was used in (Xie et al., 2017) as a method of preprocessing/denoising of projections prior to tomographic reconstruction. In contrast, we directly use it as a datafidelity 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 nonlinear 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 datafidelity 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 PoissonGaussian 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 quasilikelihood problem. On the other hand, the proposed method (RNLLS) can be interpreted as a weighted form of the wellknown 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 PoissonGaussian noise, our numerical simulations reveal that the cost function is not convex in the worst case. However, this nonconvexity did not affect the numerical results significantly.
4.7 Proposed PoissonGaussian 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 PoissonGaussian corruption of the signal. In contrast, in CT, the measured signal is a nonlinear 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 loglikelihood 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 nonconvex. 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 lowdose 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 lowdose 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 noisetosignal 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 Poissoncorrupted measurement was added to measurements. The regularization parameter was chosen omnisciently.
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 PoissonGaussian likelihood reconstructions were comparable and gave the best reconstructions for a majority of dose levels and datasets.The PoissonGaussian 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 negativevalued measurements which have to be weeded out for the Poissononly method. The nonlinear least squares method (Sec. 4.2) performed poorly. This is because the datafidelity term assumes constant variance for all signal values. In reality, the variance of Poisson noise increases as signal intensity increases. The postlog linear least squares (Sec. 4.1) failed because the linear model fails to approximate the highly nonlinear lowdose acquisition. The postlog 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 postlog 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 nonlinear least squares (RNLLS) is inbetween the performance of likelihoodbased 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 xray 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 highquality reconstructions of the same object in longitudinal studies, or highquality reconstructions of similar objects. We refer to such prior data as templates. Here, our aim is to reconstruct an object from its lowdose measurements, using templates which are previous highdose 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 ‘weightsmap’) 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 weightsmap 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 lowdose 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 weightsmap (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:

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

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

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

Let be the noisy tomographic projection of the test volume from the angle. For each , project onto , i.e., compute the eigencoefficients 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 eigenspace computed in image domain, used later in Eq. 23). Next, compute the resultant projection , i.e.,
(22) 
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 PoissonGaussian 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.

Based on the aforementioned fact, hypothesis testing is performed on to detect bins corresponding to new changes in the measurement space. We use Ztest 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 probabilitythat the sample is drawn from Normal distribution must lie in the
tailend 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. 
Once the new changes are detected in the measurement space, filtered backprojection of the vectors (containing pvalues) 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.

The final weightsmap ^{3}^{3}3An alternate method to compute a weightsmap (a simpler binary weightsmap) 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 priorsimilar regions, (b) Linear stretching: Perform linear stretching on so that the weights lie between 0 and 1.
Finally, the computed weightsmap 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 highdose 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 weightsmap . Eq. 23 is solved by alternating minimization on and until convergence is reached.
5.1 Reconstruction results
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 ^{4}^{4}4We 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 conebeam 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 conebeam measurements were simulated from fullview 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 Poissoncorrupted 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 weightsmap as seen in Fig. 6(e). The reconstructed volumes can be found in (Supplementary material, video, proofs and comparison of methods).
5.2 Reirradiation to improve reconstruction
Once the regions of new changes are detected by the weights map, this information can be used to reirradiate them with standarddose rays and further improve the quality of their reconstruction. Following are the steps of the reirradiation process:

Let the Xrays passing through the new regions have their source points denoted by , and the corresponding bins at the detector be denoted by . Let the Xrays 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 .

Block and reirradiate the object by passing standarddose 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.

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 standarddose measurements corresponding to new regions of the object and lowdose 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 Xray 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 reirradiation. The new changes within the RoI are reconstructed very well after they are reimaged with standarddose Xrays. This is also reinforced by results on the sprouts data (Fig. 10), shown in Fig. 11. The selection of bins for reirradiation and the choice of new Xray intensity can also be chosen in a supervised manner by the physician or scientist based on the particular clinical or nonclinical setting.
6 Tuning of parameters
Two parameters were used in the techniques presented in this chapter: : weight for CS term and : weight for objectprior. 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 crossvalidation. 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 PoissonGaussian noise model. The method is shown here in conjunction with the rescaled nonlinear least squares method, however in principle, it can be used with any data fidelity term. For the PoissonGaussian noise model, the cost function is given by .
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 elementwise 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) 
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 Xrays 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 reallife setting, when relative MSE cannot be computed because of absence of groundtruth, a brute force search needs to be done followed by selecting the value of that minimizes .
6.2 Selection of weightage for objectprior 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 lowdose CT imaging regime, the noise in the measurements becomes significant and needs to be accounted for during the reconstruction. Two new techniques: PoissonGaussian convolution and rescaled nonlinear least squares (RNLLS) were presented and extensively compared with many of the existing methods. RNLLS was further used in lowdose 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 reirradiating these specific regions by standarddose Xrays. 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
 Simultaneous algebraic reconstruction technique (SART): a superior implementation of the ART algorithm. Ultrasonic Imaging 6 (1), pp. 81 – 94. External Links: ISSN 01617346 Cited by: §2.
 The transformation of poisson, binomial and negativebinomial data. Biometrika 35 (3/4), pp. 246–254. External Links: ISSN 00063444 Cited by: item 5.
 A mathematical model for adaptive computed tomography sensing. IEEE Transactions on Computational Imaging 3 (4), pp. 551–565. Cited by: §2.
 A fast iterative shrinkagethresholding algorithm for linear inverse problems. Imaging Sciences 2 (1), pp. 183–202. Cited by: §4.
 Computed tomography — an increasing source of radiation exposure. New England Journal of Medicine 357 (22), pp. 2277–2284. Cited by: §1.
 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.
 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.
 CT Colonography. Note: https://idash.ucsd.edu/datacollections Cited by: Figure 2.
 On transformations used in the analysis of variance. Ann. Math. Statist. 14 (2), pp. 107–122. Cited by: item 5.
 Dynamic angle selection in xray 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.
 Image science: principles, analysis and evaluation of photographictype imaging processes. Academic Press. Cited by: §1.
 Statistical image reconstruction using mixed PoissonGaussian noise model for Xray CT. submitted in Inverse Prob. and Imaging, axXiv. Note: Submitted Cited by: §4.6.
 Practical conebeam algorithm. J. Opt. Soc. Am 1, pp. 612–619. Cited by: §4.3.
 Penalized weighted leastsquares image reconstruction for positron emission tomography. IEEE Transactions on Med. Imaging 13 (2), pp. 290–300. Cited by: §4.6.
 Object specific trajectory optimization for industrial Xray computed tomography. Scientific Reports 6 (19135). Cited by: §2.
 Comparison between prelog and postlog statistical models in ultralowdose CT reconstruction. IEEE transactions on medical imaging 36 (3), pp. 707 – 720. Cited by: §4.1.
 Lung cancer screening with ultralow dose ct using full iterative reconstruction. Japanese Journal of Radiology 35 (4), pp. 179–189. External Links: ISSN 1867108X Cited by: §2.
 Adaptive iterative reconstruction based on relative total variation for lowintensity computed tomography. Signal Processing 165, pp. 149 – 162. Cited by: §2.
 Learning from past scans: tomographic reconstruction to detect new structures. CoRR abs/1812.10998. Cited by: §5.
 Statistical learning with sparsity. CRC Press, Taylor & Francis Group. Cited by: §4.6.
 A compressed sensing approach to lowradiation CT reconstruction. 2014 9th International Symposium on Communication Systems, Networks Digital Sign (CSNDSP), pp. 793–797. Cited by: §4.1.
 An assessment of the resolution limitation due to radiationdamage in Xray diffraction microscopy. Journal of Electron Spectroscopy and Related Phenomena 170 (1), pp. 4 – 12. Cited by: §1.
 Texturepreserved lowdose CT reconstruction using region recognizable patchpriors from previous normaldose CT images. In IEEE Nuclear Science Symposium, Vol. , pp. 1–4. Cited by: §2.
 Lowdose CT reconstruction using spatially encoded nonlocal penalty. Medical physics 44 (10), pp. e376–e390. Cited by: §2.
 L1ls: simple matlab solver for l1regularized least squares problems. Note: https://stanford.edu/~boyd/l1_ls/ Cited by: §4.1.
 Iterative reconstruction technique vs filter back projection: utility for quantitative bronchial assessment on lowdose thinsection MDCT in patients with/without chronic obstructive pulmonary disease. European Radiology 24 (8), pp. 1860–1867. External Links: ISSN 14321084 Cited by: §4.3.
 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.
 Iterative reconstruction for low dose CT using plugandplay alternating direction method of multipliers (ADMM) framework. SPIE Medical Imaging 10949. Cited by: §2.
 Image restoration with noise suppression using a multiresolution support. Astronomy and Astrophysics, Suppl. Ser 112, pp. 179–189. Cited by: item 5.
 Why do commercial CT scanners still employ traditional, filtered backprojection for image reconstruction?. Inverse problems 25 (12), pp. . Cited by: §2.
 Fast compressed sensingbased CBCT reconstruction using barzilaiborwein formulation for application to online IGRT. Med. Phy. 39 (3), pp. 1207–1217. Cited by: §2.
 Philips IMR offers new capabilities to simultaneously reduce CT radiation and enhance image quality. Note: https://www.philips.com/aw/about/news/archive/standard/news/press/2013/20130617PhilipsIMRoffersnewcapabilities.html Cited by: §2.
 Chest computed tomography using iterative reconstruction vs filtered back projection (part 2): image quality of lowdose ct examinations in 80 patients. European Radiology 21 (3), pp. 636–643. External Links: ISSN 14321084 Cited by: §4.3.
 Probabilistic atlas prior for CT image reconstruction. Computer Methods and Programs in Biomedicine 128, pp. 119 – 136. Cited by: §5.
 Competitive performance of a modularized deep neural network compared to commercial algorithms for lowdose CT image reconstruction. Nature Machine Intelligence 1 (6), pp. 269–276. Cited by: §2.
 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 CT dataset. Note: https://www.unimuenster.de/Voreen/download/workspaces_and_data_sets.html Cited by: Figure 2.
 Rawdatabased iterative reconstruction versus filtered back projection: image quality of lowdose chest computed tomography examinations in 87 patients. Clin. Imaging 37 (6), pp. 1024 – 32. External Links: ISSN 08997071 Cited by: §4.3.
 Split Bregman iterative algorithm for sparse reconstruction of electrical impedance tomography. Signal Processing 92 (12), pp. 2952 – 2961. Cited by: §2.

The evolution of image reconstruction for CT—from filtered back projection to artificial intelligence
. European Radiology 29 (5), pp. 2185–2195. External Links: ISSN 14321084 Cited by: §4.3.  Iterative lowdose CT reconstruction with priors trained by artificial neural network. IEEE Transactions on Medical Imaging 36 (12), pp. 2479–2486. External Links: ISSN 02780062 Cited by: §2.

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.  Robust lowdose CT sinogram preprocessing via exploiting noisegenerating mechanism. IEEE Transactions on Medical Imaging 36 (12), pp. 2487–2498. Cited by: §4.5.
 Electronic noise modeling in statistical iterative reconstruction. IEEE Transactions on Image Processing 18 (6), pp. 1228–1238. Cited by: §1.
 Lowdose Xray CT reconstruction via dictionary learning. IEEE Transactions on Medical Imaging 31 (9), pp. 1682–1697. Cited by: §2.
 Lowdose CT reconstruction via L1 dictionary learning regularization using iteratively reweighted leastsquares. BioMedical Engineering OnLine 15 (1), pp. 66. Cited by: §2.
 Wavelet frame based poisson noise removal and image deblurring. Signal Processing 137, pp. 363 – 372. Cited by: §2.
 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.