Image Anal Stereol 2007;16:179-185 Original Research Paper LOCAL GEOMETRICAL ANALYSIS OF 3D POROUS NETWORKS BASED ON THE MEDIAL AXIS: APPLICATION TO BONE MICROARCHITECTURE MICROTOMOGRAPHY IMAGES Françoise Peyrin 1>2>3>4, Zsolt Peter 1>2>3>4, Aymeric Larrue1'2,3,4, Alexandra Bonnassie ' ' and Dominique ATTALI CNRS UMR5220, CREATIS, F-69621 Villeurbanne, France; 2Inserm U 630, CREATIS, F-69621 Villeurbanne, France; 3INSA-Lyon, CREATIS, F-69621 Villeurbanne, France; 4ESRF (European Synchrotron Radiation Facility), BP 220, 38043 Grenoble, France; 5CNRS UMR 5083, GPIPSA, 38402 St Martin d’He`res Cedex, Grenoble, France e-mail: peyrin@esrf.fr (Accepted October 17, 2007) ABSTRACT This paper proposes a new method for local evaluation the geometry of complex 3D porous networks such as bone micro-architecture. The method described here allows local quantification of the geometry of the entire volume. It is well adapted to high resolution images where the voxel size is smaller than the structure thickness. The classification of each point on the medial axis is based on a local topological analysis of a neighbourhood, the size of which is adapted locally to the studied point. This property makes this classification more robust to the small irregularities that may appear in a real biological data set. After classification, new quantitative parameters of bone micro-architecture characterizing plate and rod structures are independently defined. The method is quantitatively evaluated on numerical phantoms presenting shape irregularities. Application to experimental 3D bone microtomographic images shows that the method provides visually correct results, even on complex natural trabecular structures. The additional geometrical parameters provided by this method allow discrimination of volumes whose porosity and mean thickness are similar. Keywords: image analysis, tomography, topology. INTRODUCTION Three-dimensional analysis of complex porous networks of various materials has been generating increasing interest. One of the reasons for this is the growing accessibility of imaging systems providing 3D images of such materials. In particular microtomography systems are now commercially available and allow the examination of the bulk of materials. The availability of 3D images instead of 2D slices opens new possibilities for the quantification of material properties. Powerful stereological methods can reliably relate 2D observations to 3D statistical information such as the volume, area and length per unit volume of selected features. However, some characteristics cannot be inferred from observations of sections. Unbiased estimation of connectivity, and more generally, realistic estimation of topological properties, require three-dimensional measurements. We shall particularly focus on the quantification of complex porous networks such as bone microarchitecture. Bone micro-architecture is generally described as an interconnected network of trabeculae elements, which can be rod-like or plate-like (Parfitt, 1992). However, this plate-rod organization is subject to large variations according to the anatomical site and may evolve with bone loss. Thus, the actual structure is generally intermediate between these two ideal types. In conventional 2D histomorphometry, the calculation of morphometric parameters such as trabecular thickness typically assumes a parallel plate model (Parfitt et al., 1987). This may be a source of error in the interpretation of the results if the model is not well adapted. The parallel plate assumption can be avoided when 3D images are available by measuring a “direct” trabecular thickness (Hildebrand and Ru¨egsegger, 1997a). The structure model index (SMI) was also introduced in order to estimate if the bone structure is rod-like or plate-like (Hildebrand and Ru¨egsegger, 1997b). However, this method, only gives an estimation of the global shape of the structure. In this paper, we propose a new method of local evaluation of the geometry of complex structures. The method is based on medial axis representation followed by topological analysis. Different topology-based analysis methods have already been proposed (Saha et al., 1996; 2000; Pothuaud et al., 2000). The use of a line-skeleton estimation to estimate topological features does not allow preserving the information about plate structures (Pothuaud et al., 179 Peyrin F et al: Local geometrical analysis 2000). A complete topological analysis of bone networks allowing characterization of different types of points from their skeleton has also been proposed (Saha et al, 1996; 2000). The method consists in analyzing the topology of the 3x3x3 neighborhood of each skeleton voxel. However, these different methods were applied to MRI images with a spatial resolution of the same order as the trabecular thickness. When dealing with images at higher spatial resolution, the noise sensitivity of skeletons is a concern. Conversely to both approaches where only skeleton points were characterized, the proposed method will allow analysis of all object points. BACKGROUND Topological analysis Let X be a digital object in Z3. The Euler Number #(X) and Betti numbers ßi(X) for i=0,1,2 are conventionally used as topological characteristics of X. These quantities are related by the following relationship (Odgaard, 1997): X(X) = ß0{X) -ß1(X) + ß2{X) . (1) ß0{X) and ß2(X) respectively correspond to the number of connected components and of isolated cavities in X. ß1(X) may be interpreted as the connectivity of the structure, corresponding to the maximum number of connections that one can remove without separating a part of the structure into two connected sub-structures. The practical calculation of these topological characteristics on a discrete object requires choosing complementary adjacency systems for the object and background (typically 26-6). The Euler number #(X) can be practically computed from the discrete objectX as: X{X) =n0{X)-n1(X) + n2(X) -n3(X) , (2) where n0(X),n1 (X), n2(X) and n3(X) are respectively the number of voxels, faces, edges, and corners. Instead of performing a global topological analysis on X, we shall carry out a local analysis on a selected neighbourhood around each point on the medial axis. Medial axis The medial axis (MA) or skeletons provide simplified representations of objects. In 3D space, MA are not restricted to curves, but may be composed of curves and surfaces. They are mathematically defined as the set of the centres of maximal balls, i.e., balls that are not completely covered by any other single ball in the object. Different classes of methods allow computation of skeletons: “prairie fire” models (Blum, 1967), homotopic thinning (Lam et al, 1992) or distance transform (Di Baja and Thiel, 1994). In practice, skeletons should satisfy a number of desirable properties, such as reversibility, thinness (i.e., be one-voxel thick), and homotopy (i.e., preserve the topology of the original object). However, none of the different methods allows satisfying these properties simultaneously, so one has to relax some property when choosing a method. In our application, we privileged the reversibility property and chose a distance map approach. The MA is then defined as the ridges of the distance transform (DT) of X. The major interest of this approach in our context is that the object can be exactly recovered from its skeleton. In addition, the method is computationally efficient. Instead of the usual discrete distances d6, d18 or d26, we used chamfer distances or “weighted” discrete distances since they were shown to provide better approximations to the Euclidean distance (Borgefors, 1984). We used the chamfer distance defined by the mask (3,4,5) representing the integer weights assigned to neighbours of the current voxel. The computation of the distance map was implemented by using a sequential method in which the image is scanned in forward and backward directions with a half 3x3x3 mask. The calculation of the MA from chamfer distance has been previously described (Remy and Thiel, 2002). It consists in finding the local maxima of the distance transform after lowering some labels. Note that thanks to the reversibility property, the object may be exactly recovered from its medial axis by taking the union of the spheres centred on its points with a radius given by the distance map, which may be formalized as: X = U {B(x,DT(x))/xeMA(X)} , (3) where B(x,r) denotes the ball centred at x with radius r. PROPOSED METHOD GENERAL FRAMEWORK Some of the authors of the present paper previously described the basic principle of our method for locally analyzing the geometry in a 3D image (Bonnassie et al., 2003). The method is based on 3D topological analysis of a local region of interest around each 180 Image Anal Stereol 2007;16:179-185 point on the medial axis. The analysis process led us to introduce 4 classes of points: Boundary (0ê, Label 1), Branching (^V, Label 2), Regular (,0s, Label 3, corresponding to plates), and Arc (M, Label 4, corresponding to rods). The proposed method involves several steps which are described in the next subsections: 1. Calculation of transforms: distance map, thickness map and medial axis 2. Classification of MA points by local topological analysis 3. Classification of all object points and elimination of boundary points 4. Quantification: extraction of new quantitative geometrical parameters CALCULATION OF TRANSFORMS The distance transform and MA were computed as described above. We introduce a local thickness map TM defined on each point x of X as: TM{x) = 2 max {rjr € Biy,r),y €X} . (4) This definition is the transposition in discrete space ofthat proposed by Hildebrand and Ru¨egsegger (1997b) in continuous space. It associates to each point of the object its local thickness, defined as the size of the maximal ball including the given point. Note that DT and TM coincide on the medial axis points. According to Eq. 3, the thickness map may be built by overlapping the balls centred on the medial axis sorted by increasing radius, which may be written as: TM(x) =2 max{DTiy)/xeBiy,DTiy)),yeMAiX)}. (5) CLASSIFICATION OF MA Classification is based on the topological analysis of a region of interest (ROI) centred at each voxel of MA (X). The intersection between the surface of the object and a maximal ball enables distinguishing different types of points (regular, arc, multiple points). Because of noise sensitivity of skeletons, we do not analyze the maximal ball around the point, but a ball with an augmented diameter, constructed as follows. Let B(x) be the maximal ball centred at x. According to the definition of TM, the diameter of B(x) is T Mix). We introduce Be(x), the ball centred at x with diameter TMix) + 2e, where e is a constant which will govern the tolerance of the method. We then define a ROIVe(x) as the complement of the connected component Ce(x) containing x inside Be(x): Ve(x)=Ce(x) with Ce(x)cBe(x). (6) The ROI Ve(x) is then used for topological analysis. The Betti numbers bo(Ve(x)), b2(Ve(x)) and the Euler number c iVe (x) ) are calculated, as referred to in the first section, using a 26 connectivity. We assume that the object contains no void spaces, i.e., that bo (Ve (x) ) is always zero, then b\ iVe (x) ) is derived by using relation (1). The label L(x) is deduced from the two values bo(Ve(x)), and b2(Ve(x)) following the rules given in (Bonnassie et al, 2003). CLASSIFICATION OFX Next, the classification is propagated to the whole volume by using a process similar to the reconstruction of the object from its medial axis. The label L(x) at any point x in X is obtained from the labeling of the skeleton points by: Lix) = max {Liy)/x € B(y,DT(y)),y eMA(X)} . (7) Using this method, boundary points are not limited to the surface of the object and are thus overestimated. They do not provide useful information in the scope of the present analysis, and thus they are subsequently eliminated by reassigning their label to one of the 3 other classes. For this purpose, we use an iterative process where each boundary label is replaced by the most frequent label in its 26-neighbourhood. QUANTIFICATION: NEW QUANTITATIVE GEOMETRICAL PARAMETERS The entire classification procedure provides a partition of the 3D image into 3 classes ,jV, S? and M, being respectively the set of branching, plate and rod points such that: X=^\J,^\JM. (8) After classification, each shape may be quantified independently from each other. We compute the number of voxels in each class, ,jV, S? and M, respectively denoted as NV, PV and RV, and use their ratio to the bone volume BV or total volume TV. 181 Peyrin F et al: Local geometrical analysis The trabecular thickness, conventionally denoted Tb.Th? (in mm), may be defined as the mean value of TM calculated on the entire object X as: Tb.Th? xeX TM(x)dx x(èX dx (9) Using similar nomenclature, we introduce the mean thickness of each class, N , P and R, respectively denoted as N.Th?, P.Th?, and R.Th?. These are obtained by limiting the calculation of the mean to N , P and R respectively, so they do not require additional computation of any thickness maps. RESULTS SIMULATED IMAGES The method was evaluated on different simulated images containing plate and rod structures with various known thicknesses. As an illustration, we show results obtained on a mixed (128)3 phantom composed of the following structures: three plates (thicknesses=6, 8, 10, height and width=120), one cylinder (left-side, diameter=10, height=120), and one elliptical cylinder (right-side, half axes=10 and 14, height=120) (sizes are given in voxels). The proposed new geometrical analysis method was applied (e was set to 4). The different steps of the classification method are displayed in Fig. 1 using a color code. From a qualitative point of view, it may be seen that the classification is in good agreement with the definition of the phantom. An important point is that the method is robust with respect to the elliptical cylinder, which is seen as a rod, as expected in practice. Quantitative evaluation was also performed. Table 1 reports the relative percentages of each class in the original and 3-class classified image, as well as the thickness of each structure (expressed in voxels). For simulated images made only of cylinders or planes, the percentages were exact. When structures are interconnected, there are slight errors, which mainly come from the fact that branching points are by construction underestimated, since the algorithm privileges plates and rods. The respective thicknesses of the plates and rods are close to the theoretical values, with an error of less than one pixel. The difference in thickness between rods and plates is well estimated. (c) Fig. 1. Analysis of a mixed phantom: a): classification of MA, b) 4-class classification of the entire volume, c) 3-class classification of the entire volume. The color code is the following (yellow: boundary, red: branch, green: plate and magenta: rod). _ 182 Image Anal Stereol 2007;16:179-185 Table 1. Quantitative parameters extracted from the classified mixed phantom: 1st line: theoretical values, 2nd line: estimated values from the analysis method, 3rd line: deviation between theoretical and estimated parameters. The ratios NV/BV, PV/BV and RV/BV are expressed in percentages and the thicknesses Tb.Th?, P.Th?, and R.Th? in voxels. Phantom NVBV PVBV RVBV Tb.Th* P.Th* R.Th* Theor. 5.65% 87.67% 6.68% 9.2 9.0 12.0 Analysis 1.69% 91.87% 6.44% 8.9 8.7 11.4 Error -3.97% 4.20% -0.24% -0.3 -0.3 -0.6 BONE MICRO-CT IMAGES The proposed method was also applied to images of human bone samples acquired using synchrotron radiation microtomography (micro-CT) at the ESRF (European Synchrotron Radiation Facility) (Salome et al., 1999; Nuzzo et al., 2002). For each sample, 900 2D radiographs under different angles of view were acquired on a (1024)2 CCD-based 2D detector. The optic was set to get a pixel size on the detector of 10.13 mm. X-ray beam energy was set to 20 KeV. The 3D images were reconstructed using a 3D version of the standard filtered backprojection algorithm. A central volume of interest of (600)3 voxels was selected in each image (voxel size: 10.13 mm). The bone phase was easily selected by a simple thresholding method due to the high signal-to-noise ratio and high contrast of the 3D images. The conventional bone micro-architecture parameters (Parfitt et al., 1987; Hildebrand and Ru¨egsegger, 1997a) were evaluated: Bone Volume/Total Volume (BV/TV), Trabecular thickness Tb.Th?, Bone Surface/Bone Volume (BS/BV), and Trabecular number (Tb.N). The proposed local geometrical analysis was then performed. Fig. 2 illustrates 3D displays of bone surfaces for three images corresponding to samples with similar mean thickness (Tb.Th? 134 mm). Fig. 3a) shows a porous sample with a low BV/TV (5.6%), while the two others images correspond to samples with similar BV/TV (12.6% and 12.9%). The last two images demonstrate that although BV/TV and Tb.Th? are different characteristics calculated independently from each other, they are not sufficient to discriminate the samples. The result of the analysis of our new geometrical method is displayed in Fig. 3 using the same color code as before (red: branch, green and magenta: rod). To facilitate the visualization, only a sub-volume in each analyzed volume is displayed. Observation of Fig. 3 reveals that structures were labeled as expected. The ratio RV/PV of the Rod Volume to Plate Volume obtained for the three volumes were respectively: 1.82, 1.14 and 0.56. These values are in agreement with the appearance of the volume, which is respectively seen as rod-like, intermediate and platelike. Note that this additional geometrical parameter allows discrimination of the last two volumes. CONCLUSION We have described a new method for local characterization of the geometry in 3D images of porous network. The initial method providing a 4-class classification of the MA was improved to get a complete classification of the entire volume and eliminate border points which were overestimated. The application to simulated and experimental images shows that the method provides visually correct results, even on complex natural trabecular structures. In addition, the accuracy of the classification was tested on geometrical phantoms, establishing reduced errors. The main advantage of our method is that it is robust to surface irregularities: a structure will be recognized as a rod even if it is not purely cylindrical, conversely to other approaches. Simulations performed on elliptical cylinders confirm that the proposed method provides a large tolerance to this problem. This property is important in coping with the biological variability of trabeculae. The results effectively obtained on experimental micro-CT bone images demonstrate that rod trabeculae are properly recognized, although not ideally shaped. The proposed geometrical analysis allows characterization of the percentages of rods and plates in each sample, as well as separate quantification of the thicknesses of rods and plates. This new set of parameters is assumed to bring additional characteristics to the geometry of trabecular structures. 183 Peyrin F et al: Local geometrical analysis (b) BV/TV=12,6%, Tb.Th* =133 |im (a) RV/PV =1.82 %k. ¦%SPt-Ä. ¦^B pH !- V