Effect of data conserving respiratory motion compensation on left ventricular functional parameters assessed in gated myocardial perfusion SPECT

Background Respiratory motion compromises image quality in myocardial perfusion (MP) single-photon emission computed tomography (SPECT) imaging and may affect analysis of left ventricular (LV) functional parameters, including phase analysis-quantified mechanical dyssynchrony parameters. In this paper, we investigate the performance of two algorithms, respiratory blur modeling (RBM) and joint motion-compensated (JMC) ordered-subsets expectation maximization (OSEM), and the effects of motion compensation on cardiac-gated MP-SPECT studies. Methods Image acquisitions were carried out with a dual-detector SPECT/CT system in list-mode format. A cardiac phantom was imaged as stationary and under respiratory motion. The images were reconstructed with OSEM, RBM-OSEM, and JMC-OSEM algorithms, and compared in terms of mean squared error (MSE). Subsequently, MP-SPECT data of 19 patients were binned into dual-gated (respiratory and cardiac gating) projection images. The images of the patients were analyzed with Quantitative Gated SPECT (QGS) 2012 program (Cedars-Sinai Medical Center, USA). The parameters of interest were LV volumes, ejection fraction, wall motion, wall thickening, phase analysis, and perfusion parameters. Results In phantom experiment, compared to the stationary OSEM reconstruction, the MSE values for OSEM, RBM-OSEM, and JMC-OSEM were 8.5406·10−5,2.7190·10−5, and 2.0795·10−5, respectively. In the analysis of LV function, use of JMC had a small but statistically significant (p < 0.05) effect on several parameters: it increased LV volumes and standard deviation of phase angle histogram, and it decreased ejection fraction, global wall motion, and lateral, septal, and apical perfusion. Conclusions Compared to standard OSEM algorithm, RBM-OSEM and JMC-OSEM both improve image quality under motion. Motion compensation has a minor effect on LV functional parameters. Supplementary Information The online version contains supplementary material available at (10.1186/s40658-021-00355-w).


Background
Myocardial perfusion (MP) single-photon emission computed tomography (SPECT) examination is often performed in the diagnosis of coronary artery disease [1]. Cardiac gating of MP-SPECT examination further enables objective assessment of left ventricular (LV) functional parameters [2], which provide additional information in, for example, diagnosis of myocardial infarction [3], identification of multiple-vessel disease [4], prediction of mortality [5], and identification of patients with risk for LV remodeling [6]. In addition, evaluation of LV mechanical dyssynchrony with phase analysis [7] may help to interpret whether the patient would benefit from cardiac resynchronization therapy [8,9].
During the MP-SPECT examination, the heart undergoes quasiperiodic motion due to respiration; this movement can even exceed 20 mm [10]. As a result, the reconstructed SPECT images appear increasingly "blurred" if the motion is not accounted for in the reconstruction. Fortunately, several motion compensation methods have been introduced over the years [11][12][13][14][15]; however, as far as we know, there are only few papers where the effect of respiratory motion compensation on analysis of LV function has been investigated [16,17]. In addition, only one of these studies investigated the effect of motion compensation on phase analysis results [17]. What is more, the motion compensation method applied in reference [17] was essentially a respiratory gating method, which means that the resulting images had amplified noise levels compared to the concurrent clinical standards.
The aim of this study is to investigate the effect of respiratory motion compensation on results obtained from cardiac-gated MP-SPECT examinations. First, a cardiac phantom is imaged to validate two promising motion compensation strategies. Subsequently, these motion compensation methods are applied on clinical cardiac-gated MP-SPECT data. Finally, the clinical images reconstructed with the better motion compensation method are analyzed in terms of LV functional parameters and compared to standard clinical reconstructions. Specifically, the effect of motion compensation on phase analysis parameters is investigated.

Data acquisition and preprocessing
All image acquisition procedures were carried out on a dual-detector SPECT/CT system (Precedence; Royal Philips N.V., Netherlands) at Kuopio University Hospital. All SPECT data were acquired in list-mode format and processed in MATLAB environment (MATLAB R2015b; The MathWorks, Inc., MA, USA).

