- Email: [email protected]
3D geometrical features of anatomic structures: the example of the ulna and radius bones
Computerized Medical Imaging and Graphics, Printed in the U.S.A. All rights reserved.
Vol. 11, Nos. 415, pp. 381-386,
I
1993 copylight
0
0895-61 l/93 $6.00 + .OO 1993 Pergamon press Ltd.
3D GEOMETRICAL FEATURES OF ANATOMIC STRUCTURES: THE EXAMPLE OF THE ULNA AND RADIUS BONES Christian Roux*+,ValCrie Burdin *, Wolfgang Schiitte-Felsche *, and Christian L,&vre* *Groupe Traitement d’Images, TC16com Bretagne-Ecole Nationale Supkrieure des TCltcommunications de Bretagne, Technop6le de Brest-Iroise, B.P 832-29285 Brest Cedex-France, *Service MCary, Centre Hospitalier Rtgional Universitaire de Brest, Rue Augustin Morvan-29279 Brest Cedex-France (Received 25 March, 1993)
Abstract-This paper addresses two important issues of threedimensional (3D) natural shape analysis. The first concerns the segmentation of these shapes into 3D primitives and the second deals with their geometric characterisation by means of intrinsic features. Restricting ourselves to 3D objects formed from two-dimensional (2D) cross-sectional adjacent slices (i.e., long bones), we present a method of contextual classification and results concerning the ulnar bone. With respect to the second issue, a curvature-torsion analysis of the medulbu axis is introduced and the results are represented by means of torsion images. A first application to ulnar structures is given. It is shown that radius bones have zero torsion and that their curvatures can be studied in 2D. A preliminary statistical study of the 2D curvature of the radius is also presented. Key Words: 3D segmentation, 3D primitives, 3D geometrical analysis, Long bone structures
INTRODUCTION
Results using ulnar bone structures will show the usefulness of the method. We will show that the radius bone has a zero torsion and that it can be consequently studied in 2D, for example by way of radiological projections. A preliminary statistical study of the 2D curvatures of the radius will be presented. Finally, the conclusion will propose some comments on the method and further directions for research.
In this paper we address the problem of the description of three-dimensional (3D) shapes, restricting ourselves to 3D objects formed from two-dimensional (2D) cross-sectional adjacent slices. In two previous publications ( 1, 2 ), the reconstruction of 3D shapes using
Fourier Descriptors (FDs) has been shown on both a theoretical and practical basis. This paper presents further applications of this generic framework to the geometrical analysis of long bone structures. Two issues will be addressed, the first one concerning the problem of the segmentation of 3D shapes into geometrical primitives, the second one related to the geometrical characterisation of 3D medullar structures for prosthesis applications. In the next section, we will briefly recall the principle and the main features of the description method. A shape distance will be introduced and a contextual classification method will be presented that allows a segmentation of the 3D structures into 3D primitives. Section III is devoted to the definition and computation of pertinent 3D geometrical features when considering medullar structures for the purpose of prosthesis design.
PRIMITIVE
EXTRACTION
3 D shape description method
Because of its typical geometry, a long bone can be represented by a set of closed parallel planar contours. This is equivalent to having a cylindrical parameterization which is discrete along the Oz axis. At the height z, the curve is represented by a polar parameterization (p(s) , d(s)) with respect to its centroid, where s is the arc length (2). Ghorbel has shown in (3) that this representation can be reduced to p(s) , without information loss, under two assumptions: the arc length s is constant between two sampled points and the curve is assumed to be star-shaped with respect to its centroid. The real periodic signal Log( p(s)) is expanded into a discrete Fourier series with ak coefficients. The { (lk } are used to build a family of FDs { zk } that are
+ Correspondence should be addressed to C. Roux, Groupe Traitement d’Images, T&corn Bretagne-Ecole Nationale Supkieure des TMcommunications de Bretagne, TechnopBle de Brest-Iroise, B.P832-29285 Brest Cedex-France.
invariant with respect to scaling, rotation, translation and the starting point (2) 381
.
Computerized Medical Imaging and Graphics
382
Ik(z)
=
~k(Z)l~l(Z)lk’fork= ” 1,2,3 ,..., [G(Z)lk
and
+co p=O.l.
(1)
The set {Z,(z), k P 0 } allows us to reconstruct each contour without any significant loss of shape information because of its completeness and stability properties ( 1). We build the surface from triangular facets according to the method presented in (4). Distances between planar shapes Each contour i can be represented in the space of the FDs by the family {I:}. We have shown in ( 1) that the first r FDs are needed to represent the contour, that is to say the following vector of 43’: I’ = (I\, r:, . . . ) I:).
(2)
In this space, we define the inner product
(I’,ZZ)=
&.I/ .r;,
(3)
i=l
where Z denotes the complex conjugate of Z and the pi are related to the variances of the descriptors according to
July-October/ 1993, Volume 17, Numbers 4/5
teresting feature of this approach is that it uses the continuity property of the descriptors, but a threshold for the decision has to be defined on an arbitrary basis. This drawback vanishes when trying to group the different shapes into natural clusters according to a distance between them. However, spatial continuity has to be checked during the classification process. The method we propose here is based on a dynamic clustering approach, namely a dynamic cluster algorithm. The distance used to measure the dissimilarity between two shapes is the above-mentioned shape distance. The algorithm can be presented as follows: Let the set of samples (the slices) to be classified be {I’, 1
IL
=
u cr.
k-0
Step 0. K-l
which leads to the following norm
an initial partition II, = U Cfl is chosen, and the k=O
class centres are computed. A distance can also be defined between two contours using
Step m. K-l
partition II, = U Cr k=O
d,(Z’, Z2) = ))I’ - 1211 =
&,-)I; ,=I
-I:[‘.
(6)
The weighting factors allow us to take into account the fact that the descriptors do not have the same dynamic range. This distance is used to measure dissimilarities between planar shapes regardless of their size. For example, the distance between two contours with identical shape and different size would be equal to zero. This distance is a shape distance. Primitive extraction using slice class.ijication Because each contour is represented in the space of the IDS, the problem of segmenting the entire 3D structure into geometrical primitives can be solved using a classification approach which takes into account certain spatial constraints. In a first attempt to do this a simple decision criterion based on the concept of dominant descriptors has been proposed ( 5). One in-
Stepm+Stepm+l. for k = 1 to K - 1 let ir = I’ be the lower frontier of class C$’ First case while I’ is closer, in the sense of the shape distance, to the centre of CK, , then Zj E Ck”_: j=j+ 1 j=j+jo+ir+L
=Zj+JO
Second case while ZJ is closer, in the sense of the shape distance, to the centre of Cp, then I’ E CT” jzj-1 j=j_jo+ip+l
=IJ-Jo
Compute the centres of the new classes, End.
3D geometrical features of anatomic structures
if the centres at step mo are equal to the centres at stepmo1. In this algorithm, the centre of one class is the sample which is the closest (in the sense of the shape distance) to the computed mean value of all the samples of the class. The above algorithm has the interesting property of processing the shapes (the slices) in a predefined order, and therefore preserving spatial continuity. However, the initialisation step can be improved by taking into account the fact that the assumption of uniformly spatially distributed centre is not fulfilled. It is indeed more appropriate to use some a priori information and to have a first rough estimation of the location of these centres. This can be done using the polar diagram (2) which represents the complete structure. Such a diagram whose numbered slices are normalised between 0 and 27r is depicted in Fig. 1(a). For each slice, represented by a radius in the diagram, the moduli of the first five FDs are plotted on five concentric circles (the smallest circle for FD number 1 and the largest for FD number 5). Descriptor 0, which determines the variation of the average radius of a given slice, has not been drawn in the diagrams since its values are much greater than the values of the other FDs. Figure 1(b) depicts a result of the segmentation of an ulna, whose polar diagram is shown in Figure 1(a).
383
C. Roux et al.
??
The number of classes has been set at 12; therefore 12 centres have to be chosen. Although this choice can be made automatically, the process has been carried out interactively in the example presented here. The segmentation leads to natural primitives that are similar to those usually anatomically considered (6). GEOMETRICAL
FEATURES
Introduction
The purpose of this section is to determine a continuous mathematical representation to describe the principal axis of long bones (e.g., ulna, radius as shown in Fig. 2). Therefore, we consider the axis of the medulla to represent this course and we propose an approach for evaluating intrinsic features such as torsion and curvature. These parameters are important pieces of information for preoperative examinations in medical surgery. 3 D geometrical analysis Data related to the internal morphology of the structure are acquired using an X-ray computer&d tomography scanner which provides a set of 2D longitudinal sections of long bones. This set is reorganised in order to obtain a set of cross-sectional slices which are more suitable to the description of this family of bones. After contour extraction and segmentation into three areas (air, bone, medulla), as described in (6))
250 (a) Fig. 1. (a) Polar diagram
of the first five descriptors. (b) An ulna and its segmentation frontiers of the classes (two points of view).
into 12 classes showing the
384
Computerized Medical Imaging and Graphics
July-Octoher/l993,
Volume 17, Numbers 4/5 N
EZ =
C (PI - fh(ti))‘>
(8)
because of its better average fit, which turns out to be more suitable for noisy data. To find the optimal approximation we put at zero the derivative of E2 with respect to the coefficients qi . The equation (9 ) gives a necessary condition for E2to be a minimum: dEz=O dqi
ix1
,...,
4.
Having done this for each dimension using the arc length 1as parameter t, a continuous representation of the bone axis is available. We can now apply the equations of Serret-Frenet ( 10) to calculate its curvature cs and torsion ws: c.y =
IIS(l)x S(OIIw = IIS(OII3 s
_
(S(t) x QO)* B(f) IIS
x S(t)112 .
(lo)
where X, R and 2 denote, the first, second and third derivatives of x with respect to t .
Fig. 2. A radius bone.
we obtain in the case of an ulna or a radius about 300 cross-sections each consisting of two closed discrete contours (bone and medulla). After determination of the centre of gravity of each medulla contour, the discrete course of the bone axis is given by a set of N 3D points Pi (Xi, yi, Zi ), i = 1, . . . , N represented by vectors Pi. To calculate the torsion and curvature we look for a continuous “smooth” function S(S,, S,,, S,) possessing three derivatives that represents these points. Because of their good approximation properties, the simple way to calculate coefficients and derivatives and the fact that they minimise curvature energy (7)) we have chosen cubic spline functions, based on B-splines. These functions without interior knots are of the form: s,(t) = i siNi,4(t),
(7)
i=l
where qi( qXi, qyi, q3i) denotes the coefficients of the spline and Ni,d( t) the linearly independent cubic B-splines ( 8 ) . To measure the quality of the approximation one can either use the maximum error or the integral square error (ISE) method. Despite the fact that the approximation minimising the maximum error gives a closer fit to the points (9)) we have chosen to minimise the integral square error:
3 D application The right-handed coordinate system of Frenet is defined by the intersection of the three perpendicular vectors, tangent, normal and binormal to the curve. To examine the torsion along the calculated axis we take into account the fact that, if the torsion is zero, all binormal vectors are colinear. We use this property to visualise the changes of torsion on a grey-level image. Each pixel (m , n) corresponds to the inner product of the binormal vectors at both points P,,, and P, on the curve. This gives the images, which will be called torsion images, shown in Fig. 3. The level 0 is associated with the pixels where the binormal vectors are colinear, the level 255 when the binormal vectors are orthogonal. Thus, the diagonal of such images is black. It can be seen in Fig. 3(b) that the torsion image of the helix exhibits a periodicity which corresponds to the vertical periodicity of this curve. One can see that the medial axis of the ulna (Fig. 3(c)) is not planar, but possesses a torsion at the beginning of the last third of its length. In contrast to this the medial axis of the radius (Fig. 3 (d)) is planar, except in a small area (about 5% of the considered bone length) at one end of the bone. 2 D analysis As shown in the previous section, the estimation of curvature in the case of the radius can be reduced to a 2D problem. This enables radiological projections to be used, which are more suitable in practice because they are easier and cheaper to obtain and reduce the patient’s exposure to X-rays. In the following, it will
(4
(b)
(4
Fig. 3. Torsion image of (from left to right) (a) a 3D curve representing a third order polynomial function with an abrupt torsion at the middle of the interval, (b) a helix, (c) an ulna bone, (d) a radius bone. be assumed that the radiological projections are the projections of the 3D structure in its curvature plane. A recent study ( 11) has shown the difficulties of using X-ray photographs in medical morphology: the inner contours of the medulla cannot be segmented and are therefore not suitable for supplying a set of data points. We propose to obtain these points by detecting the border opposite the radius (see Fig. 2), because the boundary of its central part is very close to the course of the medulla. A contour in an image can be seen as the locus of the pixels where there is a significant difference in light intensity between the left and right sides. We take this into account by emphasising the border with a Prewitt kernel. Examining the picture line by line, we are looking for the location of greater changes of light intensity. The detailed method for finding these contours can be found in (12). Regarding the typical structure of the radius it seems that there are no points where the curvature direction changes. To force the still unknown function to behave accordingly, we propose to minim& the ISE
with a polynomial function of degree 2. Once the polynomial has been differentiated, the curvature can be easily computed. The histogram of curvature for 3 bones [Figure 4 (b)] shows that one value of curvature dominates. However, as shown in Figure 4(a) the curvature is not constant, but increases when approaching the elbow (on the right hand side in Figure 4(a). A study of 24 different radius bones supplies the histogram of Figure 4 (c). The calculated values lie between 25 and 72, with a maximum at about 38 cm. It is also obvious that one cannot take a single value of curvature which fits all radius bones. CONCLUSION
The segmentation of 3D natural shapes is a difficult and still open problem ( 13 ). In the framework of long bones, we have presented a solution which controls the spatial continuity of the primitives; it is based on a description space formed by smart 2D Fourier Descriptors and gives natural geometrical primitives.
60 1
0
0.2
0.4 0.6 x/' (4
0.8
1
20
30 radius
40
50
of curvature 04
60 [cm]
70
20
30 roaius
40
50
of curvature
60
70
[cm]
(c)
Fig. 4.(a)Curvature along the medullar axis of 3 radius bones; (b) histogram of curvature for the same 3 radius bones; (c) histogram of curvature for 24 radius bones.
386
Computerized Medical Imaging and Graphics
Concerning the geometrical characterisation of 3D medullar structures, we have introduced torsion images that represent the regularity of the medial axis of long bones. Classical anatomical considerations on the radius bone have been confirmed by this 3D approach which shows that this bone has zero torsion. This conclusion, which is of great importance for prosthesis design applications, has been followed by a statistical analysis of the curvature of the radius in the curvature plane. A 3D primitive extraction and a curvature-torsion analysis of the medullar axis have been presented separately. However, it is obvious that they are not independent of each other. And there is no doubt that new results will be obtained by taking into account this (almost) dual aspect in both topics. REFERENCES 1. Burdin, V.; Ghorbel, F.; de Bougrenet de la Tocnaye, J.L.; Roux,
2.
3.
4.
5.
6.
7
8 9. 10.
C. A Geometrical analysis for data compression of 3D anatomical structures. In: Laurent, P.J., Le MChaute, A., Schumaker, L.L. ( Eds), Curves and Surfaces, Academic Press, New York; 1991: 57-63. Burdin, V.; Ghorbel, F.; de Bougrenet de la Tocnaye, J.L.; Roux, C. A three-dimensional primitive extraction of long bones obtained from bi-dimensional Fourier descriptors. In: Pattern Recognition Letters, 13:213-2 17; 1992. Ghorbel, F. Vets une approche mathematique unifiee des aspects geometriques et statistiques de la reconnaissance de formes planes. Doctoral thesis; December 90; Rennes I. Burdin, V.; Lecornu, L.; de Bougrenet, J.L.; Roux, C. Description and Reconstruction of 3D Bone Structures using Fourier Descriptors. In: Proc. 13th annual International Conference of the IEEE MBS. Orlando; November 199I : 1133- I 134. Burdin, V.; Ghorbel, F.; de Bougrenet, J.L.; Roux, C. A complete and stable set of Fourier descriptors for invariant analysis and reconstruction of 3D objects. In: Proc. V EUSIPCO; Barcelona; September 1990:1655-1658. Burdin, V.; de Bougrenet de la Tocnaye, J.L.; Roux, C.; Lefevre, C. Une representation des structures osseuses tridimensionnelles par un diagramme polaire plan de descripteurs de Fourier. In: Innovation et Technologie en Biologie et Mtdecine: Vol. 12, n o 3:257-267; 1991. Gueziec, A.; Ayache, N. L&age et reconnaissance de courbes gauches bruitees. In: Proc. 8tme Congres Reconnaissance des fonnes et intelligence artificielle: Lvon (France ): 199 I :809-820. de Boor, C. A practical guide to Splines, Springer Verlag: 1976. Pavlidis, T. Algorithms for graphics and image processing, Springer Verlag: 1977. Ramis, E.; Dechamps, C.; Odoux, J. Cours de mathematiques sp&ziales. Masson: I98 I.
July-October/ 1993, Volume 17, Numbers 4/5 11. Rubin, P.J.; Leyvraz, P.F.; Heegard, J.H. Variations radiologiques
des parametres anatomiques du f&mu proximal en fonction de sa position en rotation. In: Revue de Chirurgie oF&o#&que, 75:209-2 15: 1989. 12. Schtitte-Felsche, W. Extraktion und At&se dreidimensionaler Kurven in Bildem der medizinischen D&&a&, Internatreport, GTI-ENST Bretame: A~til 1992. 13. Roux, C.; Coattie&, J.L. Some trends in 3D medical inrag@, Invited Paper. In: Proc. V EUSIPCO, Bamelo~ September 1990; 55-63.
About the Author-CHRISTIAN ROUX received the Agregation de Sciences Physiques from the Ecole Normale Sup&eure de Cachan (France) in 1978, the Ph.D. degree in Image Processing from the lnstitut National Polvtechniaue de Grenoble (France) in 1980 and the Habilitation a dihger des recherches en Sciences degree from the University of Rennes (France) in 1990. In 1978, he joined the LETICentre d’Etudes Nucleaires de Grenoble (France) as Research Assistant. He became Lecturer at the University Institute of Technology, Caracas (Venezuela)in 1981. In 1982, he was appointed Assistant Professor at the Ecole Nationale Sup&ieure des T&ommunications de Bretagne, Brest (France) and Head of the Image Processing Group. He became Professor of Image Processing and Pattern Recognition in 1987, and worked on various research projects in Medical Imaging and Computer Vision. From October 1992, he is a Visiting Professor at the Medical Image Processing Group, department of Radiology, University of Pennsylvania, Philadelphia. He served as member of the Technical Review Committee and Chairman of the track Imaging for the 14th Annual International Conference of the IEEE Enaineering in Medicine and Biology Society and as co-chairman of the-Satellite Symposium of the same Conference on 3D Advanced Image Processing in Medicine held in Rennes, France. He is also Associate Editor of the IEEE Transactions on Medical Imaging. His research interests include Image analysis and segmentation, Pattern recognition in 2D and 3D images, Image data fusion, 2D and 3D Image reconstruction. He has published more than 30 papers and holds one Patent. Christian Roux is member of SEE, AFCET, SFPT and IEEE. About the Author-vALI?RIE BURDIN received her Ph.D. degree in Signal Processing from the University of Rennes (France) in 1992. In 1988, she joined the Image Processing Group ofthe Ecole Nationale Sup&ieure des T&communications de Bretagne, Brest (France) where she works on the 3D Medical Imaging. Her research interests include the image segmentation, the analysis and the reconstruction of 3D objects. About the Author-WOLFGANG SCHCITTEFELSCHEreceived his Engineering degree from the Technical University of Munich, Germany, in 199 I. In the same year, he joined the Image Processing Group of the Ecole Nationale Sup&ieure des Tel&communications de Bretagne, Brest (France), where he worked on the 3D Medical Imaging. About the Author-CHRISTIAN LEF~VRE, M.D., works in the Department of Anatomy and Department of Orthopoedic and Traumatology Surgery of the University of Brest, France. His research interest is the locking nailing of the forearm.
Vol. 11, Nos. 415, pp. 381-386,
I
1993 copylight
0
0895-61 l/93 $6.00 + .OO 1993 Pergamon press Ltd.
3D GEOMETRICAL FEATURES OF ANATOMIC STRUCTURES: THE EXAMPLE OF THE ULNA AND RADIUS BONES Christian Roux*+,ValCrie Burdin *, Wolfgang Schiitte-Felsche *, and Christian L,&vre* *Groupe Traitement d’Images, TC16com Bretagne-Ecole Nationale Supkrieure des TCltcommunications de Bretagne, Technop6le de Brest-Iroise, B.P 832-29285 Brest Cedex-France, *Service MCary, Centre Hospitalier Rtgional Universitaire de Brest, Rue Augustin Morvan-29279 Brest Cedex-France (Received 25 March, 1993)
Abstract-This paper addresses two important issues of threedimensional (3D) natural shape analysis. The first concerns the segmentation of these shapes into 3D primitives and the second deals with their geometric characterisation by means of intrinsic features. Restricting ourselves to 3D objects formed from two-dimensional (2D) cross-sectional adjacent slices (i.e., long bones), we present a method of contextual classification and results concerning the ulnar bone. With respect to the second issue, a curvature-torsion analysis of the medulbu axis is introduced and the results are represented by means of torsion images. A first application to ulnar structures is given. It is shown that radius bones have zero torsion and that their curvatures can be studied in 2D. A preliminary statistical study of the 2D curvature of the radius is also presented. Key Words: 3D segmentation, 3D primitives, 3D geometrical analysis, Long bone structures
INTRODUCTION
Results using ulnar bone structures will show the usefulness of the method. We will show that the radius bone has a zero torsion and that it can be consequently studied in 2D, for example by way of radiological projections. A preliminary statistical study of the 2D curvatures of the radius will be presented. Finally, the conclusion will propose some comments on the method and further directions for research.
In this paper we address the problem of the description of three-dimensional (3D) shapes, restricting ourselves to 3D objects formed from two-dimensional (2D) cross-sectional adjacent slices. In two previous publications ( 1, 2 ), the reconstruction of 3D shapes using
Fourier Descriptors (FDs) has been shown on both a theoretical and practical basis. This paper presents further applications of this generic framework to the geometrical analysis of long bone structures. Two issues will be addressed, the first one concerning the problem of the segmentation of 3D shapes into geometrical primitives, the second one related to the geometrical characterisation of 3D medullar structures for prosthesis applications. In the next section, we will briefly recall the principle and the main features of the description method. A shape distance will be introduced and a contextual classification method will be presented that allows a segmentation of the 3D structures into 3D primitives. Section III is devoted to the definition and computation of pertinent 3D geometrical features when considering medullar structures for the purpose of prosthesis design.
PRIMITIVE
EXTRACTION
3 D shape description method
Because of its typical geometry, a long bone can be represented by a set of closed parallel planar contours. This is equivalent to having a cylindrical parameterization which is discrete along the Oz axis. At the height z, the curve is represented by a polar parameterization (p(s) , d(s)) with respect to its centroid, where s is the arc length (2). Ghorbel has shown in (3) that this representation can be reduced to p(s) , without information loss, under two assumptions: the arc length s is constant between two sampled points and the curve is assumed to be star-shaped with respect to its centroid. The real periodic signal Log( p(s)) is expanded into a discrete Fourier series with ak coefficients. The { (lk } are used to build a family of FDs { zk } that are
+ Correspondence should be addressed to C. Roux, Groupe Traitement d’Images, T&corn Bretagne-Ecole Nationale Supkieure des TMcommunications de Bretagne, TechnopBle de Brest-Iroise, B.P832-29285 Brest Cedex-France.
invariant with respect to scaling, rotation, translation and the starting point (2) 381
.
Computerized Medical Imaging and Graphics
382
Ik(z)
=
~k(Z)l~l(Z)lk’fork= ” 1,2,3 ,..., [G(Z)lk
and
+co p=O.l.
(1)
The set {Z,(z), k P 0 } allows us to reconstruct each contour without any significant loss of shape information because of its completeness and stability properties ( 1). We build the surface from triangular facets according to the method presented in (4). Distances between planar shapes Each contour i can be represented in the space of the FDs by the family {I:}. We have shown in ( 1) that the first r FDs are needed to represent the contour, that is to say the following vector of 43’: I’ = (I\, r:, . . . ) I:).
(2)
In this space, we define the inner product
(I’,ZZ)=
&.I/ .r;,
(3)
i=l
where Z denotes the complex conjugate of Z and the pi are related to the variances of the descriptors according to
July-October/ 1993, Volume 17, Numbers 4/5
teresting feature of this approach is that it uses the continuity property of the descriptors, but a threshold for the decision has to be defined on an arbitrary basis. This drawback vanishes when trying to group the different shapes into natural clusters according to a distance between them. However, spatial continuity has to be checked during the classification process. The method we propose here is based on a dynamic clustering approach, namely a dynamic cluster algorithm. The distance used to measure the dissimilarity between two shapes is the above-mentioned shape distance. The algorithm can be presented as follows: Let the set of samples (the slices) to be classified be {I’, 1
IL
=
u cr.
k-0
Step 0. K-l
which leads to the following norm
an initial partition II, = U Cfl is chosen, and the k=O
class centres are computed. A distance can also be defined between two contours using
Step m. K-l
partition II, = U Cr k=O
d,(Z’, Z2) = ))I’ - 1211 =
&,-)I; ,=I
-I:[‘.
(6)
The weighting factors allow us to take into account the fact that the descriptors do not have the same dynamic range. This distance is used to measure dissimilarities between planar shapes regardless of their size. For example, the distance between two contours with identical shape and different size would be equal to zero. This distance is a shape distance. Primitive extraction using slice class.ijication Because each contour is represented in the space of the IDS, the problem of segmenting the entire 3D structure into geometrical primitives can be solved using a classification approach which takes into account certain spatial constraints. In a first attempt to do this a simple decision criterion based on the concept of dominant descriptors has been proposed ( 5). One in-
Stepm+Stepm+l. for k = 1 to K - 1 let ir = I’ be the lower frontier of class C$’ First case while I’ is closer, in the sense of the shape distance, to the centre of CK, , then Zj E Ck”_: j=j+ 1 j=j+jo+ir+L
=Zj+JO
Second case while ZJ is closer, in the sense of the shape distance, to the centre of Cp, then I’ E CT” jzj-1 j=j_jo+ip+l
=IJ-Jo
Compute the centres of the new classes, End.
3D geometrical features of anatomic structures
if the centres at step mo are equal to the centres at stepmo1. In this algorithm, the centre of one class is the sample which is the closest (in the sense of the shape distance) to the computed mean value of all the samples of the class. The above algorithm has the interesting property of processing the shapes (the slices) in a predefined order, and therefore preserving spatial continuity. However, the initialisation step can be improved by taking into account the fact that the assumption of uniformly spatially distributed centre is not fulfilled. It is indeed more appropriate to use some a priori information and to have a first rough estimation of the location of these centres. This can be done using the polar diagram (2) which represents the complete structure. Such a diagram whose numbered slices are normalised between 0 and 27r is depicted in Fig. 1(a). For each slice, represented by a radius in the diagram, the moduli of the first five FDs are plotted on five concentric circles (the smallest circle for FD number 1 and the largest for FD number 5). Descriptor 0, which determines the variation of the average radius of a given slice, has not been drawn in the diagrams since its values are much greater than the values of the other FDs. Figure 1(b) depicts a result of the segmentation of an ulna, whose polar diagram is shown in Figure 1(a).
383
C. Roux et al.
??
The number of classes has been set at 12; therefore 12 centres have to be chosen. Although this choice can be made automatically, the process has been carried out interactively in the example presented here. The segmentation leads to natural primitives that are similar to those usually anatomically considered (6). GEOMETRICAL
FEATURES
Introduction
The purpose of this section is to determine a continuous mathematical representation to describe the principal axis of long bones (e.g., ulna, radius as shown in Fig. 2). Therefore, we consider the axis of the medulla to represent this course and we propose an approach for evaluating intrinsic features such as torsion and curvature. These parameters are important pieces of information for preoperative examinations in medical surgery. 3 D geometrical analysis Data related to the internal morphology of the structure are acquired using an X-ray computer&d tomography scanner which provides a set of 2D longitudinal sections of long bones. This set is reorganised in order to obtain a set of cross-sectional slices which are more suitable to the description of this family of bones. After contour extraction and segmentation into three areas (air, bone, medulla), as described in (6))
250 (a) Fig. 1. (a) Polar diagram
of the first five descriptors. (b) An ulna and its segmentation frontiers of the classes (two points of view).
into 12 classes showing the
384
Computerized Medical Imaging and Graphics
July-Octoher/l993,
Volume 17, Numbers 4/5 N
EZ =
C (PI - fh(ti))‘>
(8)
because of its better average fit, which turns out to be more suitable for noisy data. To find the optimal approximation we put at zero the derivative of E2 with respect to the coefficients qi . The equation (9 ) gives a necessary condition for E2to be a minimum: dEz=O dqi
ix1
,...,
4.
Having done this for each dimension using the arc length 1as parameter t, a continuous representation of the bone axis is available. We can now apply the equations of Serret-Frenet ( 10) to calculate its curvature cs and torsion ws: c.y =
IIS(l)x S(OIIw = IIS(OII3 s
_
(S(t) x QO)* B(f) IIS
x S(t)112 .
(lo)
where X, R and 2 denote, the first, second and third derivatives of x with respect to t .
Fig. 2. A radius bone.
we obtain in the case of an ulna or a radius about 300 cross-sections each consisting of two closed discrete contours (bone and medulla). After determination of the centre of gravity of each medulla contour, the discrete course of the bone axis is given by a set of N 3D points Pi (Xi, yi, Zi ), i = 1, . . . , N represented by vectors Pi. To calculate the torsion and curvature we look for a continuous “smooth” function S(S,, S,,, S,) possessing three derivatives that represents these points. Because of their good approximation properties, the simple way to calculate coefficients and derivatives and the fact that they minimise curvature energy (7)) we have chosen cubic spline functions, based on B-splines. These functions without interior knots are of the form: s,(t) = i siNi,4(t),
(7)
i=l
where qi( qXi, qyi, q3i) denotes the coefficients of the spline and Ni,d( t) the linearly independent cubic B-splines ( 8 ) . To measure the quality of the approximation one can either use the maximum error or the integral square error (ISE) method. Despite the fact that the approximation minimising the maximum error gives a closer fit to the points (9)) we have chosen to minimise the integral square error:
3 D application The right-handed coordinate system of Frenet is defined by the intersection of the three perpendicular vectors, tangent, normal and binormal to the curve. To examine the torsion along the calculated axis we take into account the fact that, if the torsion is zero, all binormal vectors are colinear. We use this property to visualise the changes of torsion on a grey-level image. Each pixel (m , n) corresponds to the inner product of the binormal vectors at both points P,,, and P, on the curve. This gives the images, which will be called torsion images, shown in Fig. 3. The level 0 is associated with the pixels where the binormal vectors are colinear, the level 255 when the binormal vectors are orthogonal. Thus, the diagonal of such images is black. It can be seen in Fig. 3(b) that the torsion image of the helix exhibits a periodicity which corresponds to the vertical periodicity of this curve. One can see that the medial axis of the ulna (Fig. 3(c)) is not planar, but possesses a torsion at the beginning of the last third of its length. In contrast to this the medial axis of the radius (Fig. 3 (d)) is planar, except in a small area (about 5% of the considered bone length) at one end of the bone. 2 D analysis As shown in the previous section, the estimation of curvature in the case of the radius can be reduced to a 2D problem. This enables radiological projections to be used, which are more suitable in practice because they are easier and cheaper to obtain and reduce the patient’s exposure to X-rays. In the following, it will
(4
(b)
(4
Fig. 3. Torsion image of (from left to right) (a) a 3D curve representing a third order polynomial function with an abrupt torsion at the middle of the interval, (b) a helix, (c) an ulna bone, (d) a radius bone. be assumed that the radiological projections are the projections of the 3D structure in its curvature plane. A recent study ( 11) has shown the difficulties of using X-ray photographs in medical morphology: the inner contours of the medulla cannot be segmented and are therefore not suitable for supplying a set of data points. We propose to obtain these points by detecting the border opposite the radius (see Fig. 2), because the boundary of its central part is very close to the course of the medulla. A contour in an image can be seen as the locus of the pixels where there is a significant difference in light intensity between the left and right sides. We take this into account by emphasising the border with a Prewitt kernel. Examining the picture line by line, we are looking for the location of greater changes of light intensity. The detailed method for finding these contours can be found in (12). Regarding the typical structure of the radius it seems that there are no points where the curvature direction changes. To force the still unknown function to behave accordingly, we propose to minim& the ISE
with a polynomial function of degree 2. Once the polynomial has been differentiated, the curvature can be easily computed. The histogram of curvature for 3 bones [Figure 4 (b)] shows that one value of curvature dominates. However, as shown in Figure 4(a) the curvature is not constant, but increases when approaching the elbow (on the right hand side in Figure 4(a). A study of 24 different radius bones supplies the histogram of Figure 4 (c). The calculated values lie between 25 and 72, with a maximum at about 38 cm. It is also obvious that one cannot take a single value of curvature which fits all radius bones. CONCLUSION
The segmentation of 3D natural shapes is a difficult and still open problem ( 13 ). In the framework of long bones, we have presented a solution which controls the spatial continuity of the primitives; it is based on a description space formed by smart 2D Fourier Descriptors and gives natural geometrical primitives.
60 1
0
0.2
0.4 0.6 x/' (4
0.8
1
20
30 radius
40
50
of curvature 04
60 [cm]
70
20
30 roaius
40
50
of curvature
60
70
[cm]
(c)
Fig. 4.(a)Curvature along the medullar axis of 3 radius bones; (b) histogram of curvature for the same 3 radius bones; (c) histogram of curvature for 24 radius bones.
386
Computerized Medical Imaging and Graphics
Concerning the geometrical characterisation of 3D medullar structures, we have introduced torsion images that represent the regularity of the medial axis of long bones. Classical anatomical considerations on the radius bone have been confirmed by this 3D approach which shows that this bone has zero torsion. This conclusion, which is of great importance for prosthesis design applications, has been followed by a statistical analysis of the curvature of the radius in the curvature plane. A 3D primitive extraction and a curvature-torsion analysis of the medullar axis have been presented separately. However, it is obvious that they are not independent of each other. And there is no doubt that new results will be obtained by taking into account this (almost) dual aspect in both topics. REFERENCES 1. Burdin, V.; Ghorbel, F.; de Bougrenet de la Tocnaye, J.L.; Roux,
2.
3.
4.
5.
6.
7
8 9. 10.
C. A Geometrical analysis for data compression of 3D anatomical structures. In: Laurent, P.J., Le MChaute, A., Schumaker, L.L. ( Eds), Curves and Surfaces, Academic Press, New York; 1991: 57-63. Burdin, V.; Ghorbel, F.; de Bougrenet de la Tocnaye, J.L.; Roux, C. A three-dimensional primitive extraction of long bones obtained from bi-dimensional Fourier descriptors. In: Pattern Recognition Letters, 13:213-2 17; 1992. Ghorbel, F. Vets une approche mathematique unifiee des aspects geometriques et statistiques de la reconnaissance de formes planes. Doctoral thesis; December 90; Rennes I. Burdin, V.; Lecornu, L.; de Bougrenet, J.L.; Roux, C. Description and Reconstruction of 3D Bone Structures using Fourier Descriptors. In: Proc. 13th annual International Conference of the IEEE MBS. Orlando; November 199I : 1133- I 134. Burdin, V.; Ghorbel, F.; de Bougrenet, J.L.; Roux, C. A complete and stable set of Fourier descriptors for invariant analysis and reconstruction of 3D objects. In: Proc. V EUSIPCO; Barcelona; September 1990:1655-1658. Burdin, V.; de Bougrenet de la Tocnaye, J.L.; Roux, C.; Lefevre, C. Une representation des structures osseuses tridimensionnelles par un diagramme polaire plan de descripteurs de Fourier. In: Innovation et Technologie en Biologie et Mtdecine: Vol. 12, n o 3:257-267; 1991. Gueziec, A.; Ayache, N. L&age et reconnaissance de courbes gauches bruitees. In: Proc. 8tme Congres Reconnaissance des fonnes et intelligence artificielle: Lvon (France ): 199 I :809-820. de Boor, C. A practical guide to Splines, Springer Verlag: 1976. Pavlidis, T. Algorithms for graphics and image processing, Springer Verlag: 1977. Ramis, E.; Dechamps, C.; Odoux, J. Cours de mathematiques sp&ziales. Masson: I98 I.
July-October/ 1993, Volume 17, Numbers 4/5 11. Rubin, P.J.; Leyvraz, P.F.; Heegard, J.H. Variations radiologiques
des parametres anatomiques du f&mu proximal en fonction de sa position en rotation. In: Revue de Chirurgie oF&o#&que, 75:209-2 15: 1989. 12. Schtitte-Felsche, W. Extraktion und At&se dreidimensionaler Kurven in Bildem der medizinischen D&&a&, Internatreport, GTI-ENST Bretame: A~til 1992. 13. Roux, C.; Coattie&, J.L. Some trends in 3D medical inrag@, Invited Paper. In: Proc. V EUSIPCO, Bamelo~ September 1990; 55-63.
About the Author-CHRISTIAN ROUX received the Agregation de Sciences Physiques from the Ecole Normale Sup&eure de Cachan (France) in 1978, the Ph.D. degree in Image Processing from the lnstitut National Polvtechniaue de Grenoble (France) in 1980 and the Habilitation a dihger des recherches en Sciences degree from the University of Rennes (France) in 1990. In 1978, he joined the LETICentre d’Etudes Nucleaires de Grenoble (France) as Research Assistant. He became Lecturer at the University Institute of Technology, Caracas (Venezuela)in 1981. In 1982, he was appointed Assistant Professor at the Ecole Nationale Sup&ieure des T&ommunications de Bretagne, Brest (France) and Head of the Image Processing Group. He became Professor of Image Processing and Pattern Recognition in 1987, and worked on various research projects in Medical Imaging and Computer Vision. From October 1992, he is a Visiting Professor at the Medical Image Processing Group, department of Radiology, University of Pennsylvania, Philadelphia. He served as member of the Technical Review Committee and Chairman of the track Imaging for the 14th Annual International Conference of the IEEE Enaineering in Medicine and Biology Society and as co-chairman of the-Satellite Symposium of the same Conference on 3D Advanced Image Processing in Medicine held in Rennes, France. He is also Associate Editor of the IEEE Transactions on Medical Imaging. His research interests include Image analysis and segmentation, Pattern recognition in 2D and 3D images, Image data fusion, 2D and 3D Image reconstruction. He has published more than 30 papers and holds one Patent. Christian Roux is member of SEE, AFCET, SFPT and IEEE. About the Author-vALI?RIE BURDIN received her Ph.D. degree in Signal Processing from the University of Rennes (France) in 1992. In 1988, she joined the Image Processing Group ofthe Ecole Nationale Sup&ieure des T&communications de Bretagne, Brest (France) where she works on the 3D Medical Imaging. Her research interests include the image segmentation, the analysis and the reconstruction of 3D objects. About the Author-WOLFGANG SCHCITTEFELSCHEreceived his Engineering degree from the Technical University of Munich, Germany, in 199 I. In the same year, he joined the Image Processing Group of the Ecole Nationale Sup&ieure des Tel&communications de Bretagne, Brest (France), where he worked on the 3D Medical Imaging. About the Author-CHRISTIAN LEF~VRE, M.D., works in the Department of Anatomy and Department of Orthopoedic and Traumatology Surgery of the University of Brest, France. His research interest is the locking nailing of the forearm.