Magnetic Resonance Derived Myocardial Strain Assessment Using Feature Tracking

1. Introduction

Automatic detection of borders is a fundamental issue in image analysis. In cardiac imaging, the possibility of an automatic detection of the endocardial border in the imaging of the left ventricle would give objective measurement of the ventricular volumes, and myocardial deformation (strain). This was accomplished in echocardiography with speckle tracking technique. The development of reliable methods for the automatic border detection is a challenging task that has not received a generally reliable solution in cardiac magnetic resonance (CMR). In fact, in clinical practice, borders are either drawn manually by the operator or software detects interface between myocardium and cavity ^1,2. In current article we introduce a different approach where the borders are not “detected”, rather they are “tracked”, i.e. followed in time, starting from one reliable existing instantaneous trace that is commonly -but not necessarily- manually drawn by the experienced operator over one single frame. The individual points composing such a first reliable trace are followed in time by searching the same features that are about one point in its neighborhood in the following frames. The tracked features can be the cavity-tissue boundary or anatomical elements that are different along the tissue. They are found by methods of maximum likelihood in two regions of interests between two frames.

The local frame-to-frame displacement is equivalent to evaluating the local velocity (ratio between displacement and time interval). The automatic evaluation of the velocity at a point is determined from comparison of the displacement of the image data about such point in two consecutive frames. Such methods have been used, in several different formulations, in many research fields. They fall in the general category known as Optical Flow, in advanced image analysis ^3,4. They are commonly referred as Speckle Tracking in echographic imaging when such velocities are used to follow physiological motion ^{5, 6} but also apply to any other image modality such as CMR where these methods are referred as feature tracking or border tracking.

2. Material and Methods

Feature Tracking Method

Endocardial or epicardial border of a 2D CMR cine is manually traced on one arbitrary frame (see figure 1). Mid-myocardial features can be traced as well. Such border is then defined as a sequence of N points, identified by their coordinate pairs (x_i,y_i) with i=1…N. The border tracking proceeds by tracking each single point, such a tracking is based on a hierarchical algorithm at multiple scales and by a combination of 1D tracking techniques, which guarantee higher accuracy, and 2D tracking, which is necessary to properly detect the 2D spatially extended features.

In order to first capture the large geometrical displacement of the border, the tracking is performed in the direction orthogonal to the border itself where the cavity-tissue boundary is best recognizable. The tracking along this direction is performed by using the method of transmural cuts as follow (see figure 3). A line crossing the wall, passing through the point and orthogonal to it is drawn. The pixels taken along the transmural line are placed in columns, each column corresponding to one frame of the sequence of images. In this way the evolution along a transmural cut can be represented for all instants at once in a two-dimensional representation where one axis is the distance along the line and the other axis is the time (see figure 2). This representation is similar to what is referred to as an M-mode in Echocardiography, in CMR it corresponds more to the “scout” function. To improve the quality of the analysis, in the case of poor images with a low signal to noise ratio, the space time representation is built using a line for the transmural cut with a thickness of 5 pixels. The border tracking is then performed along the space-time image.

In a second step to account for the 2D displacement of the border, a standard 2D tracking (optical flow-based) is performed, for each point independently, on a MxM moving window that is always centered on the previously estimated border point. 2D tracking is performed in two steps, where half of the first estimation is employed to center the moving windows in the second tracking passage. The window is then reduced from 32 to 16 in two additional passages.

To improve the accuracy of the motion along the border that is used to estimate rotation and torsion, the 1D tracking is performed along space-time images built from thick cuts “parallel” to the curved border (see figure 3). At each point, independently, the pixels taken along the moving border, centered at the moving border points, are placed in columns, each column corresponding to one frame of the sequence of images. To improve the quality of the analysis, and to best capture the features at the border the line is extended of 5 pixels into the tissue (sub-endocardium). The border tracking is then performed along the space-time image with the same procedure described above. To ensure the spatial coherence in the tracked border, a 3 point median filter and a 3 point Gaussian filter (of weights 0.25, 0.5, 0.25) is applied for the displacement computed at neighboring points at each step.

Tracking Along the 2D Space-Time Image

This section describes a procedure for following a border along one direction in a two-dimensional image (M-mode-like) starting from a known position at one instant.