Phantom experiment
A custom-designed phantom, which consists of a cardiac phantom and a "torso" phantom, was utilized in this part of the study [18]. The cardiac phantom was filled with 121 MBq of Tc-99m and placed inside the "torso" phantom. The "torso" phantom was then placed on a plastic plate on top of the patient table of the SPECT/CT system and attached to a respiratory phantom (AZ-733V; Anzai Medical Co.,Ltd., Japan). The respiratory phantom provided a sinusoidal motion pattern with 10 RPM frequency and 20 mm peak-topeak amplitude. The motion was directed parallel to the rotation axis of the SPECT/CT system. Before SPECT imaging, the phantom was scanned with CT. The following standard clinical settings were used: helical scan mode, tube voltage of 140 kV, tube currentexposure time product of 20 mAs, matrix size of 512 ×512, slice thickness of 5 mm, and pixel size of 1.171875 mm. The phantom was stationary during the CT scan.
The phantom was imaged twice: first with respiratory motion and then as stationary (reference scan). Standard clinical myocardial perfusion imaging acquisition protocol of Kuopio University Hospital was used. This includes 90°detector configuration, noncircular detector orbit with 64 steps from right anterior oblique to left posterior oblique angle, low-energy high-resolution collimators, energy window of 140 keV ± 10%, and acquisition time of 30 s per step.
To acquire respiratory signal for gating, the method of Sakaguchi et al. was applied [19]. Briefly, the list-mode data were binned into 50-ms frames from which the center of activity (COA) in the axial direction (z-axis) was computed. Since it was known that the phantom motion is sinusoidal, a sine wave was fit to the COA data using a Gauss-Newton scheme. This sine wave thus represented the phantom's axial position through the entire image acquisition.
In order to induce a more realistic respiratory pattern for the phantom study, distribution of acquisition time among the 64 projection angles and 10 equal-sized respiratory amplitude windows was assessed from one clinical patient data. This distribution was used to scale the phantom's acquisition time correspondingly when the list-mode data was binned into respiratory-gated projection images. Total acquisition time per projection angle applied in the binning was 4.542 s, the matrix size was 128 ×128, and the pixel size was 4.664 mm. Since the stationary scan was acquired after the moving scan, to account for decay, the acquisition time per projection angle applied in the binning was increased to 4.812 s.

Patient studies
The study population consisted of 19 patients (8 female) referred to 1-day clinical stress/rest MP-SPECT examinations. Their characteristics (mean ± standard deviation) were as follows: age, 69 ± 11 years; weight, 76 ± 13 kg; height, 169 ± 12 cm; and body mass index, 26.4 ± 2.9 kg/m 2 . The patients received 300 MBq of Tc-99m tetrofosmin prior to the stress imaging and 706 ± 11 MBq prior to the rest imaging 3 h later. For the purposes of this study, data were collected at rest only. Written informed consent was obtained from every participating patient, and the study was approved by the Research Ethics Committee of the Northern Savo Hospital District (Dno 90/2011; March 20, 2012).
Image acquisition was started 45 min after the rest injection. The patients' electrocardiogram (ECG) was monitored with a cardiac trigger monitor (CTM 3000; Ivy Biomedical Systems, Inc., USA). The patients were imaged in a supine position at each step until 30 s of R-R intervals within a 40%-wide ECG gating acceptance window had been acquired. Otherwise, the image acquisition protocol was the same as in the phantom experiment.
ECG and thoracic electrical bioimpedance (EBI) of each patient were measured with a data acquisition system equipped with amplifiers (MP150, ECG100C, and EBI100C; BIOPAC Systems, Inc., USA). Since the data acquisition system had to be started earlier than the SPECT/CT system, there was a delay t between the recorded signals and the list-mode data. To synchronize the signals with the list-mode data, the R-triggers from the cardiac trigger monitor were transferred simultaneously to the SPECT/CT system and the data acquisition system via an in-house built pulse generator. The delay t was found by maximizing the cross-correlation of these two resulting R-trigger sequences, after which the ECG and EBI signals could be shifted in time to match the list-mode data.
For cardiac gating, R waves were detected from the ECG signal using Pan-Tompkins method [20] and the R-R intervals were divided to 8 cardiac bins using forward gating approach with fixed bin length [21]. The EBI signal was low-pass filtered after which nonphysiological trends were removed with a moving average filter [17]. The filtered EBI signal was treated as the respiratory signal which was used to compute respiratory gating limits according to Dey et al. [22]. Seven equal-sized amplitude windows were used because for the majority of patients the measured respiratory motion magnitude was less than 14 mm, and a previous study suggested that respiratory blur of an individual respiratory window should be less than 2 mm [18]. Thus, ECG and EBI signals were used to bin list-mode data into a total of 3584 dual-gated projection images (8 cardiac bins ×7 respiratory windows × 64 projection angles) with a pixel size of 6.219 mm.

