- Email: [email protected]
Assessment and visualization of the curvature of the left ventricle from 3D medical images
Computerized Medical Imaging and Graphics. Printed in the U.S.A. All i-i&s reserved.
Vol. 11. Nos. 4/S, pp. 251-262,
0895-61 I l/93 $6.00 + .OO Copyright 0 1993 Pmgamon Press Ltd.
1993
ASSESSMENT AND VISUALIZATION OF THE CURVATURE OF THE LEFT VENTRICLE FROM 3D MEDICAL IMAGES Denis Friboulet*+, Isabelle E. Magnin*, Chrikophe Mathieu*, Andreas Pommert*, and Karl H. Hoehne* *URA CNRS 12 16, Institut National des Sciences appliqdes, Bit. 502, 6962 1 Villeurbanne Cedex, France *Institute of Mathematics and Computer Science in Medicine, Pav. 70, UKE, 2000 Hamburg 20, Germany (Received25 March 1993)
Abstract-We address the problem of using curvature features to assess the three-dimensional (3D) motion of the left ventricle. The adequacy of this approach depends on the actual characteristics of the curvature of the left ventricle and particularly on the spatial and temporal stability of these features. From experimental data, we compute the distribution of the Gaussian curvature over the surface of the left ventricle by using an iterative procedure. The results are visualized in 3D through a voxel-based surface rendering technique. We show that the Gaussian curvature remains stable along the cardiac cycle. This curvature feature could thus provide a reliable basis for further 3D motion analysis of the left ventricle. Key Words: Left ventricle, Curvature estimation, 3D Nonrigid motion, 3D image processing, 3D visualisation
wall is approximated by a smooth surface locally characterized through curvature features such as the Gaussian curvature, principal curvatures or strain energy. The interest of these features rely on their invariance under translation or rotation i.e. rigid motion. For instance, D. Goldgof (7) uses the Gaussian curvature as curvature feature and derives a motion field under the assumption of homothetic motion. As homothetic motion can only be a rough approximation of the actual LV motion (8,9), he then suggests to use conformal motion (i.e. a more general, preserving angles, nonrigid motion type) as a new assumption. The model of the LV used by A. Amini (6) is a thin-plate loosely constrained to stretch. He uses the potential (or strain) energy of the plate as curvature feature and proposes to recover a motion field by minimizing the bending energy and deviation from conformal motion. Minimizing the bending energy allows to favor correspondence between points representing significant shape tokens (i.e., intuitively points associated with high curvature and hence high strain energy) while minimizing the deviation from conformal motion provides correspondence in regions with less significant shape features. The adequacy of such approaches heavily depends on the actual characteristics of the curvature of the LV and on their potential variations with respect to the deformation of the LV. The selection of a proper curvature feature or the assessment of its stability could be efficiently guided by the visual inspection of the spatial distribution of the curvature of the LV along
INTRODUCIION
The quantification of the local motion of the left ventricle (LV) from medical images sequences can be a very useful diagnostic tool in assessing noninvasively the heart function. Such a tool would be for instance helpful for evaluating the severity of various cardiac diseases like infart or &hernia. Inferring the local dynamic behaviour of the LV from image sequence is a complex task because the LV experiences a nonrigid motion along the cardiac cycle and the data (i.e., the medical images) do not provide any explicit temporal correspondences information. Several studies tackle the correspondence problem by tracking surgically implanted radiopaque markers (1, 2), or opacified bifurcation points on the coronary artery tree (3). Inherent shortcoming of these methods is that they provide a sparse correspondence information and therefore do not allow to easily derive a dense motion field over the entire wall of the LV. Some authors (4, 5) proposed to overcome this problem by interpolating the surface of the LV and the displacement field from the set of tracked points. Various investigators (6, 7, 4) recently suggested to use a curvature-based approach for estimating the nonrigid motion of the LV. In this approach, the LV
+Correspondence should he addressedto Denis Friboulet, URA CNRS 12 16, lnstitut National des Sciences appliqukes, Bgt. 502, 6962 I Villeurbanne Cedex, France. 251
258
ComputerizedMedical Imaging and Graphics
the cardiac cycle. As far as we know, no experimental data nor results have been yet presented in that direction. This is precisely the topic of the present paper which deals with the computation and the three-dimensional (3D). visualization of curvature features extracted from an actual in vivo LV. The paper is organized as follows: the next section is devoted to curvature computation. We propose herein to infer the 3D curvature from the discrete representation of the surface of the LV by using an iterative relaxation scheme. The 3D display method used in this study to visualize the curvature distribution follows. It relies on voxel-based surface rendering and texture mapping techniques. The obtained experimental results are presented in the section headed, “Experimental Results.” We use here the 3D visualization tool to illustrate the behavior of the curvature computation method as well as the distribution of the Gaussian curvature. The implications of the obtained results toward the assessment of the left ventricular motion are discussed. CURVATURE
July-Cktober/l993,
Volume 17, Numbers 4/5
any direction is the same as S ( 11). The normal to this paraboloid at Q is Np = (-au - bu, -bu - cv, I)
(2)
Assuming that the original data represent a discrete approximation of S, the problem consists therefore to estimate of the above-defined paraboloid at each data point P. Let us call Qi a neighbor of P. From ( 1) each Qi = (Ui, u,, Wi)should ideally satisfy Cil = auf + 2bu,v, + cvf - 2wi = 0.
(3)
At each neighbor Qi, the gradient of the original gray-level 3D image (see “Experimental Results”) can be expressed as Gi = (x,, yi, 1) and provides an estimate of the vector normal to the surface. For each Qi, we should thus ideally have C,, = au, + bv, + x, = 0
(4)
Cl, = but + CVi+ y, = 0.
(5)
The values of the parameters a, b and c can then be determined from eqs, (3), (4), and (5) by minimizing
COMPUTATIONS (6)
The curvature of the LV was computed from the original medical images by using a scheme initially designed by Sander ( 10). This method proceeds as follows: (i) In a first stage, an initial estimate of the local differential structure is computed. This is done by fitting, in the least square sense, a paraboloid to the neighborhood of each point of the sampled surface of the LV; (ii) The second stage allows to improve the overall coherence of the curvature estimates. It is based on an iterative refinement of the initial estimates. A. Initial estimates Consider a smooth surface, S, a point P belonging to S and the unit vector, N, normal to Sat P. We define at Pa local coordinate system, Rp, based on N and two arbitrary unit vectors spanning the tangent plane. Let Q be a point of S belonging to the neighborhood of P. Approximating S through the truncated Taylor expansion of a bivariate function allows to express Q in RP as
au, VI = (UPv, w(u, w, where w(u, v) = $(a#* + 26UV+ cv2),
(1)
and a, b, c are three real coefficients. We note that eq. (1) is the equation of a paraboloid containing P and having normal N at P. It can be easily shown that the normal curvature of this paraboloid in
where the qj are weighting factors of eqs. (3), (4), and (5). The coefficient a, b and care the standard weighted least square solution of eq. (6). The principal curvatures K1, K2 and directions VK,, VK2 at P are the eigenvalues and eigenvectors of the Hessian matrix of the so-obtained local parametrization. The Gaussian curvature K, and the mean curvature, H, are given by K = K,K2; H = $(K, + K2).
B. Rejinement of the estimates
As observed by Sander (10) and as demonstrated by the experiments shown in “Experimental Results,” the initial estimates of the curvature yield poor results. The refinement step is necessary to improve the consistency of the initial estimates. At this level, each point Qi lying in the neighborhood of P is associated with a paraboloid PA(Qi). From each Qi, we may thus compute the values of the normal vector (Ni), the principal curvatures (Kli and Kzi), and the principal directions (VKli and VKzi), associated with P when P is considered as a point of PA(Qi). A new estimate of the local surface at P can then be determined by computing the best fit of this set of parameters (in the least square sense). The residual error of the fit performed at P, rp, provides a measure of the quality of the new local surface. This refinement step can ob-
3D images of the left ventricle ??D.
viously be iterated. Calling @ the overall residual associated with the Kh iteration, the refinement process stops when & - @-I < t/z, where th is the threshold defining the convergence level of the algorithm. 3D VISUALIZATION OF THE CURVATURE DISTRIBUTION It is very difficult to interprete the results only from the numerical data produced by the previously described algorithm. Our concern being the assessment of the curvature of the LV in the context of motion analysis, the understanding of these results implies an efficient representation of the computed curvature features over the surface of the LV as well as the changes of this distribution. Stated as a scientific visualization problem, this requirement means that we have to visualize the distribution of a scalar parameter over a discrete surface. Various method have been described in order to visualize the surface of organs from 3D discrete medical data sets [i.e., typically from CTor MR volumic images ( 12)]. These methods are commonly refered as voxelbased surface rendering in the computer graphics community and are well-adapted to our problem requirements. The method we actually applied proceeds as follows: (i) The surface of the LV is projected with a ray-casting algorithm and shaded through the gray level gradient method, using the Phong’s illumination model ( 13). The shading conveys the depth information and allows to render the anatomical details of the surface of the LV; and (ii) The scalar distribution to be represented (i.e., the Gaussian curvature for instance) is mapped onto the projected surface through a color scale. This operation can be done by adapting the texture mapping techniques commonly used in computer graphics ( 14) to our problem. Texture mapping consists usually in perturbating one of the parameters of the illumination model according to the two-dimensional (2D) map of an arbitrary parameter. The texture map is in our case provided by the computed 3D spatial distribution of the curvature feature and is used to compute the color of each shaded pixel. The various color maps we used for this study are based on the HLS color space (14) where H is the hue axis, L the luminance axis and S the saturation axis. All the visualizations presented in this study were performed at the Institute for Mathematics and Computer Science in Medecine (IMDM, University of Hamburg) by using the VOXEL-MAN 3D volume rendering software developed at this Institute. EXPERIMENTAL
RESULTS
A. Data acquisition and segmentation The cardiac data were obtained on the dynamic spatial reconstructor (DSR) (i.e., a high speed CT scan-
FFmOULET
et al.
259
ner) providing high resolution 3D images of the heart (15). The acquisitions were performed on a dog. The spatial resolution was about 0.9 mm (isotropic) and 18 volume images were obtained along a cardiac cycle. The contrast between the LV and the surrounding muscle was enhanced through a Roentgen contrast agent injected prior to the scan. The segmentation of LV is based on the application of a 3D gradient type edge detection filter (17). The response of the filter in the x, y and z directions provides an estimation of the 3 components of the gradient of the 3D image. The points considered as belonging to the surface of the LV can then be selected according to the following criteria: (i) They are extrema of the gradient magnitude; (ii) They form a connected component for the 26neighborhood connectivity. The 3 components of the gradient at the detected points provide the estimate of the normal to the surface of the LV required for the curvature computations. B. Curvature computation Figure 1 shows the spatial distribution of the sign of the Gaussian curvature (K) after the initial estimation stage and along the successive iterations of the refinement stage of the curvature computation algorithm. The sign of K was chosen to illustrate the behavior of the algorithm because this quantity is the product of the two principal curvatures and is therefore potentially more sensitive to noise (16). The distribution of the sign of K can be seen as a classification of a surface: the areas where K is positive correspond to elliptic regions whereas the areas related to negative values of K are associated with hyperbolic regions. The lack of coherence of the results provided by the initial step can be observed: It would be, for instance, difficult to deduce a classification of the LV surface into homogeneous and stable regions. This is probably due to the noise present in the original images (quantification and acquisition device noise ( 18)). Regarding the refinement stage, it can be observed that the variations of the distribution of Kare important for the first iterations (about 1 to 5) and that the distribution tends to stabilize as the number of iterations increases (few changes for the iterations 5 to 15). It should also be noted how the consistency of the surface classification is improved through the successive iterations. Figure 1. qualitatively demonstrates the convergence of the algorithm toward a stable and improved solution. This can be numerically confirmed by the plot of the relative variation of the overall residual depicted in Fig. 2. C. Distribution of the Gaussian curvature The distribution of the Gaussian curvature K is visualized in Fig. 3. The distribution of K at sampling
Computerized Medical Imaging and Graphics
260
Fig. 1. Spatial distribution of the sign of Gaussian curvature (K) over the surface of the LV after the initial estimation and along the successive iterations of the curvature computation algorithm.
time 1 on the free wall of the LV appears to be structurally rich: the right papillary muscle area is associated with an elliptic region (K > 0) exhibiting rather high positive values of K. This region is surrounded by an hyperbolic region (K -C0) where K takes on local negative extrema. These extrema are distributed according to a ring shaped pattern. The left papillary muscle area may be associated with the same distribution charac50
Algorithm
8
2
3
4
5
6
Convergence
7
8
9 101112131415
Number of Iterations Fig. 2. Relative variation of the overall residual as a function of the number of iterations of the curvature computation algorithm (refinement stage).
July-October/ 1993, Volume 17, Numbers 4/5
Fig. 3. Distribution of the Gaussian curvature (K) over the free wall (top) and the septal wall (bottom) of the LV at three sampling times along the cardiac cycle (times 1, 2 and 9).
teristics (elliptic region surrounded by an hyperbolic one) but without well-defined extrema. The right and left papillary muscles are separated by an elliptic bridge with relatively low values of K. Quite logically, the apical area of the LV is associated with an elliptical region where K reaches its maximum positive values. The septal wall of the LV is not structurally so rich: It corresponds to a distribution which varies smoothly between the high positive values at the apex and the high negative values associated with the junction of the main arteries. The results obtained for time 2 are very close to the results obtained for time 1 and would lead to the same description in terms of Gaussian curvature distribution. This means that the distribution of K is stable between time 1 and 2. This result is coherent with the small deformation experienced by the LV between these two times. Time 9 allows to test the behavior of the distribution of K when important motion takes place. It is interesting to observe from Fig. 3 that even though the deformation of the LV is important, the distribution of K remains structurally stable. The same structural stability may be observed for the septal wall of the LV: The distribution remains smooth and no new well-defined curvature extrema is induced by the deformation. It should also be noted that the contraction experienced by the LV during systole is reflected by the higher absolute values of K as compared with time 1 or 2.
3D images of the left ventricle ??D. FRIBCNJLET et al.
D. Discussion In the absence of any well-defined a priori deformation model of the LV, the computation of a displacement field based on curvature information will necessarily rely, to some extent, on the correspondence of points or regions exhibiting characteristic curvature properties (i.e., regions of high curvature, extrema of curvature, regions sharing the same classification properties. . .). The adequacy of this type of approach will thus rely on the assumption that the anatomical structure of the LV provides at each sampling time “enough” significant characteristic regions or points of curvature (i.e. the distribution of these shape tokens have to be sufficiently dense and well-defined to lead to consistent correspondences) and that these regions or points remain stable with respect to the deformation of the LV. As far as stability is concerned, the experimental results (observed stability of the distribution of K through time) suggest that such curvature-based approaches may be instrumental in the case of the LV. Regarding the anatomy of the LV, the results indicate that whereas the free wall of the LV provides rich and dense curvature characteristics, the structure of the septal wall appears to be poor: the derivation of correspondences in this region will be therefore difficult if only based on a shape token tracking process. CONCLUSION We described in this paper a set of tools allowing to compute and visualize the 3D distribution of the curvature of the left ventricle. The curvature computations were performed through an iterative relaxation scheme and the visualization of curvature features was obtained through voxel-based surface rendering and texture mapping techniques. The method was applied to CT volumetric data and illustrated through the visualization of the distribution of the Gaussian curvature over the surface of the left ventricle for various sampling times along the cardiac cycle. The obtained experimental results showed the stability of the Gaussian curvature spatial distribution with respect to the nonrigid motion of the left ventricle. Our current work aims at designing a method allowing to derive a motion field based on the curvature features of the left ventricle.
Acknowledgment-This work was supported by a grant of the National Institute for Health and Medical Research (INSERM). The authors wish to thank Pr R. A. Robb and D. Hanson for providing the DSR CT data.
261
SUMMARY The quantification of the local motion of the left ventricle (LV) from medical images sequences can be a very useful diagnostic tool in assessing non-invasively the heart function. Such a tool would be for instance helpful for evaluating the severity of various cardiac diseases like infart or ischemia. Inferring the local dynamic behaviour of the LV from image sequence is a complex task because the LV experiences a nonrigid motion along the cardiac cycle and the data (i.e., the medical images) do not provide any explicit temporal correspondences information. Various investigators recently suggested to use a curvature-based approach for estimating the nonrigid motion of the LV. In this approach, the LV wall is approximated by a smooth surface locally characterized through curvature features such as the Gaussian curvature, principal curvatures or strain energy. The interest of these features rely on their invariance under translation or rotation (i.e., rigid motion). The adequacy of such curvature-based approaches thus heavily depends on the actual characteristics of the curvature of the LV and on their potential variations with respect to the deformation of the LV. The selection of a proper curvature feature or the assessment of its stability could be efficiently guided by the visual inspection of the spatial distribution of the curvature of the LV along the cardiac cycle. As far as we know, no experimental data nor results have been yet presented in that direction. This is precisely the topic of the present paper which deals with the computation and the 3D visualization of curvature features extracted from an actual in vivo LV. We propose in this paper to infer the 3D curvature of the surface of the LV from its discrete representation by using an iterative relaxation scheme. The 3D display method used to visualize the curvature distribution relies on voxel-based surface rendering and texture mapping techniques. We use here the 3D visualization tool to illustrate the behavior of the curvature computation method as well as the distribution of the Gaussian curvature obtained from experimental data. The implications of the obtained experimental results toward the assessment of the left ventricular motion are discussed. REFERENCES 1. Meier, G.D.; Bove, A.A.; Santamore, W.P.; Lynch, P.R. Contractile function in canine right ventricle. Am. J. Physiol. 239: 794-804; 1980. 2. Potel, M.J.; Mac Kay, S.A.; Rubin, J.M.; Aisen, A.M.; Sayre, R.E. Three-dimensional left ventricular wall motion in man coordinate systems for representing wall movement direction. Invest. Radiol. 19:499-509; 1984.
262
Computerized Medical Imaging and Graphics
3. Kim, H.C.; Min, B.G.M.M.; Lee, L.; Seo, J.D.; Lee, Y.M.; Han, MC. Estimation of local cardiac wall deformation and regional wall stresses from biplane coronary cineangiograms. IEEE Trans. Biomed. Eng. BME-32503-512; 1985. 4. Mishm, S.K.; Goldgof, D.B. Motion analysis and modeling of epicardial surfaces from points and line correspondences. In: Proc. IEEE Workshop on visual motion, Princeton, N.J.; 1991: 300-305. 5. Chen, C.W.; Huang, T.S.; Arrott, M. Analysis and visualisation of heart motion. In: Proc. SPIE Biomedical image processing II. 1991:231-242. 6. Amini, A.A.; Owen, R.L.; Anandan, P.; Duncan, J.S. Non-rigid motion for tracking the left-ventricular wall, In: Colchester, A.C.F.; Hawkes, D.J., eds. Lecture notes in computer science: Information processing in medical imaging, Berlin: SpringerVerlag; 1991:343-357. 7. Kambhamettu, C.; Goldgof, D.B. Toward finding point corresnondences in non&id motion. In: Proc. 7th Scandinavian Conf. on Image analysis. Aalborg, Denmark, 199 I :1 I26- 1133. 8. Friboulet, D.; Magnin, I.E.; Revel, D. Assessment of a model for overall left ventricular three-dimensional motion from MRI data. Int. J. Cardiac Imaging 8:175-190; 1992. 9. Friboulet, D.; Magnin, I.E. A 3D model ofthe global deformation ofa non-rigid body. In: Prokopek, T.; Viergever, M., eds. NATO AS1 Series F, Medical images: Formation, Handling and Evaluation, Berlin: Springer-Verlag; 1991:305-324. IO. Sander. P.T.: Zucker. S.W. Inferring surface trace and differential structure from 3D images. IEEE Pattern Analysis and Machine Intelligence. 12(9):833-854; 1990. 11. Lipschutz, M.M. Schaum’s outline of theory and problems of differential geometry, New York: Mcgraw-Hill; 1969. 12. Tiede, U.; Hohne, K.H.; Bomans, M.; Pommert, A.; Riemer, M.; Wiebecke, G. Investigation of medical 3D-rendering algorithms. IEEE computer graphics and applications. 2:4 l-53: 1990.
July-Cktober/l993,
Volume 17, Numbers 4/5
13. Hohne, K.H.; Bomans, M.; Pommert, A.; Riemer, M.; Schiers, C.; Tiede, U.; Wiebecke, G. 3D visualization of tomographic volume data using the generalized voxel-model. ViiuaI computer. 6( 1):28-36; 1990. 14. Foley, J.D.; Van Dam, A.; Feiner, S.K.; Hughes, J.F. Computer graphics: Principles and Practice. Reading: Addison-Wesley; 1990. 15. Robb, R.A. X-Ray computed tomography: advanced systems and applications in biomedical research and diagnosis, In: Threedimensional Biomedical imaging, Boca Raton: CRC Press; 1985: 107-162. 16. Besl, P.J.; Jain, R.C. Invariant surface characteristics for 3D object recognition in range images. Computer vision, graphics, and image processing. 33:33-80; 1986. 17. Monga, 0.; Deriche, R.; Rocchisani, J.-M., 3D edge detection using recursive filtering. CGVIP: image understanding. 53( 1): 76-87; 1991. 18. Flynn, P.J.; Jam, A.K. On reliable curvature estimation. In: Proc. Int. Conf. on Computer Vision and Pattern Recognition. 1989: 110-116.
About the Author-DENIS FRIBOULETreceived in 1984 the engineering degree (electrical engineering) and in 1990 the PLD degree (biomedical engineering), both from the National Institute for Applied Sciences of Lyon (F). He worked from 1990 to 199 1at the Institute of Mathematics and Computer Science in Medicine at the University of Hamburg (FRG). He works since 1991 in a research unit of the National Center for Scientific Research (URA CNRS 1216) at the National Institute for Applied Sciences of Lyon, where he currently holds the position of assistant-Professor. His primary research interests are 3D cardiac data segmentation, representation and motion analysis.
Vol. 11. Nos. 4/S, pp. 251-262,
0895-61 I l/93 $6.00 + .OO Copyright 0 1993 Pmgamon Press Ltd.
1993
ASSESSMENT AND VISUALIZATION OF THE CURVATURE OF THE LEFT VENTRICLE FROM 3D MEDICAL IMAGES Denis Friboulet*+, Isabelle E. Magnin*, Chrikophe Mathieu*, Andreas Pommert*, and Karl H. Hoehne* *URA CNRS 12 16, Institut National des Sciences appliqdes, Bit. 502, 6962 1 Villeurbanne Cedex, France *Institute of Mathematics and Computer Science in Medicine, Pav. 70, UKE, 2000 Hamburg 20, Germany (Received25 March 1993)
Abstract-We address the problem of using curvature features to assess the three-dimensional (3D) motion of the left ventricle. The adequacy of this approach depends on the actual characteristics of the curvature of the left ventricle and particularly on the spatial and temporal stability of these features. From experimental data, we compute the distribution of the Gaussian curvature over the surface of the left ventricle by using an iterative procedure. The results are visualized in 3D through a voxel-based surface rendering technique. We show that the Gaussian curvature remains stable along the cardiac cycle. This curvature feature could thus provide a reliable basis for further 3D motion analysis of the left ventricle. Key Words: Left ventricle, Curvature estimation, 3D Nonrigid motion, 3D image processing, 3D visualisation
wall is approximated by a smooth surface locally characterized through curvature features such as the Gaussian curvature, principal curvatures or strain energy. The interest of these features rely on their invariance under translation or rotation i.e. rigid motion. For instance, D. Goldgof (7) uses the Gaussian curvature as curvature feature and derives a motion field under the assumption of homothetic motion. As homothetic motion can only be a rough approximation of the actual LV motion (8,9), he then suggests to use conformal motion (i.e. a more general, preserving angles, nonrigid motion type) as a new assumption. The model of the LV used by A. Amini (6) is a thin-plate loosely constrained to stretch. He uses the potential (or strain) energy of the plate as curvature feature and proposes to recover a motion field by minimizing the bending energy and deviation from conformal motion. Minimizing the bending energy allows to favor correspondence between points representing significant shape tokens (i.e., intuitively points associated with high curvature and hence high strain energy) while minimizing the deviation from conformal motion provides correspondence in regions with less significant shape features. The adequacy of such approaches heavily depends on the actual characteristics of the curvature of the LV and on their potential variations with respect to the deformation of the LV. The selection of a proper curvature feature or the assessment of its stability could be efficiently guided by the visual inspection of the spatial distribution of the curvature of the LV along
INTRODUCIION
The quantification of the local motion of the left ventricle (LV) from medical images sequences can be a very useful diagnostic tool in assessing noninvasively the heart function. Such a tool would be for instance helpful for evaluating the severity of various cardiac diseases like infart or &hernia. Inferring the local dynamic behaviour of the LV from image sequence is a complex task because the LV experiences a nonrigid motion along the cardiac cycle and the data (i.e., the medical images) do not provide any explicit temporal correspondences information. Several studies tackle the correspondence problem by tracking surgically implanted radiopaque markers (1, 2), or opacified bifurcation points on the coronary artery tree (3). Inherent shortcoming of these methods is that they provide a sparse correspondence information and therefore do not allow to easily derive a dense motion field over the entire wall of the LV. Some authors (4, 5) proposed to overcome this problem by interpolating the surface of the LV and the displacement field from the set of tracked points. Various investigators (6, 7, 4) recently suggested to use a curvature-based approach for estimating the nonrigid motion of the LV. In this approach, the LV
+Correspondence should he addressedto Denis Friboulet, URA CNRS 12 16, lnstitut National des Sciences appliqukes, Bgt. 502, 6962 I Villeurbanne Cedex, France. 251
258
ComputerizedMedical Imaging and Graphics
the cardiac cycle. As far as we know, no experimental data nor results have been yet presented in that direction. This is precisely the topic of the present paper which deals with the computation and the three-dimensional (3D). visualization of curvature features extracted from an actual in vivo LV. The paper is organized as follows: the next section is devoted to curvature computation. We propose herein to infer the 3D curvature from the discrete representation of the surface of the LV by using an iterative relaxation scheme. The 3D display method used in this study to visualize the curvature distribution follows. It relies on voxel-based surface rendering and texture mapping techniques. The obtained experimental results are presented in the section headed, “Experimental Results.” We use here the 3D visualization tool to illustrate the behavior of the curvature computation method as well as the distribution of the Gaussian curvature. The implications of the obtained results toward the assessment of the left ventricular motion are discussed. CURVATURE
July-Cktober/l993,
Volume 17, Numbers 4/5
any direction is the same as S ( 11). The normal to this paraboloid at Q is Np = (-au - bu, -bu - cv, I)
(2)
Assuming that the original data represent a discrete approximation of S, the problem consists therefore to estimate of the above-defined paraboloid at each data point P. Let us call Qi a neighbor of P. From ( 1) each Qi = (Ui, u,, Wi)should ideally satisfy Cil = auf + 2bu,v, + cvf - 2wi = 0.
(3)
At each neighbor Qi, the gradient of the original gray-level 3D image (see “Experimental Results”) can be expressed as Gi = (x,, yi, 1) and provides an estimate of the vector normal to the surface. For each Qi, we should thus ideally have C,, = au, + bv, + x, = 0
(4)
Cl, = but + CVi+ y, = 0.
(5)
The values of the parameters a, b and c can then be determined from eqs, (3), (4), and (5) by minimizing
COMPUTATIONS (6)
The curvature of the LV was computed from the original medical images by using a scheme initially designed by Sander ( 10). This method proceeds as follows: (i) In a first stage, an initial estimate of the local differential structure is computed. This is done by fitting, in the least square sense, a paraboloid to the neighborhood of each point of the sampled surface of the LV; (ii) The second stage allows to improve the overall coherence of the curvature estimates. It is based on an iterative refinement of the initial estimates. A. Initial estimates Consider a smooth surface, S, a point P belonging to S and the unit vector, N, normal to Sat P. We define at Pa local coordinate system, Rp, based on N and two arbitrary unit vectors spanning the tangent plane. Let Q be a point of S belonging to the neighborhood of P. Approximating S through the truncated Taylor expansion of a bivariate function allows to express Q in RP as
au, VI = (UPv, w(u, w, where w(u, v) = $(a#* + 26UV+ cv2),
(1)
and a, b, c are three real coefficients. We note that eq. (1) is the equation of a paraboloid containing P and having normal N at P. It can be easily shown that the normal curvature of this paraboloid in
where the qj are weighting factors of eqs. (3), (4), and (5). The coefficient a, b and care the standard weighted least square solution of eq. (6). The principal curvatures K1, K2 and directions VK,, VK2 at P are the eigenvalues and eigenvectors of the Hessian matrix of the so-obtained local parametrization. The Gaussian curvature K, and the mean curvature, H, are given by K = K,K2; H = $(K, + K2).
B. Rejinement of the estimates
As observed by Sander (10) and as demonstrated by the experiments shown in “Experimental Results,” the initial estimates of the curvature yield poor results. The refinement step is necessary to improve the consistency of the initial estimates. At this level, each point Qi lying in the neighborhood of P is associated with a paraboloid PA(Qi). From each Qi, we may thus compute the values of the normal vector (Ni), the principal curvatures (Kli and Kzi), and the principal directions (VKli and VKzi), associated with P when P is considered as a point of PA(Qi). A new estimate of the local surface at P can then be determined by computing the best fit of this set of parameters (in the least square sense). The residual error of the fit performed at P, rp, provides a measure of the quality of the new local surface. This refinement step can ob-
3D images of the left ventricle ??D.
viously be iterated. Calling @ the overall residual associated with the Kh iteration, the refinement process stops when & - @-I < t/z, where th is the threshold defining the convergence level of the algorithm. 3D VISUALIZATION OF THE CURVATURE DISTRIBUTION It is very difficult to interprete the results only from the numerical data produced by the previously described algorithm. Our concern being the assessment of the curvature of the LV in the context of motion analysis, the understanding of these results implies an efficient representation of the computed curvature features over the surface of the LV as well as the changes of this distribution. Stated as a scientific visualization problem, this requirement means that we have to visualize the distribution of a scalar parameter over a discrete surface. Various method have been described in order to visualize the surface of organs from 3D discrete medical data sets [i.e., typically from CTor MR volumic images ( 12)]. These methods are commonly refered as voxelbased surface rendering in the computer graphics community and are well-adapted to our problem requirements. The method we actually applied proceeds as follows: (i) The surface of the LV is projected with a ray-casting algorithm and shaded through the gray level gradient method, using the Phong’s illumination model ( 13). The shading conveys the depth information and allows to render the anatomical details of the surface of the LV; and (ii) The scalar distribution to be represented (i.e., the Gaussian curvature for instance) is mapped onto the projected surface through a color scale. This operation can be done by adapting the texture mapping techniques commonly used in computer graphics ( 14) to our problem. Texture mapping consists usually in perturbating one of the parameters of the illumination model according to the two-dimensional (2D) map of an arbitrary parameter. The texture map is in our case provided by the computed 3D spatial distribution of the curvature feature and is used to compute the color of each shaded pixel. The various color maps we used for this study are based on the HLS color space (14) where H is the hue axis, L the luminance axis and S the saturation axis. All the visualizations presented in this study were performed at the Institute for Mathematics and Computer Science in Medecine (IMDM, University of Hamburg) by using the VOXEL-MAN 3D volume rendering software developed at this Institute. EXPERIMENTAL
RESULTS
A. Data acquisition and segmentation The cardiac data were obtained on the dynamic spatial reconstructor (DSR) (i.e., a high speed CT scan-
FFmOULET
et al.
259
ner) providing high resolution 3D images of the heart (15). The acquisitions were performed on a dog. The spatial resolution was about 0.9 mm (isotropic) and 18 volume images were obtained along a cardiac cycle. The contrast between the LV and the surrounding muscle was enhanced through a Roentgen contrast agent injected prior to the scan. The segmentation of LV is based on the application of a 3D gradient type edge detection filter (17). The response of the filter in the x, y and z directions provides an estimation of the 3 components of the gradient of the 3D image. The points considered as belonging to the surface of the LV can then be selected according to the following criteria: (i) They are extrema of the gradient magnitude; (ii) They form a connected component for the 26neighborhood connectivity. The 3 components of the gradient at the detected points provide the estimate of the normal to the surface of the LV required for the curvature computations. B. Curvature computation Figure 1 shows the spatial distribution of the sign of the Gaussian curvature (K) after the initial estimation stage and along the successive iterations of the refinement stage of the curvature computation algorithm. The sign of K was chosen to illustrate the behavior of the algorithm because this quantity is the product of the two principal curvatures and is therefore potentially more sensitive to noise (16). The distribution of the sign of K can be seen as a classification of a surface: the areas where K is positive correspond to elliptic regions whereas the areas related to negative values of K are associated with hyperbolic regions. The lack of coherence of the results provided by the initial step can be observed: It would be, for instance, difficult to deduce a classification of the LV surface into homogeneous and stable regions. This is probably due to the noise present in the original images (quantification and acquisition device noise ( 18)). Regarding the refinement stage, it can be observed that the variations of the distribution of Kare important for the first iterations (about 1 to 5) and that the distribution tends to stabilize as the number of iterations increases (few changes for the iterations 5 to 15). It should also be noted how the consistency of the surface classification is improved through the successive iterations. Figure 1. qualitatively demonstrates the convergence of the algorithm toward a stable and improved solution. This can be numerically confirmed by the plot of the relative variation of the overall residual depicted in Fig. 2. C. Distribution of the Gaussian curvature The distribution of the Gaussian curvature K is visualized in Fig. 3. The distribution of K at sampling
Computerized Medical Imaging and Graphics
260
Fig. 1. Spatial distribution of the sign of Gaussian curvature (K) over the surface of the LV after the initial estimation and along the successive iterations of the curvature computation algorithm.
time 1 on the free wall of the LV appears to be structurally rich: the right papillary muscle area is associated with an elliptic region (K > 0) exhibiting rather high positive values of K. This region is surrounded by an hyperbolic region (K -C0) where K takes on local negative extrema. These extrema are distributed according to a ring shaped pattern. The left papillary muscle area may be associated with the same distribution charac50
Algorithm
8
2
3
4
5
6
Convergence
7
8
9 101112131415
Number of Iterations Fig. 2. Relative variation of the overall residual as a function of the number of iterations of the curvature computation algorithm (refinement stage).
July-October/ 1993, Volume 17, Numbers 4/5
Fig. 3. Distribution of the Gaussian curvature (K) over the free wall (top) and the septal wall (bottom) of the LV at three sampling times along the cardiac cycle (times 1, 2 and 9).
teristics (elliptic region surrounded by an hyperbolic one) but without well-defined extrema. The right and left papillary muscles are separated by an elliptic bridge with relatively low values of K. Quite logically, the apical area of the LV is associated with an elliptical region where K reaches its maximum positive values. The septal wall of the LV is not structurally so rich: It corresponds to a distribution which varies smoothly between the high positive values at the apex and the high negative values associated with the junction of the main arteries. The results obtained for time 2 are very close to the results obtained for time 1 and would lead to the same description in terms of Gaussian curvature distribution. This means that the distribution of K is stable between time 1 and 2. This result is coherent with the small deformation experienced by the LV between these two times. Time 9 allows to test the behavior of the distribution of K when important motion takes place. It is interesting to observe from Fig. 3 that even though the deformation of the LV is important, the distribution of K remains structurally stable. The same structural stability may be observed for the septal wall of the LV: The distribution remains smooth and no new well-defined curvature extrema is induced by the deformation. It should also be noted that the contraction experienced by the LV during systole is reflected by the higher absolute values of K as compared with time 1 or 2.
3D images of the left ventricle ??D. FRIBCNJLET et al.
D. Discussion In the absence of any well-defined a priori deformation model of the LV, the computation of a displacement field based on curvature information will necessarily rely, to some extent, on the correspondence of points or regions exhibiting characteristic curvature properties (i.e., regions of high curvature, extrema of curvature, regions sharing the same classification properties. . .). The adequacy of this type of approach will thus rely on the assumption that the anatomical structure of the LV provides at each sampling time “enough” significant characteristic regions or points of curvature (i.e. the distribution of these shape tokens have to be sufficiently dense and well-defined to lead to consistent correspondences) and that these regions or points remain stable with respect to the deformation of the LV. As far as stability is concerned, the experimental results (observed stability of the distribution of K through time) suggest that such curvature-based approaches may be instrumental in the case of the LV. Regarding the anatomy of the LV, the results indicate that whereas the free wall of the LV provides rich and dense curvature characteristics, the structure of the septal wall appears to be poor: the derivation of correspondences in this region will be therefore difficult if only based on a shape token tracking process. CONCLUSION We described in this paper a set of tools allowing to compute and visualize the 3D distribution of the curvature of the left ventricle. The curvature computations were performed through an iterative relaxation scheme and the visualization of curvature features was obtained through voxel-based surface rendering and texture mapping techniques. The method was applied to CT volumetric data and illustrated through the visualization of the distribution of the Gaussian curvature over the surface of the left ventricle for various sampling times along the cardiac cycle. The obtained experimental results showed the stability of the Gaussian curvature spatial distribution with respect to the nonrigid motion of the left ventricle. Our current work aims at designing a method allowing to derive a motion field based on the curvature features of the left ventricle.
Acknowledgment-This work was supported by a grant of the National Institute for Health and Medical Research (INSERM). The authors wish to thank Pr R. A. Robb and D. Hanson for providing the DSR CT data.
261
SUMMARY The quantification of the local motion of the left ventricle (LV) from medical images sequences can be a very useful diagnostic tool in assessing non-invasively the heart function. Such a tool would be for instance helpful for evaluating the severity of various cardiac diseases like infart or ischemia. Inferring the local dynamic behaviour of the LV from image sequence is a complex task because the LV experiences a nonrigid motion along the cardiac cycle and the data (i.e., the medical images) do not provide any explicit temporal correspondences information. Various investigators recently suggested to use a curvature-based approach for estimating the nonrigid motion of the LV. In this approach, the LV wall is approximated by a smooth surface locally characterized through curvature features such as the Gaussian curvature, principal curvatures or strain energy. The interest of these features rely on their invariance under translation or rotation (i.e., rigid motion). The adequacy of such curvature-based approaches thus heavily depends on the actual characteristics of the curvature of the LV and on their potential variations with respect to the deformation of the LV. The selection of a proper curvature feature or the assessment of its stability could be efficiently guided by the visual inspection of the spatial distribution of the curvature of the LV along the cardiac cycle. As far as we know, no experimental data nor results have been yet presented in that direction. This is precisely the topic of the present paper which deals with the computation and the 3D visualization of curvature features extracted from an actual in vivo LV. We propose in this paper to infer the 3D curvature of the surface of the LV from its discrete representation by using an iterative relaxation scheme. The 3D display method used to visualize the curvature distribution relies on voxel-based surface rendering and texture mapping techniques. We use here the 3D visualization tool to illustrate the behavior of the curvature computation method as well as the distribution of the Gaussian curvature obtained from experimental data. The implications of the obtained experimental results toward the assessment of the left ventricular motion are discussed. REFERENCES 1. Meier, G.D.; Bove, A.A.; Santamore, W.P.; Lynch, P.R. Contractile function in canine right ventricle. Am. J. Physiol. 239: 794-804; 1980. 2. Potel, M.J.; Mac Kay, S.A.; Rubin, J.M.; Aisen, A.M.; Sayre, R.E. Three-dimensional left ventricular wall motion in man coordinate systems for representing wall movement direction. Invest. Radiol. 19:499-509; 1984.
262
Computerized Medical Imaging and Graphics
3. Kim, H.C.; Min, B.G.M.M.; Lee, L.; Seo, J.D.; Lee, Y.M.; Han, MC. Estimation of local cardiac wall deformation and regional wall stresses from biplane coronary cineangiograms. IEEE Trans. Biomed. Eng. BME-32503-512; 1985. 4. Mishm, S.K.; Goldgof, D.B. Motion analysis and modeling of epicardial surfaces from points and line correspondences. In: Proc. IEEE Workshop on visual motion, Princeton, N.J.; 1991: 300-305. 5. Chen, C.W.; Huang, T.S.; Arrott, M. Analysis and visualisation of heart motion. In: Proc. SPIE Biomedical image processing II. 1991:231-242. 6. Amini, A.A.; Owen, R.L.; Anandan, P.; Duncan, J.S. Non-rigid motion for tracking the left-ventricular wall, In: Colchester, A.C.F.; Hawkes, D.J., eds. Lecture notes in computer science: Information processing in medical imaging, Berlin: SpringerVerlag; 1991:343-357. 7. Kambhamettu, C.; Goldgof, D.B. Toward finding point corresnondences in non&id motion. In: Proc. 7th Scandinavian Conf. on Image analysis. Aalborg, Denmark, 199 I :1 I26- 1133. 8. Friboulet, D.; Magnin, I.E.; Revel, D. Assessment of a model for overall left ventricular three-dimensional motion from MRI data. Int. J. Cardiac Imaging 8:175-190; 1992. 9. Friboulet, D.; Magnin, I.E. A 3D model ofthe global deformation ofa non-rigid body. In: Prokopek, T.; Viergever, M., eds. NATO AS1 Series F, Medical images: Formation, Handling and Evaluation, Berlin: Springer-Verlag; 1991:305-324. IO. Sander. P.T.: Zucker. S.W. Inferring surface trace and differential structure from 3D images. IEEE Pattern Analysis and Machine Intelligence. 12(9):833-854; 1990. 11. Lipschutz, M.M. Schaum’s outline of theory and problems of differential geometry, New York: Mcgraw-Hill; 1969. 12. Tiede, U.; Hohne, K.H.; Bomans, M.; Pommert, A.; Riemer, M.; Wiebecke, G. Investigation of medical 3D-rendering algorithms. IEEE computer graphics and applications. 2:4 l-53: 1990.
July-Cktober/l993,
Volume 17, Numbers 4/5
13. Hohne, K.H.; Bomans, M.; Pommert, A.; Riemer, M.; Schiers, C.; Tiede, U.; Wiebecke, G. 3D visualization of tomographic volume data using the generalized voxel-model. ViiuaI computer. 6( 1):28-36; 1990. 14. Foley, J.D.; Van Dam, A.; Feiner, S.K.; Hughes, J.F. Computer graphics: Principles and Practice. Reading: Addison-Wesley; 1990. 15. Robb, R.A. X-Ray computed tomography: advanced systems and applications in biomedical research and diagnosis, In: Threedimensional Biomedical imaging, Boca Raton: CRC Press; 1985: 107-162. 16. Besl, P.J.; Jain, R.C. Invariant surface characteristics for 3D object recognition in range images. Computer vision, graphics, and image processing. 33:33-80; 1986. 17. Monga, 0.; Deriche, R.; Rocchisani, J.-M., 3D edge detection using recursive filtering. CGVIP: image understanding. 53( 1): 76-87; 1991. 18. Flynn, P.J.; Jam, A.K. On reliable curvature estimation. In: Proc. Int. Conf. on Computer Vision and Pattern Recognition. 1989: 110-116.
About the Author-DENIS FRIBOULETreceived in 1984 the engineering degree (electrical engineering) and in 1990 the PLD degree (biomedical engineering), both from the National Institute for Applied Sciences of Lyon (F). He worked from 1990 to 199 1at the Institute of Mathematics and Computer Science in Medicine at the University of Hamburg (FRG). He works since 1991 in a research unit of the National Center for Scientific Research (URA CNRS 1216) at the National Institute for Applied Sciences of Lyon, where he currently holds the position of assistant-Professor. His primary research interests are 3D cardiac data segmentation, representation and motion analysis.