X is defined as the horizontal direction and y the vertical direction. Columns are annotated x_i, i=1…M, where M is the number of columns in the image. The tracking is given by determination of a discrete sequence of real numbers y_i=y(x_i), starting from a known point y_k corresponding to the columns x_k.

The displacement from the known point y_k to the point y_k+1 is estimated by evaluating the cross-correlation between the entire column at x_k with the entire column at x_k+1. The cross-correlation function will present a maximum, the position of the maximum gives the value of the vertical displacement required to maximize the similarity between the two columns, therefore y_k+1 is estimated by adding such a displacement to y_k. This procedure is repeated between all pairs of nearby columns and the result is an estimate of the entire border y_i, i=1…M. The cross-correlation is here computed using a Fast Fourier Transform algorithm to reduce calculation time.

The first estimate y_i is further refined iteratively. To accomplish this aim a subset of the image is extracted by taking a few points above and below the previous estimate y_i and a new image whose center corresponds to the sequence y_i is generated and used for the correction tracking. This refinement is repeated until no correction is found.

An improved and more natural result is then achieved by a final snake procedure [5] to follow, in the space-time image, the image brightness level that passes through the fixed point y_k. The entire process makes use of the time periodicity to ensure a periodic result and avoid the drift effect.

Technical Limitation of Feature Tracking

The border tracking technique, like any speckle tracking method, is based on quantification of changes on pixel brightness from one frame to the other. This gives a lower limit to velocity related to the need to see a speckle that is one pixel at one frame, moving to the neighboring pixel in the next frame. This limit is therefore

Equation 1
where Δx is the pixel size and Δt is the time interval between the two frames. The coefficient k depends on the quality of the tracking algorithm and on its ability to evaluate dynamic sub-pixel variations. This limit means that velocities that are well above this limit are estimated with great accuracy, such accuracy is reduced when velocity values approach and fall below such a limit.

This limitation also implies that an increase in the acquisition frame-rate (reduction of Δt) on one side allows an easier evaluation of large velocities and their rapid variations (like during the isovolumic phases). On the other side, the increase of the frame rate (reduction of Δt) increase this limit and imply a reduced accuracy in the evaluation of lower velocities until it is not accompanied by a similar increase of spatial resolution (reduction of Δx).

Phantom Image Preparation

A series of artificial computer-generated loops has been prepared to allow testing of the image analysis procedure in simple and perfectly controlled conditions. For this, a phantom in a short axis projection of an ideal left ventricle was prepared as follows.

The endocardial and epicardial borders are represented by two concentric circles with radius R₀(t) and R₁(t), respectively. The image is prepared by making the annulus, which represents the tissue between the two borders, as uniformly colored gray on a black background. Then an 8×8 top-hat linear filter is applied to avoid unphysical discontinuities.

The epicardium movement is taken, in [mm], as R₀(t)=10+5cos(2πt/T) where T is the heartbeat period taken as T=1s. The theoretical endocardial kinematics is constant along the border and depends on time only, velocity is only radial and given by V0(t)= dR₀/dt=-π sin(2πt/T), in [cm/s]. Percentage strain, computed relative to the length the border has at time zero, is St₀(t)=100x(R₀(t)-R₀(0))/R₀(0)=100(cos(2πt/T) -1)/3, and strain rate follows from (1) as SR₀(t)=10 V0/R₀, in [s^-1]. The epicardium is assumed either as moving accordingly to a constant thickness, R₁(t)= R₀(t)+5mm, or as still R₁(t)= R₀(0)+5mm.

Each image is square of size of 48mm, centered on the tissue annulus, and has a resolution NxN. Example images are shown in figure 4, plates a and b; the strain and strain rate time profiles are shown in figure 4, plates c and d. The loops are prepared by varying the resolution N, the frame-rate FR, and the epicardial type of motion.
The endocardial tracking method is applied to such images by taking on the first frame a number N_p of points uniformly spaced along the circular endocardium.

3. Representative Results

Phantom Study

The application of the image analysis method to the computer-generated phantom images is here analyzed. A global measure of the eventual error is computed by the root mean square percentage difference. The root mean square, average and maximum errors in the endocardium strain are defined as

Equation 2
where St₀(t) is the exact value, St(t) is the value computed by the image analysis, and the summations extend over all frames N_F=FRxT. The same definition is used for the radius, velocity, and strain rate. The tracking is about independent from the position along the endocardium, the differences between the different points is well below 1%.

