Learning Personalized Representation for Inverse Problems in Medical Imaging Using Deep Neural Network
Abstract
Recently deep neural networks have been widely and successfully applied in computer vision tasks and attracted growing interests in medical imaging. One barrier for the application of deep neural networks to medical imaging is the need of large amounts of prior training pairs, which is not always feasible in clinical practice. In this work we propose a personalized representation learning framework where no prior training pairs are needed, but only the patient’s own prior images. The representation is expressed using a deep neural network with the patient’s prior images as network input. We then applied this novel image representation to inverse problems in medical imaging in which the original inverse problem was formulated as a constraint optimization problem and solved using the alternating direction method of multipliers (ADMM) algorithm. Anatomically guided brain positron emission tomography (PET) image reconstruction and image denoising were employed as examples to demonstrate the effectiveness of the proposed framework. Quantification results based on simulation and real datasets show that the proposed personalized representation framework outperform other widely adopted methods.
1 Introduction
Over the past several years, deep neural networks have been widely and successfully applied to computer vision tasks such as image segmentation, object detection, and image super resolution, by demonstrating better performance than stateofart methods when large amounts of datasets are available. For medical imaging tasks such as lesion detection and regionofinterest (ROI) quantification, obtaining highquality diagnostic images is essential. Recently the neural network method has been applied to inverse problems in medical imagingWang et al. (2016); Sun et al. (2016); Kang et al. (2017); Chen et al. (2017); Hammernik et al. (2018); Wu et al. (2017); Adler and Öktem (2018); Schlemper et al. (2018); Gong et al. (2017a); Kim et al. (2018).
During network training, images reconstructed from high dose or longscanned duration are needed as training labels. However, collecting large amounts of training labels is not an easy task: highdose computed tomography (CT) has potential safety concerns; longscanned dynamic PET is not employed in routine clinical practice; in cardiac magnetic resonance imaging (MRI), it is impossible to acquire breathhold and fully sampled 3D images. With limited amounts of highquality patient datasets available, overfitting can be a potential pitfall: if a new patient dataset does not lie in the training space due to population difference, the trained network cannot accurately recover unseen structures. In addition, lowquality images are often simulated by artificially downsampling the fulldose/highcount data, which may not reflect the real physical conditions of lowdose imaging. This mismatch between training and the real clinical environment can reduce the network performance.
Apart from using training pairs to perform supervised learning, a lot of prior arts focus on exploiting prior images acquired from the same patient to improve the image quality. The priors can come from temporal information (Chen et al., 2008; Wang and Qi, 2015), different physics settings (Gnahm and Nagel, 2015), or even other imaging modalities (Gong et al., 2017b). They are included into the maximum posterior estimation or sparse representation framework using predefined analytical expressions or prelearning steps. The predefined expressions might not be able to extract all the useful information, and the prelearnt model might not be optimal for the later inverseproblem optimization as no dataconsistency constraint is enforced during prelearning. Ideally the learning process should be included inside the inverseproblem solving process.
Recently, the deep image prior framework proposed in (Ulyanov et al., 2017) shows that convolutional neural networks (CNNs) can learn intrinsic structure information from the corrupted images based on random noise input. No prior training pairs or learning process are needed in this framework. The ability of neural networks to learn the structure information is also revealed in generative adversarial networks (GANs) (Goodfellow et al., 2014; Radford et al., 2015), where various distributions can be generated based on random noise input. These prior arts show that the network itself can be treated as a regularizer, and it is possible to train a network without highquality reference images. Furthermore, it has been shown in conditional GAN works (Mirza and Osindero, 2014; Isola et al., 2017; Zhu et al., 2017) that when the input is not random noise, but the associated prior information, the prediction results can be improved.
Inspired by the prior arts, we propose a new framework to solve inverse problems in medical imaging based on learning personalized representation. The personalized representation is learnt through a deep neural network with prior images of the same patient as network input. ADMM algorithm was employed to decouple the inverseproblem optimization and network training steps. Modified 3D Unet (Çiçek et al., 2016) was employed as the network structure. The total number of trainable parameters for the network is about 1.46 million, which is less than the unknowns in the original inverse problem.
PET is a molecular imaging modality widely used in neurology studies Gong et al. (2016). The image resolution of current PET scanners is still limited by various physical degradation factors Gong et al. (2017c). Improving PET image resolution is essential for a lot of applications, such as dopamine neurotransmitter imaging, brain tumor staging and early diagnosis of Alzheimer’s disease. For the past decades, various efforts are focusing on using MR or CT to improve PET image quality Bai et al. (2013). In this work, anatomically guided brain PET image reconstruction and denoising are employed as two examples to demonstrate the effectiveness of the proposed framework.
The main contributions of this work include:(1) proposing to employ the patient’s prior images as network input to construct a personalized representation;(2) embedding this personalized representation in complicated inverse problems to accurately estimate unknown images; (3) proposing a ADMM framework to decouple the whole optimization into a penalized inverse problem and a network training problem; (4) demonstrating the effectiveness of the proposed framework using clinical brain imaging applications.
2 Methodology
2.1 Related works
In inverse problems such as image reconstruction and denoising, the measured data can be modelled as a collection of independent random variables and its mean is related to the original image through an affine transform
(1) 
where is the transformation matrix and is the known additive item. For the past decades, a lot of efforts focus on representing the original image by different basis functions (Do and Vetterli, 2005; Elad and Aharon, 2006; Wang and Qi, 2015). For example, in dictionary learning proposed by Elad and Aharon (2006) , , the patch of , can be represented as
(2) 
where is the prelearnt overcomplete dictionary based on highquality reference images, and is the coefficient for the patch. constraint is added on to enforce sparse representation. In the kernel method proposed by Wang and Qi (2015), the original image is represented as
(3) 
where is the kernel matrix calculated using the radial basis function, and is the kernel coefficient. Both the dictionary learning and the kernel method assumes the local/global of the original unknown image can be linearly represented using prelearnt/predefined basis vectors. In a recent work by Gong et al. (2017a) for image reconstruction tasks, the original unknown image is represented using CNN as
(4) 
where represents the neural network, are the parameters of the neural network, and denotes the network input. In this CNN representation, the network is prelearnt using large number of training pairs from different patient datasets, network parameters are thus fixed, and the network input is keeping updated during the reconstruction process. Recently Ulyanov et al. (2017) uses the same representation as shown in (4) for image superresolution and denoising tasks. Different from Gong et al.’s approach, the network input is a predefined random noise, and the network parameters are keeping updated in the image restoration process. One biggest advantage of Ulyanov et al.’s framework is that pretraining steps and large number of training pairs are not needed.
2.2 Personalized representation using deep neural network
Inspired by the prior arts, here we propose a personalized representation framework for inverse problems in medical imaging. In this proposed framework, the original unknown image is represented by a deep neural network as shown in (4). The prior image of the same patient is employed as network input to construct a personalized representation. The network parameters are updated in the inverseproblem optimization process. Compared to employing random noise as network input as presented in (Ulyanov et al., 2017), using prior images of the patient can make the representation more accurate. To illustrate this, Fig. 2 presents the comparison of network outputs with different network inputs using the brain dataset introduced in Section 3.2. Clearly, when the patient’ own MR image is employed as network input, the cortices regions are clearer and noise in white matter regions is reduced. After substituting with the neural network representation (4), the original data model shown in (1) can be rewritten as
(5) 
Supposing the measured random variable follows a distribution of , the log likelihood for the measured data can be written as
(6) 
The maximum likelihood estimate of the unknown image can be calculated in two steps as
(7) 
2.2.1 Optimization
The objective function in (7) is difficult to solve due to the coupling between the likelihood function and the neural network. Here we transfer it to the constrained format as below
max  (8)  
s.t. 
We use the augmented Lagrangian format for the constrained optimization problem in (8) as
(9) 
which can be solved by the ADMM algorithm iteratively in three steps
(10)  
(11)  
(12) 
After using the ADMM algorithm, we decoupled the constraint optimization problem into a network training problem with L2 norm as loss function (10), and a penalized inverse problem (11). For subproblem (10), currently network training is mostly based on first order methods, such as the Adam algorithm (Kingma and Ba, 2014) and the Nesterov’s accelerated gradient (NAG) algorithm (Nesterov, 1983). The LBFGS algorithm (Liu and Nocedal, 1989) is a quasinewton method, combining a history of updates to approximate the Hessian matrix. It is not widely used in network training as it requires large batch size to accurately calculate the descent direction, which is less effective than first order methods for largescale applications. In this proposed framework, as only the patient’s own prior images are employed as network input, the data size is much smaller than the traditional network training. In our case, the LBFGS method is preferred to solve subproblem (10) due to its stability and better performance observed in the experiments (results shown in Section. 2.2.2). For the penalized inverse problem shown in subproblem (11), it is application dependent and is discussed as below.
PET image reconstruction
In PET image reconstruction, is the detection probability matrix, with denoting the probability of photons originating from voxel being detected by detector . denotes the expectation of scattered and random events. is the number of lines of response (LOR). The measured photon coincidences follow Poisson distribution, and the loglikelihood function can be explicitly written as
(13) 
Subproblem (11) thus corresponds to a penalized reconstruction problem. The optimization transfer method Lange et al. (2000) is used to solve it. As in is coupled together, we first construct a surrogate function for to decouple the image pixels so that each pixel can be optimized independently. is constructed as follows
(14) 
where and is calculated by
(15) 
It can be verified that the constructed surrogate function fulfills the following two conditions:
(16)  
(17) 
After getting this surrogate function, subproblem (10) can be optimized pixel by pixel. For pixel , the surrogate objective function for subproblem (10) is
(18) 
The final update equation for pixel after maximizing (18) is
(19) 
PET image restoration
For PET image denoising problems, is an identity matrix and the noisy images follow the Gaussian distribution. Thus objective function (7) can be seen as a network training problem, which can be directly solved using the LBFGS algorithm with L2 norm as the loss function. For PET image deblurring problems, is a blurring matrix. Solving the optimization problem follows the same flowchart as shown in the previous mentioned PET image reconstruction problem.


2.2.2 Network structure and implementation details
The network structure employed in this work is based on the modified 3D Unet (Çiçek et al., 2016) as shown in Fig. 4. In our implementation, there are several modifications compared to the original 3D Unet: (1) using convolutional layer with stride 2 to downsample the image instead of using max pooling, to construct a fully convolutional network; (2) directly adding the left side feature to the right side instead of concatenating, to reduce the number of training parameters, thus the risk of overfitting; (3) using the bilinear interpolation instead of the deconvolution upsampling, to reduce the checkboard artifact; (4) using leaky ReLU instead of ReLU. We compared the behaviors of different network training algorithms for subproblem (10). The Adam, NAG, and LBFGS algorithms were compared after running 300 epochs using the real brain data mentioned in Section. 3.2. When comparing different algorithms, we computed the normalized cost value, which is defined as ,where and is the cost value after running Adam for 700 iterations and 1 iteration, respectively. Fig. 7(a) plots the normalized cost value curves for different algorithms. The LBFGS algorithm is monotonic decreasing while the Adam algorithm is not due to the adaptive learning rate implemented. The NAG algorithm is slower than the other two algorithms. Due to the monotonic property, the network output using the LBFGS algorithm is more stable and less influenced by the image noise when running multiple realizations. This is the reason we chose LBFGS algorithm to solve subproblem (10). All the neural network related training was implemented using TensorFlow 1.4. As the penalty parameter used in the ADMM framework has a large impact on the convergence speed, we examined the loglikelihood to determine the penalty parameter used in practice. As an example, Fig. 7 (b) shows the loglikelihood curves using different penalty parameters for the real brain data mentioned in Section. 3.2. Considering the convergence speed and stability of the likelihood, = 3e3 was chosen.
3 Experiment
3.1 Brain PET image reconstruction: simulation study
A 3D brain phantom from BrainWeb AubertBroche et al. (2006) was used in the simulation. Corresponding T1 weighted MR image was used as the prior image. The voxel size is 222 and the phantom image size is 128128105. To simulate mismatches between the MR and PET images, twelve hot spheres of diameter 16 mm were inserted into the PET image as tumor regions, which are not visible in the MR image. In this experiment, the last 5 min frame of a onehour FDG scan was used as the groundtruth image. The computer simulation modeled the geometry of a Siemens mCT scanner. Noisefree sinogram data were generated by forwardprojecting the groundtruth image using the system matrix and the attenuation map. Poisson noise was then introduced to the noisefree data by setting the total count level to be equivalent to last 5 min scan with 5 mCi injection. Gaussian postfiltering method (denoted as EM+filter) and the kernel method (denoted as KMRI) (Wang and Qi, 2015) were employed as comparison methods. Fig. 11 shows three orthogonal views of the reconstructed images using different methods. The kernel method and the proposed method both reveal more cortex structures compared to the EMplusfilter method due to the boundary information provided by the MR priors. Compared to the kernel method, the proposed method can recover even more details of the cortices and the white matter regions are cleaner. Besides, the tumor uptake using the proposed method is higher and the tumor shape is closer to the ground truth. This means that even when there are mismatches between the PET and MR images, the proposed method can still recover the true PET intensities and shapes. Fig. 22(a,b) shows the contrast recovery coefficient(CRC). vs standard deviation (STD) curves for different methods. For both the gray matter region and the tumor region, the proposed method outperforms other methods.



