Improved low-count quantitative PET reconstruction with an iterative neural network

1 Impro v ed lo w-count quantitati v e PET reconstruction with an iterati ve neural netw ork Hongki Lim, Student Member , IEEE, Il Y ong Chun, Member , IEEE, Y uni K. De waraja, Member , IEEE, and Jef frey A. Fessler , F ellow , IEEE Abstract —Image r econstruction in lo w-count PET is particu- larly challenging because gammas from natural radioacti vity in Lu-based crystals cause high random fractions that lower the measurement signal-to-noise-ratio (SNR). In model-based image reconstruction (MBIR), using mor e iterations of an unregularized method may increase the noise, so incorporating regularization into the image reconstruction is desirable to control the noise. New regularization methods based on lear ned con volutional operators are emerging in MBIR. W e modify the architecture of an iterative neural network, BCD-Net , for PET MBIR, and demonstrate the efficacy of the trained BCD-Net using XCA T phantom data that simulates the low true coincidence count- rates with high random fractions typical for Y -90 PET patient imaging after Y -90 microsphere radioembolization. Numerical results show that the proposed BCD-Net significantly impro ves CNR and RMSE of the reconstructed images compared to MBIR methods using non-trained regularizers, total v ariation (TV) and non-local means (NLM). Mor eover , BCD-Net successfully generalizes to test data that differs from the training data. Impro vements wer e also demonstrated for the clinically relev ant phantom measurement data where we used training and testing datasets ha ving very different activity distributions and count- levels. Index T erms —Iterative neural network, Regularized model- based image reconstruction, Low-count quantitative PET , Y -90 I . I N T R O D U C T I O N Image reconstruction in low-count PET is particularly chal- lenging because dominant gammas from natural radioactivity in Lu-based crystals cause high random fractions, lowering the measurement signal-to-noise-ratio (SNR) [1]. T o accurately reconstruct images in low-count PET , regularized model-based image reconstruction (MBIR) solves the following optimiza- tion problem consisting of 1) a data fidelity term f ( x ) that models the physical PET imaging system, and 2) a regulariza- tion term R ( x ) that penalizes image roughness and controls noise [2]: ˆ x = arg min x ≥ 0 f ( x ) + R ( x ) (1) f ( x ) := 1 T ( Ax + ¯ r ) − y T log( Ax + ¯ r ) . This work was supported by grant R01 EB022075, awarded by the National Institute of Biomedical Imaging and Bioengineering, National Institute of Health, U.S. Department of Health and Human Services. Hongki Lim and Jeffre y A. Fessler are with the Department of Electrical Engineering and Computer Science, University of Michigan, Ann Arbor , MI 48109 USA (email: { hongki, fessler } @umich.edu). Il Y ong Chun is with the Department of Electrical Engineering, University of Hawai’i–M ¯ anoa, HI 96822 USA (iychun@hawaii.edu). Y uni K. Dewaraja is with the Department of Radiology , University of Michigan, Ann Arbor, MI 48109 USA (yuni@med.umich.edu). Here, f ( x ) is the Poisson neg ativ e log-likelihood for measure- ment y and estimated measurement means ¯ y ( x ) = Ax + ¯ r , the matrix A denotes the system model, and ¯ r denotes the mean background ev ents such as scatter and random coincidences. Recently , applying learned regularizers to R ( x ) is emerging for MBIR [3]. While there is much ongoing research on machine learn- ing or deep-learning techniques applied to CT [4]–[8] and MRI [9]–[13] reconstruction problems, fewer studies have applied these techniques to PET . Most past PET studies used deep learning in image space without exploiting the physical imaging model in (1). For example, [14] applied a deep neural network (NN) mapping between reconstructed PET images with normal dose and reduced dose and [15] applied a multilayer perceptron mapping between reconstructed images using maximum a posteriori algorithm and a reference (true) image, and their frame work uses the acquisition data only to form the initial image. Recently , [16] trained a NN to reconstruct a 2D image directly from PET sinogram and [17], [18] proposed a PET MBIR framework using a deep- learning based regularizer . Our proposed MBIR framew ork, BCD-Net , also uses a regularizer that penalizes differences between the unkno wn image and “denoised” images gi ven by a regression neural network in an iterati ve manner . In particular , whereas [17], [18] trained only a single image denoising NN, the proposed method is an iterativ e framew ork that includes multiple trained NNs. This iterative frame work enables the NNs in the later stages to learn how to recover fine details. Our proposed BCD-Net also differs from [17], [18] in that our denoising NNs are defined by an optimization formulation with a mathematical motiv ation (whereas, for the trained reg- ularizer , [17], [18] brought U-Net [19] and DnCNN that were [20] dev eloped for medical image segmentation and general Gaussian denoising, respectively) and characterized by fewer parameters, thereby avoiding o ver -fitting and generalizing well to unseen data especially when training samples are limited. Iterativ e NNs [8]–[11], [21]–[24] are a broad family of methods that originate from an unrolling algorithm for solving an optimization problem and BCD-Net [25] is a specific example of an iterative NN. BCD-Net is constructed by unfolding a block coordinate descent (BCD) MBIR algorithm using “learned” con v olutional analysis operators [26]–[28], leading to significantly improv ed image recovery accuracy in extreme imaging applications, e.g., low-dose CT [29], dual- energy CT [30], highly undersampled MRI [25], denoising low-SNR images [25], etc. A preliminary version of this paper was presented at the 2018 Nuclear Science Symposium and Copyright c  2019 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works. 2 Medical Imaging Conference [31]. W e significantly extended this work by applying our proposed method to measured PET data with ne wly dev eloped techniques. W e also added detailed analysis of our proposed method as well as comparisons to related works. T o show the ef ficacy of our proposed BCD-Net method in low-count PET imaging, we performed both digital phantom simulation and e xperimental measurement studies with acti vity distributions and count-rates that are relev ant to clinical Y -90 PET imaging after li ver radioembolization. No vel therapeutic applications have sparked gro wing interest in quantitative imaging of Y -90, an almost pure beta emitter that is widely used in internal radionuclide therapy . In addition to the FDA approv ed Y -90 microsphere radioembolization and Y -90 ibritu- momab radioimmunotherapy , there are 50 active clinical trials for Y -90 labeled therapies (www .clinicaltrials.go v). Ho we ver , the lack of gamma photons complicates imaging of Y -90; it inv olves SPECT via bremsstrahlung photons produced by the betas [32] or PET via a very low abundance positron in the presence of bremsstrahlung that leads to lo w signal-to- noise [33]. This paper applies a BCD-Net that is trained for realistic lo w-count PET imaging environments and compares its performance with those of non-trained re gularizers. Our proposed BCD-Net applies to PET imaging in general, partic- ularly in other imaging situations that also have lo w counts. Using shorter scan times and lower tracer acti vity in diagnostic PET has cost benefits and reduces radiation exposure, but at the expense of reduced counts that makes traditional iterative reconstruction challenging. Section II develops the proposed BCD-Net architecture for PET MBIR. Section II also explains the simulation studies in the setting of Y -90 radioembolization and provides details on how we perform the physical phantom measurement. Section III presents how the different reconstruction methods perform with the simulation and measurement data. Section IV discusses what training and imaging factors most affect gener- alization performance of BCD-Net. Section V concludes with future works. I I . M E T H O D S This section presents the problem formulation of the BCD- Net and giv es a detailed deri vation that inspires the final form of BCD-Net. W e also provide several techniques for BCD- Net that we specifically de vised for PET data where each measurement has dif ferent count-le vel (and noise-lev el). Then we re view the related works that we compare with BCD- Net such as MBIR methods using con v entional non-trained regularizers. This section also describes the simulation setting and details on the measurement data and what e valuation metrics are used to assess the efficac y of each reconstruction algorithm. A. BCD algorithm for MBIR using “learned” con volutional r e gularization Con ventional PET regularizers penalize differences between neighboring pixels [34]. That approach is equiv alent to assum- ing that conv olving the image with the [1,-1] finite difference Fig. 1. Architecture of the proposed BCD-Net for PET . The proposed BCD- Net has an iterativ e NN architecture: each BCD-Net iteration uses three inputs – fixed measurement and mean background { y , ¯ r } , and the image x ( n − 1) reconstructed at the previous BCD-Net iteration – and provides the reconstructed image x ( n ) . A circular arrow abov e MAP EM update indicates inner iterations. g 1 ( · ) and g 2 ( · ) are the normalization and scaling functions described in Section II-C. filter along different directions produces sparse outputs. Using such “hand-crafted” filters is unlikely to be the best approach. A recent trend is to use training data to learn filters c k that produce sparse outputs when con v olved with images of interest [26], [27], [35], [36]. Such learned filters can be used to define a re gularizer that prefers images ha ving sparse outputs, as follows [37]: R ( x ) = min { z k } β K X k =1 1 2 k c k ∗ x − z k k 2 2 + α k k z k k 1 ! , (2) where β is regularization parameter , { c k ∈ R R : k = 1 , . . . , K } is a set of conv olutional filters, { z k ∈ R n p : k = 1 , . . . , K } is a set of sparse codes, { α k ∈ R : k = 1 , . . . , K } is a set of thresholding parameters controlling the sparsity of { z k } , n p is the number of image vox els, and R and K denote the size and number of learned filters, respecti vely . BCD-Net is inspired by this type of “learned” regularizer . Ultimately , we hope that the learned regularizer can better separate true signal from noisy components compared to hand-crafted filters [29]. A natural BCD algorithm solves (1) with regularizer (2) by alternativ ely updating { z k } and x : { z ( n +1) k } = arg min { z k } 1 2   c k ∗ x ( n ) − z k   2 2 + α k k z k k 1 = T ( c k ∗ x ( n ) , α k ) (3) x ( n +1) = arg min x ≥ 0 f ( x ) + β 2 K X k =1    c k ∗ x − z ( n +1) k    2 2 ! , (4) 3 where T ( · , · ) is the element-wise soft thresholding operator: T ( t , q ) j := sign( t j ) max( | t j | − q, 0) . Assuming that learned filters { c k } satisfy the tight-frame condition, P K k =1 k c k ∗ x k 2 2 = k x k 2 2 ∀ x [26], we rewrite the updates in (3)-(4) as follows: u ( n +1) = K X k =1 ˜ c k ∗  T  c k ∗ x ( n ) , α k  (5) x ( n +1) = arg min x ≥ 0 f ( x ) + β 2    x − u ( n +1)    2 2 , (6) where ˜ c k denotes a rotated version of c k . The operations of con volution, soft thresholding and then filtering again with summation typically hav e the effect of denoising the image x ( n ) . For efficient image reconstruction (6) in PET , we use the standard EM-surrogate of Poisson log-likelihood function [38]: f ( x ) + β 2    x − u ( n +1)    2 2 = n d X i =1 [ Ax ] i + ¯ r i − y i log([ Ax ] i + ¯ r i ) + β 2 n p X j =1 ( x j − u ( n +1) j ) 2 ≤ n p X j =1  − e j ( x ( n 0 ) )( x ( n 0 ) j ) log( x j ) + a j x j + β 2 ( x j − u ( n +1) j ) 2  = n p X j =1 Q j ( x j ) where n 0 denotes n 0 th inner-iteration in (6), e j ( x ( n 0 ) ) = P n d i =1 a ij y i ¯ y i ( x ( n 0 ) ) , a ij denotes an element of the system model at i th ro w and j th column, and n d is the number of rays. Equating ∂ Q j ( x j ) ∂ x j to zero is equiv alent to finding the root of the following quadratic formula: β x 2 j +  a j − β u ( n +1) j  x j − e j ( x ( n 0 ) ) x ( n 0 ) j = 0 , and finding the root [39] leads to the minimizer: x ( n 0 +1) j =    √ λ 2 + β ν − λ β , λ < 0 ν √ λ 2 + β ν + λ , λ ≥ 0 , where λ = 1 2 ( a j − β u ( n +1) j ) , ν = e j ( x ( n 0 ) ) x ( n 0 ) j , a j = P n d i =1 a ij . B. BCD-Net for PET MBIR and training its denoising module T o further improve denoising capability by providing more trainable parameters, we extend the conv olutional image de- noiser (CID) in (5) [25], by replacing { ˜ c k } with separate decoding filters { d k } . W e define BCD-Net to use the follo wing updates for each iteration: u ( n +1) = K X k =1 d ( n +1) k ∗  T  c ( n +1) k ∗ x ( n ) , α ( n +1) k  (7) x ( n +1) = arg min x ≥ 0 f ( x ) + β 2    x − u ( n +1)    2 2 , (8) Algorithm 1 BCD-Net for PET MBIR Require: { c ( n ) k , d ( n ) k , α ( n ) k : n = 1 , . . . , T } , y , ¯ r , A , c Initialize: x (0) using EM algorithm Calculate a j = P n d i =1 a ij for n = 0 , . . . , T − 1 do u ( n +1) = P K k =1 d ( n +1) k ∗  T  c ( n +1) k ∗ g 1  x ( n )  , α ( n +1) k  β ( n +1) =     a j − P n d i =1 a ij y i ¯ y i ( x ( n ) )     2 k x ( n ) − g 2 ( u ( n +1) ) k 2 · c for n 0 = 0 , . . . , T 0 − 1 do λ = 1 2  a j − β ( n +1) g 2  u ( n +1) j  ν = x ( n 0 ) j  P n d i =1 a ij y i ¯ y i ( x ( n 0 ) )  x ( n 0 +1) j =    √ λ 2 + β ( n +1) ν − λ β ( n +1) , λ < 0 ν √ λ 2 + β ( n +1) ν + λ , λ ≥ 0 end for x ( n +1) = x ( T 0 ) end for where separate encoding and decoding filters { c k } and { d k } are learned for each iteration. Fig. 1 shows the corresponding BCD-Net architecture. W e refer to the u and x updates in (7)-(8) as two modules : 1) image denoising module and 2) image reconstruction module. The final output image is from the reconstruction module. The image denoising module consists of encoding and decoding filters { c ( n +1) k } , { d ( n +1) k } and thresholding values { α ( n +1) k } . W e train these parameters to “best map” from noisy images into high-quality reference images (e.g., true images if av ailable) in the sense of mean squared error: arg min { c k } , { d k } , { α k } L X l =1      x true ,l − K X k =1 d k ∗  T  c k ∗ x ( n ) l , α k       2 2 , (9) where L is the total number of training samples, { x true ,l ∈ R n p : l = 1 , . . . , L } is a set of true images and { x ( n ) l ∈ R n p : l = 1 , . . . , L } is a set of images estimated by image reconstruction module in the n th iteration. W e train the set of filters and thresholding values iteration-by-iteration and do not include the system matrix or sinograms for training as shown in (9). Moreover , we do not enforce the tight-frame condition when training the filters. One can further extend the CID in (7) to a general regression NN, e.g., a deep U-Net [19]. W e inv estigated if the iterative BCD-Net combined with U-Net denoisers (by replacing the denoising module in (7) with a U-Net) performs better than the proposed BCD-Net using CID (7). Section II-G2 giv es the details of the U-Net implementation. C. Adaptive BCD-Net generalizing to various count-levels 1) Normalization and scaling scheme: Different PET im- ages can hav e very different intensity v alues due to v ariations 4 in scan time and acti vity , and it is important for trained meth- ods to be able to generalize to a wide range of count lev els. T ow ards this end, we implemented normalization and scaling techniques in BCD-Net. [18] extended [17] by implementing “local linear fitting” to ensure that the denoising NN output has similar intensity as the input patch from the current estimated image. Our approach is different in that we normalize and scale the image with a global approach, not a patch-based approach. In particular , we modify the architecture in (7)-(8) as: u ( n +1) = K X k =1 d ( n +1) k ∗  T α ( n +1) k  c ( n +1) k ∗ g 1 ( x ( n ) )  (10) x ( n +1) = arg min x ≥ 0 f ( x ) + β 2    x − g 2 ( u ( n +1) )    2 2 , (11) where the normalization function g 1 ( · ) is defined by g 1 ( v ) := 1 P j v j v to ensure that 1 T g 1 ( v ) = 1 , and the scaling function g 2 ( · ) is defined by g 2 ( v ) := { arg min s f ( s · v ) } v . W e solve the optimization problem over s using Newton’ s method: s ( n +1) = s ( n ) − ∇ s f ( s ( n ) · v ) ∇ 2 s f ( s ( n ) · v ) = s ( n ) − P n d i =1 [ Av ] i − y i [ Av ] i s ( n ) [ Av ] i + ¯ r i P n d i =1 y i  [ Av ] i ( s ( n ) [ Av ] i + ¯ r i )  2 . (12) T o be consistent with the modified CID in (10), we also apply this image-based normalization technique when training the con volutional filters and thresholding values: arg min { c k } , { d k } , { α k } L X l =1      g 1 ( x true ,l ) − K X k =1 d k ∗  T α k  c k ∗ g 1 ( x ( n ) l )       2 2 . 2) Adaptive r e gularization parameter scheme: The best regularization parameter value can also vary greatly between scans, depending on the count le vel. Therefore, instead of choosing one specific value for the regularization parameter, we set the β value for each iteration based on ev aluation on current gradients of data-fidelity term and re gularization term: β ( n +1) =   ∇ x f ( x ( n ) )   2   ∇ x R ( x ( n ) )   2 · c =   a j − e j ( x ( n ) )   2   x ( n ) − g 2 ( u ( n +1) )   2 · c, n = 0 , . . . , T − 1 , (13) where c is a constant specifying ho w we balance between the data-fidelity term and regularization term and n denotes n th outer-iteration. Algorithm (1) gives detailed pseudocode of the proposed method. T denotes the total number of outer- iterations and T 0 denotes the number of inner iterations used for (8). W e use x ( n ) as the initial image when solving (11). D. Conventional MBIR methods: Non-trained re gularizer s W e compared the proposed BCD-Net with two MBIR methods that use standard non-trained regularizers. T ABLE I D E T A IL S O N X C A T SI M U L A TI O N DA TA : V A R IATI O N S BE T W E EN T R AI N I N G A N D T E S TI N G DA TA . T raining data T esting data Concentration ratio (hot:warm) 9:1 4:1 T otal net trues 200 K 500 K Random fraction (%) 90.9 87.5 T ABLE II D E T A IL S O N P H A N TO M M EA S U RE M E N T D A T A : AC TI V I T Y C O N CE N T R A TI O N R A T IO B E T WE E N HO T A N D W AR M R E GI O N S A N D RA N D O MS F R AC T IO N S F O R T W O P H A N TO M ST U D I ES . Sphere Liv er-torso T otal activity (GBq) 0.65 1.9 Concentration ratio (hot:warm) 8.9:1 5.4:1 T otal prompts 3.2 - 6.3 M 2.3 M T otal randoms 2.9 - 5.7 M 2.1 M T otal net trues 308 - 599 K 220 K Random fraction(%) 90.3 - 90.5 90.7 1) T otal-variation (TV): TV regularization penalizes the sum of absolute v alue of differences between adjacent voxels: R ( x ) = β k C x k 1 , where C is finite differencing matrix. Recent work [40] applied Primal-Dual Hybrid Gradient (PDHG) [41] for PET MBIR using TV regularization and demonstrated that PDHG- TV is superior than clinical reconstruction (e.g., OS-EM) for low-count datasets in terms of sev eral image quality ev aluation metrics such as contrast recovery and variability . 2) Non-local means (NLM): NLM regularization penalizes the differences between nearby patches in image: R ( x ) = β X i,j ∈ S i p  k N i x − N j x k 2 2  , where p ( t ) is a potential function of a scalar variable t , S i is the search neighborhood around the i th vox el, and N i is a patch extraction operator at the i th vox el. W e used the Fair potential function for p ( t ) : p ( t ) = σ 2 f s t σ 2 f N f + log 1 + s t σ 2 f N f !! , where σ f is a design parameter and N f is the number of vox els in the patch N i x . Unlike con ventional local filters that assume similarity between only adjacent vox els, NLM filters can av erage image intensities ov er distant voxels. As in [42], we used ADMM to accelerate algorithmic con ver gence with an adaptiv e penalty parameter selection method [43]. E. Experimental setup: Digital phantom simulation and ex- perimental measur ement 1) Y -90 PET/CT XCA T simulations: W e used the XCA T [44] phantom (Fig. 2) to simulate Y -90 PET follo wing ra- dioembolization. W e set the image size to 128 × 128 × 100 with a vox el size 4.0 × 4.0 × 4.0 (mm 3 ) and chose 100 slices ranging from lung to liver . T o simulate extremely low count scans with high random fractions, typical for Y -90 PET , we set total true coincidences and random fractions based on numbers from 5 T ABLE III D E T A IL S O N T Y P I CA L PA TI E N T M E A S UR E M E NT DAT A : TOTA L TR UE S A ND R A ND O M S F R AC T I ON S . Patient A T otal activity (GBq) 2.55 T otal prompts 2.7 M T otal randoms 2.3 M T otal net trues 380 K Random fraction(%) 85.8 patient PET imaging performed after radioembolization [45]. T o test the generalization capability of the trained BCD-Net, we changed all imaging factors between training and testing dataset. Here, imaging factors include activity distribution (shape and size of tumor and liv er background, concentration ratio between hot and warm region) and count-lev el (total true coincidences and random fraction). Fig. 2 and T able I provide details on ho w we changed the testing dataset from the training dataset. W e trained BCD-Net using five pairs ( L = 5 ) of 3D true images and estimated images at each iteration (1 true image, 5 realizations). W e generated multiple (5) realizations to train the denoising NN to deal with the Poisson noise. W e also generated 5 realizations (1 true image, 5 realizations) as a testing dataset to ev aluate the noise across realizations. 2) Y90 PET/CT physical phantom measur ements and pa- tient scan: For training BCD-Net, we used PET measure- ments of a sphere phantom (Fig. 4) where six ‘hot’ spheres (2,4,8,16,30 and 113 mL, 0.5 MBq/ml) are placed in a ‘warm’ background (0.057 MBq/ml) with total activity of 0.65 GBq. The phantom was scanned for 40 (3 acquisitions) - 80 (1 acquisition) ( L = 4) minutes on a Siemens Biograph mCT PET/CT . F or testing BCD-Net and other reconstruction algo- rithms, we used an anthropomorphic liv er/lung torso phantom (Fig. 4) with total activity and distribution that is clinically realistic for imaging following radioembolization with Y -90 microspheres: 5% lung shunt, 1.17 MBq/mL in liver , 3 hepatic lesions (4 and 16 mL spheres, 29 mL ovoid) of 6.6 MBq/ml. The phantom with total activity of 1.9 GBq was scanned 5 times (each 30 minutes) on a Siemens Biograph mCT PET/CT . Fig. 4 and T able II provide details on the count-level (random fraction) and activity distrib ution dif ferences between training (sphere phantom) and testing (li ver phantom) dataset. W e also tested BCD-Net with an actual Y -90 patient scan and T able III provides count-lev el information. W e acquired all measurement data with time of flight TOF information. The measurement data size is 200 × 168 × 621 × 13. The last dimension of measurement indicates the number of time bin. The reconstructed image size is 200 × 200 × 112 with a voxel size 4.07 × 4.07 × 2.03 (mm 3 ). T o reconstruct the image with measurement data, we used a SIEMENS TOF system model ( A in (1)) along with manufacturer giv en attenuation/normalization correction, PSF modelling, and randoms/scatters estimation. F . Evaluation metrics For the XCA T phantom simulation, we ev aluated each reconstruction with contrast recov ery (CR) (volume-of-interest (V OI): cold spot indicated in Fig. 2), noise across realizations, root mean squared error (RMSE), and contrast to noise ratio (CNR). For the physical phantom measurement, we used CR (V OI: hot spheres) and CNR averaged over multiple hot spheres. W e define each V OI’ s mask based on attenuation map interpolated to PET vox el size. For the patient measurement, we used CNR and the field of view (FO V) activity bias since the total acti vity in FO V is known (equal to the injected activity because the microspheres are trapped) wheareas the activity distribution is unknown: CR (V OI: cold spot) =  1 − C V OI C BKG  × 100 (%) CR (V OI: hot sphere) = C VOI C BKG − 1 R True − 1 × 100 (%) Noise = v u u u u u t 1 J Liver P j ∈ Li ver  1 M − 1 P M m =1 ( ˆ x m [ j ] − 1 M P M m 0 =1 ˆ x m 0 [ j ]) 2  1 J Liver P j ∈ Li ver x true [ j ] × 100 % RMSE = s P j ( x true [ j ] − ˆ x [ j ]) 2 J FO V × 100 (%) CNR = C Lesion − C BKG STD BKG FO V bias = P j ˆ x [ j ] − x true [ j ] P j x true [ j ] × 100 (%) , where C V OI is mean counts in the V OI, R True is true ratio between hot and warm region, x [ j ] denotes the j th vox el of an image x , M is the number of realizations ( M = 5 in both XCA T phantom simulation and physical phantom measurement) and J Liv er is the number of vox els in the volume of li ver , STD BKG is standard deviation between vox el v alues in uniform background liv er (indicated in Fig. 2), and J FO V is the total number of voxels in the FO V . As the background region when calculating the patient CNR, we used a part of liv er region that has relatively uniform activity distribution. G. T raining details W e trained the denoising network in each iteration with a stochastic gradient descent method using the PyT orch [46] deep-learning library . 1) BCD-Net with CID: W e trained a set of CID for the denoising module in BCD-Net where each iteration has 78 sets of thresholding values and con volutional encoding/decoding filters ( K = 78) . W e set the size of each filter as 3 × 3 × 3 ( R = 3 3 ) , and set the initial thresholding values by sorting the initial estimate of image and getting a 10% largest value of sorted initial image. W e used the Adam optimization method [47] to train the NN. W e applied the learning rate decay scheme. Due to the large size of 3D input, we set the batch size as 1. 2) BCD-Net with U-Net: W e implemented a 3-D version of U-Net by modifying a shared code 1 (implemented for denoising 2-D MRI images) for fastMRI challenge [48]. W e used a batch normalization layer instead of the instance normalization layer used in the baseline code. The ‘encoder’ part of U-Net consists of multiple sets of 1) max pooling layer , 2) 3 × 3 × 3 conv olutional layer, 3) batch normalization (BN) 1 https://github .com/facebookresearch/fastMRI 6 Attenuation map (coronal) Attenuation map (axial) T rue acti vity (training) T rue acti vity (testing) Zoomed in 1 128 1 100 0 0.0192 Attenuation Coefficients 1 128 1 128 0 0.0192 Attenuation Coefficients 1 128 1 128 0 9.2 Relative Activity 1 128 1 128 0 4.2 Relative Activity 1 51 1 51 0 4.2 Relative Activity Background liver Cold spot EM TV NLM BCD-Net-CID BCD-Net-UNet BCD-Net-UNet (params: 4K) (params: 4K) (params: 1.4M) 1 51 1 51 0 4.2 Relative Activity 1 51 1 51 0 4.2 Relative Activity 1 51 1 51 0 4.2 Relative Activity 1 51 1 51 0 4.2 Relative Activity 1 51 1 51 0 4.2 Relative Activity 1 51 1 51 0 4.2 Relative Activity Fig. 2. XCA T phantom simulation: (First ro w) coronal and axial view of attenuation map and true relative activity distribution corresponding to axial attenuation map. (Second row) reconstructed images of one slice from different reconstruction methods. BCD-Net-CID/UNet is the BCD-Net with CID/UNet and params indicates the number of trainable parameters. 0 70 Cold Spot Contrast Recovery (%) 0 20 40 Noise EM PDHG-TV ADMM-NLM BCD-Net-CID (params: 4K) BCD-Net-UNET (params: 4K) BCD-Net-UNET (params: 1.4M) 0 20 40 Iterations 6 9 12 RMSE EM PDHG-TV ADMM-NLM BCD-Net-CID (params: 4K) BCD-Net-UNET (params: 4K) BCD-Net-UNET (params: 1.4M) 0 20 40 Iterations 5 10 CNR EM PDHG-TV ADMM-NLM BCD-Net-CID (params: 4K) BCD-Net-UNET (params: 4K) BCD-Net-UNET (params: 1.4M) (a) (b) (c) Fig. 3. (a) Plot of noise in background liv er vs contrast recovery in cold spot (b) RMSE vs iteration (c) Contrast to noise ratio vs iteration. W e initialized regularized methods with the 10th iterate of EM reconstruction. layer , 4) ReLU layer and the ‘decoder’ part of U-Net consists of multiple sets of 1) upsampling with trilinear interpolation [17], 2) 3 × 3 × 3 con v olutional layer , 3) BN layer, 4) ReLU layer . For training the U-Net, we used the same training dataset that we used for training the CID. W e also used the Adam optimization method and identical settings (number of epochs, learning rate decay , batch size) as those of the CID. W e trained and tested two different U-Nets sizes. At each BCD- Net iteration, the U-Net has either about 4 K (similar size to the CID) or 1.4 M trainable parameters. W e set the number of conv olutional filter channels of the first encoder layer as 12 with 4 times of contraction/expansion for the U-Net with 1.4 M parameters and 5 with 1 time of contraction/expansion for the U-Net with 4 K parameters. I I I . R E S U LT S A. Reconstruction setup W e compared the proposed BCD-Net method to the stan- dard EM (1 subset), TV -based MBIR with PDHG algorithm (PDHG-TV), and NLM-based MBIR with ADMM algorithm (ADMM-NLM). For regularized MBIR methods including BCD-Net, we used 10 EM algorithm iterations to get the initial image x (0) . For each regularization method, we finely tuned the regularization parameter β (within range [2 − 15 , 2 15 ] ) by considering the recovery accuracy and noise. For NLM, we additionally tuned the window and search sizes. For the XCA T simulation data, we used 40 iterations for EM and 30 iterations ( T = 30) for PDHG-TV , ADMM-NLM, and BCD-Net. W e used 1 inner-iteration ( T 0 = 1) for the reconstruction module (11) for each outer-iteration of BCD-Net. For the measured data, we used 20 iterations for EM and 10 iterations ( T = 10) for PDHG-TV , ADMM-NLM, and BCD-Net. W e used 1 inner- iteration ( T 0 = 1) for the reconstruction module (11). W e set c = 0 . 01 in (13) in the XCA T simulation study and c = 0 . 005 in both the phantom measurement and patient studies. B. Results: Reconstruction (testing) on simulation data Fig. 2-3 shows that the proposed iterativ e NN, BCD- Net, significantly improv es ov erall reconstruction performance ov er the other non-trained regularized MBIR methods. Fig. 3 reports av eraged ev aluation metrics ov er realizations. Fig. 3 shows that BCD-Net with a trained CID achie ves the best results in most ev aluation metrics. In particular , BCD-Net with a CID improv es CNR and RMSE compared to PDHG-TV and ADMM-NLM. BCD-Net also improved contrast recov ery in the cold region while not increasing noise compared to the 7 Sphere phantom attenuation map (coronal) Attenuation map (axial) T rue acti vity image x (0) of regularized methods (EM) 0 0.01 Attenuation Coefficient 0 0.01 Attenuation Coefficient Liv er phantom attenuation map (coronal) Attenuation map (axial) T rue acti vity image x (0) of regularized methods (EM) 0 0.01 Attenuation Coefficient 0 0.01 Attenuation Coefficient TV NLM BCD-Net-CID BCD-Net-UNet BCD-Net-UNet (params: 4K) (params: 4K) (params: 1.4M) Fig. 4. Y90 PET/CT ph ysical phantom measurement: (First row: training data, Second ro w: testing data) Attenuation map, true acti vity , and x (0) of regularized methods of sphere and liv er phantom used for training and testing BCD-Net. (Third row) Reconstructed images of one slice from different reconstruction methods. 0 50 100 Hot Spot Contrast Recovery (%) 0 35 70 Noise EM PDHG-TV ADMM-NLM BCD-Net-CID (params: 4K) BCD-Net-UNET (params: 4K) BCD-Net-UNET (params: 1.4M) 0 10 20 Iterations 0 10 20 CNR EM PDHG-TV ADMM-NLM BCD-Net-CID (params: 4K) BCD-Net-UNET (params: 4K) BCD-Net-UNET (params: 1.4M) (a) (b) Fig. 5. Liv er phantom measurement: (a) Plot of noise in background liv er vs contrast recov ery in hot spheres (b) Contrast to noise ratio vs iteration. W e initialized regularized methods with the 10th iterate of EM reconstruction. initial EM reconstruction, whereas PHDG-TV and ADMM- NLM improv ed noise while de grading the CR. For Fig. 2, we selected the iteration number for EM to obtain the highest CNR and the last iteration number for other methods. Fig. 2 shows that BCD-Net’ s reconstructed image with a CID is closest to the true image whereas PHDG-TV and ADMM- NLM exceedingly blur the cold region. BCD-Net with the U-Net denoiser shows good recovery for the cold region, howe v er , it blurs the hot region. Moreover , the larger sized U-Net (params: 1.4 M) denoiser worsens the performance of BCD-Net possibly due to over -fitting the training dataset. C. Results: Reconstruction (testing) on measurement data 1) Phantom study: Similar to the simulation results, Fig. 4-5 shows that, BCD-Net improved overall reconstruction per- formance o ver the other reconstruction methods. Fig. 4 shows that reconstructed images using PHDG-TV and ADMM-NLM show uniform texture in background liv er compared to EM, howe v er , those exceedingly blur around hot spheres. The blurred hot region is more evident in the quantification results in Fig. 5. BCD-Net giv es more visibility for hot spheres with noisier texture in uniform liv er region. Fig. 5 shows that BCD-Net with a CID improves CNR compared to PDHG- TV and ADMM-NLM. BCD-Net with CID also impro ved contrast recovery in hot spheres while slightly increasing noise compared to the initial EM reconstruction. In Fig. 5 (a), BCD- 8 Patient attenuation map (coronal) Attenuation map (axial) OSEM w/ filter (coronal) OSEM w/ filter (axial) 0 0.0127 Attenuation Coefficient 0 0.0133 Attenuation Coefficient TV (coronal) TV (axial) NLM (coronal) NLM (axial) BCD-Net-CID (coronal) BCD-Net-CID (axial) BCD-Net-UNet (coronal) BCD-Net-UNet (axial) Fig. 6. Y90 PET/CT patient measurement: Attenuation map and reconstructed images of one slice (coronal and axial view) using OSEM, TV , NLM, and BCD-Net. W e visualized the reconstructed image of BCD-Net-UNet with 4 K parameters 0 10 20 Iterations 0 50 100 FOV Bias EM PDHG-TV ADMM-NLM BCD-Net-CID (params: 4K) BCD-Net-UNET (params: 4K) BCD-Net-UNET (params: 1.4M) 0 10 20 Iterations 5 25 CNR EM PDHG-TV ADMM-NLM BCD-Net-CID (params: 4K) BCD-Net-UNET (params: 4K) BCD-Net-UNET (params: 1.4M) Fig. 7. Patient scan: (a) Field of view bias vs iteration. BCD-Net shows similar results compared to other methods. (b) Contrast to noise ratio vs iteration. Net with U-Net denoiser shows a fluctuation with iterations, howe v er , the plot trend is similar to that of BCD-Net with CID. 2) P atient study: Because of the unkno wn true activity distribution, we quantitatively ev aluated each reconstruction method with FO V acti vity bias. In this quantitati ve ev aluation, BCD-Net showed similar results compared to other methods. See Fig. 6-7. Fig. 6 sho ws that the quality of image using different methods in patient study is similar to that of phantom measurement study shown in Fig. 4. Fig. 7 (b) shows that the CNR trend in the patient study is similar to that of the XCA T simulation and the liv er phantom measurement. I V . D I S C U S S I O N In this study we showed the efficac y of trained BCD-Net on both qualitativ e and quantitativ e Y -90 PET/CT imaging and compared between con ventional non-trained regularizers. The proposed approach uses learned denoising NNs to lift estimated signals and thresholding operations to remove un- wanted signals. In particular , the iterative framework of BCD- Net enables one to train the filters and thresholding values to deal with the different image roughness at its each iteration. W e experimentally demonstrate its generalization capabilities with simulation and measurement data. In the XCA T PET/CT simulation with activity distrib utions and count-rates mimick- ing Y -90 PET imaging, total counts in the cold spot were ov erestimated with standard reconstruction and other MBIR methods using non-trained regularization, yet approached the 9 (a) Impact of number of filter (b) Impact of size of filter (Fixed R = 3 3 ) (Fixed K = 2 5 ) 0 15 30 Iterations 5 7 9 RMSE Training dataset K = 2 4 K = 2 5 K = 2 6 K = 2 7 K = 2 7 (l1 loss) 0 15 30 Iterations 5 7 9 RMSE Testing dataset 0 15 30 Iterations 5 7 9 RMSE Training dataset R = 3 3 R = 5 3 R = 7 3 R = 9 3 0 15 30 Iterations 5 7 9 RMSE Testing dataset (c) T rained with l 1 loss (d) K = 2 4 (e) R = 9 3 1 51 1 51 0 4.2 Relative Activity 1 51 1 51 0 4.2 Relative Activity 1 51 1 51 0 4.2 Relative Activity Fig. 8. (a)-(b) Impact of number/size of filter and training loss on testing dataset RMSE. (c) Reconstructed image from BCD-Net-CID with filters and thresholding values trained with l 1 -loss. β in training and testing case when c = 0 . 005 5 10 Iterations 29 32 Testing dataset: Liver phantom (Total prompts: 2.3 M) Testing dataset: Patient A (Total prompts: 2.7 M) Testing dataset: Patient B (Total prompts: 3.8 M) Fig. 9. Efficacy of adaptive selection of regularization parameter β . true value with the proposed approach. Impro vements were also demonstrated for the measurement data where we used training and testing datasets ha ving very different acti vity distribution and count-lev els. The architecture and size of denoising NN significantly affect the performance of BCD- Net. In both simulation and measurement experiments, the CID outperformed the U-Net architectures. Using a U-Net with more trainable parameters degraded the performance, especially in the simulation study , due to the small size of dataset. Size of the denoising NN should be set with consideration of training dataset size. W e tested which imaging variable most affects the gener- T ABLE IV I M P AC T O F I M AG I NG V A R IA B L E O N G EN E R A LI Z A T I ON C A P A B IL I T Y OF B C D- N E T - C ID . Changed imaging variable T raining T esting RMSE Drop (%) Identical - 4.74 - Shape and size See Fig. 2 5.49 15.9 Concentration ratio 9:1 4:1 5.55 17.1 Concentration ratio 1.7:1 4:1 5.81 22.5 T rues Count-lev el 2 × 10 5 5 × 10 5 5.01 5.7 T rues Count-lev el 11 × 10 5 5 × 10 5 5.71 20.5 T ABLE V C O MPA R IS O N BE T W E EN P O ST - R EC O N ST R UC T I O N P RO C E SS I N G AN D I T ER ATI V E NN . Method g 2  u (1)  RMSE x (30) RMSE ∆ (%) CID (params: 4 k) 6.87 6.36 7.52 U-Net (params: 4 k) 7.39 6.70 9.30 U-Net (params: 1.4 M) 7.67 7.04 8.15 (a) CID (b) U-Net (c) U-Net (params: 4 K) (params: 4 K) (params: 1.4 M) 1 51 1 51 0 4.2 Relative Activity 1 51 1 51 0 4.2 Relative Activity 1 51 1 51 0 4.2 Relative Activity Fig. 10. g 2 ( u (1) ) generated by (a) CID and (b)-(c) U-Net. Using more parameters degraded the visual quality and RMSE value as in the iterativ e NN. alization performance of the proposed BCD-Net. T able IV shows how BCD-Net performs when training and testing data had the same activity distribution and count-lev el (only difference is Poisson noise) and how the performance of BCD-Net is degraded when each imaging variable is changed between training and testing dataset. W e changed one of three factors (shape and size of tumor and liv er , concentration ratio, count-lev el) in training dataset compared to testing dataset. The result shows that generalization performance of the pro- posed BCD-Net depends largely on all imaging variables. Howe v er , training with higher contrast and lower count-level dataset (compared to testing dataset) gave less degradation of performance compared to the opposite cases. This result suggests that it is better to hav e noisier data in training dataset than testing dataset. In other words, training for extra noise reduction than needed is better than less noise reduction than needed. W e also in vestigated ho w each factor in training of de- noising module (7) impacts the generalization capability of BCD-Net. Fig. 8(a)-(b) show the impact of number and size of filters on performance. Plots show that the proposed BCD- Net achieved lower training RMSE when using larger number and size of filters; howe ver , it did not decrease testing RMSE compared to smaller number and size of filters and BCD- Net with lar ger size of filter exceedingly blurs image thereby resulting in higher RMSE. See Fig. 8(e). W e also tested l 1 10 training loss to see if it improves the performance over the l 2 loss (MSE) in (9). Howe ver , it led to unnaturally piece- wise constant images and details in small cold regions were ignored. Fig. 9 sho ws ho w the re gularization parameter β in (13) changes with iterations in training and testing datasets. The β value in each iteration con v erges to different limits in training and testing cases. The adapti ve scheme automatically increases the β value when the count-lev el decreases. This behavior concurs with the general knowledge that more regularization is needed when the noise-lev el increases. These empirical results underscore the importance of such adaptiv e regularization parameter selection schemes proposed in Section II-C2 in PET imaging. Many related works [5], [14], [15], [49] use single image denoising (deep) NN (e.g., U-Net) as a post-reconstruction processing and we in vestigated how the denoising NN de- tached from the data-fit term performs compared to iterati ve NN. Fig. 10 illustrates u (1) generated by CID and U-Net. As in the iterative NN, using more trainable parameters degraded the visual quality and RMSE value in U-Net case and CID achiev ed better result than U-Net. In all cases, iterativ e NN achiev ed lower RMSE compared to those post-reconstruction processed images as shown in T able V. BCD-Net is trained for a specific number of iterations and its practical use would be akin to how ML-EM is used with a fixed number of iterations in clinical systems. If one is interested in con vergence guarantees with running more iterations, then one can extend the sequence conv ergence guarantee of BCD-Net in [23] by setting the n th adaptiv e denoiser as e D ( n ) = g 2 ( D ( n ) ( g 1 ( x ( n − 1) ))) with some n th denoiser D ( n ) (e.g., CID (7) and U-Net), ∀ n , using sufficient T 0 (so MAP EM finds a critical point), and additionally assuming that β ( n ) con ver ges. W e empirically observed that the β ( n ) tends to con ver ge to some constant in this Y -90 PET as well as another application of Lu-177 SPECT . T o more practically guarantee the con ver gence, one could use training and testing dataset having similar count-lev el and a fixed re gularization parameter value across iterations using an initial estimated image and a corresponding denoised image as follows: β =   ∇ x f ( x (0) )   2   ∇ x R ( x (0) )   2 · c =   a j − e j ( x (0) )   2   x (0) − g 2 ( u (1) )   2 · c. The con ver gence properties depend on additional technical assumptions detailed in [23]. V . C O N C L U S I O N It is important for a “learned” regularizer to hav e gener- alization capability to help ensure good performance when applying it to an unseen dataset. For low-count PET recon- struction, the proposed iterativ e NN, BCD-Net, showed reli- able generalization capability even when the training dataset is small. The proposed BCD-Net achie ved significant qualitativ e and quantitative improvements over the conv entional MBIR methods using “hand-crafted” non-trained regularizers: TV and NLM. In particular , these con v entional MBIR methods hav e a trade-off between noise and recov ery accuracy , whereas the proposed BCD-Net improv es CR for hot regions while not increasing the noise when the regularization parameter is appropriately set. V isual comparisons of the reconstructed images also sho w that the proposed BCD-Net significantly improv es PET image reconstruction performance compared to MBIR methods using non-trained regularizers. Future work includes in vestigating performance of BCD- Net trained with end-to-end training principles and adaptiv e selection of trainable parameter numbers depending on the size of training dataset. V I . A C K N OW L E D G M E N T W e acknowledge Se Y oung Chun (UNIST) for providing NLM re gularization codes. W e also ackno wledge Maurizio Conti and Deepak Bharkhada (SIEMENS Healthcare Molecu- lar Imaging) for providing the forward/back projector for TOF measurement data. This work was supported by NIH-NIBIB grant R01EB022075. R E F E R E N C E S [1] T . Carlier, K. P . Willo wson, E. Fourkal, D. L. Bailey , M. Doss, and M. Conti, “Y90-PET imaging: exploring limitations and accuracy under conditions of low counts and high random fraction, ” Med. Phys. , v ol. 42, no. 7, pp. 4295–309, Jun. 2015. [2] S. Ahn, S. G. Ross, E. Asma, J. Miao, X. Jin, L. Cheng, S. D. W ollenweber, and R. M. Manjeshwar, “Quantitati ve comparison of OSEM and penalized likelihood image reconstruction using relativ e difference penalties for clinical PET, ” Physics in Medicine & Biology , vol. 60, no. 15, p. 5733, 2015. [3] G. W ang, J. C. Y e, K. Mueller , and J. A. Fessler, “Image reconstruction is a ne w frontier of machine learning, ” IEEE T rans. Med. Ima g. , vol. 37, no. 6, pp. 1289–96, Jun. 2018. [4] H. Chen, Y . Zhang, M. K. Kalra, F . Lin, Y . Chen, P . Liao, J. Zhou, and G. W ang, “Low-dose CT with a residual encoder -decoder con volutional neural network, ” IEEE T rans. Med. Imag. , vol. 36, no. 12, pp. 2524– 2535, 2017. [5] K. H. Jin, M. T . McCann, E. Froustey , and M. Unser, “Deep conv o- lutional neural network for in verse problems in imaging, ” IEEE T rans. Image Pr ocess. , vol. 26, no. 9, pp. 4509–4522, 2017. [6] J. C. Y e, Y . Han, and E. Cha, “Deep con v olutional framelets: A general deep learning framework for in verse problems, ” SIAM J ournal on Imaging Sciences , vol. 11, no. 2, pp. 991–1048, 2018. [7] H. Gupta, K. H. Jin, H. Q. Nguyen, M. T . McCann, and M. Unser , “CNN-based projected gradient descent for consistent CT image recon- struction, ” IEEE T r ans. Med. Imag. , vol. 37, no. 6, pp. 1440–1453, 2018. [8] I. Y . Chun, H. Lim, Z. Huang, and J. A. Fessler , “Fast and conv ergent iterativ e signal recovery using trained conv olutional neural networkss, ” in Pr oc. Allerton Conf. on Commun., Contr ol, and Comput. , Allerton, IL, Oct. 2018, pp. 155–159. [9] H. K. Aggarwal, M. P . Mani, and M. Jacob, “MoDL: model-based deep learning architecture for inv erse problems, ” IEEE T r ans. Med. Imag. , vol. 38, no. 2, pp. 394–405, Feb. 2019. [10] K. Hammernik, T . Klatzer , E. Kobler , M. P . Recht, D. K. Sodickson, T . Pock, and F . Knoll, “Learning a v ariational network for reconstruction of accelerated MRI data, ” Magn. Reson. Imaging , vol. 79, no. 6, pp. 3055–3071, 2018. [11] J. Sun, H. Li, Z. Xu et al. , “Deep ADMM-Net for compressive sensing MRI, ” in Pr oc. NIPS , 2016, pp. 10–18. [12] M. Mardani, E. Gong, J. Y . Cheng, S. S. V asanawala, G. Zaharchuk, L. Xing, and J. M. Pauly , “Deep generativ e adversarial neural networks for compressiv e sensing MRI, ” IEEE T rans. Med. Imag. , vol. 38, no. 1, pp. 167–179, 2019. [13] G. Y ang, S. Y u, H. Dong, G. Slabaugh, P . L. Dragotti, X. Y e, F . Liu, S. Arridge, J. Kee gan, Y . Guo et al. , “D A GAN: deep de-aliasing generativ e adversarial networks for fast compressed sensing MRI re- construction, ” IEEE T r ans. Med. Imag. , vol. 37, no. 6, pp. 1310–1321, 2018. 11 [14] J. Xu, E. Gong, J. Pauly , and G. Zaharchuk, “200x low-dose PET reconstruction using deep learning, ” arXiv preprint , 2017. [15] B. Y ang, L. Y ing, and J. T ang, “Artificial Neural Network Enhanced Bayesian PET Image Reconstruction, ” IEEE T rans. Med. Imag. , vol. 37, no. 6, pp. 1297–1309, June 2018. [16] I. Haggstrom, C. R. Schmidtlein, G. Campanella, and T . J. Fuchs, “DeepPET : A deep encoder-decoder network for directly solving the PET image reconstruction inverse problem, ” Med. Im. Anal. , vol. 54, pp. 253–62, May 2019. [17] K. Gong, J. Guan, K. Kim, X. Zhang, J. Y ang, Y . Seo, G. El Fakhri, J. Qi, and Q. Li, “Iterati ve PET image reconstruction using con volutional neural network representation, ” IEEE T rans. Med. Imag . , vol. 38, no. 3, pp. 675–685, 2019. [18] K. Kim, D. W u, K. Gong, J. Dutta, J. H. Kim, Y . D. Son, H. K. Kim, G. El Fakhri, and Q. Li, “Penalized PET reconstruction using deep learning prior and local linear fitting, ” IEEE T r ans. Med. Ima g. , v ol. 37, no. 6, pp. 1478–1487, 2018. [19] O. Ronneberger , P . Fischer , and T . Brox, “U-net: Conv olutional networks for biomedical image segmentation, ” in Proc. Med. Image Compt. and Computer Assist. Interven. (MICCAI) . Springer , 2015, pp. 234–241. [20] K. Zhang, W . Zuo, Y . Chen, D. Meng, and L. Zhang, “Beyond a Gaussian denoiser: Residual learning of deep CNN for image denoising, ” IEEE Tr ans. Image Process. , vol. 26, no. 7, pp. 3142–3155, July 2017. [21] K. Gregor and Y . LeCun, “Learning fast approximations of sparse coding, ” in Pr oc. ICML , 2010, pp. 399–406. [22] Y . Chen and T . Pock, “T rainable nonlinear reaction diffusion: A flexible framew ork for fast and effecti ve image restoration, ” IEEE T rans. P attern Anal. Mach. Intell. , vol. 39, no. 6, pp. 1256–1272, 2017. [23] I. Y . Chun, Z. Huang, H. Lim, and J. A. Fessler, “Momentum-net: Fast and con vergent iterative neural network for in verse problems, ” arXiv pr eprint arXiv:1907.11818 , Jul. 2019. [24] S. Y e, Y . Long, and I. Y . Chun, “Momentum-Net for low-dose CT image reconstruction, ” arXiv pr eprint arXiv:2002.12018 , Feb . 2020. [25] I. Y . Chun and J. A. Fessler, “Deep BCD-net using identical encoding- decoding CNN structures for iterative image recovery , ” in Pr oc. Im- age, V ideo, and Multidim. Signal Pr ocess. (IVMSP) W orkshop , Zagori, Greece, Apr . 2018, pp. 1–5. [26] ——, “Con volutional analysis operator learning: acceleration and con- ver gence, ” IEEE T rans. Image Process. , vol. 29, no. 1, pp. 2108–2122, 2020. [27] I. Y . Chun, D. Hong, B. Adcock, and J. A. Fessler , “Conv olutional analysis operator learning: Dependence on training data, ” IEEE Signal Pr ocess. Lett. , vol. 26, no. 8, pp. 1137–1141, 2019. [28] C. Crockett, D. Hong, I. Y . Chun, and J. A. Fessler, “Incorporating handcrafted filters in con volutional analysis operator learning for ill- posed inv erse problems, ” in Pr oc. IEEE Intl. W orkshop on Compt. Adv . in Multi-Sensor Adaptive Process. (CAMSAP) , Guadeloupe, W est Indies, Dec. 2019, pp. 316–320. [29] I. Y . Chun, X. Zheng, Y . Long, and J. A. Fessler, “BCD-Net for low- dose CT reconstruction: Acceleration, con ver gence, and generalization, ” in Pr oc. Med. Image Compt. and Computer Assist. Interven. (MICCAI) , Shenzhen, China, Oct. 2019, pp. 31–40. [30] Z. Li, I. Y . Chun, and Y . Long, “Image-domain material decomposition using an iterative neural network for dual-energy CT, ” in Pr oc. IEEE Intl. Symp. Biomed. Imag . (ISBI) (to appear) , Io wa City , IA, Apr . 2020. [31] H. Lim, Z. Huang, J. A. Fessler , Y . K. Dewaraja, and I. Y . Chun, “ Application of trained deep BCD-Net to iterativ e lo w-count PET image reconstruction, ” in Proc. IEEE Nuclear Science Symposium and Medical Imaging Conference (NSS-MIC) , Sydney , Australia, No v . 2018, pp. 1–4. [32] M. Elschot, M. G. Lam, M. A. van den Bosch, M. A. Vier gev er , and H. W . de Jong, “Quantitati ve Monte Carlo-based 90Y SPECT reconstruction, ” J. Nucl. Sci. , vol. 54, no. 9, pp. 1557–1563, 2013. [33] A. S. Pasciak, A. C. Bourgeois, J. M. McKinney , T . T . Chang, D. R. Osborne, S. N. Acuf f, and Y . C. Bradley , “Radioembolization and the dynamic role of 90Y PET/CT, ” F r ontiers in oncology , vol. 4, p. 38, 2014. [34] J. Nuyts, D. Beque, P . Dupont, and L. Mortelmans, “ A concav e prior penalizing relativ e differences for maximum-a-posteriori reconstruction in emission tomography , ” IEEE T r ans. Nucl. Sci. , vol. 49, no. 1, pp. 56–60, 2002. [35] L. Pfister and Y . Bresler, “Learning sparsifying filter banks, ” in W avelets and Sparsity XVI , vol. 9597. International Society for Optics and Photonics, 2015, p. 959703. [36] J.-F . Cai, H. Ji, Z. Shen, and G.-B. Y e, “Data-dri ven tight frame con- struction and image denoising, ” Applied and Computational Harmonic Analysis , vol. 37, no. 1, pp. 89–105, 2014. [37] I. Y . Chun and J. A. Fessler , “Conv olutional analysis operator learning: Application to sparse-vie w CT, ” in Pr oc. Asilomar Conf. on Signals, Syst., and Comput. , Pacific Grove, CA, Oct. 2018, pp. 1631–1635. [38] A. R. De Pierro, “ A modified expectation maximization algorithm for penalized likelihood estimation in emission tomography , ” IEEE T r ans. Med. Imag. , vol. 14, no. 1, pp. 132–7, Mar. 1995. [39] W . H. Press, B. P . Flannery , S. A. T eukolsky , and W . T . V etterling, Numerical r ecipes in C . Ne w Y ork: Cambridge Univ . Press, 1988. [40] Z. Zhang, S. Rose, J. Y e, A. E. Perkins, B. Chen, C.-M. Kao, E. Y . Sidky , C.-H. T ung, and X. P an, “Optimization-based image reconstruction from low-count, list-mode T OF-PET data, ” IEEE T ransactions on Biomedical Engineering , vol. 65, no. 4, pp. 936–946, 2018. [41] A. Chambolle and T . Pock, “ An introduction to continuous optimization for imaging, ” Acta Numerica , vol. 25, pp. 161–319, 2016. [42] S. Y . Chun, Y . K. Dewaraja, and J. A. Fessler, “ Alternating direction method of multiplier for tomography with nonlocal regularizers, ” IEEE T rans. Med. Imag. , vol. 33, no. 10, pp. 1960–1968, 2014. [43] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers, ” F ound. & T r ends in Machine Learning , vol. 3, no. 1, pp. 1–122, 2010. [44] W . Segars, G. Sturgeon, S. Mendonca, J. Grimes, and B. M. Tsui, “4D XCA T phantom for multimodality imaging research, ” Medical Physics , vol. 37, no. 9, pp. 4902–4915, 2010. [45] H. Lim, Y . K. Dewaraja, and J. A. Fessler, “ A PET reconstruction formulation that enforces non-negativity in projection space for bias reduction in Y -90 imaging, ” Phys. Med. Biol. , vol. 63, no. 3, p. 035042, Feb . 2018. [46] A. P aszke, S. Gross, S. Chintala, G. Chanan, E. Y ang, Z. DeV ito, Z. Lin, A. Desmaison, L. Antiga, and A. Lerer, “ Automatic differentiation in PyT orch, ” in NIPS-W , 2017. [47] D. P . Kingma and J. Ba, “Adam: A method for stochastic optimization, ” arXiv pr eprint arXiv:1412.6980 , 2014. [48] J. Zbontar, F . Knoll, A. Sriram, M. J. Muckley , M. Bruno, A. Defazio, M. Parente, K. J. Geras, J. Katsnelson, H. Chandarana et al. , “fastMRI: An open dataset and benchmarks for accelerated MRI, ” arXiv pr eprint arXiv:1811.08839 , 2018. [49] E. Kang, J. Min, and J. C. Y e, “ A deep conv olutional neural network us- ing directional w avelets for low-dose X-ray CT reconstruction, ” Medical Physics , vol. 44, no. 10, pp. e360–e375, 2017.