Results are summarized in Table I for 15 phantoms with varying spatial resolution, frame-rate, and epicardial border type of motion; the effect of varying the number of points used to track the endocardial border is also shown.

Errors are in all cases very small for the integral quantities (radius and strain) and slightly larger for the differential quantities (velocity and strain rate) that are related to the derivative of the former. This was expected because the derivative operator amplifies errors. The quality of results is degraded when the resolution is reduced; in fact, the accuracy is related to the pixel size that represents (in a loose sense) the minimum displacement readable from one frame to the other. The time resolution does not affect the results significantly until the frame-rate is sufficient, at very high frame-rate result do not improve because frame by frame displacements becomes lower than the pixel size. This shows that an increase in the frame-rate is of little or no utility when it is not accompanied by an increase in the spatial resolution.

However, the simple sinusoidal motion here considered does not require an extreme time resolution. Similarly, the use of as few as 8 points is sufficient to follow the simple, circular, endocardial shape. Endocardial results are not appreciably influenced by the type of motion that the epicardium undergoes. We have also verified that the results are not significantly affected by the adopted image filtering.

A visual presentation of results is given in Figure 4 where the computed endocardial border at two instants is reported over the phantom images (plates a and b). The strain and strain rate are reported in (plates c and d) for the case #1 and the small resolution case #8. The strain and strain rate in case #1 (squares) presents an excellent agreement with the theoretical value, the mean error being equal to 0.6% and 3%, respectively. The agreement is only a little worse in case #8, where image resolution is halved, with errors 0.9% and 4.5% for strain and strain rate, respectively.

Clinical Validation ¹.

We compared mid-LV whole slice circumferential myocardial strain (ε_cc) by the Harmonic Phase Imaging (HARP) and FT techniques in 191 Duchene Muscular Dystrophy patients grouped according to age and severity of cardiac dysfunction and 42 age-matched, control subjects . Retrospective, off line analysis was performed on matched tagged and SSFP slices. For the entire study population (n=233), mean FT ε_cc (-13.3 ± 3.8 %) were highly correlated with HARP ε_cc (-13.6±3.4 %) with a Pearson correlation coefficient of 0.899. The mean ε_cc of DMD patients determined by HARP (-12.52 ± 2.69 %) and FT (-12.16 ± 3.12 %) were not significantly different (p=NS). Similarly, the mean ε_cc of the control subjects by determined HARP (-18.85 ± 1.86) and FT (-18.81 ± 1.83) were not significantly different (p=NS). We concluded that FT-based assessment of ε_cc correlates highly with ε_cc derived from tagged images in a large DMD patient population with a wide range of cardiac dysfunction.

Table 1. Phantom analysis of endocardial border tracking: root mean square and maximum percentual errors [%] are computed for the main quantities in correspondence of different phantom parameters. The parameters marked bold indicate the variations from the Phantom #1. The dependence on frame-rate, resolution, and number of tracked points is considered. The influence of the type of epicardial motion is considered for the two limiting cases when the epicardial border does not move (no motion) or moved with endo (no thickening). The last phantom (*) is built without filtering the basic stepwise images with brightness changing abruptly in one pixel. Errors above 10% are marked bold.

Figure 1. cMR image of the left ventricle, in long axis view (left picture) and in short axis view (right picture), with a traced endocardial border drawn on top.

Figure 2. Space-time representation, where space is along a transmural cut, of the image sequence. The transmural cut is taken as for the starting point in figure 3. The time evolution of the starting point, tracked automatically, is reported.

Figure 3. Image of the left ventricle, in long axis view, with transmural cuts and cuts parallel to curved border.

Figure 4. Phantom study. Two images (case #2) at maximum expansion (plate a) and contraction (plate b), the computed endocardial border points are overlapped. The strain (plate c) and the strain rate (plate d) computed with two different phantoms (cases #1 and #8) are shown in comparison to the effective values.

Figure 5. Examples of global circumferential (black curve) and segmental strain (color curves) in normal patients (a). Example of global circumferential (black curve) and segmental circumferental starain (color curves) in patients with depressed left ventricular function and left bundle branch block (b). Note different timing of peak circumferential strain indicative of presence of left ventricular dyssynchrony