**Motion Analysis of the Right Ventricle From**MRI Images Edith Haber 1, Dimitris N. Metaxas 2, and Leon Axel 3 1 Department of Bioengineering 2 VAST Lab, Department of Computer and Information Science 3 Department of Radiology University of Pennsylvania, Philadelphia, PA 19104-6389, USA edith@seas, upenn, edu http://~, cis. upenn, edu/~edith/ Abstract. Both normal and abnormal right ventricular (RV) wall mo- tion is not well understood. In this paper, we use data from tagged MRI images to perform the first 3D motion study of the entire right ven- tricle to date. Our technique is an adaptation of a physics-based de- formable modeling methodology that was successfully used on the left ventricle(LV). As opposed to the previous approach, currently we use segmented contours to generate the geometry, 1D tags for our input data (due to the thinner RV), and localized degrees of freedom (DOFs) with finite elements. Although we build a biventricular model, our results focus on method validation and visualizing clinically useful parameters that describe RV wall motion. 1 Abnormal motion of the RV serves as an indicator of several types of heart disease, such as RV ischemia and hypertrophy [3]. However, there is no in-depth knowledge of the RV motion and its correlation to the various diseases. One reason that researchers do not agree about the exact pumping mechanism of the normal RV may be regional variations in contraction patterns occuring at different phases of systole. Since studies have found that RV contraction varies with increases in pressure and volume, and ischemia [3], a method for accurately assessing RV wall motion can both answer questions into its normal function and can be used as an indicator of various diseases. The RV receives blood from the right atrium and pumps it into the pul- monary artery. It appears crescent-like in a cross-sectional view and, unlike the LV, is difficult to define with any parameterized 3D shape. The RV shares a septum with the LV, while its outer free wall is in mechanical contact with the pericardium and the lungs. The RV cavity can be conceptually separated into an inflow tract, a highly-trabeculated apical portion, and a relatively smooth out- flow tract [6] (Fig. 1). The 3ram-thick free wall is thin relative to the 7ram-thick LV free wall. This relative thinness and complex geometry make it difficult to capture RV wall motion using regular imaging modalities. This paper presents a new methodology for modeling and analyzing the RV shape and motion from MRI-tagged data. We make significant modifications to Introduction178 RV Inflow Apex Fig. 1. Schematic of right ventricular compartments the physics-based deformable modeling framework of [8, 12], which was applied to MRI tag intersections for LV wall motion reconstruction. At the currently available image resolution, a small number of tag point intersections fall within the normally thin-walled RV. In order to capture the full 3D motion, the new method uses different constraints to incorporate input data consisting of the entire lengths of 1D tags extracted from multi-view MRI-SPAMM images. Also, due to the complex shape of the RV, we cannot construct our initial model by fitting a deformable geometric primitive to the initial contour data. Instead, we build a deformable model of the biventricular geometry directly from the initial contours. Although this model includes both the right and left ventricles, we focus on the analysis of the RV in this paper. The initial geometric model is a volumetric finite element model (FEM) mesh that we fit to subsequent contour and tag data. Due to the complex RV geometry and motion, the analysis consists of using the maximum possible DOFs. Instead of lumping them into parameter functions as in [12], we use stiffness from finite elements to smooth the sparse input data which result from imaging the thin RV wall. However, in order to make our analysis clinically useful we convert the model's DOFs into equivalent lumped parameters such as scaling and twisting. In addition to the quantitative analysis, we use computer graphics techniques to visualize the reconstructed RV shape and motion by color mapping the distri- bution of derived parameters onto our model. Our preliminary results quantify previously documented knowledge of the RV motion and provide a more detailed motion analysis of the RV compared to previous methods. In addition, having a biventricular model opens the way for the quantitative study of inter-ventricular dependence and relative motion. 2 Right ventricular wall motion has been studied invasively by implanting markers into small portions of the RV wall and using an imaging modality to track the motion of the markers [16, 4,18, 14]. While results from these techniques showed regional variations in contraction, questions remain about whether the invasive nature of these techniques affect contraction patterns [4]. Although, the MRI tagging method has provided a non-invasive method for studying heart wall motion, most of these studies have centered on the LV. Since the two-dimensional Related Work

179 (2D) images do not capture important through-plane motion, data from multiple views must be integrated by fitting a 3D model to the data [12, 19, 9]. To the best of our knowledge, the only work that has used entire tag lines (as opposed to tag intersections) is that of O'Dell,et al. [11]. These researchers developed a displacement field fitting technique which was used to fit an ellipsoidal finite element model to LV data. Up to now the only MRI tagging study of the RV free wall was performed by Young, et al. [20], who fit a surface model to intersections of mid-wall contours with the tags. In the past, RV geometry has been reconstructed from both CT and MR images as these technologies developed [10]. Integrated models for the RV-LV geometry have also been developed using tools from constructive solid geometry [13]. In addition, the right ventricutar free wall was modeled as a thin shell by fitting biquadratic surface patches [15]. 3 Image Acquisition and Data Extraction Fig. 2. Mid-wall short-axis 1D tagged MR images at (a) beginning systole (b) mid- systole (c) end~systole. In order to capture the motion in the ventricular walls using MRI images, we use a 1D SPAMM (SPAtial Modulation of Magnetization) tagging technique [2]. This non-invasive technique produces a family of parallel tagging planes by using a sequence of non-selective radio frequency excitation separated by intervals of magnetic field gradient in the direction perpendicular to these planes prior to imaging. The intersection of these tagging planes with the imaging plane pro- duces dark stripes (Fig. 2) which represent actual tissue. Therefore, the initially straight lines deform as the heart contracts. Due to through-plane motion, the portion of the tag plane seen in an initial image may move out of the imaging plane and be replaced by the intersection of a different portion of the same tag plane with the imaging plane. Thus, the stripes only provide information about motion in the direction perpendicular to the initial tag plane and we need 3 different views to capture the full 3D motion, e.g., 2 short axis views (one with horizontal and one with vertical stripes), and a long-axis view (Fig. 3). The fig-

180 ure depicts how each tag plane/image plane combination provides information about motion in mutually perpendicular directions. z z x X Y Y X / Fig. 3. 1st row: Image and tag (dark) planes. 2nd row: 2D image planes with example of tag motion from initial (dark lines) to final (dashed lines) time. Highlighted arrows indicate for which direction we get motion information. We extract tags and contours from the images using SPAMMVU, a program developed in the Department of Radiology at our institution [1]. Tags from multiple time phases (between end-diastole and end-systole) were tracked by adapting the deformable mesh scheme of [21] for 1D stripes. The extracted stripes are approximated as a series of points by sampling the tag lines at 2ram intervals. 4 4.1 Biventricular Model Geometry Deformable Model Although the RV cannot be defined by any simple parameterized geometric primitive, an accurate geometric model is important for a model-based fitting technique. Since the septum plays a role in the function of both ventricles, our approach is to create a biventricular geometric model. Due to the complex geome- tries involved, we build a discretized mesh of volumetric finite elements directly from contours extracted from end-diastolic images [5]. The short-axis contours were used to generate the finite element mesh, with the insertion points of the RV free wall into the septum as guide points. Points sampled for the contours then became the nodes of the finite elements. This geometric model could also be used for volume measurement, and, given cardiac material properties, for stress analysis. Our deformable model fitting technique results in the displacement of over 500 nodes through time. In order to make these results clinically meaningful we convert nodal displacements into scaling and twisting parameters. Since it would be less relevant to describe the twisting of the RV with respect to the long axis of the LV (or vice versa), we define a local coordinate system for

181 Fig. 4. Initial finite element mesh. (a) Anterior view of with sections cut away to show local coordinate systems. Note that the RV apex is above that of the LV apex. (b) Cross sectional view of a layer elements which connects contours from two image planes each ventricle. The axes of the coordinate systems are the eigenvectors of the matrix of central moments [8], with the eigenvector corresponding to the smallest eigenvalue is the one with the most inertia being labeled the z-axis (see Fig. 4). The transformation of position x with respect to the global, non-inertial frame to the local position s, is x = c + Rs (1) where c is the origin of the local frame and R is the rotation matrix which takes a point from the global to the local coordinate system. The RV free wall is in the RV coordinate system, while the LV free wall and septum are both in the LV coordinate system. Some of the finite elements connect nodes across the model so that two systems remain connected. The finite element mesh, with a shaded RV and LV endocardium, is also shown in Fig. 4. 4.2 Model Dynamics Our basic approach is to fit the model to the data with physics-based deforma- tions as in [12]. Unlike the previous work, our model deforms with different kinds of DOFs, i.e., the position, q~, of each node i of the finite element mesh. We use a non-inertial model and integrate the following equation of motion to solve for qi: r = fi,intcrnal + fi,c=ternat (2) where fi,inte~nat is the finite element stiffness force, and fi,~=te~,~at is the force derived from the image data. The forces are only used to deform our model and are not meant to replicate the actual forces in the heart wall muscle. The large number of DOFs requires we have internal stiffness forces in order to main- tain continuity in the motion. As mentioned earlier, we later lump those DOFs to parameters which are more recognizable to physicians. In the following, we

182 describe the calculation of the internal force and two types of external images forces, which we call contour and SPAMM forces. External Forces From Contour Data We use the extracted contours to maintain the shape of the deformable model. These contour forces fi,c, are cal- culated in the same way as done by [12] and applied to the appropriate node(s). In order to distribute these forces as uniformly as possible, the contours were first sampled at 5ram intervals. External Forces from SPAMM data In order to capture within-wall motion, we apply 'SPAMM forces' to our model that will mimic the motion of the tag stripes. We can accomplish this because we have a time correspondence between stripes and the finite elements allow us to register the initial stripes locations to the model. We register each point on the stripe to the non-deforming, stationary local coordinate system of the appropriate finite element. The transformation from the local position, (e,n, s), to the global position of a point, (x, y, z), is written in terms of the finite element shape functions: n n ~. (3) x = ~ Nj(e,n,s)xj, y = ~ Nj(e,n,s)yj, z = Z Nj(e,n,s)zj j=l j=l j=l where xj, yj, zj is the position of the jth node in the element numbering system. The shape functions, Nj, can be seen as weighting the global coordinate of a node according to where a point lies in the (e, n, s) system. In this paper, we use six-noded wedge and eight-noded parallelepiped elements, whose linear shape functions are given in [22]. Since we only know the positions of the nodes (xj,yj,zj) and the position of a tag point on the stripe (x,y, z), we solve the set of three equations (Eq. 3) for the local coordinates (e, n, s) using Newton- Raphson. Since each finite element is defined to be a 2 by 2 cube centered at the origin of its local coordinate system, the point falls within a particular element if: -1 <= e,n,s <= 1. f f ~ ~ normal time I time 2 i/ eO -21 Fig. 5. Application of SPAMM forces

183 For each material point from an extracted tag, we define the spring-like SPAMM force, fs, to be a scalar multiple of the distance vector, e : fs -- ae, where a is the strength, e is a vector in the direction of the tag plane normal which extends from a material point to its corresponding tag stripe at the next time phase (see Fig. 5). After time 0, the material point will move out of the plane due to forces applied on other points in other directions. For these cases, i.e. time 1, the projected point (P) of the material point (M(1)) onto the original image plane, rather than the material point itself, is used to calculate e. The SPAMM force, fs, is applied within the element and must be distributed to the nodes of that element. The force on each node is weighted by the shape functions: fj,s = Njfs where the subscript j refers to the node number within an element. This is the manner in which concentrated loads are distributed to the nodes in finite element theory, so that the collection of material points registered to the model have a 'line of forces' applied to them. If the same node receives forces from different material points, the final force is the average of all fj,s. The total external force on each global node, i, is the sum of the contour and the SPAMM forces, i.e., fi,ezternat = fi,c + fi,s. (4) Internal forces due to stiffness In our formulation, we model stiffness as an internal force, fi,internal. We use stiffness to impose a continuity and smoothing constraint on our many DOFs, and not to model the actual cardiac material properties (a topic of future research). As a result, we consider each element to consist of an isotropic, linear, incompressible material. The stiffness force, fi,internal, on a node from a particular element can be computed from fi,internat = Kd (5) where d are the nodal displacements. The stiffness matrix, K, incorporates the geometry and material properties of the element and is computed from where D is the stress-strain matrix and B is the strain-displacement matrix, whose formulations are given in [22]. 4.3 We convert the reconstructed displacement results into a twisting parameter, 0, of a point (s~,Sy,Sz) about the z-axis of the local coordinate system. The equations which define twisting are s~ = sx cos 0 - sy sin 0 and su = sx sin 0 + sy cos 0, where an uppercase subscript denotes the point at its initial position. Since simply subtracting the angles will give erroneous results when a point moves across the 0 = 0 line, we solve the two equations for sin ~ and cos 0, and solve for 0: 0 = arctan SXSy - 8ySx 8y8y ~- 8X8x Model Parameters: Twisting (7)

184 5 Results 5.1 In order to validate the fitting method, we built a computational phantom in the shape of a thick-walled, annular cylinder. We deform the model using equations in cylindrical coordinates [19, 7]: Validation z = Az + r (8) e = O + - + 7(Z - + (9) (10) where (R, O, Z) = initial position, (r, 0, z) = deformed position, Ri = inter- nal radius, and Zmin ~- minimum height. Twisting within the cylinder is the superposition of a minimum twist (Omin), twisting between short-axis layers (controlled by ~2), and twisting between the inner and outer walls (controlled by 7)- Axial deformation is varied in the axial and radial directions with the use of the )~ and r parameters, respectively. Finally, an incompressibility constraint leads to radial deformation being a function of the initial radius. We generate parallel tag planes with a 7ram separation with tag points sam- pled every 2ram along a stripe. Simulated short-axis image plane (5 total) were separated at 7.2ram intervals, and simulated rotating long-axis image planes (as in our actual imaging protocol) were separated at 20 degree intervals. As done in [19] , we use Eq. (8 -10 ) to solve for unknown r, 0, Z tag and contour po- sitions from known R, O, z for the short axis, and to solve for unknown r, O, z tag and contour positions from known R, 0, Z for the long axis. Geometric and deformation parameters are shown in Table 1. Phantom Parameters Geometric mm Ri Rou~r Zmi,~ Height Deformation 0.82 0.1 ~ 0.79~ 23 30 10 45 A r "~ 0.1~ Omi,~ ri 6 ~ 16mm Table 1. Geometric and deformation parameters of phantom. We use the same element and fitting technique as the real data between initial and deformed times. Our fitting criterion, the tag error, is the magnitude of the vector e, which was used to calculate the SPAMM forces. Fitting can be controlled manually or automatically by increasing the SPAMM strength and boundary strength as the model deforms in order to overcome the stiffness of the elements. All RMS tag errors were reduced to < 0.2ram while the average errors in nodal positions were all reduced to less than 5%. Results are shown in Fig. 6. The color map, seen here and in other figures as grey-scale without any loss of information, is that of the twisting parameter.

185 Fig. 6. Validation results: fitting method was applied to a computational, cylindrical phantom. (a) Initial configuration and (b) deformed shape with a color plot of twisting parameter on inner and outer walls. 5.2 Reconstructed Normal RV Wall Motion We applied our method to a set of 3 orthogonal 1D tagged MRI images of normal ventricles. The tag errors during the model fitting stages are shown in Table 2. Although we capture the motion of the LV in addition to the RV, we concentrate on results for RV deformation which concur with previous studies and show the possibility of providing more detail. For example, one group of researchers measured systolic long axis displacement of tag points in MRI images [17]. Unlike their reporting of average values in certain regions, we can display the spatial and temporal distribution with our model. These researchers reported that the maximum motion (about 25ram) occured in the basal region of the RV free wall. The color map in Fig. 7 is a plot of absolute axial displacement along the LV long-axis which was prescribed during imaging). This figure shows maximum displacement of about 19mm at the RV base and in the RV outflow tract. Similar to observations by [3], we found significant displacement of the free wall towards the septum and septal wall thickening which contributed to a smaller RV cavity. In a cross-sectional view, the most significant contraction is that of the posterior basal portion towards the outflow portion. Fig. 7 also shows the paths of cardiac material points located at the centers of the elements. Note that these points are different from those used during model fitting. It is apparent that the initial displacement of these points is greater than their final displacement. We also attempted to plot twisting on the RV wall. The apparently large concentration of twisting in the free wall shown in Fig. 8 results from a large forward displacement in a portion of the wall that is near the model axes rather than an actual twist. Thus, global twisting may not be a helpful parameter to use for the RV because of its non-circular cross-section. However, differences in twisting (i.e., between walls or between levels) may be good indicators of regional deformation.

186 Fig. 7. Color plot of axial displacement on RV as it deforms through 4 time intervals. The paths of RV mid-wall points are shown from time 1 (white) to time 5 (red) It is too early to conduct a more thorough comparison of our results with those found in literature. At this point we demonstrate the capability of the fitting technique to capture major trends in the deformation of the normal RV and identify areas needing improvement. We plan to validate the model with data from several patients and to identify which parameters derived from the reconstruction are clinically useful. 6 Conclusion We have developed a novel approach which fits a geometrically complex biven- tricular model to data from 1D tagged MR] images. In this paper, we validated the fitting technique on a computational phantom with predefined geometry and deformations. The results presented here show possibly useful clinical informa- tion about RV contraction. Through this work, we learned that the complicated RV wall motion may be best described with regional deformation parameters, i.e., strains. An in-depth study of these deformations should lead to greater knowledge of normal and abnormal heart wall motion. Our preliminary results of reconstructed normal RV wall motion demonstrate the plausibility of using this method in a clinical setting.

187 RMS Errors (ram) During Model Fitting View and tag orientation Time 1-2 Time 2-3 Time 3-4 Time 4-5 SA-Vertical 0.235 SA-Horizontal 0.244 LA-Horizontal 0.269 0.270 0.274 0.278 0.364 0.351 0.281 0.393 0.361 0.335 Table 2. Error between model material points and SPAMM stripes during model fitting for each time interval (SA = Short - axis, LA = Long - axis) Fig. 8. Color plot of twisting parameter on RV free wall 7 Acknowledgements This research has been funded by grants from the NSF (Career Award 96-24604), the Whitaker Foundation, and the NIH/NLM (LM06638-01). References 1. L. Axel, D. Bloomgarden, C.N. Chang, D. Kraitchman, and A.A. Young. SPAM- MVU: a program for the analysis of dynamic tagged MRI. In Proceedings of the Soc. of Magnetic Resonance in Medicine 12th Annual Meeting, page 724, 1993. 2. L. Axel and L. Dougherty. Heart wall motion: Improved method of spatial modu- lation of magnetization for MR imaging. Radiology, 272:349-50, 1989. 3. D. Barnard and J.S. Alpert. Right ventricular function in health and disease. Curt Probl Cardiol, 12:422-29, 1987. 4. C.J. Chuong, M.S. Sacks, G. Templeton, F. Schwiep, and Jr. R.L. Johnson. Re- gional deformation and contractile function in canine right ventricular free wall. Am J Physiol, 260:H1224-1235, 1991. 5. E. Haber, D.N. Metaxas, and L. Axel. Three-dimensional geometric modeling of cardiac right and left ventricles. In Biomedical Engineering Soe Annual Meeting, San Diego, California, 1997. 6. J.W. Hurst. Atlas of the Heart. Gower Medical, New York, 1988. 7. A.D. McCulloch and J.H. Omens. Non-homogeneous analysis of three-dimensional transmural finite deformation in canine ventricular myocardium. Y Biomechanics, 24:539-48, 1991. 8. D.N. Metaxas. Physics-based deformable models: applications to computer vision, graphics, and medical imaging. Kluwer Academic Publishers, Cambridge, 1996. 9. C.C. Moore, W.G. O'Dell, E.R. McVeigh, and E.A. Zerhouni. Calculation of three dimensional left ventricular strains from biplanar tagged MR images. Journal of Magnetic Resonance Imaging, 2:165-75, 1992.

188 10. P.M.F. Nielsen, I.J. Le Grice, and B.H. Smaill P.J. Hunter. Mathematical model of geometry and fibrous structure of the heart. J Appl Physiol, 260:H1365-H1378, 1991. 11. W.G. O'Dell, C.C. Moore, W.C. Hunter, E.A. Zerhouni, and E.R. McVeigh. 3- dimensional myocardial deformations: calculation with displacement field fitting to tagged MR images. Radiology, 195:829-35, 1996. 12. J. Park, D. Metaxas, and L. Axel. Analysis of left ventricular wall motion based on volumetric deformable models and MRI-SPAMM. Med Image Analysis J, 1(1):53- 71, 1996. 13. J.S. Pirolo, S.J. Bresina, L.L. Creswell, K.W. Myers, B.A. Szabo, M.W. Vannier, and M.K. Pasque. Mathematical three-dimensional solid modeling of biventricular geometry. Annals of Biomed Eng, 21:199-219, 1993. 14. R.A. Raines, M.M. LeWinter, and J.W. Covell. Regional shortening patterns in canine right ventricle. Am J Physiol, 231:H1395-1400, 1976. 15. M.S. Sacks, C.J. Chuong, G.H. Templeton, and R. Peshock. In vivo 3-d recon- struction and geometric characterization of the right ventricular free wall. Annals of Biomed Eng, 21:263-275, 1993. 16. W.P. Santamore, G.D. Meier, and A.A. Bove. Contractile function in canine right ventricle. Am J Physiol, 239:H794-804, 1980. 17. M. Stuber, E. Nagel, S.E. Fischer, M.B. Scheidegger, and P. Boesiger. Systolic long axis contraction of the human myocardium. In Proceedings of the SMR, 1995. 18. L.K. Waldman, J.J. Allen, R.S. Pavelec, and A.D. McCulloch. Distributed me- chanics of the canine right ventricle: effects of varying preload. J Biomechanics, 29(3):373-81, 1996. 19. A.A. Young and L. Axel. Non-rigid heart wall motion using MR tagging. In Proceedings of the IEEE Computer Society Conference on computer vision and pattern recognition, pages 399-404, Campaign, Illinois, 1992. 20. A.A. Young, Z.A. Fayad, and L. Axel. Right ventricular mid-wall surface motion and deformation using magnetic resonance tagging. J Appl Physiol, 271:2677-88, 1996. 21. A.A. Young, D. L. Kraitchman, L. Dougherty, and L. Axel. Trackiag and finite element analysis of stripe deformation in magnetic resonance tagging. IEEE Trans on Med Imaging, 14(3):413-21, September 1995. 22. O.C. Zienkiewicz and R.L. Taylor. The Finite Element Method. McGraw-Hill, New York, fourth edition, 1989.