Automatic Classification of Left Ventricular Regional Wall Motion Abnormalities in Echocardiography Images Using Nonrigid Image Registration

Data Description

The study was approved by the Regional Committee for Medical Research Ethics and performed according to the Helsinki Declaration. Written informed consent was obtained from each participant. The 2-D image sequences of ten healthy volunteers and 14 patients (ischemic heart disease) from two apical views (apical four-chamber (A4C) and apical two-chamber (A2C)) were acquired. All recordings were made using a General Electric ultrasound machine with a frame rate of 46–74 per second, including the electrocardiogram (ECG) display. In the present study, for each view, one cardiac cycle was stored, with care being taken to ensure that there was no respiratory movement or probe displacement during the data acquisition. One cardiac cycle was identified by selecting the two consecutive R-wave of the ECG signal which is synchrone with the end of diastole phase. Echocardiography images were acquired during standard clinical examinations and no study was excluded because of image quality.

The American Heart Association (AHA) proposed a standard LV division of 17 segments that had a correspondence with the irrigation areas of the main coronary arteries (Figs. 1 and 2) [2]. In this study, RWM of LV myocardium was evaluated visually by examining the motion of the interior and exterior of the LV myocardium and also thickening of each segment. Then, each segment was scored according to a three-point standard scale—1: normal (normal wall thickening and motion), 2: hypokinetic (reduced wall thickening and motion), 3: akinetic (absence of wall thickening and motion) [2, 3]. The segments were analyzed and scored independently by two highly experienced echocardiographers, and then the consensual visual RWM scores between the two echocardiographers were used as reference scores (gold standard) for all comparisons. Currently, for each subject, we only used A2C and A4C views. Consequently, according to the reference scoring, among the 336 analyzed segments (48 sequences), 188 (56 %), 44 (13 %), and 104 (31 %) were normal, hypokinetic, and akinetic, respectively.

LV Standardized Division

The approximate region of LV was extracted from echocardiography images (Fig. 3a) in each apical view. This was done by manually defining a rectangular region of interest (ROI) around LV on end-diastole image. This image was chosen because it has maximum LV volume at one cycle of heart and easily identified from the R-wave on the ECG signal. Then, the coordinates of the extracted rectangular ROI were automatically applied on all images of one cardiac cycle. These images are used for the following analysis in this study.

The LV in each apical view was divided up into seven regions in agreement with the LV standardized division by the AHA. To do this, for each view, three anatomical landmarks (the apex P1 and each side of the mitral valve P2, P3) and a distance d were manually located on end-diastole image by echocardiographer (Fig. 3a). Then, the LV division was automatically performed in three steps. First, the apex point P1 was connected to the mid-point P4 of the mitral valves points by a line, defined as the LV main axis. Second, the LV main axis was divided into three equal sections named apical, mid-, and basal using two orthogonal lines with length 2d. Third, apical section was divided into three equal regions by two radial lines with length d. Therefore, we defined a mask with seven regions located on the LV excluding the mitral valve (Fig. 3b). It should be noted that this LV division had been previously used for echocardiography images [27–31].

LV Displacement Field Estimation Using Image Registration

We consider an image sequence (images of one cardiac cycle) f(x, y, t) with t = t₀ (t = 0) ,… , t_I–₁ (t = I–1) and 0 ≤ x < X, 0 ≤ y < Y, where f(x, y, t) is the intensity of each pixel at position x, y and time t. Moreover, I is the total number of images of one cardiac cycle; X and Y denote the domain of the image.

Our goal is to estimate a displacement field, DF(x, y, t), for all pixels of end-diastole image (t = t₀) as a reference image over the cardiac cycle, that represents the position at time t of a pixel that was at position x, y at time t₀. In other words, we consider the end-diastole image as a spatial reference and we want to estimate the displacement field for all pixels of this image over the cardiac cycle. To do this, we register all images in a cardiac cycle to end-diastole image using a hierarchical transformation model [36]. The goal of image registration is to find the optimal parameters of the transformation model T(x, y). In conclusion, the displacement field for all pixels of the end-diastole image over the cardiac cycle is represented by the resultant transformation models as follow:

1

T_t(x, y) is the estimated hierarchical transformation model using the image registration at time t to end-diastole image.

The motion of the LV is nonrigid. Moreover, it has a global motion and a local deformation. Therefore, in this paper, the transformation is represented by a nonrigid hierarchical model which consists of a global and a local model as follow [36]:

2

The global motion model describes the overall motion of the LV over the whole sequence and is modeled by an affine transformation as follow:

3

The coefficients (a₁₁, a₁₂, a₂₁, a₂₂, a₁₃, a₂₃) parameterize the 6 degrees of freedom of the transformation, describing the rotation and translation of the LV.

The local deformation of the LV is modeled by a FFD transformation based on B-spline basis functions [37]. This nonrigid transformation is a powerful tool for modeling 2-D deformable objects and has been previously applied to the LV motion analysis in echocardiography images [38, 39]. FFD transformation deforms the shape of a 2-D object by manipulating an underlying mesh of control points [37]. The FFD is written as the 2-D tensor product of standard 1-D cubic B-spline:

4

Where i = ⌊x/n_x⌋ − 1, j = ⌊y/n_y⌋ − 1, u = x/n_x − ⌊x/n_x⌋, v = y/n_y − ⌊y/n_y⌋. Φ denotes a n_x × n_y mesh of control points (∅_i,j) with uniform spacing δ. Moreover, B_m(u) and B_n(v) are the basis functions of the B-spline evaluated at u and v, respectively as follow [37]:

The control points Φ act as parameters of the FFD transformation and the resolution of the control point mesh determines, simultaneously, the number of degrees of freedom of nonrigid deformation and consequently, the computational complexity. In the initial configuration, the control point mesh lies at their initial positions (Φ_initial) and by displacing or manipulating the control points, a desired deformation of the image is derived.

The result of image registration as an optimization problem is to estimate the optimal parameters of the hierarchical transformation (,Φ) so that a cost function defined in Eq. 5 associated with the image similarity and the smoothness of transformation is minimized.

5

We have used the sum of the squared intensities difference between homologous pixels in the two images as an image similarity criterion (Eq. 6). We chose to use this criterion because of its simplicity, fast computation time and good result even in the presence of decorrelated speckle noise. The similarity criterion is written as follow:

6

Where N is the number of pixels in the image. Moreover, the transformation should be constrained to be smooth to prevent artifacts such as folding [40]. A common smoothing term which regularizes the transformation [40] is written as follow:

7

In Eq. 5, λ is the weighting parameter which defines the tradeoff between the smoothness of the transformation and the image similarity criterion.

For computational efficiency, the optimization process proceeds in several stages. First, the affine transformation parameters are optimized by minimizing the cost function defined in Eq. 5. It is noted that the smoothness term of the cost function is zero for any affine transformation [40]. Then, the parameters of B-spline FFD transformation are optimized by minimizing the cost function defined in Eq. 5. In this step, we employ an iterative multiresolution optimization approach in which the resolution of the mesh of control points is increased, along with the image resolution, in a coarse to fine fashion to obtain a higher speed and robustness [41]. In addition, in each step, we apply an iterative gradient descent method [42] which steps in the direction of the gradient vector with a certain step size μ. The algorithm stops if the ║∇║ ≤ ε for some small positive threshold, where ∇ is the gradient vector of the cost function. Consequently, the final resultant positions of the control points of B-spline FFD (Φ_final) are the solution of our image registration problem.

A New Quantitative Regional Index

To quantify the RWM of LV, a new quantitative regional index which is related to final resultant position of control points of B-spline FFD model is proposed. To do this, after registering each image in one cycle of heart to the reference image, the sum of the squared Euclidean distance between homologous positions of initial (Φ_initial) and the final resultant (Φ_final) of the control points of B-spline FFD model for each region is computed. This value is divided by the number of control points in each region to normalize LV size and allow comparison between different cases.

