Image acquisition
A retrospective dataset consisting of 270 clinical perfusion CMR studies was analyzed. All studies were performed under procedures and protocols approved by the institutional review board of the National Heart, Lung and Blood Institute (NHLBI), and all subjects gave written informed consent (ClinicalTrials.gov Identifier: NCT00027170). Each study contained two perfusion imaging series, one at rest and one during vasodilator stress. Gadolinium-DTPA (Magnevist, Berlex Laboratories, Wayne, NJ, USA) was administered (0.05 mmol/kg) at 5 ml/s during stress and rest perfusion imaging followed by a saline flush.
Perfusion CMR was performed on a 1.5 T scanner (Siemens Healthcare, Erlangen, Germany) with a saturation recovery dual-sequence technique [28] at every R-R interval over 60 heart beats. A steady-state free precession (SSFP) sequence was used in 245 studies; while a fast low-angle shot (FLASH) sequence was used in 25 studies. Typical imaging parameters for the myocardial perfusion image series included: 90° composite saturation preparation pulse, 50° (SSFP) or 12° (FLASH) flip angle, 90 ms (SSFP) or 100 ms (FLASH) inversion time, 1.2 ms echo time, 2.3 ms repetition time, 8 mm slice thickness, 360 × 270 mm field of view, 128 × 80 acquisition matrix, 256 × 192 image matrix after interpolation, and parallel imaging factor of 2 [29]. For each perfusion imaging, a FLASH low-resolution dedicated AIF image series was also acquired with a separate saturation pulse and typical parameters of 8° flip angle, 5.0 ms inversion time, 0.7 ms echo time, 1.3 ms repetitive time, 10 mm slice thickness, and 64 × 48 acquisition and image matrix size. Examples of both the low-resolution dedicated AIF image series and the myocardial series are shown in Fig. 1. At the start of each perfusion acquisition, two proton density (PD) weighted images were also acquired with no saturation preparation pulse for surface coil intensity normalization.
Image processing algorithm
Figure 2 shows the image processing pipeline of the proposed automated algorithm to measure AIFs from the image series. Note that the same algorithm is applicable to either dedicated AIF or standard myocardial image series.
Motion correction
Displacement of the heart can occur during image acquisition due to the patient’s cardiac or respiratory motion. A non-rigid body image registration technique based on optical flow computations is used to correct for anatomical structure motion invariably present in the perfusion image series [30]. This motion correction is applied to each series before subsequent image processing.
Heart region detection
The location of the heart region, in the form of a bounding box, is determined in order to aid further processing. This is achieved by identifying regions most indicative of the heart ventricles. Candidate ventricle regions are dynamically thresholded from a pixel-wise standard deviation map of the time series images. This standard deviation map highlights pixels with large intensity changes, such as those caused by contrast agent perfusion, while removing regions that remain constantly bright, such as the chest wall. This map is thresholded at one standard deviation above the mean for the AIF series, and at two standard deviations above the mean for the myocardial series to obtain the candidate ventricle regions. After binarization, regions that do not match the temporal signal characteristics of the ventricles are identified and removed. Specifically, regions whose intensity increase is less than twice their baseline intensity, indicating minimal contrast enhancement, are removed. Regions whose peak intensity occurs within the first or last 3 frames of the time series are also removed. Next, a similarity check is performed to examine whether each region represents a unique ventricular candidate. Similar regions are grouped as a single ventricular region that has been split by papillary muscles, image artifacts or slice placement. Similar regions are identified as having average time-signal intensity curves with a cross-correlation coefficient of more than 0.75, and whose minimum Euclidean distance is less than the sum of each region’s average radius. The final regions are subject to a linear voting scheme to iteratively determine which two candidate regions are most characteristic of the RV and LV cavities based on their time-signal features. Features used for voting classification include: distance to the image center, distance to previously selected candidate regions, size of each region, signal intensity upslope, peak value (PV), time to peak (TTP), full width at half maximum (FWHM), and an M-value [13, 14] which combines the previous three features as shown in Eq. 1.
$$ M=\frac{PV}{\left(TTP\ast FWHM\right)} $$
(1)
For each feature, the candidate ventricles are ranked by how well they match typical ventricle characteristics; that is those with larger region size, PV, upslope, M-value, smaller TTP, FWHM, and shorter distances to image center and to the previously selected ventricles. The ranks are converted to a score of 1 to N, 1 being the lowest rank and N, the number of candidate regions, the highest. The feature scores for each region are totaled, and the region with the highest total score is selected as a ventricle region. The second ventricle is selected similarly. The two selected ventricle regions are used to create a bounding box around the heart for subsequent processing (see Fig. 2).
Intensity correction
Surface coil intensity correction was based on the corresponding PD weighted images to improve the signal intensity uniformity of the heart region in the perfusion series. A slightly modified and automated version of the surface coil intensity correction algorithm presented in [31] is applied to reduce such inhomogeneity. Based on a hierarchical weighting scheme for foreground and background regions, this method estimates a polynomial signal intensity surface profile from the PD images which is used to adjust the signal intensity in the image series. After the surface coil intensity correction, the image series is further adjusted to remove baseline intensity based on pre-contrast perfusion images.
Ventricular pixel detection
To further the detection of ventricular blood pool pixels in the images, an ICA algorithm [27] is first used to obtain representative time-signal intensity curves from the previously identified ventricular regions. Assuming a mixture of two independent sources of signal (the RV and LV) in all ventricular regions, ICA separates and extracts the two primary time-signals that represent the dynamic contrast of the two ventricles. All of the pixels in the bounding box are classified to the RV, LV or background regions by computing their cross correlation to the RV and LV time-signals after the ICA process. Pixels with a cross correlation greater than a statistically determined value of 0.7 are assigned to the matching ventricle; the remaining pixels are then classified as the background region. The RV is identified as being the first region to reach peak intensity, which is followed by the LV region (see Figs. 2 and 3).
AIF extraction
The final step of the algorithm is to select LV pixels that are brighter than a default threshold to compute an average intensity value of the blood pool in each image. This step excludes any papillary muscle pixels that may have been included in the LV region from the previous steps because they do not receive any contrast agent and remain dark. This step also excludes pixels that contain potential partial volume errors, as these pixels will also be darker than the average LV pixel. This closely replicates manual analysis, where the LV cavity is a relatively small but bright region within the heart. The default threshold was computed from the maximal intensity projection image as the 75th percentile of the maximal intensity range of the LV region. Example results are shown in Figs. 2 and 4, where the red pixels indicate the pixels selected for the AIF measurement.
Finally, the AIF curve is linearly re-sampled at every half second to convert the time unit from image frames to seconds, since the perfusion imaging is performed on every R-R interval. We refer to this re-sampled curve as the AIF time-signal intensity curve for the remainder of the paper and use it in our statistical comparisons.
AIF timing point calculation
In order to calculate important time-signal features for quantitative perfusion analysis as well as for candidate ventricle region selection, three contrast enhancement time points: baseline, start of, and peak contrast enhancement are derived from the AIF time-signal intensity curve automatically. First, the peak time is detected simply from the peak value of the time-signal intensity curve. The baseline time is determined next as the point of the curve with minimal intensity variation with its neighbors (the immediately adjacent points) between the beginning of the series and the rising peak (indicated by the point of maximal intensity change before the peak time). Finally the start time is detected by fitting a line to the rising peak and selecting the point of the curve geometrically closest to the intersection of this fitted line and the baseline intensity. As an example, automatically detected AIF timing points are shown in Fig. 5.
Quantitative evaluation
Since many of the measurements made were not normally distributed, the median and 95 % confidence interval (CI) of the median were calculated throughout the analysis using SPSS Statistics software (International Business Machines Corp, Armonk, New York). A non-parametric Mann-Whitney U test [32] was used to determine if there were significant differences between the automated and manual AIF measurements, with p < 0.05 considered statistically significant.
AIF contrast enchantment characteristics comparisons
The performance of the proposed automated method was evaluated on both rest and stress image series from all clinical perfusion studies using custom image analysis software developed in Interactive Data Language (IDL, Exelis Visual Information Solutions, Boulder, Colorado). The results of the automated AIF detection were compared to a reference created by manually drawing an ROI in the LV blood cavity throughout the entire dataset. Pearson’s correlation coefficient and normalized root mean square error (RMSE) were used to evaluate the similarity between the automated and manual AIF time-signal intensity curves. Descriptive time-signal statistics of the AIF were compared which included TTP, FWHM, PV, signal intensity upslope, and M-value, as described in the previous sections. As shown in Fig. 6, accurate AIFs are characterized by high values of PV, upslope, M-value, and low values of TTP and FWHM [13, 14]. Such improvements will likely be due to selection of the brightest LV pixels, exclusion of papillary muscles, and exclusion of pixels containing partial volume errors.
AIF timing point and myocardial blood flow comparisons
Further evaluation was performed on a subset of 21 clinical studies consisting of dedicated AIF image series of normal healthy volunteers from the SSFP dataset. First, the automatically detected first-pass contrast timing point, the start of peak and peak time, from the automated AIF curve were compared against manually selected ones. Second, fully quantitative MBF estimates using the automatically and manually measured AIFs from the dedicated AIF image series were compared. To facilitate these comparisons, the endocardial and epicardial borders of the myocardium were manually traced by a trained CMR expert onto the myocardial image series using the custom software. The defined myocardial regions were subdivided into six equiangular sectors to derive myocardial time-signal intensity curves. The MBF calculation was based on a model constrained deconvolution technique [33]. To observe only the effect of the AIF upon the MBF, the manually selected timing points were used in all MBF calculations.
Next, MBF estimates were evaluated again with differently sized AIF ROI selections. Here we generated three differently sized ROIs by selecting three different percentile thresholds in the AIF extraction part of the automated algorithm (see “AIF extraction” section): the default ROI using the 75th percentile, a larger ROI using the 50th percentile, and the largest ROI using the 25th percentile. MBF estimates were calculated using each of these AIFs as described previously.
Finally, the execution time required by the automated algorithm to extract the AIF was measured on these 21 studies using an Intel Core i7-3770 3.4GHz central processing unit (CPU). For comparison, the time required to manually perform the equivalent processing of the same input image series to measure the AIF using our custom software was also measured. The execution time for both automatic and manual processing included the time required for image loading, motion correction or registration, LV ROI selection, and AIF signal intensity calculation from the ROI. It does not include the myocardial contouring time. It should be noted that the custom image analysis software was specifically designed to facilitate quick and efficient manual ROI drawing and registration.