3.2 Brain PET image reconstruction: clinical data study









A 70minutes dynamic PET scan of a human subject acquired on a Siemens Brain MRPET scanner after 5 mCi FDG injection was employed in the real data evaluation. The data were reconstructed with an image array of 256256153 and a voxel size of 1.25 1.25 1.25 . A simultaneous acquired T1weighted MR image has the same image array and voxel size as the PET images. Correction factors for attenuation, randoms, scatters were estimated using the standard software provided by the manufacturer and included during reconstruction. The motion correction was performed in the LOR space based on the simultaneously acquired MR navigator signal(Catana et al., 2011). To generate multiple realizations for quantitative analysis, the last 40 minutes PET data were binned together and resampled with a 1/8 ratio to obtain 20 i.i.d. datasets that mimic 5minutes frames. As the ground truth of the regional uptake is unknown, a hot sphere with diameter 12.5 mm, mimicking a tumor, was added to the PET data (invisible in the MRI image). For tumor quantification, images with and without the inserted tumor were reconstructed and the difference was taken to obtain the tumor only image and compared with the ground truth. The tumor contrast recovery (CR) was calculated as , where is the mean tumor uptake inside the tumor ROI, is the ground truth of the tumor uptake, and is the number of the realizations. For the background, 11 circular ROIs with a diameter of 12.5 mm were drawn in the white matter to calculate the standard deviation. Fig. 18 shows the reconstructed images and the corresponding MR prior images of the real brain dataset using different methods. For the methods with MR information included, more cortex details are recovered and the image noise in the white matter is much reduced. The cortex shape in the proposed method is clearer than the kernel method. For the tumor region which is unobserved in the MR image, the uptake is higher in the proposed method compared with the kernel method. Fig. 22(c) shows the CRSTD curves for different methods. Clearly the proposed method has the best CRSTD tradeoff compared with other methods.
3.3 Clinical PET image denoising
The patient images were acquired with a Siemens Biograph mCT PET/CT system. Patients were injected with 111 MBq (3 mCi) of 68GaPRGD2. PET image at 60 min post injection was acquired with 5 bed positions. For the proposed method, the corresponding CT image is employed as the network input. The Gaussian denoising, and the nonlocal mean (NLM) filtering guided by corresponding CT images Chan et al. (2014) were employed as the reference methods. To evaluate the performance of different methods quantitatively, the contrasttonoise ratio (CNR) regarding the lesion and muscle regions was used as the figure of merit, which is defined as , where and denote the mean intensity inside the lesion ROI and muscle ROIs, respectively, and is the pixeltopixel standard deviation inside the chosen muscle ROIs. Fig. 24 shows the coronal view of one patient data processed using different methods given the optimal parameter for each method. It can be observed that the proposed method has the best visual appearance, which smooths out the noise but also keeps the detailed organ and tumor structures. Table 1 shows the CNRs for all the patient datasets using different methods. The proposed method has the best performance for all the patient datasets. In this example, clinical wholebody PET denoising were employed to show the superior performance of the proposed framework quantitatively. When applied to brain imaging denoising applications, the proposed framework can also perform well. Fig. 2 is a qualitative example. More quantitative evaluations are needed for brain imaging restoration applications.
Patient index  1  2  3  4  5  6  7  8  9  10 

Gaussian  
NLM  
Proposed 
4 Discussion and Conclusion
In this work, we propose a personalized learning framework for inverse problems in medical imaging using deep neural network. Prior training pairs are not needed in this process, but only the patient’s own prior images. An ADMM framework is proposed to decouple the whole optimization into a penalized inverse problem and a network training problem. Brain PET image reconstruction and denoising were employed as examples to demonstrate the effectiveness of this proposed framework. Quantitative results show that the proposed personalized representation framework performs better than other widely adopted methods. One advantage of decoupling the whole optimization is that all current algorithms for inverse problem and network training can be employed, which makes the proposed framework more adoptable. One weakness of the proposed framework is its lack of theoretical convergence guarantee though monotonic convergence is observed in Fig. 7(b). Unet structure is employed in the framework due to its strong representation ability, so that the number of trainable network parameters is less than the original inverse problem. Finding a network structure with better representation ability will improve the performance of the proposed framework. Apart from PET imaging, this proposed framework can also be extended to other inverse problems in medical imaging where personalized priors are available, such as dynamic brain imaging, longitudinal imaging, and multiparametric imaging scenarios. The superior performance of the proposed method is essentially due to a better representation of image. Future work will focus on theoretical convergence analysis and more thorough clinical evaluations.
References
 Wang et al. (2016) Shanshan Wang, Zhenghang Su, Leslie Ying, Xi Peng, Shun Zhu, Feng Liang, Dagan Feng, and Dong Liang. Accelerating magnetic resonance imaging via deep learning. In Biomedical Imaging (ISBI), 2016 IEEE 13th International Symposium on, pages 514–517. IEEE, 2016.
 Sun et al. (2016) Jian Sun, Huibin Li, Zongben Xu, et al. Deep admmnet for compressive sensing mri. In Advances in Neural Information Processing Systems, pages 10–18, 2016.
 Kang et al. (2017) Eunhee Kang, Junhong Min, and Jong Chul Ye. A deep convolutional neural network using directional wavelets for lowdose xray ct reconstruction. Medical physics, 44(10), 2017.
 Chen et al. (2017) Hu Chen, Yi Zhang, Weihua Zhang, Peixi Liao, Ke Li, Jiliu Zhou, and Ge Wang. Lowdose ct via convolutional neural network. Biomedical optics express, 8(2):679–694, 2017.
 Hammernik et al. (2018) Kerstin Hammernik, Teresa Klatzer, Erich Kobler, Michael P Recht, Daniel K Sodickson, Thomas Pock, and Florian Knoll. Learning a variational network for reconstruction of accelerated mri data. Magnetic resonance in medicine, 79(6):3055–3071, 2018.
 Wu et al. (2017) Dufan Wu, Kyungsang Kim, Georges El Fakhri, and Quanzheng Li. Iterative lowdose ct reconstruction with priors trained by artificial neural network. IEEE transactions on medical imaging, 36(12):2479–2486, 2017.
 Adler and Öktem (2018) Jonas Adler and Ozan Öktem. Learned primaldual reconstruction. IEEE transactions on medical imaging, 37(6):1322–1332, 2018.
 Schlemper et al. (2018) Jo Schlemper, Jose Caballero, Joseph V Hajnal, Anthony N Price, and Daniel Rueckert. A deep cascade of convolutional neural networks for dynamic mr image reconstruction. IEEE transactions on Medical Imaging, 37(2):491–503, 2018.
 Gong et al. (2017a) Kuang Gong, Jiahui Guan, Kyungsang Kim, Xuezhu Zhang, Georges El Fakhri, Jinyi Qi, and Quanzheng Li. Iterative pet image reconstruction using convolutional neural network representation. arXiv preprint arXiv:1710.03344, 2017a.
 Kim et al. (2018) Kyungsang Kim, Dufan Wu, Kuang Gong, Joyita Dutta, Jong Hoon Kim, Young Don Son, Hang Keun Kim, Georges El Fakhri, and Quanzheng Li. Penalized pet reconstruction using deep learning prior and local linear fitting. IEEE transactions on medical imaging, 37(6):1478–1487, 2018.
 Chen et al. (2008) GuangHong Chen, Jie Tang, and Shuai 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):660–663, 2008.
 Wang and Qi (2015) Guobao Wang and Jinyi Qi. PET image reconstruction using kernel method. IEEE Transactions on Medical Imaging, 34(1):61–71, 2015.
 Gnahm and Nagel (2015) Christine Gnahm and Armin M Nagel. Anatomically weighted secondorder total variation reconstruction of 23na mri using prior information from 1h mri. NeuroImage, 105:452–461, 2015.
 Gong et al. (2017b) Kuang Gong, Jinxiu ChengLiao, Guobao Wang, Kevin T Chen, Ciprian Catana, and Jinyi Qi. Direct patlak reconstruction from dynamic pet data using the kernel method with mri information based on structural similarity. IEEE Transactions on Medical Imaging, 2017b.
 Ulyanov et al. (2017) Dmitry Ulyanov, Andrea Vedaldi, and Victor Lempitsky. Deep image prior. arXiv preprint arXiv:1711.10925, 2017.
 Goodfellow et al. (2014) Ian Goodfellow, Jean PougetAbadie, Mehdi Mirza, Bing Xu, David WardeFarley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio. Generative adversarial nets. In Advances in neural information processing systems, pages 2672–2680, 2014.
 Radford et al. (2015) Alec Radford, Luke Metz, and Soumith Chintala. Unsupervised representation learning with deep convolutional generative adversarial networks. arXiv preprint arXiv:1511.06434, 2015.
 Mirza and Osindero (2014) Mehdi Mirza and Simon Osindero. Conditional generative adversarial nets. arXiv preprint arXiv:1411.1784, 2014.
 Isola et al. (2017) Phillip Isola, JunYan Zhu, Tinghui Zhou, and Alexei A Efros. Imagetoimage translation with conditional adversarial networks. arXiv preprint, 2017.
 Zhu et al. (2017) JunYan Zhu, Taesung Park, Phillip Isola, and Alexei A Efros. Unpaired imagetoimage translation using cycleconsistent adversarial networks. arXiv preprint arXiv:1703.10593, 2017.
 Çiçek et al. (2016) Özgün Çiçek, Ahmed Abdulkadir, Soeren S Lienkamp, Thomas Brox, and Olaf Ronneberger. 3d unet: learning dense volumetric segmentation from sparse annotation. In International Conference on Medical Image Computing and ComputerAssisted Intervention, pages 424–432. Springer, 2016.
 Gong et al. (2016) Kuang Gong, Stan Majewski, Paul E Kinahan, Robert L Harrison, Brian F Elston, Ravindra Manjeshwar, Sergei Dolinsky, Alexander V Stolin, Julie A BrefczynskiLewis, and Jinyi Qi. Designing a compact high performance brain pet scanner—simulation study. Physics in Medicine & Biology, 61(10):3681, 2016.
 Gong et al. (2017c) Kuang Gong, Jian Zhou, Michel Tohme, Martin Judenhofer, Yongfeng Yang, and Jinyi Qi. Sinogram blurring matrix estimation from point sources measurements with rankone approximation for fully 3D PET. IEEE Transactions on Medical Imaging, 2017c.
 Bai et al. (2013) Bing Bai, Quanzheng Li, and Richard M Leahy. Magnetic resonanceguided positron emission tomography image reconstruction. Seminars in nuclear medicine, 43(1):30–44, 2013.
 Do and Vetterli (2005) Minh N Do and Martin Vetterli. The contourlet transform: an efficient directional multiresolution image representation. IEEE Transactions on image processing, 14(12):2091–2106, 2005.
 Elad and Aharon (2006) Michael Elad and Michal Aharon. Image denoising via sparse and redundant representations over learned dictionaries. IEEE Transactions on Image processing, 15(12):3736–3745, 2006.
 Kingma and Ba (2014) Diederik Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014.
 Nesterov (1983) Yurii Nesterov. A method of solving a convex programming problem with convergence rate o (1/k2). In Soviet Mathematics Doklady, volume 27, pages 372–376, 1983.
 Liu and Nocedal (1989) Dong C Liu and Jorge Nocedal. On the limited memory bfgs method for large scale optimization. Mathematical programming, 45(13):503–528, 1989.
 Lange et al. (2000) Kenneth Lange, David R Hunter, and Ilsoon Yang. Optimization transfer using surrogate objective functions. Journal of computational and graphical statistics, 9(1):1–20, 2000.
 AubertBroche et al. (2006) Berengere AubertBroche, Mark Griffin, G Bruce Pike, Alan C Evans, and D Louis Collins. Twenty new digital brain phantoms for creation of validation image data bases. IEEE transactions on medical imaging, 25(11):1410–1416, 2006.
 Catana et al. (2011) Ciprian Catana, Thomas Benner, Andre van der Kouwe, Larry Byars, Michael Hamm, Daniel B Chonde, Christian J Michel, Georges El Fakhri, Matthias Schmand, and A Gregory Sorensen. Mriassisted pet motion correction for neurologic studies in an integrated mrpet scanner. Journal of Nuclear Medicine, 52(1):154–161, 2011.
 Chan et al. (2014) Chung Chan, Roger Fulton, Robert Barnett, David Dagan Feng, and Steven Meikle. Postreconstruction nonlocal means filtering of wholebody pet with an anatomical prior. IEEE Transactions on medical imaging, 33(3):636–650, 2014.