Longitudinal strain from velocity encoded cardiovascular magnetic resonance: a validation study

Background Regional myocardial function is typically evaluated by visual assessment by experienced users, or by methods requiring substantial post processing time. Visual assessment is subjective and not quantitative. Therefore, the purpose of this study is to develop and validate a simple method to derive quantitative measures of regional wall function from velocity encoded Cardiovascular Magnetic Resonance (CMR), and provide associated normal values for longitudinal strain. Method Both fast field echo (FFE) and turbo field echo (TFE) velocity encoded CMR images were acquired in three long axis planes in 36 healthy volunteers (13 women, 23 men), age 35±12 years. Strain was also quantified in 10 patients within one week after myocardial infarction. The user manually delineated myocardium in one time frame and strain was calculated as the myocardium was tracked throughout the cardiac cycle using an optimization formulation and mechanical a priori assumptions. A phantom experiment was performed to validate the method with optical tracking of deformation as an independent gold standard. Results There was an excellent agreement between longitudinal strain measured by optical tracking and longitudinal strain measured with TFE velocity encoding. Difference between the two methods was 0.0025 ± 0.085 (ns). Mean global longitudinal strain in the 36 healthy volunteers was −0.18 ± 0.10 (TFE imaging). Intra-observer variability for all segments was 0.00 ± 0.06. Inter-observer variability was −0.02 ± 0.07 (TFE imaging). The intra-observer variability for radial strain was high limiting the applicability of radial strain. Mean longitudinal strain in patients was significantly lower (−0.15± 0.12) compared to healthy volunteers (p<0.05). Strain (expressed as percentage of normal strain) in infarcted regions was lower compared to remote areas (p<0.01). Conclusion In conclusion, we have developed and validated a robust and clinically applicable technique that can quantify longitudinal strain and regional myocardial wall function and present the associated normal values for longitudinal strain.

This appendix describes the motion tracking and strain calculation algorithm in detail. Starting from velocity encoded data the myocardium is tracked throughout the cardiac cycle.
The myocardium is represented as a spatiotemporal finite-element mesh that undergoes a nonrigid cyclic motion.
In an ideal case without noise the displacement of mesh-point could be calculated as follows: where   , u x t is the displacement of a node over time t, ˆ( ) v x is the measured velocity, and   , 0 0 u x  . The goal is now to construct a robust method of solving (A1) also in the presence of noise and artefacts. To do this we solve an optimization problem where we minimize the difference with the measured velocities and the time derivative of the tracked nodes as: where   , u x t is the displacement of a node over time t, ˆ( ) v x is the measured velocity.
Instead of solving Equation A2, directly for the deformation u, we instead represent the deformation on the boundary with local basis functions. This will introduce spatial regularity and smoothness. The deformation in the interior of the myocardium is based on the deformation at the boundary and assuming a material model. Furthermore, the whole computation process needs to be able to handle constraints in form of user supplied manual corrections of the deformation when necessary for an accurate tracking.

A2
A flowchart of the entire algorithm is given in Figure A1. Figure A1. Flowchart of the complete motion tracking and strain calculation algorithm The optimization problem is solved iteratively. The geometry of myocardium is approximated by manual delineation of myocardium boundary on the long axis view in the first time frame. The region is divided into small triangular elements using Delaunay triangulation algorithm. An example of the triangulated mesh is illustrated in Figure A2. The calculated deformation provides the mesh configuration at the subsequent time frames. In the following algorithm description we will differentiate between boundary nodes, which are finite element nodes that reside on the boundary of the myocardium. Internal nodes are finite element nodes that do not lie on the boundary of the myocardium.
A3 Figure A2. The image shows the finite element mesh in the first timeframe of a typical healthy normal. The red box indicates location of a zoom to show typical mesh density.
Step 1: Set up boundary function The deformation of each boundary node is calculated from a set of local basis functions. Seven local basis functions are used and their centrums are located equally distributed along the center line of the myocardium. For each local basis function the following four modes of motion are allowed; radial expansion, shear, and outward motion, longitudinal contraction (including passive motion). These four modes of motion have been chosen to be physiologically plausible.
The longitudinal motion deserves a special comment. Consider a case where the there is only contraction in a midventricular section, no contraction either basally or apically. Then we would expect that the basal portion of the ventricle is passively dragged towards the apex (no A4 longitudinal strain). This is the purpose of including the passive motion in longitudinal contraction mode. In order to describe the motion modes algebraically we define α i is the local angle of the node i, d i as the distance of the node i to the centre line (defined as positive on the endocardial side of the centre line, and negative on the epicardial centre line), a i the coordinate of the node i along the arc length, where the arc length is normalized to -1..1, a j is the coordinate o the centre of the local basis function j along the arc length, w ij as a weighting function for node i with respect to local basis j, p ij is a weighting function that includes passive motion of longitudinal contraction for node i and basis function j. The index i goes from 1 ... number of boundary points, and index j goes from 1..7. The weighting function is calculated as: where σ = 0.05. The weighting function that includes passive motion of longitudinal contraction was defined as: where sign(a j ) is the sign function defined as -1 for a j <0, and 1 for a j >0, iexp is the integral of the exponential function.
Radial motion for node i and basis function j was defined as: Shear motion for node i and basis function j was defined as: Outward motion for node i and basis function j was defined as: Finally the longitudinal motion including passive motion was defined as: The basis function are formulated in a matrix notion. This matrix is denoted as P and is a sparse matrix with 4 j columns and one row for each of the boundary node i.

