Cardiorespiratory motion-compensated micro-CT image reconstruction using an artifact model-based motion estimation Marcus Brehm Varian Medical System Imaging Laboratory, Täfernstrasse 7, Baden–Dättwil 5405, Switzerland and German Cancer Research Center (DKFZ), Im Neuenheimer Feld 280, Heidelberg 69120, Germany
Stefan Sawall, Joscha Maier, Sebastian Sauppe, and Marc Kachelrießa) German Cancer Research Center (DKFZ), Im Neuenheimer Feld 280, Heidelberg 69120, Germany
(Received 3 November 2014; revised 28 January 2015; accepted for publication 26 February 2015; published 26 March 2015) Purpose: Cardiac in vivo micro-CT imaging of small animals typically requires double gating due to long scan times and high respiratory rates. The simultaneous respiratory and cardiac gating can either be done prospectively or retrospectively. In any case, for true 5D imaging, i.e., three spatial dimensions plus one respiratory-temporal dimension plus one cardiac temporal dimension, the amount of information corresponding to a given respiratory and cardiac phase is orders of magnitude lower than the total amount of information acquired. Achieving similar image quality for 5D than for usual 3D investigations would require increasing the amount of data and thus the applied dose to the animal. Therefore, the goal is phase-correlated imaging with high image quality but without increasing the dose level. Methods: To achieve this, the authors propose a new image reconstruction algorithm that makes use of all available projection data, also of that corresponding to other motion windows. In particular, the authors apply a motion-compensated image reconstruction approach that sequentially compensates for respiratory and cardiac motion to decrease the impact of sparsification. In that process, all projection data are used no matter which motion phase they were acquired in. Respiratory and cardiac motion are compensated for by using motion vector fields. These motion vector fields are estimated from initial phase-correlated reconstructions based on a deformable registration approach. To decrease the sensitivity of the registration to sparse-view artifacts, an artifact model-based approach is used including a cyclic consistent nonrigid registration algorithm. Results: The preliminary results indicate that the authors’ approach removes the sparse-view artifacts of conventional phase-correlated reconstructions while maintaining temporal resolution. In addition, it achieves noise levels and spatial resolution comparable to that of nongated reconstructions due to the improved dose usage. By using the proposed motion estimation, no sensitivity to streaking artifacts has been observed. Conclusions: Using sequential double gating combined with artifact model-based motion estimation allows to accurately estimate respiratory and cardiac motion from highly undersampled data. No sensitivity to streaking artifacts introduced by sparse angular sampling has been observed for the investigated dose levels. The motion-compensated image reconstruction was able to correct for both, respiratory and cardiac motion, by applying the estimated motion vector fields. The administered dose per animal can thus be reduced for 5D imaging allowing for longitudinal studies at the highest image quality. C 2015 American Association of Physicists in Medicine. [http://dx.doi.org/10.1118/1.4916083] Key words: micro-CT, phase correlation, motion compensation
1. INTRODUCTION In vivo cone-beam micro-CT imaging of small animals is an appreciated technology in preclinical research due to its ability to provide 3D information of living objects with very high spatial resolution. An active field of research is cardiac imaging, as it has already been for clinical CT the last decade. However, small animals have higher respiratory rates and faster heart beats than humans, e.g., a mouse breathes between 90 and 200 times/min and a mouse heart beats between 300 and 600 times/min. For this very reason, contrast-enhanced microCT imaging of animal hearts is a particular challenge. 1948
Med. Phys. 42 (4), April 2015
In case of cardiac imaging, also respiratory motion has to be considered due to breathing-induced translation of the heart. One or more temporal windows of interest are described by the centers of the respiratory phase r and the cardiac phase c and by the respective widths ∆r and ∆c. The required double gating, a simultaneous respiratory and cardiac gating, can either be done prospectively1 or retrospectively.2 In any case, retrieving true 5D information, i.e., three spatial dimensions plus one respiratory-temporal dimension plus one cardiac temporal dimension, means that just a fraction (∆r = 10% and ∆c = 20% ⇒ 2% total) of the total amount of projection data acquired corresponds to a given respiratory and cardiac phase.
0094-2405/2015/42(4)/1948/11/$30.00
© 2015 Am. Assoc. Phys. Med.
1948
1949
Brehm et al.: 5D motion-compensated micro CT imaging
The consequences are severe sparse-view artifacts when using conventional reconstruction algorithms (see Fig. 1). The possibilities on the acquisition side are quite low to overcome the problem of limited projection data in each motion phase. In particular, for long term studies, dose issues become important in small animal imaging. For this reason, the limitation in amount of projection data has to be considered by dedicated reconstruction algorithms. The standard reconstruction approach for temporal-correlated volumes is to sort the data into phase bins and reconstruct each bin separately. However, the resulting volumes suffer from streak artifacts induced by the increased angular spacing of projections used for reconstruction. This problem is also known from C-arm CBCT systems in interventional radiology3 and onboard CBCT imaging in radiation therapy.4 One way to deal with the problem of missing projection data in each motion phase is to utilize data from other phases by explicitly incorporating the temporal dimension into the reconstruction process.5–7 For cardiac micro-CT imaging, a variant of the McKinnon–Bates (MKB) image reconstruction algorithm in combination with an edge-preserving noisesuppression filter has been proposed.8 Another way for regularization is the usage of iterative reconstruction algorithms that incorporate some a priori knowledge into the reconstruction process.9,10 This approach has been applied to micro-CT data without and with incorporating temporal dimension.11,12 Both non-motion-compensated approaches, filter-based 5D noise suppression and iterative reconstruction with 5D regularization, considerably reduce the image noise and artifact levels. Furthermore, those are straightforward and need only minor changes to existing algorithms. But full dose usage, i.e., contribution of each projection to every phase volume, is not reached by any of those two approaches. The reason for that is that the prior knowledge imposed by a motion compensation (MoCo) approach is much stronger: There motion is modeled and voxels may change their position.With the more conventional approaches, this knowledge is neither explicitly nor implicitly used. In this work, we want to extend the idea of incorporating cardiac- and respiratory-temporal dimension one step further and apply all projection data for each phase reconstruction by motion compensation.13–16 Here, all acquired projections contribute to each volume, independent of the desired cardiac and respiratory phase, without a loss in temporal resolution
1949
due to the incorporation of motion vector fields (MVFs). In contrast to Ref. 16, we do not require an additional prospectivegated scan for motion estimation which is time-consuming and increases dose. Instead, we apply an artifact-insensitive estimation method that is made robust by incorporating two main ideas: The consideration of the impact of image artifacts by an artifact model that we successfully applied for onboard CBCT imaging in radiation therapy17 and the consideration of the presence of two separate motions influencing each other by a sequential double gating.
2. MATERIALS AND METHODS 2.A. MoCo based on phase-correlated Feldkamp
We use the well-known Feldkamp–Davis–Kress (FDK) filtered backprojection18 as our standard reconstruction algorithm. The FDK algorithm is denoted by f Std = X−1 p, where X is the x-ray transform operator, X−1 is the reconstruction operator, p are the projection data, and f Std the reconstructed image. For phase-correlated Feldkamp (PCF) reconstruction, the projection images are retrospectively sorted into several phase bins. Subsequently, each phase bin is reconstructed separately using the FDK algorithm. For an arbitrary phase −1 bin n, the corresponding operator is denoted as XPCF(n) and the respective PCF image f PCF(n) is given by −1 f PCF(n) = XPCF(n) p.
However, if just a limited number of projections are acquired, the binning induces a large angular spacing of projections used for a single reconstruction. This results in viewaliasing artifacts that deteriorate image quality of these phasecorrelated images. Using a motion compensation approach, projection data from other phases are additionally considered for reconstruction to improve the angular sampling rate. An application of all projection data without considering motion, as it is the case for a standard FDK reconstruction, obviously causes inconsistencies in reconstruction data and ends up in motion blurring. But if motion is modeled accurately, a MoCo volume can be obtained without motion blurring as we have shown previously for respiratory motion in radiation therapy onboard CBCT scans.17,19 The compensation is done by warping the PCF reconstructions according to the motion model before superimposition. For arbitrary phase bins n and m, the MVF describing the motion from phase m to phase n is denoted as Tnm and the respective MoCo image f MoCo(n) is given by ( ) −1 f MoCo(n) = XPCF(m) p ◦ Tnm . m
2.B. Sequential motion compensation for respiratory and cardiac motion
F. 1. Ungated 3D reconstruction (left) using all measured data and doublegated 5D reconstruction (right) using only 2% of the data, corresponding to a given respiratory and cardiac motion phase (C = 200 HU, W = 1200 HU). Medical Physics, Vol. 42, No. 4, April 2015
For cardiac imaging, two different motion patterns are superimposed and have to be taken into account. Besides the cardiac motion, also the respiration-induced displacement of the heart is of importance. We consider both motion components
1950
Brehm et al.: 5D motion-compensated micro CT imaging
separately within a sequential motion compensation approach. In a first step, only respiratory motion is considered while cardiac motion is completely ignored, i.e., the respiratory phasecorrelated images (see Fig. 2), given by gating of respiratory motion only, are used for the first step’s motion estimation and compensation. In consequence, for each pair of respiratory phase bins (r n ,r m ), the corresponding motion vector field Trr nm is given such that ( ) −1 f MoCo(r n ) = p ◦ Trr nm . XPCF(r ) m rm
In a second step, cardiac motion is considered while respiratory motion is already compensated for a given respiratory phase, i.e., the phase-correlated images used for the second step’s motion estimation and compensation are given by gating of cardiac motion and by compensating of respiratory motion according to the motion model estimated in step one (see Fig. 2). For each pair of cardiac phase bins (cn ,cm ), the corresponding motion vector field Tcc nm is derived from those phasecorrelated images. The double motion compensated image f MoCo(r n,c n ) is then given by f MoCo(r n,c n ) =
( ) −1 * XPCF(r p ◦ Trr nm + ◦ Tcc nm . m,c m ) cm , rm -
1950
2.C. Artifact model-based motion estimation
As no other data are available, motion estimation is performed on the PCF reconstructions. Unfortunately, the presence of view-aliasing artifacts within the PCF images and their change in appearance for different phases influences, the registration process such that conventional approaches fail. To decrease the sensitivity on these artifacts, we are incorporating knowledge from motion patterns and an artifact model into the estimation process, an approach that has been successfully applied to respiratory motion in radiation therapy onboard CBCT scans.17 The applied artifact model contains information about location and intensities, e.g., of view-aliasing artifacts. While containing equivalent artifacts as the PCF reconstruction, the model is designed to be free of subject motion. A separation of estimated motion induced by subject motion and induced by artifacts is thus possible later on. In particular, the model is a second series of 3D images, called 4D artifact images in the following. Those 4D artifact images are generated as follows. In order to have a subject-specific model, we use a prior image similar to that used in Ref. 20. The prior image is generated by an FDK reconstruction of all measured data followed by a segmentation. By applying segmentation, all motion artifacts like blurring and streaking can be removed from the FDK
F. 2. Four different phase-correlated reconstructions of a selected mouse dataset: The double-gated reconstructions (left column) suffer from severe sparse-view artifacts. In contrast, the images used for the sequential gating and motion compensation approach, the respiratory-gated images (second column left) and the respiratory-compensated cardiac-gated images (second column right), include less artifacts. This allows for a more robust motion estimation and thus for motion-compensated images (right column) of high image quality (C = 200 HU, W = 1200 HU). Medical Physics, Vol. 42, No. 4, April 2015
1951
Brehm et al.: 5D motion-compensated micro CT imaging
1951
F. 3. Prior image derived from the FDK image by simple multilevel thresholding using automatically chosen values via cluster analysis (C = 200 HU, W = 1200 HU).
image. Thus, motion cannot be a longer part of a forwardprojected dataset in order to obtain a motion-free model.21 The segmentation consists of an identification for regions like air, lung tissue, fat tissue, water tissue, and bone by simple
thresholding with automatically chosen values via cluster analysis.22 By this classification, the cardiovascular system which is enhanced by the use of a contrast agent will be assigned to the bone region. Subsequently, the segmented regions are set
F. 4. Illustration of the proposed motion estimation method using an artifact model: For both, measured projection data as well as simulated projection data obtained by forward projection of a prior image, a PCF reconstruction is performed taking respiratory motion into account only. Subsequently, MVFs are estimated for both image series using deformable registration and then combined to separate respiratory motion from artifact motion. The procedure is repeated for cardiac motion while respiratory motion is already compensated. Medical Physics, Vol. 42, No. 4, April 2015
1952
Brehm et al.: 5D motion-compensated micro CT imaging
to the mean value of the particular region except for voxels of the bone region that retain their CT-values, as they vary too much to properly model them with just one value (see Fig. 3). In a second step, the measurement and reconstruction process is simulated. A forward projection is performed from the resulting prior image (see Fig. 3) in the same geometry the measured projection data were acquired. These simulated projection data undergo the same phase-correlated image reconstruction using the same respiratory signal as the measured projection data. The resulting 4D CBCT images model the desired artifacts while being free of subject motion. Two independent registrations are performed for both image series, based on measured and on simulated projection data, such that two sets of MVFs are obtained (see Fig. 4). While the first MVF contains components induced by subject motion and induced by image artifacts, the second one contains components induced by image artifacts only. In a final step, the second MVF is subtracted from the first one in order to reduce the artifact-induced MVF components of that MVF. Thus, a MVF is obtained that consists mainly of subject motion. The registration method used for estimating the MVFs of an image series is a deformable registration algorithm in combination with a suited strategy we used previously.19 Basis of this motion estimation is the demons algorithm,23–25 a nonrigid intensity-based registration algorithm, and the strategy to just estimate motion between adjacent phases instead of doing so for each pair of motion phases. Here, MVFs of nonadjacent phases are given by concatenation of MVFs from adjacent phases. This strategy allows for a cyclic consistent registration by incorporating temporal constraints like the quasiperiodic
1952
breathing or cardiac motion pattern. Thus, the sensitivity of the registration process to image artifacts is reduced. For the evaluation, we perform motion-compensated image reconstructions without and with our approach. Within a standard motion compensation, denoted by sMoCo, each required MVF Tnm is estimated independently by a standard deformable intensity-based 3D-3D image registration method matching the corresponding two PCF reconstructions f PCF(m) and f PCF(n). In our case, we used the demons algorithm.23–25 Here, the MVFs originate from the PCF reconstructions only and are not corrected using the artifact model. In accordance with Fig. 4, artifact model-based motion compensation, denoted by acMoCo, is performed applying MVFs that are obtained by cyclic registration and correction based on the artifact model. 2.D. Low-dose phase-correlated (LDPC) reconstruction
For evaluation purposes, we additionally perform image reconstructions using another reconstruction paradigm allowing for the reconstruction of double phase-correlated images. Here, we apply the LDPC algorithm proposed by Sawall et al. in Ref. 8 that incorporates the temporal dimension into the reconstruction process and consists of two components. The first component of the LDPC algorithm is the MKB reconstruction method5,6 that was introduced to reduce streak artifacts in the PCF reconstructions of the human heart. The idea of the MKB method is to add an image f corr to the PCF reconstruction that corrects for the streak artifacts. To obtain
F. 5. Three orthogonal slices are shown for images originating from a standard (FDK, left), a conventional phase-correlated (PCF, middle), and a motion-compensated reconstruction (acMoCo, right) at end-expiration (r = 0%) and end-systole (c = 60%) (C = 200 HU, W = 1200 HU). Medical Physics, Vol. 42, No. 4, April 2015
1953
Brehm et al.: 5D motion-compensated micro CT imaging
1953
the correction image, a forward projection of a standard reconstruction f Std is performed. A phase-correlated reconstruction of that forward projected projection data yields the standard image with similar streak artifacts as the PCF reconstruction. The final correction image f corr that exhibits only the streak artifacts can be received by subtracting this streaky image from f Std, i.e., for an arbitrary phase bin n the resulting McKinnon–Bates image f MKB(n) is given by −1 (p − X f Std). f MKB(n) = f Std + XPCF(n)
The drawback of the MKB method is that image noise is increased according to Gaussian error propagation. Thus, the variances of image pixels in the MKB reconstruction Var( f MKB) are the sum of the variance in the PCF reconstruction Var( f PCF) and the variance in the correction image Var( f corr). The second component of the LDPC algorithm is thus a 5D edge-preserving bilateral filter26 in order to reduce the noise level. These five dimensions sum up as three spatial dimensions, a respiratory-temporal, and a cardiac-temporal dimension. An extensive evaluation of this method can be found in Ref. 8. 2.E. Measurements
More than ten in vivo scans of mice have been acquired and the results for two mice are shown in particular. The data of all mice were measured with a dedicated in vivo cone-beam micro-CT scanner (TomoScope Synergy Twin, CT Imaging GmbH, Erlangen, Germany) in single source mode installed at the Institute of Medical Physics, Erlangen, Germany. The focus to isocenter and the focus to detector distance are RF = 170 mm and RFD = 209 mm, respectively. Operating in 2 × 2 pixel binning mode, the flat detector generates images with a matrix size of 516 × 506 pixels, with a pixel spacing of 100 × 100 µm, and with a data acquisition frame rate of 25 frames/s. The scans were conducted at 65 kV tube voltage with a tube current of 0.3 mA. Per mouse, a total of 7200 projections were acquired within 288 s in a circular trajectory over ten rotations, i.e., an angular range of 3600◦. The respective tube current time product was 86.4 mA s and the absorbed dose was measured as 500 mGy in each case using an ionization chamber (PTW–UNIDOS, PTW, Freiburg, Germany). The detector integration time of 40 ms limits the number of acquired projections during one respiration cycle to 11 and during one cardiac cycle to 6. This, in turn, implies that the respiratory-temporal resolution is limited to about 10% of the respiratory cycle and that the cardiac-temporal resolution is limited to about 20% of the cardiac cycle. The mice were anesthetized with a combination of ketamine and rompun. Some mice were administered with the blood pool agent Exitron nano 12000 (Miltenyi BioTec GmbH, Bergisch Gladbach, Germany and nanoPET Pharma GmbH, Berlin, Germany). For the others, the blood pool agent eXIA 160XL (Binitio Biomedical, Inc., Ottawa, Canada) was used. The mean respiratory rate was between 120 and 150 rpm (respirations per minute) and the mean heart rate was between 260 and 300 bpm (beats per minute). All animal studies Medical Physics, Vol. 42, No. 4, April 2015
F. 6. For acMoCo reconstructions from Fig. 5, the respective difference images are shown to FDK reconstructions (left) and to PCF reconstructions (right) (C = 0 HU, W = 500 HU).
were approved by the ethical committee at the University of Erlangen-Nürnberg. 2.F. Reconstruction parameters and processing time
For each reconstruction algorithm, the images were reconstructed in a 512 × 512 × 512 voxel volume with a 80 × 80 × 80 µm voxel size. A kymogram-based motion detection method27 has been applied to the rawdata to retrieve the respective respiratory and cardiac motion phase information used for retrospective gating. Respiratory-cardiac-correlated image reconstructions were carried out with a respiratory window of ∆r = 10% and a cardiac window width of ∆c = 20%. Volumes were generated with a step size of 5% in the respiratory cycle and a step size of 10% in the cardiac cycle, i.e., the reconstructed volumes overlap in respiratory-temporal and cardiactemporal directions. The total reconstruction time is about one hour. Here, the main part is related to the motion estimation part which consists mainly of nonoptimized code. The cyclic consistent registration required 16.1±1.2 min for the 20 respiratory-gated volumes and 7.4±0.5 min for the 10 cardiac-gated volumes. Since we perform registration a second time on the artifact images, the total registration time is about 47 min. All computations were performed on a 64 bit version of Windows 7 system equipped with two Intel Xeon X5680 processors. However, by optimizing the code and by using GPU hardware, considerably lower timings should be achievable for acMoCo.
1954
Brehm et al.: 5D motion-compensated micro CT imaging
1954
F. 7. Two profiles, as depicted in subfigure (a), were evaluated for FDK (dotted line, circles), PCF (dashed-dotted line, triangles), and acMoCo (dashed line, diamonds). (b) Edges between a contrast-enhanced part of the heart to unenhanced soft tissue and (c) edges between vertebral body (bone) and surrounding soft tissue. (a) ROI positions for edge profiles. (b) Contrast-enhanced part of the heart. (c) Vertebral body.
In addition to the case where all measured projection data corresponding to an absorbed dose of 500 mGy were used for reconstruction, we also investigated cases with less dose by an additional sparsification of the projection data. In more detail, we had a look on cases corresponding to 250, 100, 50, and 25 mGy. Therefore, we considered just a part of the measured projection data, i.e., only projections from the first five rotations out of a total of ten rotations in case of 250 mGy, Medical Physics, Vol. 42, No. 4, April 2015
from the first two rotations in case of 100 mGy, from the first rotation in case of 50 mGy, and only every other projection from the first rotation in case of 25 mGy. 3. RESULTS Figure 5 compares the results of FDK (left), PCF (middle), and acMoCo (right). By applying all available projection data
1955
Brehm et al.: 5D motion-compensated micro CT imaging
1955
F. 8. Phase-correlated reconstructions centered at r = 0% and c = 50%. The upper row shows the PCF reconstructions while the bottom row contains the acMoCo images. Dose levels ranging from 25 to 250 mGy were obtained by using only fractions of the data available (C = 200 HU, W = 1200 HU).
but taking motion not into account, the FDK reconstructions show a low noise level but moving structures appear blurred. Due to the chosen gating parameters, a respiratory window of ∆r = 10% width and a cardiac window of ∆c = 20% width, only about 2% of the projection data correspond to the desired motion phase. Hence, the conventional simple gating approach of PCF suffers from highly increased image noise and streaking artifacts. By acMoCo, we remove all streaking artifacts of PCF while maintaining the temporal resolution. In addition, the noise level is comparable to that of the FDK reconstruction. The respective difference images are shown in Fig. 6. Edge profiles [see Fig. 7(a)] were measured within the volumes reconstructed by FDK, PCF, and acMoCo to further quantify that the proposed motion compensation does not decrease spatial resolution of FDK and that it does maintain temporal resolution of PCF. Figure 7(b) shows the change in
CT values between contrast-enhanced structure of the heart and the unenhanced soft tissue nearby. In Fig. 7(c), edges between vertebral body and surrounding soft tissue are shown. All profiles indicate that there is no loss of edge information due to acMoCo compared to a PCF reconstruction in terms of temporal resolution and compared to a FDK reconstruction in terms of spatial resolution. We further investigated the robustness of motion compensation against decreasing dose levels which lead to increased image noise and streak artifacts within gated reconstructions used for motion estimation (see Figs. 8 and 9). Image quality of conventional phase-correlated reconstruction dramatically decreases and anatomical details vanish due to the severe artifacts and high noise. In comparison, acMoCo almost maintains high image quality of a high dose scan. Due to the decreased amount of projection data contributing to the reconstruction
F. 9. Phase-correlated reconstructions centered at r = 0% and c = 50% for another mouse. The upper row shows the PCF reconstructions while the bottom row contains the acMoCo images. Dose levels ranging from 25 to 250 mGy were obtained by using only fractions of the data available (C = 200 HU, W = 1200 HU). Medical Physics, Vol. 42, No. 4, April 2015
1956
Brehm et al.: 5D motion-compensated micro CT imaging
1956
F. 10. Phase-correlated reconstructions centered at r = 0% and c = 20% applying LDPC and acMoCo are shown for two dose levels, 500 mGy (two columns left) and 100 mGy (two columns right) (C = 200 HU, W = 1200 HU).
at decreased dose levels, the image noise is increasing as it is the case for a FDK reconstruction. But the streaking artifacts are still removed, the temporal resolution of the heart is maintained, and anatomical details remain being well recovered. While this holds for all investigated dose levels in case of the first data set shown in Fig. 8, it is slightly different for the 25 mGy dose level of the second data set. In case of 25 mGy, no projection data remained after double gating for certain phase windows. While a PCF is not possible in that case as shown in Fig. 9, acMoCo overcomes this issue by its sequential motion compensation approach. And even in that case, the temporal resolution was maintained for the left and the right ventricle. However, the severe data sparsification results in a loss of temporal resolution, as can be clearly seen near the left atrium and the diaphragm.
The evaluation also includes a comparison to the results of the LDPC algorithm, another recently published reconstruction algorithm that is incorporating the temporal dimension into the reconstruction process. Figure 10 shows the reconstruction results of LDPC and acMoCo for two different dose levels. LDPC is capable to completely remove the streaking artifacts for the higher dose level, but for the lower dose levels, severe streak artifacts remain. In comparison, acMoCo completely removes the streaking artifacts for both dose levels. For acMoCo, the noise level, measured in a homogeneous soft tissue region, is 24 HU in the 500 mGy case and is 37 HU in the 100 mGy case. In comparison, the noise level for LDPC is 41 HU in the 500 mGy case and 67 HU in the 100 mGy case. Due to the increased dose usage, the noise level of acMoCo is lower than for LDPC and more anatomical details are visible.
F. 11. Importance of sequential gating and artifact model-based motion estimation: Results for sMoCo are shown with concurrent gating (left) and with sequential gating (middle) when applying data corresponding to 50 mGy. For comparison, results for acMoCo are shown with sequential gating (right) (C = 200 HU, W = 1200 HU). Medical Physics, Vol. 42, No. 4, April 2015
1957
Brehm et al.: 5D motion-compensated micro CT imaging
1957
F. 12. Motion-compensated reconstructions acMoCo are shown in different cardiac motion phases (C = 200 HU, W = 1200 HU).
The quality, however, depends on the success of motion estimation to retrieve the motion vector fields required for motion compensation. To illustrate how important the proposed sequential double gating and the artifact model-based motion estimation is, we prepared Fig. 11. If we apply a standard motion compensation, sMoCo with simultaneous double gating streaking artifacts partially remain. When using the sequential double gating approach for sMoCo, the results improve but the robustness against streaking artifacts is lower than for proposed acMoCo. These results apply not only for the shown respiratory and cardiac window but also for all others. Figure 12 shows axial slices of acMoCo from different cardiac phase bins for a fixed respiratory phase bin.
4. SUMMARY AND DISCUSSION We proposed a motion-compensated image reconstruction algorithm acMoCo using a sequential double gating approach combined with an artifact model-based motion estimation including a cyclic consistent nonrigid registration algorithm. This technique improves phase-correlated imaging of highly undersampled data. We demonstrated that dose savings of about an order of magnitude are possible compared to the conventional phase-correlated reconstruction approach. Both, sequential double gating and artifact model-based motion estimation, are essential parts of successful motion compensation for decreasing dose levels. Concurrent double gating increases the impact of sparse-view artifacts, such that a standard motion compensation sMoCo fails in this case. But also when using sequential double gating, sMoCo shows sensitivity to the remained streaking artifacts whose magnitude increases with decreasing dose levels. The artifact modelbased motion estimation including cyclic consistent nonrigid registration shows no sensitivity to these streaking artifacts even for cases with a one magnitude lower dose level. In addition, by the sequential double gating approach, higher Medical Physics, Vol. 42, No. 4, April 2015
sparsifications of the projection data are possible for acMoCo. While conventional phase-correlated reconstruction with its concurrent double gating fails because no projection is left for some desired phase window, by sequential double gating, the missing data can be borrowed from other respiratory phases. Furthermore, acMoCo outperforms the LDPC algorithm in terms of noise level and robustness with decreasing dose level. During the evaluation process, we used cardiac window width of 20% of the RR interval for image reconstruction. This is the minimum reasonable window size with respect to the detector integration time of 40 ms and the mean RR interval of about 200 ms. For accurate cardiac function estimation, e.g., ejection fraction, this window size could be insufficient as accuracy is affected by temporal blurring, in particular in systole which is the shortest cardiac phase. We want to emphasize that this limitation is only due to the detector integration time of the imaging system and not due to the image reconstruction approach. If shorter detector integration times are available and thus smaller cardiac window sizes are reasonable, a reconstruction algorithm has to deal with a reduced dose fraction per cardiac phase at constant total dose level. And we have demonstrated that acMoCo can produce reliable results with decreasing dose level. With the proposed motion compensation acMoCo, it is possible to perform 4D or 5D phase-correlated imaging at the same dose level as required for conventional 3D studies. As a consequence acMoCo allows for longitudinal in vivo studies of the rodent heart as well as for long term studies in preclinical research.
ACKNOWLEDGMENTS This study was supported by Varian Medical Systems, Palo Alto, CA and in parts by the Deutsche Forschungsgemeinschaft (DFG) under Grant No. FOR 661. The highspeed image reconstruction software RayConStruct–IR was provided by RayConStruct GmbH, Nürnberg, Germany. The
1958
Brehm et al.: 5D motion-compensated micro CT imaging
authors thank PD Dr. Andreas Hess and Sandra Strobelt for help with the mouse measurements. The authors thank Miltenyi BioTec GmbH, Bergisch Gladbach, Germany; nanoPET Pharma GmbH, Berlin, Germany; and Binitio Biomedical, Inc., Ottawa, Canada for providing their latest contrast agents. a)Author
to whom correspondence should be addressed. Electronic mail: [email protected]; Telephone: +49 (6221) 42 3067; Fax: +49 (6221) 42 2585. 1C. T. Badea, B. Fubara, L. W. Hedlund, and G. A. Johnson, “4–D micro–CT of the mouse heart,” Mol. Imaging 4, 110–116 (2005). 2M. Drangova, N. L. Ford, S. A. Detombe, A. R. Wheatley, and D. W. Holdsworth, “Fast retrospectively gated quantitative four–dimensional (4D) cardiac micro computed tomography imaging of free–breathing mice,” Invest. Radiol. 42, 85–94 (2007). 3S. Kriminski, M. Mitschke, S. Sorensen, N. M. Wink, P. E. Chow, S. Tenn, and T. D. Solberg, “Respiratory correlated cone–beam computed tomography on an isocentric C–arm,” Phys. Med. Biol. 50, 5263–5280 (2005). 4J.-J. Sonke, L. Zijp, P. Remeijer, and M. van Herk, “Respiratory correlated cone beam CT,” Med. Phys. 32, 1176–1186 (2005). 5G. C. McKinnon and R. Bates, “Towards imaging the beating heart usefully with a conventional CT scanner,” IEEE Trans. Biomed. Eng. BME-28, 123–127 (1981). 6K. L. Garden and R. A. Robb, “3–D reconstruction of the heart from few projections: A pratical implementation of the McKinnon–Bates algorithm,” IEEE Trans. Med. Imaging 5, 233–239 (1986). 7S. Leng, J. Zambelli, B. Nett, R. Tolakanahalli, P. Munro, J. Star-Lack, B. Paliwal, and G.-H. Chen, “Streaking artifacts reduction in four–dimensional cone–beam computed tomography,” Med. Phys. 35, 4649–4659 (2008). 8S. Sawall, F. Bergner, R. Lapp, M. Mronz, M. Karolczak, A. Hess, and M. Kachelrieß, “Low–dose cardio–respiratory phase–correlated cone–beam micro–CT of small animals,” Med. Phys. 38, 1416–1424 (2011). 9M. Persson, D. Bone, and H. Elmqvist, “Total variation norm for three– dimensional iterative reconstruction in limited view angle tomography,” Phys. Med. Biol. 46, 853–866 (2001). 10E. Y. Sidky and X. Pan, “Image reconstruction in circular cone–beam computed tomography by constrained, total–variation minimization,” Phys. Med. Biol. 53, 4777–4807 (2008). 11J. Song, Q. H. Liu, G. A. Johnson, and C. T. Badea, “Sparseness prior based iterative image reconstruction for retrospectively gated cardiac micro–CT,” Med. Phys. 34, 4476–4483 (2007). 12L. Ritschl, S. Sawall, M. Knaup, A. Hess, and M. Kachelrieß, “Iterative 4D cardiac micro–CT image reconstruction using an adaptive spatio–temporal sparsity prior,” Phys. Med. Biol. 57, 1517–1525 (2012).
Medical Physics, Vol. 42, No. 4, April 2015
1958 13C.
R. Crawford, K. King, C. Ritchie, and J. D. Godwin, “Respiratory compensation in projection imaging using a magnification and displacement model,” IEEE Trans. Med. Imaging 15, 327–332 (1996). 14C. Ritchie, C. Crawford, J. D. Godwin, K. King, and Y. Kim, “Correction of computed tomography motion artifacts using pixel–specific back–projection,” IEEE Trans. Med. Imaging 15, 333–342 (1996). 15T. Li, E. Schreibmann, Y. Yang, and L. Xing, “Motion correction for improved target localization with on–board cone–beam computed tomography,” Phys. Med. Biol. 51, 253–267 (2006). 16C. T. Badea, E. Schreibmann, and T. Fox, “A registration based approach for 4D cardiac micro–CT using combined prospective and retrospective gating,” Med. Phys. 35, 1170–1179 (2008). 17M. Brehm, P. Paysan, M. Oelhafen, and M. Kachelrieß, “Artifact–resistant motion estimation with a patient–specific artifact model for motion– compensated cone–beam CT,” Med. Phys. 40, 101913 (13pp.) (2013). 18L. Feldkamp, L. Davis, and J. Kress, “Practical cone–beam algorithm,” J. Opt. Soc. Am. 1, 612–619 (1984). 19M. Brehm, P. Paysan, M. Oehlhafen, P. Kunz, and M. Kachelrieß, “Self– adapting cyclic registration for motion–compensated cone–beam CT in image–guided radiation therapy,” Med. Phys. 39, 7603–7618 (2012). 20E. Meyer, R. Raupach, M. Lell, B. Schmidt, and M. Kachelrieß, “Normalized metal artifact reduction (NMAR) in computed tomography,” Med. Phys. 37, 5482–5493 (2010). 21M. Sun, M. Oelhafen, L. Carvalho, J. Pavkovich, T. Berkus, and J. Star-Lack, “Improving motion accuracy for the McKinnon–Bates 4D– CBCT reconstruction algorithm,” in Proccedings of 2nd International Meeting on Image Formation in X-ray CT (Salt Lake City, Utah, 2012), pp. 131–134. 22P.-S. Liao, T.-S. Chen, and P.-C. Chung, “A fast algorithm for multilevel thresholding,” J. Inf. Sci. Eng. 17, 713–727 (2001). 23J.-P. Thirion, “Image matching as a diffusion process: An analogy with Maxwell’s demons,” Med. Image Anal. 2, 243–260 (1998). 24H. Wang, L. Dong, J. O’Daniel, R. Mohan, A. S. Garden, K. K. Ang, D. A. Kuban, M. Bonnen, J. Y. Chang, and R. Cheung, “Validation of an accelerated ‘demons’ algorithm for deformable image registration in radiation therapy,” Phys. Med. Biol. 50, 2887–2905 (2005). 25T. Vercauteren, X. Pennec, A. Perchant, and N. Ayache, “Diffeomorphic demons: Efficient non–parametric image registration,” NeuroImage 45, 61–72 (2009). 26C. Tomasi and R. Manduchi, “Bilateral filtering for gray and color images,” in Proceedings of the Sixth International Conference on Computer Vision (IEEE Computer Society, Washington, 1998), pp. 839–846. 27M. Kachelrieß, D.-A. Sennst, W. Maxlmoser, and W. A. Kalender, “Kymogram detection and kymogram–correlated image reconstruction from subsecond spiral computed tomography scans of the heart,” Med. Phys. 29, 1489–1503 (2002).
Stefan Sawall, Joscha Maier, Sebastian Sauppe, and Marc Kachelrießa) German Cancer Research Center (DKFZ), Im Neuenheimer Feld 280, Heidelberg 69120, Germany
(Received 3 November 2014; revised 28 January 2015; accepted for publication 26 February 2015; published 26 March 2015) Purpose: Cardiac in vivo micro-CT imaging of small animals typically requires double gating due to long scan times and high respiratory rates. The simultaneous respiratory and cardiac gating can either be done prospectively or retrospectively. In any case, for true 5D imaging, i.e., three spatial dimensions plus one respiratory-temporal dimension plus one cardiac temporal dimension, the amount of information corresponding to a given respiratory and cardiac phase is orders of magnitude lower than the total amount of information acquired. Achieving similar image quality for 5D than for usual 3D investigations would require increasing the amount of data and thus the applied dose to the animal. Therefore, the goal is phase-correlated imaging with high image quality but without increasing the dose level. Methods: To achieve this, the authors propose a new image reconstruction algorithm that makes use of all available projection data, also of that corresponding to other motion windows. In particular, the authors apply a motion-compensated image reconstruction approach that sequentially compensates for respiratory and cardiac motion to decrease the impact of sparsification. In that process, all projection data are used no matter which motion phase they were acquired in. Respiratory and cardiac motion are compensated for by using motion vector fields. These motion vector fields are estimated from initial phase-correlated reconstructions based on a deformable registration approach. To decrease the sensitivity of the registration to sparse-view artifacts, an artifact model-based approach is used including a cyclic consistent nonrigid registration algorithm. Results: The preliminary results indicate that the authors’ approach removes the sparse-view artifacts of conventional phase-correlated reconstructions while maintaining temporal resolution. In addition, it achieves noise levels and spatial resolution comparable to that of nongated reconstructions due to the improved dose usage. By using the proposed motion estimation, no sensitivity to streaking artifacts has been observed. Conclusions: Using sequential double gating combined with artifact model-based motion estimation allows to accurately estimate respiratory and cardiac motion from highly undersampled data. No sensitivity to streaking artifacts introduced by sparse angular sampling has been observed for the investigated dose levels. The motion-compensated image reconstruction was able to correct for both, respiratory and cardiac motion, by applying the estimated motion vector fields. The administered dose per animal can thus be reduced for 5D imaging allowing for longitudinal studies at the highest image quality. C 2015 American Association of Physicists in Medicine. [http://dx.doi.org/10.1118/1.4916083] Key words: micro-CT, phase correlation, motion compensation
1. INTRODUCTION In vivo cone-beam micro-CT imaging of small animals is an appreciated technology in preclinical research due to its ability to provide 3D information of living objects with very high spatial resolution. An active field of research is cardiac imaging, as it has already been for clinical CT the last decade. However, small animals have higher respiratory rates and faster heart beats than humans, e.g., a mouse breathes between 90 and 200 times/min and a mouse heart beats between 300 and 600 times/min. For this very reason, contrast-enhanced microCT imaging of animal hearts is a particular challenge. 1948
Med. Phys. 42 (4), April 2015
In case of cardiac imaging, also respiratory motion has to be considered due to breathing-induced translation of the heart. One or more temporal windows of interest are described by the centers of the respiratory phase r and the cardiac phase c and by the respective widths ∆r and ∆c. The required double gating, a simultaneous respiratory and cardiac gating, can either be done prospectively1 or retrospectively.2 In any case, retrieving true 5D information, i.e., three spatial dimensions plus one respiratory-temporal dimension plus one cardiac temporal dimension, means that just a fraction (∆r = 10% and ∆c = 20% ⇒ 2% total) of the total amount of projection data acquired corresponds to a given respiratory and cardiac phase.
0094-2405/2015/42(4)/1948/11/$30.00
© 2015 Am. Assoc. Phys. Med.
1948
1949
Brehm et al.: 5D motion-compensated micro CT imaging
The consequences are severe sparse-view artifacts when using conventional reconstruction algorithms (see Fig. 1). The possibilities on the acquisition side are quite low to overcome the problem of limited projection data in each motion phase. In particular, for long term studies, dose issues become important in small animal imaging. For this reason, the limitation in amount of projection data has to be considered by dedicated reconstruction algorithms. The standard reconstruction approach for temporal-correlated volumes is to sort the data into phase bins and reconstruct each bin separately. However, the resulting volumes suffer from streak artifacts induced by the increased angular spacing of projections used for reconstruction. This problem is also known from C-arm CBCT systems in interventional radiology3 and onboard CBCT imaging in radiation therapy.4 One way to deal with the problem of missing projection data in each motion phase is to utilize data from other phases by explicitly incorporating the temporal dimension into the reconstruction process.5–7 For cardiac micro-CT imaging, a variant of the McKinnon–Bates (MKB) image reconstruction algorithm in combination with an edge-preserving noisesuppression filter has been proposed.8 Another way for regularization is the usage of iterative reconstruction algorithms that incorporate some a priori knowledge into the reconstruction process.9,10 This approach has been applied to micro-CT data without and with incorporating temporal dimension.11,12 Both non-motion-compensated approaches, filter-based 5D noise suppression and iterative reconstruction with 5D regularization, considerably reduce the image noise and artifact levels. Furthermore, those are straightforward and need only minor changes to existing algorithms. But full dose usage, i.e., contribution of each projection to every phase volume, is not reached by any of those two approaches. The reason for that is that the prior knowledge imposed by a motion compensation (MoCo) approach is much stronger: There motion is modeled and voxels may change their position.With the more conventional approaches, this knowledge is neither explicitly nor implicitly used. In this work, we want to extend the idea of incorporating cardiac- and respiratory-temporal dimension one step further and apply all projection data for each phase reconstruction by motion compensation.13–16 Here, all acquired projections contribute to each volume, independent of the desired cardiac and respiratory phase, without a loss in temporal resolution
1949
due to the incorporation of motion vector fields (MVFs). In contrast to Ref. 16, we do not require an additional prospectivegated scan for motion estimation which is time-consuming and increases dose. Instead, we apply an artifact-insensitive estimation method that is made robust by incorporating two main ideas: The consideration of the impact of image artifacts by an artifact model that we successfully applied for onboard CBCT imaging in radiation therapy17 and the consideration of the presence of two separate motions influencing each other by a sequential double gating.
2. MATERIALS AND METHODS 2.A. MoCo based on phase-correlated Feldkamp
We use the well-known Feldkamp–Davis–Kress (FDK) filtered backprojection18 as our standard reconstruction algorithm. The FDK algorithm is denoted by f Std = X−1 p, where X is the x-ray transform operator, X−1 is the reconstruction operator, p are the projection data, and f Std the reconstructed image. For phase-correlated Feldkamp (PCF) reconstruction, the projection images are retrospectively sorted into several phase bins. Subsequently, each phase bin is reconstructed separately using the FDK algorithm. For an arbitrary phase −1 bin n, the corresponding operator is denoted as XPCF(n) and the respective PCF image f PCF(n) is given by −1 f PCF(n) = XPCF(n) p.
However, if just a limited number of projections are acquired, the binning induces a large angular spacing of projections used for a single reconstruction. This results in viewaliasing artifacts that deteriorate image quality of these phasecorrelated images. Using a motion compensation approach, projection data from other phases are additionally considered for reconstruction to improve the angular sampling rate. An application of all projection data without considering motion, as it is the case for a standard FDK reconstruction, obviously causes inconsistencies in reconstruction data and ends up in motion blurring. But if motion is modeled accurately, a MoCo volume can be obtained without motion blurring as we have shown previously for respiratory motion in radiation therapy onboard CBCT scans.17,19 The compensation is done by warping the PCF reconstructions according to the motion model before superimposition. For arbitrary phase bins n and m, the MVF describing the motion from phase m to phase n is denoted as Tnm and the respective MoCo image f MoCo(n) is given by ( ) −1 f MoCo(n) = XPCF(m) p ◦ Tnm . m
2.B. Sequential motion compensation for respiratory and cardiac motion
F. 1. Ungated 3D reconstruction (left) using all measured data and doublegated 5D reconstruction (right) using only 2% of the data, corresponding to a given respiratory and cardiac motion phase (C = 200 HU, W = 1200 HU). Medical Physics, Vol. 42, No. 4, April 2015
For cardiac imaging, two different motion patterns are superimposed and have to be taken into account. Besides the cardiac motion, also the respiration-induced displacement of the heart is of importance. We consider both motion components
1950
Brehm et al.: 5D motion-compensated micro CT imaging
separately within a sequential motion compensation approach. In a first step, only respiratory motion is considered while cardiac motion is completely ignored, i.e., the respiratory phasecorrelated images (see Fig. 2), given by gating of respiratory motion only, are used for the first step’s motion estimation and compensation. In consequence, for each pair of respiratory phase bins (r n ,r m ), the corresponding motion vector field Trr nm is given such that ( ) −1 f MoCo(r n ) = p ◦ Trr nm . XPCF(r ) m rm
In a second step, cardiac motion is considered while respiratory motion is already compensated for a given respiratory phase, i.e., the phase-correlated images used for the second step’s motion estimation and compensation are given by gating of cardiac motion and by compensating of respiratory motion according to the motion model estimated in step one (see Fig. 2). For each pair of cardiac phase bins (cn ,cm ), the corresponding motion vector field Tcc nm is derived from those phasecorrelated images. The double motion compensated image f MoCo(r n,c n ) is then given by f MoCo(r n,c n ) =
( ) −1 * XPCF(r p ◦ Trr nm + ◦ Tcc nm . m,c m ) cm , rm -
1950
2.C. Artifact model-based motion estimation
As no other data are available, motion estimation is performed on the PCF reconstructions. Unfortunately, the presence of view-aliasing artifacts within the PCF images and their change in appearance for different phases influences, the registration process such that conventional approaches fail. To decrease the sensitivity on these artifacts, we are incorporating knowledge from motion patterns and an artifact model into the estimation process, an approach that has been successfully applied to respiratory motion in radiation therapy onboard CBCT scans.17 The applied artifact model contains information about location and intensities, e.g., of view-aliasing artifacts. While containing equivalent artifacts as the PCF reconstruction, the model is designed to be free of subject motion. A separation of estimated motion induced by subject motion and induced by artifacts is thus possible later on. In particular, the model is a second series of 3D images, called 4D artifact images in the following. Those 4D artifact images are generated as follows. In order to have a subject-specific model, we use a prior image similar to that used in Ref. 20. The prior image is generated by an FDK reconstruction of all measured data followed by a segmentation. By applying segmentation, all motion artifacts like blurring and streaking can be removed from the FDK
F. 2. Four different phase-correlated reconstructions of a selected mouse dataset: The double-gated reconstructions (left column) suffer from severe sparse-view artifacts. In contrast, the images used for the sequential gating and motion compensation approach, the respiratory-gated images (second column left) and the respiratory-compensated cardiac-gated images (second column right), include less artifacts. This allows for a more robust motion estimation and thus for motion-compensated images (right column) of high image quality (C = 200 HU, W = 1200 HU). Medical Physics, Vol. 42, No. 4, April 2015
1951
Brehm et al.: 5D motion-compensated micro CT imaging
1951
F. 3. Prior image derived from the FDK image by simple multilevel thresholding using automatically chosen values via cluster analysis (C = 200 HU, W = 1200 HU).
image. Thus, motion cannot be a longer part of a forwardprojected dataset in order to obtain a motion-free model.21 The segmentation consists of an identification for regions like air, lung tissue, fat tissue, water tissue, and bone by simple
thresholding with automatically chosen values via cluster analysis.22 By this classification, the cardiovascular system which is enhanced by the use of a contrast agent will be assigned to the bone region. Subsequently, the segmented regions are set
F. 4. Illustration of the proposed motion estimation method using an artifact model: For both, measured projection data as well as simulated projection data obtained by forward projection of a prior image, a PCF reconstruction is performed taking respiratory motion into account only. Subsequently, MVFs are estimated for both image series using deformable registration and then combined to separate respiratory motion from artifact motion. The procedure is repeated for cardiac motion while respiratory motion is already compensated. Medical Physics, Vol. 42, No. 4, April 2015
1952
Brehm et al.: 5D motion-compensated micro CT imaging
to the mean value of the particular region except for voxels of the bone region that retain their CT-values, as they vary too much to properly model them with just one value (see Fig. 3). In a second step, the measurement and reconstruction process is simulated. A forward projection is performed from the resulting prior image (see Fig. 3) in the same geometry the measured projection data were acquired. These simulated projection data undergo the same phase-correlated image reconstruction using the same respiratory signal as the measured projection data. The resulting 4D CBCT images model the desired artifacts while being free of subject motion. Two independent registrations are performed for both image series, based on measured and on simulated projection data, such that two sets of MVFs are obtained (see Fig. 4). While the first MVF contains components induced by subject motion and induced by image artifacts, the second one contains components induced by image artifacts only. In a final step, the second MVF is subtracted from the first one in order to reduce the artifact-induced MVF components of that MVF. Thus, a MVF is obtained that consists mainly of subject motion. The registration method used for estimating the MVFs of an image series is a deformable registration algorithm in combination with a suited strategy we used previously.19 Basis of this motion estimation is the demons algorithm,23–25 a nonrigid intensity-based registration algorithm, and the strategy to just estimate motion between adjacent phases instead of doing so for each pair of motion phases. Here, MVFs of nonadjacent phases are given by concatenation of MVFs from adjacent phases. This strategy allows for a cyclic consistent registration by incorporating temporal constraints like the quasiperiodic
1952
breathing or cardiac motion pattern. Thus, the sensitivity of the registration process to image artifacts is reduced. For the evaluation, we perform motion-compensated image reconstructions without and with our approach. Within a standard motion compensation, denoted by sMoCo, each required MVF Tnm is estimated independently by a standard deformable intensity-based 3D-3D image registration method matching the corresponding two PCF reconstructions f PCF(m) and f PCF(n). In our case, we used the demons algorithm.23–25 Here, the MVFs originate from the PCF reconstructions only and are not corrected using the artifact model. In accordance with Fig. 4, artifact model-based motion compensation, denoted by acMoCo, is performed applying MVFs that are obtained by cyclic registration and correction based on the artifact model. 2.D. Low-dose phase-correlated (LDPC) reconstruction
For evaluation purposes, we additionally perform image reconstructions using another reconstruction paradigm allowing for the reconstruction of double phase-correlated images. Here, we apply the LDPC algorithm proposed by Sawall et al. in Ref. 8 that incorporates the temporal dimension into the reconstruction process and consists of two components. The first component of the LDPC algorithm is the MKB reconstruction method5,6 that was introduced to reduce streak artifacts in the PCF reconstructions of the human heart. The idea of the MKB method is to add an image f corr to the PCF reconstruction that corrects for the streak artifacts. To obtain
F. 5. Three orthogonal slices are shown for images originating from a standard (FDK, left), a conventional phase-correlated (PCF, middle), and a motion-compensated reconstruction (acMoCo, right) at end-expiration (r = 0%) and end-systole (c = 60%) (C = 200 HU, W = 1200 HU). Medical Physics, Vol. 42, No. 4, April 2015
1953
Brehm et al.: 5D motion-compensated micro CT imaging
1953
the correction image, a forward projection of a standard reconstruction f Std is performed. A phase-correlated reconstruction of that forward projected projection data yields the standard image with similar streak artifacts as the PCF reconstruction. The final correction image f corr that exhibits only the streak artifacts can be received by subtracting this streaky image from f Std, i.e., for an arbitrary phase bin n the resulting McKinnon–Bates image f MKB(n) is given by −1 (p − X f Std). f MKB(n) = f Std + XPCF(n)
The drawback of the MKB method is that image noise is increased according to Gaussian error propagation. Thus, the variances of image pixels in the MKB reconstruction Var( f MKB) are the sum of the variance in the PCF reconstruction Var( f PCF) and the variance in the correction image Var( f corr). The second component of the LDPC algorithm is thus a 5D edge-preserving bilateral filter26 in order to reduce the noise level. These five dimensions sum up as three spatial dimensions, a respiratory-temporal, and a cardiac-temporal dimension. An extensive evaluation of this method can be found in Ref. 8. 2.E. Measurements
More than ten in vivo scans of mice have been acquired and the results for two mice are shown in particular. The data of all mice were measured with a dedicated in vivo cone-beam micro-CT scanner (TomoScope Synergy Twin, CT Imaging GmbH, Erlangen, Germany) in single source mode installed at the Institute of Medical Physics, Erlangen, Germany. The focus to isocenter and the focus to detector distance are RF = 170 mm and RFD = 209 mm, respectively. Operating in 2 × 2 pixel binning mode, the flat detector generates images with a matrix size of 516 × 506 pixels, with a pixel spacing of 100 × 100 µm, and with a data acquisition frame rate of 25 frames/s. The scans were conducted at 65 kV tube voltage with a tube current of 0.3 mA. Per mouse, a total of 7200 projections were acquired within 288 s in a circular trajectory over ten rotations, i.e., an angular range of 3600◦. The respective tube current time product was 86.4 mA s and the absorbed dose was measured as 500 mGy in each case using an ionization chamber (PTW–UNIDOS, PTW, Freiburg, Germany). The detector integration time of 40 ms limits the number of acquired projections during one respiration cycle to 11 and during one cardiac cycle to 6. This, in turn, implies that the respiratory-temporal resolution is limited to about 10% of the respiratory cycle and that the cardiac-temporal resolution is limited to about 20% of the cardiac cycle. The mice were anesthetized with a combination of ketamine and rompun. Some mice were administered with the blood pool agent Exitron nano 12000 (Miltenyi BioTec GmbH, Bergisch Gladbach, Germany and nanoPET Pharma GmbH, Berlin, Germany). For the others, the blood pool agent eXIA 160XL (Binitio Biomedical, Inc., Ottawa, Canada) was used. The mean respiratory rate was between 120 and 150 rpm (respirations per minute) and the mean heart rate was between 260 and 300 bpm (beats per minute). All animal studies Medical Physics, Vol. 42, No. 4, April 2015
F. 6. For acMoCo reconstructions from Fig. 5, the respective difference images are shown to FDK reconstructions (left) and to PCF reconstructions (right) (C = 0 HU, W = 500 HU).
were approved by the ethical committee at the University of Erlangen-Nürnberg. 2.F. Reconstruction parameters and processing time
For each reconstruction algorithm, the images were reconstructed in a 512 × 512 × 512 voxel volume with a 80 × 80 × 80 µm voxel size. A kymogram-based motion detection method27 has been applied to the rawdata to retrieve the respective respiratory and cardiac motion phase information used for retrospective gating. Respiratory-cardiac-correlated image reconstructions were carried out with a respiratory window of ∆r = 10% and a cardiac window width of ∆c = 20%. Volumes were generated with a step size of 5% in the respiratory cycle and a step size of 10% in the cardiac cycle, i.e., the reconstructed volumes overlap in respiratory-temporal and cardiactemporal directions. The total reconstruction time is about one hour. Here, the main part is related to the motion estimation part which consists mainly of nonoptimized code. The cyclic consistent registration required 16.1±1.2 min for the 20 respiratory-gated volumes and 7.4±0.5 min for the 10 cardiac-gated volumes. Since we perform registration a second time on the artifact images, the total registration time is about 47 min. All computations were performed on a 64 bit version of Windows 7 system equipped with two Intel Xeon X5680 processors. However, by optimizing the code and by using GPU hardware, considerably lower timings should be achievable for acMoCo.
1954
Brehm et al.: 5D motion-compensated micro CT imaging
1954
F. 7. Two profiles, as depicted in subfigure (a), were evaluated for FDK (dotted line, circles), PCF (dashed-dotted line, triangles), and acMoCo (dashed line, diamonds). (b) Edges between a contrast-enhanced part of the heart to unenhanced soft tissue and (c) edges between vertebral body (bone) and surrounding soft tissue. (a) ROI positions for edge profiles. (b) Contrast-enhanced part of the heart. (c) Vertebral body.
In addition to the case where all measured projection data corresponding to an absorbed dose of 500 mGy were used for reconstruction, we also investigated cases with less dose by an additional sparsification of the projection data. In more detail, we had a look on cases corresponding to 250, 100, 50, and 25 mGy. Therefore, we considered just a part of the measured projection data, i.e., only projections from the first five rotations out of a total of ten rotations in case of 250 mGy, Medical Physics, Vol. 42, No. 4, April 2015
from the first two rotations in case of 100 mGy, from the first rotation in case of 50 mGy, and only every other projection from the first rotation in case of 25 mGy. 3. RESULTS Figure 5 compares the results of FDK (left), PCF (middle), and acMoCo (right). By applying all available projection data
1955
Brehm et al.: 5D motion-compensated micro CT imaging
1955
F. 8. Phase-correlated reconstructions centered at r = 0% and c = 50%. The upper row shows the PCF reconstructions while the bottom row contains the acMoCo images. Dose levels ranging from 25 to 250 mGy were obtained by using only fractions of the data available (C = 200 HU, W = 1200 HU).
but taking motion not into account, the FDK reconstructions show a low noise level but moving structures appear blurred. Due to the chosen gating parameters, a respiratory window of ∆r = 10% width and a cardiac window of ∆c = 20% width, only about 2% of the projection data correspond to the desired motion phase. Hence, the conventional simple gating approach of PCF suffers from highly increased image noise and streaking artifacts. By acMoCo, we remove all streaking artifacts of PCF while maintaining the temporal resolution. In addition, the noise level is comparable to that of the FDK reconstruction. The respective difference images are shown in Fig. 6. Edge profiles [see Fig. 7(a)] were measured within the volumes reconstructed by FDK, PCF, and acMoCo to further quantify that the proposed motion compensation does not decrease spatial resolution of FDK and that it does maintain temporal resolution of PCF. Figure 7(b) shows the change in
CT values between contrast-enhanced structure of the heart and the unenhanced soft tissue nearby. In Fig. 7(c), edges between vertebral body and surrounding soft tissue are shown. All profiles indicate that there is no loss of edge information due to acMoCo compared to a PCF reconstruction in terms of temporal resolution and compared to a FDK reconstruction in terms of spatial resolution. We further investigated the robustness of motion compensation against decreasing dose levels which lead to increased image noise and streak artifacts within gated reconstructions used for motion estimation (see Figs. 8 and 9). Image quality of conventional phase-correlated reconstruction dramatically decreases and anatomical details vanish due to the severe artifacts and high noise. In comparison, acMoCo almost maintains high image quality of a high dose scan. Due to the decreased amount of projection data contributing to the reconstruction
F. 9. Phase-correlated reconstructions centered at r = 0% and c = 50% for another mouse. The upper row shows the PCF reconstructions while the bottom row contains the acMoCo images. Dose levels ranging from 25 to 250 mGy were obtained by using only fractions of the data available (C = 200 HU, W = 1200 HU). Medical Physics, Vol. 42, No. 4, April 2015
1956
Brehm et al.: 5D motion-compensated micro CT imaging
1956
F. 10. Phase-correlated reconstructions centered at r = 0% and c = 20% applying LDPC and acMoCo are shown for two dose levels, 500 mGy (two columns left) and 100 mGy (two columns right) (C = 200 HU, W = 1200 HU).
at decreased dose levels, the image noise is increasing as it is the case for a FDK reconstruction. But the streaking artifacts are still removed, the temporal resolution of the heart is maintained, and anatomical details remain being well recovered. While this holds for all investigated dose levels in case of the first data set shown in Fig. 8, it is slightly different for the 25 mGy dose level of the second data set. In case of 25 mGy, no projection data remained after double gating for certain phase windows. While a PCF is not possible in that case as shown in Fig. 9, acMoCo overcomes this issue by its sequential motion compensation approach. And even in that case, the temporal resolution was maintained for the left and the right ventricle. However, the severe data sparsification results in a loss of temporal resolution, as can be clearly seen near the left atrium and the diaphragm.
The evaluation also includes a comparison to the results of the LDPC algorithm, another recently published reconstruction algorithm that is incorporating the temporal dimension into the reconstruction process. Figure 10 shows the reconstruction results of LDPC and acMoCo for two different dose levels. LDPC is capable to completely remove the streaking artifacts for the higher dose level, but for the lower dose levels, severe streak artifacts remain. In comparison, acMoCo completely removes the streaking artifacts for both dose levels. For acMoCo, the noise level, measured in a homogeneous soft tissue region, is 24 HU in the 500 mGy case and is 37 HU in the 100 mGy case. In comparison, the noise level for LDPC is 41 HU in the 500 mGy case and 67 HU in the 100 mGy case. Due to the increased dose usage, the noise level of acMoCo is lower than for LDPC and more anatomical details are visible.
F. 11. Importance of sequential gating and artifact model-based motion estimation: Results for sMoCo are shown with concurrent gating (left) and with sequential gating (middle) when applying data corresponding to 50 mGy. For comparison, results for acMoCo are shown with sequential gating (right) (C = 200 HU, W = 1200 HU). Medical Physics, Vol. 42, No. 4, April 2015
1957
Brehm et al.: 5D motion-compensated micro CT imaging
1957
F. 12. Motion-compensated reconstructions acMoCo are shown in different cardiac motion phases (C = 200 HU, W = 1200 HU).
The quality, however, depends on the success of motion estimation to retrieve the motion vector fields required for motion compensation. To illustrate how important the proposed sequential double gating and the artifact model-based motion estimation is, we prepared Fig. 11. If we apply a standard motion compensation, sMoCo with simultaneous double gating streaking artifacts partially remain. When using the sequential double gating approach for sMoCo, the results improve but the robustness against streaking artifacts is lower than for proposed acMoCo. These results apply not only for the shown respiratory and cardiac window but also for all others. Figure 12 shows axial slices of acMoCo from different cardiac phase bins for a fixed respiratory phase bin.
4. SUMMARY AND DISCUSSION We proposed a motion-compensated image reconstruction algorithm acMoCo using a sequential double gating approach combined with an artifact model-based motion estimation including a cyclic consistent nonrigid registration algorithm. This technique improves phase-correlated imaging of highly undersampled data. We demonstrated that dose savings of about an order of magnitude are possible compared to the conventional phase-correlated reconstruction approach. Both, sequential double gating and artifact model-based motion estimation, are essential parts of successful motion compensation for decreasing dose levels. Concurrent double gating increases the impact of sparse-view artifacts, such that a standard motion compensation sMoCo fails in this case. But also when using sequential double gating, sMoCo shows sensitivity to the remained streaking artifacts whose magnitude increases with decreasing dose levels. The artifact modelbased motion estimation including cyclic consistent nonrigid registration shows no sensitivity to these streaking artifacts even for cases with a one magnitude lower dose level. In addition, by the sequential double gating approach, higher Medical Physics, Vol. 42, No. 4, April 2015
sparsifications of the projection data are possible for acMoCo. While conventional phase-correlated reconstruction with its concurrent double gating fails because no projection is left for some desired phase window, by sequential double gating, the missing data can be borrowed from other respiratory phases. Furthermore, acMoCo outperforms the LDPC algorithm in terms of noise level and robustness with decreasing dose level. During the evaluation process, we used cardiac window width of 20% of the RR interval for image reconstruction. This is the minimum reasonable window size with respect to the detector integration time of 40 ms and the mean RR interval of about 200 ms. For accurate cardiac function estimation, e.g., ejection fraction, this window size could be insufficient as accuracy is affected by temporal blurring, in particular in systole which is the shortest cardiac phase. We want to emphasize that this limitation is only due to the detector integration time of the imaging system and not due to the image reconstruction approach. If shorter detector integration times are available and thus smaller cardiac window sizes are reasonable, a reconstruction algorithm has to deal with a reduced dose fraction per cardiac phase at constant total dose level. And we have demonstrated that acMoCo can produce reliable results with decreasing dose level. With the proposed motion compensation acMoCo, it is possible to perform 4D or 5D phase-correlated imaging at the same dose level as required for conventional 3D studies. As a consequence acMoCo allows for longitudinal in vivo studies of the rodent heart as well as for long term studies in preclinical research.
ACKNOWLEDGMENTS This study was supported by Varian Medical Systems, Palo Alto, CA and in parts by the Deutsche Forschungsgemeinschaft (DFG) under Grant No. FOR 661. The highspeed image reconstruction software RayConStruct–IR was provided by RayConStruct GmbH, Nürnberg, Germany. The
1958
Brehm et al.: 5D motion-compensated micro CT imaging
authors thank PD Dr. Andreas Hess and Sandra Strobelt for help with the mouse measurements. The authors thank Miltenyi BioTec GmbH, Bergisch Gladbach, Germany; nanoPET Pharma GmbH, Berlin, Germany; and Binitio Biomedical, Inc., Ottawa, Canada for providing their latest contrast agents. a)Author
to whom correspondence should be addressed. Electronic mail: [email protected]; Telephone: +49 (6221) 42 3067; Fax: +49 (6221) 42 2585. 1C. T. Badea, B. Fubara, L. W. Hedlund, and G. A. Johnson, “4–D micro–CT of the mouse heart,” Mol. Imaging 4, 110–116 (2005). 2M. Drangova, N. L. Ford, S. A. Detombe, A. R. Wheatley, and D. W. Holdsworth, “Fast retrospectively gated quantitative four–dimensional (4D) cardiac micro computed tomography imaging of free–breathing mice,” Invest. Radiol. 42, 85–94 (2007). 3S. Kriminski, M. Mitschke, S. Sorensen, N. M. Wink, P. E. Chow, S. Tenn, and T. D. Solberg, “Respiratory correlated cone–beam computed tomography on an isocentric C–arm,” Phys. Med. Biol. 50, 5263–5280 (2005). 4J.-J. Sonke, L. Zijp, P. Remeijer, and M. van Herk, “Respiratory correlated cone beam CT,” Med. Phys. 32, 1176–1186 (2005). 5G. C. McKinnon and R. Bates, “Towards imaging the beating heart usefully with a conventional CT scanner,” IEEE Trans. Biomed. Eng. BME-28, 123–127 (1981). 6K. L. Garden and R. A. Robb, “3–D reconstruction of the heart from few projections: A pratical implementation of the McKinnon–Bates algorithm,” IEEE Trans. Med. Imaging 5, 233–239 (1986). 7S. Leng, J. Zambelli, B. Nett, R. Tolakanahalli, P. Munro, J. Star-Lack, B. Paliwal, and G.-H. Chen, “Streaking artifacts reduction in four–dimensional cone–beam computed tomography,” Med. Phys. 35, 4649–4659 (2008). 8S. Sawall, F. Bergner, R. Lapp, M. Mronz, M. Karolczak, A. Hess, and M. Kachelrieß, “Low–dose cardio–respiratory phase–correlated cone–beam micro–CT of small animals,” Med. Phys. 38, 1416–1424 (2011). 9M. Persson, D. Bone, and H. Elmqvist, “Total variation norm for three– dimensional iterative reconstruction in limited view angle tomography,” Phys. Med. Biol. 46, 853–866 (2001). 10E. Y. Sidky and X. Pan, “Image reconstruction in circular cone–beam computed tomography by constrained, total–variation minimization,” Phys. Med. Biol. 53, 4777–4807 (2008). 11J. Song, Q. H. Liu, G. A. Johnson, and C. T. Badea, “Sparseness prior based iterative image reconstruction for retrospectively gated cardiac micro–CT,” Med. Phys. 34, 4476–4483 (2007). 12L. Ritschl, S. Sawall, M. Knaup, A. Hess, and M. Kachelrieß, “Iterative 4D cardiac micro–CT image reconstruction using an adaptive spatio–temporal sparsity prior,” Phys. Med. Biol. 57, 1517–1525 (2012).
Medical Physics, Vol. 42, No. 4, April 2015
1958 13C.
R. Crawford, K. King, C. Ritchie, and J. D. Godwin, “Respiratory compensation in projection imaging using a magnification and displacement model,” IEEE Trans. Med. Imaging 15, 327–332 (1996). 14C. Ritchie, C. Crawford, J. D. Godwin, K. King, and Y. Kim, “Correction of computed tomography motion artifacts using pixel–specific back–projection,” IEEE Trans. Med. Imaging 15, 333–342 (1996). 15T. Li, E. Schreibmann, Y. Yang, and L. Xing, “Motion correction for improved target localization with on–board cone–beam computed tomography,” Phys. Med. Biol. 51, 253–267 (2006). 16C. T. Badea, E. Schreibmann, and T. Fox, “A registration based approach for 4D cardiac micro–CT using combined prospective and retrospective gating,” Med. Phys. 35, 1170–1179 (2008). 17M. Brehm, P. Paysan, M. Oelhafen, and M. Kachelrieß, “Artifact–resistant motion estimation with a patient–specific artifact model for motion– compensated cone–beam CT,” Med. Phys. 40, 101913 (13pp.) (2013). 18L. Feldkamp, L. Davis, and J. Kress, “Practical cone–beam algorithm,” J. Opt. Soc. Am. 1, 612–619 (1984). 19M. Brehm, P. Paysan, M. Oehlhafen, P. Kunz, and M. Kachelrieß, “Self– adapting cyclic registration for motion–compensated cone–beam CT in image–guided radiation therapy,” Med. Phys. 39, 7603–7618 (2012). 20E. Meyer, R. Raupach, M. Lell, B. Schmidt, and M. Kachelrieß, “Normalized metal artifact reduction (NMAR) in computed tomography,” Med. Phys. 37, 5482–5493 (2010). 21M. Sun, M. Oelhafen, L. Carvalho, J. Pavkovich, T. Berkus, and J. Star-Lack, “Improving motion accuracy for the McKinnon–Bates 4D– CBCT reconstruction algorithm,” in Proccedings of 2nd International Meeting on Image Formation in X-ray CT (Salt Lake City, Utah, 2012), pp. 131–134. 22P.-S. Liao, T.-S. Chen, and P.-C. Chung, “A fast algorithm for multilevel thresholding,” J. Inf. Sci. Eng. 17, 713–727 (2001). 23J.-P. Thirion, “Image matching as a diffusion process: An analogy with Maxwell’s demons,” Med. Image Anal. 2, 243–260 (1998). 24H. Wang, L. Dong, J. O’Daniel, R. Mohan, A. S. Garden, K. K. Ang, D. A. Kuban, M. Bonnen, J. Y. Chang, and R. Cheung, “Validation of an accelerated ‘demons’ algorithm for deformable image registration in radiation therapy,” Phys. Med. Biol. 50, 2887–2905 (2005). 25T. Vercauteren, X. Pennec, A. Perchant, and N. Ayache, “Diffeomorphic demons: Efficient non–parametric image registration,” NeuroImage 45, 61–72 (2009). 26C. Tomasi and R. Manduchi, “Bilateral filtering for gray and color images,” in Proceedings of the Sixth International Conference on Computer Vision (IEEE Computer Society, Washington, 1998), pp. 839–846. 27M. Kachelrieß, D.-A. Sennst, W. Maxlmoser, and W. A. Kalender, “Kymogram detection and kymogram–correlated image reconstruction from subsecond spiral computed tomography scans of the heart,” Med. Phys. 29, 1489–1503 (2002).
