4D X-Ray CT Reconstruction using Multi-Slice Fusion

4D X-Ray CT Reconstruction using Multi-Slice Fusion Soumendu Majee 1 , Thilo Balke 1 , Craig A. J. K emp 2 , Gregery T . Buzzard 3 , and Charles A. Bouman 1 1 School of Electrical and Computer Engineering, Purdue Univ ersity , W est Lafayette, IN, USA 2 Eli Lilly and Company , Indianapolis, IN, USA 3 Department of Mathematics, Purdue Univ ersity , W est Lafayette, IN, USA There is an increasing need to reconstruct objects in four or more dimensions corresponding to space, time and other independent parameters. The best 4D reconstruction algorithms use regularized iterative reconstruction approaches such as model based iterative reconstruction (MBIR), which depends critically on the quality of the prior modeling. Recently , Plug-and-Play methods ha ve been shown to be an effective way to incorporate advanced prior models using state-of-the-art denoising algorithms designed to remove additive white Gaussian noise (A WGN). However , state-of-the-art denoising algorithms such as BM4D and deep con volutional neural networks (CNNs) are primarily available for 2D and sometimes 3D images. In particular , CNNs are difficult and computationally expensive to implement in four or more dimensions, and training may be impossible if there is no associated high-dimensional training data. In this paper , we present Multi-Slice Fusion , a novel algorithm for 4D and higher -dimensional reconstruction, based on the fusion of multiple low-dimensional denoisers. Our approach uses multi-agent consensus equilibrium (MA CE), an extension of Plug-and-Play , as a framework f or integrating the multiple lower -dimensional prior models. W e apply our method to the problem of 4D cone-beam X-ray CT reconstruction for Non Destructive Evaluation (NDE) of moving parts. This is done by solving the MA CE equations using lower -dimensional CNN denoisers implemented in parallel on a heterogeneous cluster . Results on experimental CT data demonstrate that Multi-Slice Fusion can substantially impro ve the quality of reconstructions relative to traditional 4D priors, while also being practical to implement and train. Index T erms —In verse Problems, Plug-and-Play , Consensus Equilibrium I . I N T RO D U C T I O N Improv ements in imaging sensors and computing power hav e made it possible to solve increasingly dif ficult recon- struction problems. In particular , the dimensionality of recon- struction problems has increased from the traditional 2D and 3D problems representing space to more difficult 4D or e ven 5D problems representing space-time and, for example, heart or respiratory phase [1], [2]. These higher-dimensional reconstruction problems pose sur- prisingly dif ficult challenges both computationally and perhaps more importantly , in terms of algorithmic design and training due to the curse of dimensionality [3]. Howe ver , the high dimensionality of the reconstruction also presents important opportunities to improv e reconstruction quality by exploiting regularity . In particular , it has been sho wn that the temporal resolution of 4D CT imaging can be increased by an order of magnitude by exploiting the time space regularity of objects being imaged [2], [4]–[6]. These approaches use model-based iterativ e reconstruction (MBIR) [7], [8] to enforce regularity using simple space-time prior models. Recently , it has been demonstrated that Plug-and-Play (P&P) priors [9]–[12] can dramatically improve reconstruction quality by enabling the use of state-of-the-art denoisers as prior models in MBIR. So P&P has great potential to improv e reconstruction quality in 4D CT imaging problems. Howe ver , state-of-the-art denoisers such as deep Con volutional Neural Networks (CNNs) and BM4D are primarily a vailable for 2D and sometimes 3D images, and it is difficult to extend them to higher dimensions [3], [13], [14]. In particular , extending CNNs to 4D requires v ery computationally , and memory in- tensiv e 4D conv olution applied to 5D feature tensor structures. This problem is further compounded by the lack of GPU accelerated routines for 4D conv olution from major Deep- Learning frameworks such as T ensorflow , Keras, PyT orch. Currently only 1D, 2D, and 3D con volutions are supported with GPU acceleration. Furthermore, 4D CNNs require 4D ground truth data to train the P&P denoisers, which might be difficult or impossible to obtain. Fig. 1. Illustration of our Multi-Slice Fusion approach. W e fuse lower di- mensional denoisers with the measurement model to produce a 4D regularized reconstruction. In this paper, we present a nov el 4D X-ray CT reconstruc- tion algorithm that combines multiple low-dimensional CNN denoisers to implement a highly ef fectiv e 4D prior model. Our approach, which we call Multi-Slice Fusion , integrates the multiple low-dimensional priors using multi-agent consensus equilibrium (MA CE) [15]. MACE is an extension of the P&P framew ork that formulates the inv ersion problem using an equilibrium equation—as opposed to an optimization—and allows for the use of multiple prior models and agents. Figure 1 illustrates the basic concept of our approach. Multi- Slice Fusion integrates together three distinct denoisers each of which is designed to remove noise along lower dimensional slices of the 4D object. When MA CE fuses the denoisers it simultaneously enforces the constraints of each denoising agent, so the reconstructions are constrained to be smooth in all four dimensions. Consequently , Multi-Slice Fusion results in high-quality reconstructions that are practical to train and compute even when the dimensionality of the reconstruction is high. In our implementation, one MACE agent estimates the cone-beam tomographic in version. The remaining 3 agents are CNN denoisers trained to remov e additive white Gaussian noise (A WGN) along 2.5D slices of the reconstruction [3], [16], [17]. The MA CE solution may be computed using a v ariety of algorithms, including v ariants of the Plug-and-Play algorithm based on ADMM or other approaches [9], [10], [18], [19]. W e implement our method for the problem of 4D cone- beam X-ray CT reconstruction [20] of moving parts, and the MA CE solution is computed on a heterogeneous cluster in which dif ferent agent updates are computed on different cluster nodes. In particular, the algorithm runs in parallel on heterogeneous cluster nodes, with the in version portion running on CPU nodes and the prior portion running on GPU nodes. Experimental results on both simulated data and real data obtained from an North Star Imaging X50 X-ray cone-beam scanner indicate that the Multi-Slice Fusion reconstruction using CNN priors can substantially improve reconstructed image quality as compared to MBIR reconstruction using traditional 4D MRF priors. I I . P RO B L E M F O R M U L A T I O N In X-Ray CT imaging, an object is rotated and several 2D projections (radiographs) of the object are measured for dif fer- ent angles. The problem is then to reconstruct the 4D array of X-ray attenuation coef ficients from these measurements, where 3 dimensions correspond to the spatial dimensions and the fourth dimension corresponds to time. Let N t be the number of time-points, and N s be the number of voxels at each time point of the 4D image. For each time- point, n ∈ { 1 , . . . , N t } , define y n ∈ R M n to be the vector of sinogram measurements at time n , and x n ∈ R N s to be the 3D volume of X-ray attenuation coef ficients for that time- point. Let us stack all the measurements to form a vector Y = [ y > 1 , .., y > N t ] > ∈ R M and stack all the 3D v olumes to form a vector X = [ x > 1 , . . . , x > N t ] > ∈ R N , where N = N t N s and M = P N t n =1 M n . The 4D reconstruction problem then becomes the task of recovering the 4D image of attenuation coefficients, X , from the series of sinogram measurements, Y . The image X can be reconstructed as a maximum a poste- riori (MAP) estimate as X ∗ = arg min X { f ( X ) + R ( X ) } , (1) where f ( X ) is the data-fidelity or log-likelihood term and R ( X ) is the 4D regularizer or prior model for the image, X . The data-fidelity term of the 4D image, X , can be written in a separable fashion as f ( X ) = N t X n =1 f n ( x n ) , (2) where f n ( x n ) = 1 2 α k y n − A n x n k 2 Λ n (3) is the data fidelity term for time n , where A n is the system matrix, and Λ n is the weight matrix, and α is a normalizing weight scalar [20]. The values in the diagonal weight matrix, Λ n , account for the non-uniform noise variance for each measurement. If the prior model, R ( X ) , can be expressed analytically like a 4D Markov Random Field (MRF) as in [2], then the expres- sion in Equation (1) can be minimized iteratively . Howe ver , in practice, it can be difficult to represent an advanced prior model in the form of a simple cost function R ( X ) that can be minimized. Consequently , P&P methods hav e been created as a method for representing prior models as denoising operations [9], [10]. More recently , P&P methods have been generalized to the consesus equilibrium framework as a way to integrate multiple criteria [15]. Fig. 2. Architecture of CNN denoiser. Different sizes of input and output necessitate a selection operator for the residual connection. Each rectangular block denotes a tensor, and each ellipsoidal block denotes an operation. Blue ellipses specify the shape of the con volution kernel. I I I . F U S I N G D E N O I S E R S V I A C O N S E N S U S E Q U I L I B R I U M W e use the MA CE framew ork to fuse multiple denoisers to form a prior model for reconstruction. This allo ws us to make a prior model out of denoisers that are not optimization based [15]. Let H k for k ∈ { 1 , . . . , K } be a lower -dimensional denoiser that denoises a full 4D v olume by processing the volume “slice by slice”. Furthermore, each denoiser , H k , is trained to remov e additiv e white Gaussian noise (A WGN) of standard deviation σ from the image. The denoisers serve as prior models ev en though a corresponding regularization term, R ( X ) , does not necessarily exist. The design of these denoisers is described in more detail in section IV. Each of the denoisers is a separate agent in the MA CE framework. The data-fidelity is achieved using a proximal map opera- tion, F , as another agent. The F -operator is defined as the proximal map of f as F ( X ) = arg min Z ∈ R N  f ( Z ) + β σ 2 k X − Z k 2 2  , (4) where σ is the same as above and the unit-less parameter β is used to vary the relati ve weight of the F -operator with with respect to the denoisers and in turn the amount of regularization in the reconstruction can be adjusted. The operator F enforces fidelity to the measurements while each of the denoisers H 1 , . . . , H K enforces regularity of the image in orthogonal image planes. Our method imposes a consensus between the operators F , H 1 , . . . , H K to achiev e a balanced reconstruction that lies at the intersection of the solution space of the measurement model and each of the prior models. In order to formally state the consensus equilibrium criteria and present the algorithm for computing the solution, we will need to introduce some additional notation. First, define the column vector X = [ X > 1 , . . . , X > K +1 ] > denotes a tall vector ∈ R ( K +1) N as the concatination of K + 1 state vectors X k ∈ R N where each state vector represents the input to one of the K + 1 agents. Furthermore, let L ( X ) denote the operator formed by the application of all K + 1 agents on each corresponding state vector . More precisely , we define the agent operator , L , as L ( X ) =    H 1 ( X 1 ) . . . H K ( X K ) F ( X K +1 )    . (5) So L ( X ) simply denotes the parallel application of each agent to its corresponding state vector . Next we define the averaging operator , G , as G ( X ) =    X X . . . X    , (6) where X denotes the weighted av erage of the state vectors for each agent giv en by X = 1 2 X K +1 + 1 2 K K X k =1 X k . (7) Notice that we choose the weighting so that there is equal weighting of the data-fidelity (i.e., agent K + 1 ) and regular - ization (i.e., agents 1 to K ), irrespective of ho w many denoiser we use to form our prior model. Using, the above notation, it can be shown [15], [21] that when the agents are at consensus, the equilibrium equation L ( X ) = G ( X ) (8) holds. The equilibrium condition of equation (8) implies that at equilibrium each agent with different inputs produces the same output. This output is then the consensus solution that represents a balance point between the forces of each agent. From equation (8) it follo ws that the consensus solution is the fixed point of the map T = (2 G − I )(2 L − I ) and can found efficiently by using partial update Mann iterations [15], [21], [22]. Partial update Mann iterations use approximate updates for the proximal map operator , F . W e define ˜ L ( X ; V ) to be the corresponding partial update approximation of L ( X ) ˜ L ( X ; V ) =     H 1 ( X 1 ) . . . H K ( X K ) ˜ F ( X K +1 ; V K +1 )     , (9) where ˜ L ( X ; V ) is the same as L ( X ) but uses only an approximation proximal mapping formed by a single update of a iterativ e optimization method ˜ F ( X K +1 ; V K +1 ) initialized by V K +1 . That is, the iterati ve optimization in the F -operator in equation (4) is initialized by V K +1 and terminated early . The the partial update Mann iterations for fusing denoisers can be simplified to form Algorithm 1. The weighing param- eter ρ ∈ (0 , 1) is used to control speed of con ver gence. The initial reconstruction used should be close to the final solution for faster conv ergence. Output of a fast conv entional algorithm can be used as the initial reconstruction. W e continue the partial update Mann iterations until the differ ences between state v ectors X k become smaller than a fixed threshold. W e use three passes of Iterativ e Coordinate Descent (ICD) for com- puting the approximate proximal mapping, ˜ F ( X K +1 ; V K +1 ) . Algorithm 1: Fusing denoisers using MACE Input: Initial Reconstruction: X (0) ∈ R N Output: Final Reconstruction: X ∗ 1 X ← W ← X (0) . . . X (0) ! 2 while not conver ged do 3 X ← ˜ L ( W ; X ) 4 Z ← G (2 X − W ) 5 W ← W + 2 ρ ( Z − X ) 6 X ∗ ← X 1 Fig. 3. Illustration of our training volume. W e extract 3D patches from the volume to train our network. I V . C N N S PAC E - T I M E D E N O I S E R Since we iterativ ely apply multiple denoisers to multiple 4D volumes to compute the reconstruction, it is essential for each denoiser to be fast. W e use Deep Con volutional Neural Networks since they provide state-of-the-art quality in denoising while at the same time being fast due to optimized GPU libraries. The architecture of the CNN denoiser is shown in Figure 2. Each of the CNN denoisers used in the Multi-Slice Fusion is designed to perform 2.5D denoising of a 3D volume using a 2D moving windo w . By 2.5D we mean that as an input the neural network takes a thin 3D patch consisting of five consecutiv e 2D slices while the output is the denoised center slice. The other slices are being denoised by shifting the 5- slice moving window . The 2D slices are processed using 2D con volution while the adjacent slices are treated as separate input channels. The neural network architecture is modeled after the residual network architecture presented in [17] but extended so that it works on 3D rather than 2D images. In [3], [16], it has been shown that this type of 2.5D processing is a computationally efficient way of performing ef fectiv e 3D denoising with CNNs. W e note that it is similar to the method presented in [17] where the 5 adjacent slices are analogous to the 3 color channels used in denoising of color images. Our 3D CNN denoiser can theoretically be trained with space-time data to get dif ferent time and space regularization. Howe ver , in the absence of space-time training data, we trained the network on 3D-space volumes. For training, we used a CT reconstruction of a bottle and screw cap made from different plastics that has low noise. A slice from the training volume is illustrated in Figure 3. While our training and testing volumes are different, all of them are representati ve of a variety of NDE problems in which objects to be imaged are constructed from a relativ ely small number of distinct materials. W e extracted 63 , 000 partially overlapping 3D patches of size 40 × 40 × 5 from our training volume of size 534 × 546 × 100 . As a data augmentation technique random rotation, flipping, and mirroring were applied during training. Pseudo- random additiv e white Gaussian noise (A WGN) was applied to the patches and the network was trained to remove the noise. The A WGN noise used does not represent CT measurement noise b ut is a result of using the MA CE framew ork. More specifically , the A WGN noise that the denoisers removes in MA CE is an ’artificial’ noise induced from the quadratic norm used in equation (4). V . 4 D P R I O R A S M U LT I - S L I C E F U S I O N In this section we describe ho w we fuse the space-time denoisers described in section IV using the approach in section III to form a 4D reconstruction. Although, the approach described in section III can fuse arbitrary number of denoisers, we fuse three 3D CNN denois- ers to form our 4D prior . W e refer to denoisers H 1 , H 2 , H 3 in equation (5) as H xy ,t , H y z ,t , and H z x,t , respectiv ely . The denoiser H xy ,t is a CNN space-time denoiser that does con volution in the xy plane and uses fiv e input channels to Scanner Model North Star Imaging (Roger s, MN) X50 Source-Detector Distance 839 mm Magnification 5.57 Number of V iews per Rotation 150 Cropped Detector Array 731 × 91 , (0 . 25 mm) 2 V ox el Size (0 . 0456 mm) 3 Reconstruction Size ( x, y , z , t ) 731 × 731 × 91 × 16 T ABLE I E X PE R I M EN TA L S ET U P F O R R E A L X - R AY D AT A input slices from neighboring time-points. The denoisers H y z ,t and H z x,t are analogous to H xy ,t but are applied along the yz and zx plane, respectiv ely . Each of the denoisers, despite being 3D denoisers, denoise the 4D volume by processing “slice by slice”. Each denoiser incorporates information along a spatial slice and uses a moving temporal windo w of 5 slices to denoise a single slice that is in the center . Howe ver , our Multi-Slice Fusion fuses all the denoisers H xy ,t , H y z ,t , and H z x,t to incorporate information along all four dimensions. V I . E X P E R I M E N T A L R E S U LT S T o demonstrate the improved reconstruction quality of our method, we present results on both real and simulated 4D X-ray CT data for NDE applications. In both cases we compare Multi-Slice Fusion reconstruction with 4D Markov Random Field based MBIR reconstruction (MBIR+4D-MRF) and conv entional 3D Filtered Back Projection reconstruction (FBP). The MBIR+4D-MRF used was a variant of the QG- GMRF model described in [2] with 26 spatial and 2 temporal neighbors. W e also compare with reconstructions using single CNNs to illustrate the effect of model fusion. W e refer to reconstructions using single CNN denoisers H xy ,t , H y z ,t , and H z x,t as MBIR+ H xy ,t , MBIR+ H y z ,t , and MBIR+ H z x,t , respectiv ely . In all our experiments the cone-beam in version within Multi-Slice Fusion was executed on two CPU cluster nodes, each with 20 Kaby Lake CPU cores. It was parallelized over time-points trivially and within time points using multiple threads similar to [20]. W e perform three equits per each Mann iteration for the cone-beam inv ersion. The CNN denoisers ran in parallel with the cone-beam in version in a separate cluster with a Nvidia T esla P100 GPU. A. Real Dataset In this section we present results on real data to ev aluate our method. The data is from a dynamic cone-beam X-ray scan of a glass vial, with elastomeric stopper and aluminum crimp-seal, using a North Star Imaging (Rogers, MN) X50 X- ray CT system. The vial is undergoing dynamic compression during the scan, to capture the mechanical response of the components. Of particular interest is the moment when the aluminum seal is no longer in contact with the underside of the glass neck finish. This indicates the moment when the force applied exceeds that ex erted by the rubber on the glass; this is kno wn as the “residual seal force” [23]. During the scan the vial was held in place by fixtures going in and out of the Field of V ie w thus resulting in artifacts. All our 4D results incorporate preprocessing to correct these artifacts. The experimental setup is summarized in T able I. (a) FBP (3D) (b) MIBIR+4D-MRF (c) Multi-Slice Fusion (d) Plot of cross-section Fig. 4. Plot of cross-section through the vial at a time when the aluminum and glass have physically separated. Multi-Slice Fusion is able to resolve the junction between materials better while simultaneously producing a smoother reconstruction within materials compared to MBIR+4D-MRF and FBP . Fig. 5. Illustration of temporal resolution. Multi-Slice Fusion results in improved space-time resolution of the separation of aluminum and glass. Figure 6 compares Multi-Slice Fusion with sev eral other methods. Each image is a slice through the reconstructed vial Source-Detector Distance 839 mm Magnification 5.57 Number of V iews per Rotation 75 Cropped Detector Array 240 × 28 , (0 . 25 mm) 2 V ox el Size (0 . 0456 mm) 3 Reconstruction Size ( x, y , z , t ) 240 × 240 × 28 × 8 T ABLE II E X PE R I M EN TA L S ET U P F O R S I M UL ATE D X - R A Y D AT A for one time point along the spatial xy plane. Both FBP and MBIR+4D-MRF suf fer from ob vious windmill artif acts, higher noise and hav e blurred contours. In contrast to that, the Multi- Slice Fusion reconstruction has smooth and uniform textures while preserving edge definition. Figure 6 also illustrates the effect of model fusion by comparing Multi-Slice Fusion with MBIR+ H xy ,t , MBIR+ H y z ,t , and MBIR+ H z x,t . MBIR+ H y z ,t and MBIR+ H z x,t suffer from horizontal and vertical streaks respectiv ely since the denoisers were applied in those planes. MBIR+ H xy ,t does not suffer from streaks in the figure since we are viewing a slice along the xy plane, but it suffers from other artifacts. MBIR+ H xy ,t cannot reconstruct the outer ring since the slice displayed is at the edge of the aluminum seal and the xy plane does not contain sufficient information. In contrast, Multi-Slice Fusion can resolve the edges of the rings better than either of MBIR+ H xy ,t , MBIR+ H y z ,t , and MBIR+ H z x,t since it uses information from all the spatial coordinates. Next, we plot a cross-section through the object for Multi- Slice Fusion, MBIR+4D-MRF and FBP in Figure 4. For this, we choose a time-point where we know the aluminum and glass ha ve separated spatially . The Multi-Slice Fusion reconstruction has a steeper and more defined slope in the junction of aluminum and glass. This supports that Multi-Slice Fusion is able to preserve fine details in spite of producing a smooth regularized image. Finally in Figure 5 we plot a cross-section through the object with time to show the improv ed space-time resolution of our method. W e do this for both the 4D methods: MBIR+4D-MRF and Multi-Slice Fusion. Multi-Slice Fusion results in improved space-time resolution of the separation of aluminum and glass. B. Simulated Data In this section we present results on simulated data to ev aluate our method. W e take a lo w-noise CT reconstruction of a bottle and screw cap and denoise it further using BM4D [14] to generate a clean 3D phantom. W e then move the 3D volume vertically for each time point to generate a 4D phantom. W e forward project the phantom to generate sinogram data and use that to reconstruct the object. The simulated experimental setup is summarized in T able II. Figure 8 compares Multi-Slice Fusion with sev eral other methods. Each image is a slice through the reconstructed object for one time point along the spatial xy plane. Both FBP and MBIR+4D-MRF suffer from high noise and jagged edges and fail to recov er the small hole in the bottom of the image. MBIR+ H y z ,t and MBIR+ H z x,t suffer from horizontal and vertical streaks respectively since the denoisers were applied (a) FBP (3D) (b) MBIR+4D-MRF (c) Multi-Slice Fusion (d) MBIR+ H xy,t (Missing Feature) (e) MBIR+ H yz ,t (Horizontal Streaks) (f) MBIR+ H zx,t (V ertical Streaks) Fig. 6. Comparison of different methods. Each image is a slice through the reconstructed vial for one time point along the spatial xy plane. Both FBP and MBIR+4D-MRF suffer from obvious windmill artifacts, higher noise and have blurred contours. In contrast to that, the Multi-Slice Fusion reconstruction has smooth and uniform textures while preserving edge definition. MBIR+ H yz ,t and MBIR+ H zx,t suffer from horizontal and vertical streaks. MBIR+ H xy,t cannot reconstruct the outer ring since the slice displayed is at the edge of the aluminum seal and the xy plane does not contain sufficient information. Multi-Slice Fusion can resolve the edges of the rings better than either of MBIR+ H xy,t , MBIR+ H yz ,t , and MBIR+ H zx,t since it has information from all the spatial coordinates. in those planes. MBIR+ H xy ,t does not suffer from streaks in the figure since we are viewing a slice along the xy plane, but it suf fers from other artifacts. MBIR+ H xy ,t cannot reconstruct the small hole in the bottom of the image since the xy plane does not contain sufficient information. Next we plot a cross-section through the object for Multi- Slice Fusion, MBIR+4D-MRF , FBP , and the phantom Fig- ure 7. Multi-Slice Fusion results in the most accurate recon- struction of the gap between materials. (a) Phantom (b) FBP (3D) (c) MIBIR+4D- MRF (d) Multi-Slice Fusion (e) Plot of cross-section Fig. 7. Plot of cross-section through the phantom and reconstructions from simulated data. Multi-Slice Fusion results in the most accurate reconstruction of the gap between materials. Finally we report the Peak Signal to Noise Ratio (PSNR) and Structural Similarity Index (SSIM) with respect to the phantom for each method in T able III to objectiv ely measure image quality . W e define the PSNR for a giv en reconstruction X with a phantom X 0 as P S N R ( X ) = 20 log 10  Rang e ( X 0 ) RM S E ( X , X 0 )  , (10) where range is computed from the 0 . 1 and 99 . 9 percentiles of the phantom. As can be seen from T able III, Multi-Slice Fusion results in the highest PSNR and SSIM scores. V I I . C O N C L U S I O N In this paper , we proposed a no vel 4D X-Ray CT reconstruc- tion algorithm, Multi-Slice Fusion, that combines multiple low dimensional denoisers to form a 4D prior . Our method allows the formation of an adv anced 4D prior using state-of- the-art CNN denoisers without needing to train on 4D data. Furthermore, it allo ws for multiple le vels of parallelism; thus enabling reconstruction of large volumes in a reasonable time. Method PSNR(dB) SSIM FBP 19.690 0.609 MBIR+4D-MRF 25.837 0.787 Multi-Slice Fusion 29.071 0.943 MBIR+ H xy,t 29.026 0.922 MBIR+ H yz ,t 28.040 0.932 MBIR+ H zx,t 28.312 0.926 T ABLE III Q UA N TI TA T I V E E V A L UA T I O N . M U LTI - S L IC E F U S I ON H A S T H E H IG H E S T P S NR A N D S S IM M E T RI C A M O N G A LL T H E M E TH O D S . Although we focused on 4D X-Ray CT reconstruction for Non Destructi ve Evaluation (NDE) applications, our method can be used for any reconstruction problem in volving multiple dimensions. A C K N O W L E D G M E N T The authors would lik e to acknowledge support from Eli Lilly and Company under research project funding agreement 17099289. Charles A. Bouman and Gregery T . Buzzard were supported in part by NSF grant CCF-1763896. W e also thank M. Cory V ictor and Dr . Coralie Richard from Eli Lilly for their assistance and guidance in setting up the residual seal force test experiment. R E F E R E N C E S [1] C. Huang, J. L. Ackerman, Y . Petibon, T . J. Brady , G. El Fakhri, and J. Ouyang, “MR-based motion correction for PET imaging using wired active MR microcoils in simultaneous PET -MR: Phantom study , ” Medical physics , vol. 41, no. 4, 2014. [2] K. A. Mohan, S. V enkatakrishnan, J. W . Gibbs, E. B. Gulsoy , X. Xiao, M. De Graef, P . W . V oorhees, and C. A. Bouman, “TIMBIR: A method for time-space reconstruction from interlaced views. ” IEEE T rans. Computational Imaging , vol. 1, no. 2, pp. 96–111, 2015. [3] A. Ziabari, D. H. Y e, K. D. Sauer, J. Thibault, and C. A. Bouman, “2.5D deep learning for CT image reconstruction using a multi-GPU imple- mentation, ” in Signals, Systems, and Computers, 2018 52nd Asilomar Confer ence on . IEEE, 2018. [4] J. Gibbs, K. A. Mohan, E. Gulsoy , A. Shahani, X. Xiao, C. Bouman, M. De Graef, and P . V oorhees, “The three-dimensional morphology of growing dendrites, ” Scientific r eports , vol. 5, p. 11824, 2015. [5] K. A. Mohan, S. V enkatakrishnan, J. W . Gibbs, E. B. Gulsoy , X. Xiao, M. De Graef, P . W . V oorhees, and C. A. Bouman, “4D model-based iterativ e reconstruction from interlaced views. ” in ICASSP , 2015, pp. 783–787. [6] X. W ang, K. A. Mohan, S. J. Kisner , C. Bouman, and S. Midkiff, “Fast vox el line update for time-space image reconstruction, ” in 2016 IEEE International Conference on Acoustics, Speech and Signal Pr ocessing (ICASSP) . IEEE, 2016, pp. 1209–1213. [7] S. J. Kisner, E. Haneda, C. A. Bouman, S. Skatter, M. K ourinny , and S. Bedford, “Model-based CT reconstruction from sparse views, ” in Second International Confer ence on Image F ormation in X-Ray Computed T omography , 2012, pp. 444–447. [8] K. Sauer and C. Bouman, “ A local update strategy for iterativ e recon- struction from projections, ” IEEE T ransactions on Signal Pr ocessing , vol. 41, no. 2, pp. 534–548, 1993. [9] S. Sreehari, S. V . V enkatakrishnan, B. W ohlberg, G. T . Buzzard, L. F . Drummy , J. P . Simmons, and C. A. Bouman, “Plug-and-play priors for bright field electron tomography and sparse interpolation, ” IEEE T ransactions on Computational Imaging , vol. 2, no. 4, pp. 408–423, 2016. [10] S. V . V enkatakrishnan, C. A. Bouman, and B. W ohlberg, “Plug-and-play priors for model based reconstruction, ” in Global Confer ence on Signal and Information Processing (GlobalSIP), 2013 IEEE . IEEE, 2013, pp. 945–948. [11] Y . Sun, B. W ohlberg, and U. S. Kamilov , “ An online plug-and- play algorithm for regularized image reconstruction, ” arXiv preprint arXiv:1809.04693 , 2018. (a) Phantom (b) FBP (3D) (c) MBIR+4D-MRF (d) Multi-Slice Fu- sion (e) MBIR+ H xy,t (f) MBIR+ H yz ,t (g) MBIR+ H zx,t Fig. 8. Comparison of different methods. Each image is a slice through the reconstructed object for one time point along the spatial xy plane. Both FBP and MBIR+4D-MRF suffer from high noise and jagged edges and fail to recover the small hole in the bottom of the image. MBIR+ H yz ,t and MBIR+ H zx,t suffer from horizontal and vertical streaks respectively since the denoisers were applied in those planes. MBIR+ H xy,t cannot reconstruct the small hole in the bottom of the image since the xy plane does not contain sufficient information. [12] U. S. Kamilov , H. Mansour, and B. W ohlberg, “ A plug-and-play priors approach for solving nonlinear imaging in verse problems, ” IEEE Signal Pr ocessing Letters , vol. 24, no. 12, pp. 1872–1876, 2017. [13] K. Dabov , A. Foi, V . Katko vnik, and K. Egiazarian, “Image denoising by sparse 3D transform-domain collaborativ e filtering, ” IEEE Tr ansactions on image pr ocessing , vol. 16, no. 8, pp. 2080–2095, 2007. [14] M. Maggioni, G. Boracchi, A. Foi, and K. Egiazarian, “V ideo denoising using separable 4D nonlocal spatiotemporal transforms, ” in Image Pr o- cessing: Algorithms and Systems IX , vol. 7870. International Society for Optics and Photonics, 2011, p. 787003. [15] G. T . Buzzard, S. H. Chan, S. Sreehari, and C. A. Bouman, “Plug- and-play unplugged: Optimization-free reconstruction using consensus equilibrium, ” SIAM J ournal on Imaging Sciences , vol. 11, no. 3, pp. 2001–2020, 2018. [16] D. Jiang, W . Dou, L. V osters, X. Xu, Y . Sun, and T . T an, “Denoising of 3D magnetic resonance images with multi-channel residual learning of con volutional neural network, ” Japanese journal of radiology , vol. 36, no. 9, pp. 566–574, 2018. [17] K. Zhang, W . Zuo, Y . Chen, D. Meng, and L. Zhang, “Beyond a gaussian denoiser: Residual learning of deep CNN for image denoising, ” IEEE T ransactions on Image Processing , vol. 26, no. 7, pp. 3142–3155, 2017. [18] Y . Sun, B. W ohlberg, and U. S. Kamilov , “Plug-in stochastic gradient method, ” arXiv preprint , 2018. [19] Y . Sun, S. Xu, Y . Li, L. T ian, B. W ohlberg, and U. S. Kamilov , “Regu- larized fourier ptychography using an online plug-and-play algorithm, ” arXiv pr eprint arXiv:1811.00120 , 2018. [20] T . Balke, S. Majee, G. T . Buzzard, S. Poveromo, P . Howard, M. A. Groeber , J. McClure, and C. A. Bouman, “Separable models for cone- beam MBIR reconstruction, ” Electr onic Imaging , vol. 2018, no. 15, pp. 181–1, 2018. [21] V . Sridhar, G. T . Buzzard, and C. A. Bouman, “Distributed framework for fast iterative CT reconstruction from view-subsets, ” Electronic Imag- ing , vol. 2018, no. 15, pp. 102–1, 2018. [22] V . Sridhar, X. W ang, G. T . Buzzard, and C. A. Bouman, “Distributed memory framework for fast iterativ e CT reconstruction from view- subsets using multi-agent consensus equilibrium, ” Manuscript in pr epa- ration for IEEE T ransactions on Computational Imaging , 2018. [23] R. Mathaes, H.-C. Mahler, J.-P . Buettiker , H. Roehl, P . Lam, H. Brown, J. Luemkemann, M. Adler , J. Huwyler, A. Streubel et al. , “The phar- maceutical vial capping process: Container closure systems, capping equipment, regulatory framework, and seal quality tests, ” Eur opean Journal of Pharmaceutics and Biopharmaceutics , vol. 99, pp. 54–64, 2016.