Study population
The proposed method was evaluated in a group of 21 subjects composed of 11 healthy volunteers (10 females, 1 male) with no history of prior or current cardiovascular disease or cardiac medication, mean age 67.3±3.6 years, range 60–72; and 10 patients (2 females, 8 males) with heart failure of different etiologies (ischemic cardiomyopathy and idiopathic dilated cardiomyopathy), mean age 61.1±13.4 years, range 32–78.
The patients were enrolled from the Department of Cardiology, Linköping University, Sweden. Exclusion criteria for the patients were: significant ventricular arrhythmia, heart rate below 40 bpm or over 100 bpm, cardiovascular shunt, or more than mild to moderate valvular disorder. The research was performed in line with the declaration of Helsinki and was approved by the Linköping ethics board, reference number 2010/273-31. All subjects gave written informed consent.
CMR examinations
CMR examinations were performed on a clinical 3T Philips Ingenia scanner (Philips Healthcare, Best, the Netherlands). All subjects received a gadolinium contrast agent (Magnevist, Bayer Schering Pharma AG) prior to the acquisition for a late gadolinium enhancement (LGE) study.
4D flow CMR datasets were acquired during free-breathing, using a navigator gated gradient-echo pulse sequence with interleaved three-directional flow-encoding and retrospective vector cardiogram controlled cardiac gating. Scan parameters included: Candy cane view adjusted to cover both ventricles, velocity encoding (VENC) 120 cm/s, flip angle 10°, echo time 2.6 ms, repetition time 4.4 ms, parallel imaging (SENSE) speed up factor 3 (AP direction), k-space segmentation factor 3, acquired temporal resolution of 52.8 ms, spatial resolution 2.7 × 2.7 × 2.8 mm3, elliptical k-space acquisition, scan time: 7–8 min excluding and 10–15 min including the navigator efficiency at heart rate 60 bpm.
Two-dimensional cine through-plane phase-contrast-CMR velocity data were acquired in a slice perpendicular to the main flow direction in the ascending aorta just downstream from the aortic valve, above the coronary arteries. The velocity data were acquired in a breath hold (duration approximately 10–15 sec) and the acquisition was retrospectively gated to the ECG using the following settings: slice thickness 8 mm, field of view (FOV) 350 mm × 300 mm, sensitivity encoding (SENSE) factor 2, velocity encoding range (VENC) 200 cm/s, echo time (TE) 2.2 ms, repetition time (TR) 3.8 ms, flip angle 10°. Five lines of k-space were acquired per heartbeat, resulting in a temporal resolution of approximately 20 ms.
The 4D flow CMR data were corrected for concomitant gradient fields on the scanner. Offline processing corrected for phase wraps using a temporal phase unwrapping method and background phase errors using a weighted 2nd order polynomial fit to static tissue [23, 24]. The 2D velocity data were corrected for concomitant gradient fields and background phase errors on the scanner. Following these processing steps, all datasets were screened for large residual background velocity offset errors that would affect the evaluation.
Atlas-based segmentation and flow analysis
The proposed method utilizes atlas-based segmentation and is schematically illustrated in Fig. 1. Atlas-based segmentation enables automatic segmentation of specific regions-of-interest by registering an image volume in which these regions have been delineated i.e. an atlas, to an unsegmented image volume (hereafter termed input dataset). The deformed atlas will constitute a segmentation mask in the input dataset.
The atlas was created from one dataset of a healthy volunteer. The aorta, pulmonary artery, and vena cava were segmented semi-automatically. Additionally, the atlas includes the location of 2D planes to be used for flow analysis. The remaining datasets where analyzed using this atlas.
The method can be summarized as follows:
-
1.
Atlas-based segmentation: The atlas is fitted to the input dataset by means of registration.
-
2.
4D vessel segmentation: Segmentation of the great thoracic vessels in 3D over the entire cardiac cycle.
-
3.
Flow volume assessment: Extraction of 2D vessel planes and net flow volume calculation in each plane.
The method was implemented in MATLAB (The MathWorks, Inc., Natick, Massachusetts, United States), using a registration toolkit for MATLAB that relies on non-parametric methods for non-rigid registration [25]. The analyses were performed on an HP Z820 Workstation with an Intel Xeon 6 Core 2.5 GHz processor and 64 GB RAM. Each step of the approach is described in more detail below.
Atlas creation
An atlas representing the great thoracic vessels was constructed by delineating each major vessel in one dataset belonging to a healthy subject. The segmentation needed to create the binary masks was performed in a systolic 3D phase-contrast Magnetic Resonance Angiography (PC-MRA) computed from 4D flow MRI data according to (1).
$$ \text{PC\_MRA} = M \sqrt{{V_{x}^{2}} + {V_{y}^{2}} + {V_{z}^{2}}} $$
((1))
where V
x
, V
y
and V
z
are the averages of the components of the three directional blood flow velocity over the systolic time frames, and M is the average magnitude of the signals acquired over the systolic time frames, which acts as a noise suppressor in areas of very low signal such as the lungs. The resulting PC-MRA is exemplified in Fig. 2a.
Additionally, 2D planes were defined at multiple locations of interest in the atlas, in this case: perpendicular to the ascending aorta, proximal descending aorta, distal descending aorta, pulmonary trunk, left and right branches of the pulmonary artery, inferior and superior vena cava. These planes were added in order to calculate the net flow volume passing through them over a heartbeat; and they also serve as a non-subjective, numerical way of evaluating the segmentation method.
The atlas creation process is demonstrated in Fig. 2. Note that the atlas only needs to be generated once, and it can be used subsequently to segment an arbitrary number of input datasets. Furthermore, the specific locations of the 2D planes described were selected for evaluation purposes; however, the method allows for any number of planes at any location close to the major vessels.
Atlas-based segmentation
For each input dataset, a 3D PC-MRA was calculated from 4D flow CMR data. The atlas’ PC-MRA was registered to each input dataset’s PC-MRA using affine registration followed by non-rigid registration. In this way, the robustness of the affine registration was used to obtain an initial rough deformation, while the accuracy of non-rigid registration was used for fine-tuning.
Both registration methods were performed using three scales, five iterations per scale, and linear interpolation. The non-rigid registration method used during this study was the Morphon algorithm, which uses local displacement estimations iteratively in order to fit a source image to a target image [26–28].
Diffeomorphic field accumulation was used in all registration instances, since it enables compression and deformation of the data, while avoiding tearing or folding. Both fluid and elastic regularization of the displacement field were also used in order to limit the amount of deformation permitted during the registration. Elastic regularization limits the deformation to model an elastic material, while fluid regularization models the deformation as a viscous fluid. A representative example of the appearance of the images before and after the registration process is shown in Fig. 3.
4D vessel segmentation
After registering the atlas to the input 4D flow CMR dataset at systole, the atlas was adapted to all other time frames in order to account for the motion of cardiovascular structures over the cardiac cycle. This was achieved by applying non-rigid registration between the magnitude data of one specific systolic time frame and each of the remaining time frames. The systolic time frame chosen was the time frame with the highest signal magnitude, as this was the one that best resembled the fitted atlas obtained in the previous step. The resulting deformations were applied to the fitted atlas, generating a time-resolved segmentation that follows the motion of the vessels in the input dataset during the cardiac cycle.
Non-rigid registration during this step was performed using five scales, five iterations per scale, and linear interpolation. The amount of scales was increased in order to make the registration more sensitive to smaller differences between the images. Both fluid and elastic regularization of the displacement field were used. Additionally, this registration method uses tensor magnitudes to calculate the edges of the shapes in the images and bases the process on these edges. While this approach is usually used when executing multimodal registration, our initial tests indicated that it was helpful in detecting the small differences between time frames in our data. The process described in this section is schematically presented in Fig. 4.
Flow volume quantification
The deformation obtained when fitting the atlas to the input dataset was also used to locate the 2D planes positioned in the major vessels. As the planes were deformed by this process, a flat plane was fitted into the points of each deformed plane using principal component analysis (PCA) [29]. The principal components resulting from the PCA were used as the vectors that define the location of the plane.
The fitted atlas was further improved in each vessel plane by using a 2D circular averaging filter with a nine pixel radius, thus making the vessel region rounder, smoother, and also a bit larger. Having a segmentation that includes a region slightly larger than the vessel doesn’t considerably affect the flow volume obtained since this value is not dependent on the vessel area, and the pixels that directly surround the vessel usually have very low velocities. This strategy guaranteed that none of the pixels with higher velocities in the vessel region were left out of the calculation.
The velocity data were extracted from each of the 2D vessel planes using linear interpolation in the 4D flow CMR data masked by the time-resolved atlas created in the previous steps. An example of the appearance of the mask in 2D and how it follows the vessels during the heartbeat can be seen in Fig. 5. Finally, the flow volume in each vessel plane was calculated by integrating the volumetric flow rates over the cardiac cycle.
Evaluation
The accuracy of the results was evaluated visually and quantitatively for every input dataset. For the visual evaluation, each major vessel was divided into sections which were assessed and scored independently, as shown on Fig. 6. A numerical value between one and four was given as the score according to the following scale: 1 = poor segmentation, large errors are clearly visible; 2 = fair segmentation, intermediate errors visible; 3 = good segmentation, small errors visible; 4 = very good segmentation, no visible errors.
In the datasets where the segmentation was successful, net flow volumes obtained in the proximal ascending aorta with the proposed method were compared to those obtained in semi-automatically placed and segmented 2D planes. In order to determine the effect of the automatic positioning of the planes, the semi-automatically placed 2D planes were also segmented using the proposed automatic approach.
Additionally, flow volumes were calculated automatically at the following locations:
The net flow volumes at these locations are supposed to be closely related in subjects without cardiovascular shunts, which was also used to evaluate the proposed approach.
Furthermore, flow volumes in the proximal descending aorta and distal descending aorta were also calculated and compared between each other. The net flow volume is expected to only differ slightly between these two locations, as there are only minor branches in this section of a normal aorta.
Continuous variables are presented as mean ± standard deviation (SD). Simple linear regression analysis was used in order to evaluate the accuracy of the automatically calculated flow volumes versus the semi-automatically calculated ones in the ascending aorta. The same method was used to qualify the relationship between the flow volumes obtained at different vessel locations that are expected to be closely related to each other. A p-value <0.05 was considered statistically significant.