8

Where n_r is the number of control points in region r and (Φ_final,r,t) is the final resultant positions of control points of B-spline FFD for image at time t and region r. Finally, for each myocardial region, the maximal value of the I_t,r through the images of one cardiac cycle is defined as the quantitative index for this region Ind_r.

9

Quantification of RWM of LV

After registering each image in one cycle of heart to the reference image for each case, we have computed the quantitative index Ind_r for each of seven regions according to Eq. 9 to quantify the RWM of LV. This regional index Ind_r which is related to final resultant position of control points of B-spline FFD model over a cardiac cycle, reflects wall motion and thickening of each myocardial segment of the LV. Therefore, it decreases with the severity of the RWM abnormality, being the highest values in normal segments and the lowest values in akinetic segments. The curves obtained from the values of the I_t,r (Eq. 8) for the images of one cardiac cycle in the apical, mid-, and basal segments of interventricular septum (left wall of LV in A4C view), separately are demonstrated in Fig. 5 for the two above cases. According to the reference visual scoring, these segments for the healthy case and the patient case are scored normal and akinetic, respectively. The maximal value of each curve is defined the quantitative index for this region Ind_r. The mean values and standard deviations of the quantitative indices Ind_r for the apical, mid-, and basal segments (A4C and A2C views) independently, in the three reference visual scores (normal, hypokinetic, and akinetic) are calculated and summarized in Table 1. These values indicate that the proposed quantitative index decreases gradually according to the severity of RWM abnormality. Moreover, results show that the basal segments have slightly greater values than the mid- and apical segments. Besides, the mid-segments have slightly greater values than the apical segments. These results confirm the fact that the wall motion increases from apex to base.

Classification of RWM of LV

Since the proposed regional quantitative index Ind_r decreases according to the severity of RWM abnormality, the estimation of two threshold values TH₁ and TH₂ separating akinetic from hypokinetic and respectively hypokinetic from normal is required to classify RWM of LV in a three-point scale. This process is done for apical, mid-, and basal segments of the LV, independently. Therefore, we must determine six threshold values (TH_1a, TH_2a, TH_1m, TH_2m, TH_1b, TH_2b).

To estimate the threshold values, a cross-validation procedure [43] is achieved and repeated ten times. For each trial, 2/3 of the cases (A2C and A4C views) are randomly selected for training and the remaining cases are used for validation. The distribution of the three reference RWM scores in randomly selected training subsets and the entire data sets is similar approximately to have an unbiased estimation of inter-scores thresholds during the training step. For each training subset, the thresholds are estimated by maximizing the weighted kappa coefficient [44] calculated from the 3 × 3 contingency table obtained from the comparison between the reference visual scoring and the proposed automated scoring. As mentioned above, this process is done for apical, mid-, and basal segments of the LV, independently. Consequently, for each trial, the six thresholds are estimated using this process. These thresholds are used for the quantitative scoring of wall motion of each segment of LV as normal, hypokinetic, or akinetic. The resultant scoring obtained from the proposed method is compared against those assigned by the reference visual analysis using two parameters, the absolute and the relative agreement. The absolute agreement is defined as the percentage of segments for which the scores obtained by the two methods are equal and the relative agreement is defined as the percentage of segments for which the scores obtained by the two methods differed is no more than 1. The values of agreement between two methods averaged over the ten subsets in the validation and training subsets for the apical, mid-, and basal segments, independently are demonstrated in Table 2. Results show that the apical segments which are known to be more complex to analyze have fewer absolute agreement values than the mid- and basal segments. Moreover, results show that the difference between the automatic and the reference visual scores in only apical segments is more than 1. In addition, there is absolute agreement of RWM scores between two methods averaged over the ten subsets in 88 ± 1.1 % of normal, 64 ± 2.2 % of hypokinetic, and 82 ± 1.3 % of akinetic segments as classified by the reference visual scoring. Finally, the same performances are obtained on both the validation and the training subsets indicating that the cross-validation procedure performs correctly without ‘over learning’.