Reconstructions
We present the reconstruction-related equations for the dual-gated data, noting that the phantom data are also dual-gated data-with a single cardiac bin.

Parameters
For the sake of comprehensibility, we first present all the parameters required in the following sections. Let C, R, and L denote the number of cardiac bins, respiratory windows, and projection angles, respectively. Let M and N denote the number of pixels in a projection image and number of voxels in the unknown/estimated image, respectively. Let τ c,r denote the image acquisition time at cardiac bin c, respiratory window r, and projection angle . Let g c,r ∈ R M×1 , c = 1, . . . , C, r = 1, . . . , R, = 1, . . . , L denote the dual-gated projection image at cardiac bin c, respiratory window r, and projection angle . Let denote the unknown, true image at cardiac bin c and let denote the observation matrix block at cardiac bin c, respiratory window r, and projection angle . Let g c,r denote a concatenated vector such that and let K c,r denote a concatenated matrix such that Let P ∈ R M×N , = 1, . . . , L denote the forward projection matrix at projection angle and let R ∈ R N×N , = 1, . . . , L denote a rotation matrix at projection angle . Finally, let M r →r ∈ R N×N denote a transformation matrix from a reference respiratory window r to respiratory window r and let 1 ∈ R ML×1 denote a vector of ones.

Non-motion-compensated reconstruction
For reference, it is assumed that no respiratory motion occurs during the image acquisition and the data is not respiratory-gated. In this case, neglecting attenuation, scattered photons, and measurement errors, our imaging model is where the observation matrix blocks are [23] K c,r = τ c,r P R .
Imagesf c can now be reconstructed with maximum likelihood expectation maximization (MLEM) algorithm [24] where superscripts (new) and (old) denote new and old estimates off c , respectively.

Motion-compensated reconstruction
Respiratory blur modeling The first motion compensation strategy that we apply has previously been used in positron emission tomography [25,26]. It bases on an assumption that the data is motion-blurred (not respiratory-gated) and this is taken into account in the imaging model; therefore, we call it respiratory blur modeling (RBM). In practice, this means that imaging model (1) applies, and the observation matrix blocks are modified such that where B c is a blurring matrix such that The reconstruction is carried out with Eq. (2) using the modified observation matrix blocks (Eq. (3)).

Joint motion-compensated reconstruction
The second motion compensation strategy that we apply was presented by Qi et al. [13] and is called joint motion-compensated (JMC) reconstruction. This method bases on a common gating assumption that the object to be imaged moves instantaneously from one position to another, but stays still while occupying a position. This means that the imaging model is where the observation matrix blocks are Derivation of JMC-MLEM algorithm follows the same path as that of the original MLEM algorithm [24]. That is, the solution is searched by formulating the likelihood function p X c |f c for the so-called complete data tensor X c where X c i, j, r is the number of photons emitted by image voxel j and collected at projection image pixel i at respiratory window r and cardiac bin c. Taking the conditional expectation of log-likelihood with respect to data vector g c,r and current estimate of f c (expectation step), and setting its derivatives with respect to the unknown values of f c (j) to zero (maximization step), one has derived the JMC-MLEM algorithm. In matrix form, this becomeŝ

Implementation
Projection and rotation matrices Projection matrices P were implemented such that they also model collimator-detector response as a distance-dependent Gaussian function [27]. Rotation matrices R were implemented using Gaussian interpolation [28]. Acquisition times τ c,r were calculated during list-mode data processing.

Transformation matrices
To estimate the motion between respiratory windows r and r , a sequence of respiratory-gated images was reconstructed based on data vectors C c=1 g c,r using the standard MLEM algorithm. We ran the algorithm for 15 iterations and applied post-reconstruction Gaussian smoothing with a standard deviation of 1 voxel. We chose r = 0.5R as the reference respiratory window from which the myocardium-containing region of interest (ROI) was segmented using an approach similar to the one introduced by Germano et al. [29]. Voxels outside the ROI were set to zero, and the other respiratory windows were registered with the reference window by assuming rigid-body translation and using the sum of squared differences as the term to be minimized with Gauss-Newton optimization [30]. The motivation for segmenting myocardium from the reference window was to focus the image registration on the myocardial region because respiratory motion of the heart may be different than that of extracardiac structures [22]. After computing the translation vectors, their directions were reversed to finally form the matrices M r →r applying linear interpolation.

Reconstruction details
In both phantom and patient studies, the reconstructions with Eqs. (2) and (6) were accelerated by using the ordered subsets variant of the EM algorithm (OSEM) [31]. The data were arranged to 8 subsets, and the algorithm was run for 10 iterations. In the phantom experiment, the stationary reference data were reconstructed with standard OSEM algorithm, and the moving data were reconstructed with OSEM, RBM-OSEM, and JMC-OSEM algorithms. The computation times for each method were recorded. In the patient studies, the data were reconstructed with OSEM, RBM-OSEM, and JMC-OSEM algorithms.
In patient studies, the reconstructed transaxial images were manually rotated to shortaxis images using Gaussian interpolation. Extracardiac activity was masked, and threedimensional post-reconstruction Gaussian filtering with full width at half maximum of 13.3 mm was applied. In phantom experiment, no rotation was necessary as the coronal plane was parallel to the short-axis plane.

Image analysis: phantom experiment
In order to quantitatively study the visibility of cardiac defects, we computed the contrast between each defect region and a corresponding normal myocardial region according to the study of Kortelainen et al. [18]. Briefly, the CT image of the stationary phantom was converted to μ-map using piecewise linear conversion [32]. This μ-map was resampled to SPECT matrix size and coregistered with SPECT image applying normalized gradient fields method [33]. This μ-map was then converted back to Hounsfield units (HUs). This was done to identify the defect regions, which had higher HU values than the surrounding water. An HU value of 50 was used as the threshold to binarize the converted μ-map, which permitted the selection of defect-containing regions. The normal myocardial regions (references) were determined as regions of similar shape and size, located adjacent to the defect regions (Fig. 1).
Having computed the coordinates of defect and reference regions, OSEM, RBM-OSEM, and JMC-OSEM reconstructed moving images were registered with the stationary image and post-reconstruction Gaussian filtering with a standard deviation of 1.0 was applied. For each image, contrast for each defect was computed as [18] where A reference and A defect are average voxel values on reference ROI and defect ROI, respectively. In addition, mean squared error (MSE) of each of the moving images was computed, considering the stationary image as reference.

Image analysis: patient studies
Based on the results from the phantom experiment, we opted to analyze only the OSEM and JMC-OSEM reconstructed patient images, excluding the RBM-OSEM images. The analysis was performed with Quantitative Gated SPECT (QGS) 2012 program (Cedars-Sinai Medical Center, USA) [34]. Studied parameters included LV end-diastolic volume (EDV), end-systolic volume (ESV), stroke volume (SV), ejection fraction (EF), and wall motion (WM) and wall thickening (WT) assessed in anterior, lateral, inferior, septal, and apical wall segments. In addition, we computed global wall motion (WM GLB ) and global wall thickening (WT GLB ) as average values over these five segments.
In phase analysis, QGS fits a first Fourier harmonic curve to the time-activity data of each myocardial sampling point and records amplitude and phase of each curve [35]. Before progressing, the program discards 5% of the phase angles whose corresponding amplitude values are the lowest. The remaining phase angles are binned into a histogram using 360 bins, and histogram bandwidth (BW), standard deviation (StD), and entropy (ENT) are computed to describe the phase angle distribution. In order to evaluate the effect of motion compensation on myocardial perfusion, we applied the "Freeze" command of QGS. This command warps the myocardial activity from all cardiac bins c to the end-diastolic bin through landmark-based image registration and sums them together, thus eliminating cardiac contractile motion [36]. Perfusion of these "motion-frozen" images was assessed using the same five-segment polar map as in the assessment of wall motion and wall thickening.
Statistical significance of differences between motion-compensated and non-motioncompensated data was evaluated with SPSS Statistics 23 (IBM Corporation, USA). If the data were normally distributed based on the Shapiro-Wilk test (p > 0.05), we used pairedsamples t test for comparisons. If the data were not normally distributed (p < 0.05), we used Wilcoxon signed-rank test for comparisons.

Phantom experiment
The computation times for each method are given in Table 1. Figures 2 and 3 present images of the phantom as stationary and reconstructed with OSEM, moving and reconstructed with OSEM, moving and reconstructed with RBM-OSEM, and moving and reconstructed with JMC-OSEM. Contrasts for individual defects and the MSE values of images are given in Table 2. While RBM-OSEM successfully reduced respiratory blur compared to the standard OSEM as seen from the lower MSE value, it could only partially recover the visibility of cardiac defects. The visibility of defects was further improved with JMC-OSEM.

Patient studies
Respiratory motion magnitude observed in each patient study was computed as Euclidean distance between respiratory windows 1 and 7. For this study population, the (mean ± standard deviation) respiratory motion magnitude was 10.9 ± 4.9 mm. Figures 4, 5, and 6 present "motion-frozen" images of three patient cases, with the images reconstructed with OSEM, RBM-OSEM, and JMC-OSEM algorithms. Motion compensation either retained or improved image quality of the "motion-frozen" images, depending on the measured respiratory motion magnitude. The improvements could be seen as higher contrast between the myocardium and LV cavity and better visibility of hypoperfused regions. The best visibility of hypoperfused regions was achieved with JMC-OSEM.
All QGS analysis results are given in Table 3. Use of JMC increased EDV and ESV and decreased EF compared to standard reconstructions (p < 0.001) and left SV unaffected (p > 0.05). Use of JMC also decreased wall motion in inferior (p < 0.05) and lateral (p < 0.01) walls, and also globally (p < 0.05). In addition, it decreased wall thickening in septal and inferior walls (p < 0.05). Use of JMC increased phase analysis parameter StD (p < 0.05). Use of JMC also decreased myocardial perfusion in lateral and septal walls (p < 0.01) and in apex (p < 0.05), but left anterior and inferior walls unaffected (p > 0.05).

Discussion
In this work, we studied the performance of two respiratory motion compensation strategies. Both methods, RBM-OSEM and JMC-OSEM, reduced respiratory motion blur and improved visibility of hypoperfused regions compared to the standard OSEM algorithm; however, JMC-OSEM proved to be superior to RBM-OSEM. Thus, we chose the JMC-OSEM algorithm to study the effect of respiratory motion compensation on quantification of LV functional parameters.
These motion compensation algorithms are just a few of several methods that have been presented over the years, and they were chosen for the present work because they fit well in the MLEM framework which is currently the recommended way to  reconstruct myocardial perfusion images [37]. Another appropriate method is the popular reconstruct-transform-average (RTA) algorithm [12]. The reason why RTA was not included in this study was because we concluded earlier using phantom data that, in terms of MSE, RBM-OSEM outperformed it in case where the patient's breathing pattern included a baseline shift (Additional file 1: Supplementary material). It was not obvious for us whether RBM-OSEM could outperform JMC-OSEM, so we included both algorithms in this study. While JMC-OSEM was superior to RBM-OSEM in terms of image quality (MSE), one must decide whether it is worth the computational cost; as one can observe from Eq. (6), each iteration of JMC requires R forward projections and 2R backprojections, whereas RBM requires just 1 forward projection and 2 backprojections, as shown in Eq. (2). Appropriately, in this work, reconstruction of phantom data of 10 respiratory windows with JMC-OSEM lasted nearly 7 times longer than reconstruction with RBM-OSEM (Table 1). This apparent discrepancy of 7 times longer computation time instead of 10 is related to multiplication with sparse matrices M r →r in JMC-OSEM and B c in RBM-OSEM. B c is a weighted average of R different M r →r matrices and thus contains more non-zero elements than a single M r →r matrix, which increases the computation time of RBM-OSEM algorithm.
Respiratory motion has been demonstrated to decrease the apparent myocardial activity especially in the anterior and inferior walls in several studies [11,18,38,39], and the same observation was made in the present phantom experiment (Figs. 2 and 3). However, we did not obtain similar results from the patient studies (Table 3)-some JMC reconstructions displayed higher perfusion values in anterior and inferior walls, but overall, the differences were statistically non-significant. One explanation could be the way QGS presents the perfusion results: they are scaled to the maximal myocardial count value. Unfortunately, we cannot obtain absolute "motion-frozen" perfusion values for comparisons from QGS. On the other hand, the observed perfusion decreases in lateral, septal, and apical walls are concordant with the results of Ko et al. [39]. Use of JMC increased LV EDV and ESV, as could be expected based on previous studies [17,39,40]. In addition, the decrease of global wall motion as a result of JMC is concordant with the simulation results of Bitarafan-Rajabi et al. [40], although one could have expected a larger decrease. Although the change was statistically non-significant, we observed a small decrease in global wall thickening, which was also suggested by Bitarafan-Rajabi et al. [40].
Interestingly, we found that JMC caused a small but statistically significant increase in phase analysis parameter StD. This finding indicates that respiratory motion is a mechanism that partially conceals mechanical dyssynchrony. As far as we know, this finding has not been reported before. In the previous study by Kortelainen et al., no differences in phase analysis parameters were found between respiratory-gated and non-respiratorygated images [17]. This might be due to the fact that the authors used only 15-s-long time window in projection image binning, as opposed to the present study where we used all the acquired emission photons that occurred during accepted R-R intervals (∼ 30 s). Therefore, the effect of their respiratory motion compensation method may have been overshadowed by the effect of increased noise. This study had some limitations. First, the number of patients was low. Unfortunately, the camera used in the current study was upgraded to a new model which prevented acquisition of a larger study population. Second, in our current method, we estimate the respiratory motion before the final reconstruction by pre-reconstructing the respiratory windows. However, it was demonstrated by Song et al. that simultaneous reconstruction and motion estimation could be more effective in compensation of respiratory motion [41]. Third, besides the phantom experiment, we are unable to provide a reference "ground truth" data for the patient studies. However, as pointed out in the previous paragraphs, a part of our results are concordant with the ones obtained in previous studies, which enhances the credibility of our work. Fourth, the applied motion model (translation only) is simplistic to accurately describe real cardiac respiratory motion. However, since the pre-reconstructions of respiratory windows are likely to suffer from limitedangle errors due to irregular breathing [22], we opted for the most robust motion model (only three parameters to be estimated).
The degree of visual improvement obtained using respiratory motion compensation likely depends on the respiratory motion magnitude. Previously, Pretorius et al. established linear regression models for segmental count differences as a function of respiratory motion in a large study population (n = 1103) [42]. The larger the respiratory motion was, the larger were the count differences between respiratory motion-compensated and non-compensated reconstructions [42]. Although our study population is too small to definitely establish such trends, our results agree with those of Pretorius et al. For patients whose respiratory motion magnitude was smaller than 5 mm, the visual changes due to motion compensation were virtually non-existent. For the three patients presented in Figs. 4, 5, and 6, use of RBM and JMC improved the visibility of hypoperfused myocardial regions and the contrast between myocardium and LV cavity. These promising results warrant further studies and suggest that respiratory motion compensation could be performed routinely in near future.

Conclusions
We studied two respiratory motion compensation expectation maximization algorithms, RBM-OSEM and JMC-OSEM, in the context of myocardial perfusion SPECT imaging. The algorithms were validated with a phantom experiment. The JMC-OSEM algorithm was found to be superior in terms of image quality and was thus subsequently chosen to reconstruct clinical cardiac-gated myocardial perfusion images. Use of motion compensation slightly increased phase analysis-quantified mechanical dyssynchrony. It also improved image quality and had minor effects on left ventricular volumes, wall motion, and wall thickening.