(A9)
Step 2: Mapping of boundary functions onto the interior nodes As previously mentioned we are only explicitly calculating the deformation of the boundary nodes. The displacement of the interior nodes is introduced in the calculation by using a matrix G, which is constructed by adjoining zeroes to matrix P (from Step 1). The process is illustrated in the following example: Suppose the index of boundary nodes are 1,3,5,6,7,9,10,11 and index of interior nodes are 2, 4, and 8. G is constructed such that: A6 Rows 2, 4, and 8 are zero rows.
The matrix G will be used in step 4.

Step 3: Discretization of the material operator
The Cauchy-Navier operator N is used to ensure spatial smoothness and spatial continuity of the deformation. Since we are tracking a periodically moving and deforming object, we consider a finite element mesh whose node points are attached to a set of material points of the underlying myocardium, such that the mesh moves and deforms along with the material region. The use of the Cauchy-Navier operator can be seen as way of adding linear material properties in the tracking formulation. The material operator N is discretized by finite element operators. This is formulated in matrix notation as: where k are 1/(1-2ν) where ν is a material property constant and was set to 0.45, , D x and are derivative matrices constructed from the derivative values at each node of all the triangular elements with respect to x and y respectively, and M is a diagonal matrix with diagonal entries as the areas of the triangular elements. N is defined for both boundary nodes and interior nodes.

4: Construction of projection matrices
The displacements of all the nodes (both boundary and interior) are calculated as a matrix multiplication.
where c is a coefficient vector of basis functions. A is a mapping matrix for all nodes, and is constructed as where N is the Cauchy-Navier operator, P is the basis function matrix, and G is the mapping matrix onto the interior nodes. N, P and G are defined in the previous steps. When taking the temporal derivative of Equation A13 twice we get the following This can be expressed in matrix notation as where H is a temporal second derivative operator written as .
The Equation A16 can be re-arranged into Step 5 (the temporal second derivative of displacement) is calculated using forward differentiation of the measure velocity data.
Step 6: Iteratively optimizing the coefficient vector ). This iterative process for calculating the coefficient vector is repeated 10 times.
Step 7: Tracking after manual corrections The user interface for manual interaction allows the user to pre-specify the displacement of specific boundary node(s) at specific time frame(s). The nearest boundary node to each manual correction is stored. To maintain local dependency of manual corrections, the boundary functions are reset by including local first order polynomial basis into the matrix P (for details see Step 8) which give rise to matrix new P . Using new P , matrices A and K are recalculated as given in Equation A20 .
Each manual correction introduces a constraint (u = Ac) into the system of linear equations (Hc = J) of coefficients c. We are forming a new matrices H new and J new , where we add one displacement constraint one for each manual correction. Now we can again iteratively solve for the coefficient vector c, but now instead of using Equation A19, we use the following equation This is repeated 10 times.
Step 9: Strain tensor calculations As a last step the strain tensor is calculated. First the displacements of all nodes are calculated according to Equation A13. If x u and y u are the displacement components in x and y directions, then 2 dimensional strain components x  , y  and xy  on the faces of a strain element aligned with the xy coordinate system are obtained as follows: , , If the motion of myocardium is described in lateral (x') and axial (y') directions, then the strain components ' x  (longitudinal strain), ' y  (radial strain) and ' ' x y  (shear strain) with respect to an element aligned with the ' ' x y coordinate system, rotated by angle from xy are: