- Email: [email protected]
A robust method for measuring trabecular bone orientation anisotropy at in vivo resolution using tensor scale
Pattern Recognition 37 (2004) 1935 – 1944 www.elsevier.com/locate/patcog
A robust method for measuring trabecular bone orientation anisotropy at in vivo resolution using tensor scale Punam Kumar Saha∗ , Felix W. Wehrli Department of Radiology, University of Pennsylvania, 4th Floor, Blockley Hall, 423 Guardian Drive, Philadelphia, PA 19104, USA Received 17 November 2003; accepted 19 December 2003
Abstract Trabecular bone (TB) is a network of interconnected struts and plates that constantly remodels to adapt dynamically to the stresses to which it is subjected in such a manner that the trabeculae are oriented along the major stress lines (Wol3’s Law). Next to bone density, TB anisotropy has been found to be one of the most signi6cant determinants of the bone’s biomechanical behavior. Typically, orientational anisotropy of TB networks is expressed in terms of the fabric tensor, obtained by measuring the mean intercept length between structure elements along test lines. This method, however, can provide only a global statistical average of TB orientation anisotropy and, in general, requires a large sampling volume. Here, we present a new method, based on the recently conceived notion of “tensor scale”, which provides regional information of TB orientation anisotropy. Regional structure is represented by local best 6t ellipsoid (ellipse in 2D) and the structural orientation is determined from the eigenvectors along the semi-axes. The method is found to be remarkably robust over a wide range of resolution regimes and image rotation as shown with micro-CT images from specimens of the human distal radius, the latter showing the characteristic di3erences in structural anisotropy for transverse and longitudinal sections. Finally, the method’s applicability to in vivo MR imaging is demonstrated with data from the distal tibia. ? 2004 Pattern Recognition Society. Published by Elsevier Ltd. All rights reserved. Keywords: Tensor scale; Fuzzy subset; Trabecular bone; Orientation; Anisotropy; -MRI-CT;
1. Introduction Analysis of regional architecture is of paramount interest in many medical imaging applications. This paper illustrates an application of the recently developed “tensor scale” concept [1] for quantitatively measuring tissue architectural parameters. Although tensor scale here is applied to the analysis of local structure morphology, the notion of scale has a long history in applications to computer vision and image processing. “Scale” may be thought of as the spatial resolution, or, more generally, a range of resolutions needed to ensure a suDcient yet compact object representation. Scale plays an important role in determining ∗ Corresponding author. Tel.: +1-215-662-6780; fax: +1-215898-9145. E-mail address: [email protected] (P.K. Saha).
the optimum trade-o3 between noise smoothing and perception of structures. Also, scale is helpful in breaking a computer vision and image-processing task into a hierarchy of tasks where those at higher levels deal with larger structures. Marr [2] pointed out many bene6ts of multi-scale image analysis and computer vision. Although the concept of scale originally evolved from scale-space theory [2–5], recently, the concept of “local scale” [6–8] has attracted signi6cant interest as it allows space-variant process control 6tting with local structures. The notion of local morphometric scale was presented in [7,8] using the spherical model. A major limitation of this model is its ignoring orientation and anisotropy, which makes it unsuited for many applications, especially in biomedical imaging where structures are inherently anisotropic. The local scale using a tensor model [1] (an ellipse in 2D and an ellipsoid in 3D) provides a uni6ed and parametric representation of size, orientation,
0031-3203/$30.00 ? 2004 Pattern Recognition Society. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.patcog.2003.12.019
1936
P.K. Saha, F.W. Wehrli / Pattern Recognition 37 (2004) 1935 – 1944
and anisotropy of local structures. Although tensor scale is a concept with wide-ranging applications in computer vision and image processing, here we focus on its application in quantitative analysis of regional morphology only; a few other applications of tensor scale have been reported in [9,10]. The focus of this paper is to explore the potential of tensor scale to analyze global and regional trabecular bone (TB) morphology, although the method may also be applicable to other tissues, e.g., the characterization of vascular micro architecture in during angiogenesis [11,12]. TB consists of a meshwork of interconnected struts and plates providing the bone high mechanical strength at minimum weight. Quantitative analysis [13] of trabecular bone structure by in vivo MR imaging [14] will be useful for fracture risk evaluation, early assessment of disease status, therapy response monitoring and screening of high risk subjects. The most common bone disorder is osteoporosis, which leads to enhanced bone fragility and increase the risk of fracture caused by both bone loss and micro architectural deterioration of bone tissue [15]. It is generally agreed [16] that 50–60% of the bone’s mechanical competence can be explained by variations in the apparent density (bone mass/tissue volume). However, clinical and experimental evidence for an independent contribution from trabecular architecture is compelling [17,18]. A few studies [19–22] involving histomorphometry in patients with and without vertebral fractures, matched for gender and bone mineral density (BMD), found the two groups of subjects to di3er in histomorphometric indices. Trabecular bone remodels, i.e., it adapts to the stresses to which it is subjected (Wol3’s Law [23]) in such a manner that the trabeculae follow the major stress lines. Thus, TB at most skeletal locations is highly anisotropic. Structural anisotropy (fabric) can be expressed in terms of a tensor and, although intuitively obvious, it has recently been shown that the stress tensor is aligned with the fabric tensor [24]. Next to material density (i.e., bone volume fraction) structural anisotropy most critically determines TB’s mechanical behavior [25]. Classically, the fabric tensor can be obtained from images by determining the mean intercept length (MIL) [26] of test lines along all spatial directions. This method, however, requires a suDciently large number of intercepts, a requirement that is only met for suDciently large sampling volume. Also, the MIL method provides only a statistical average of orientation anisotropy. This paper introduces a new method for computation of global and regional TB orientation anisotropy using the recently developed tensor scale theory and presents results of preliminary evaluative experiments. Speci6cally, the robustness of the method under varying resolutions and image rotations is studied and the characteristic differences between TB anisotropy computed from transverse and longitudinal sections are examined. Finally, potential applications of the method to in vivo TB MR images are demonstrated.
Fig. 1. Illustration of tensor scale computation at an image point (the red dot). Finitely many radially outward directed straight-line segments, called sample lines (blue), distributed at approximately uniform angular spacing, are traced and the optimum edge location (orange and pink dots) is identi6ed on each sample line. Exploiting the axial symmetry of the ellipse, the edge points on every pair of radially opposite sample lines are repositioned (orange and yellow dots) by selecting the one closer to the candidate point and then reMecting it onto the opposite sample line. Finally, tensor scale is computed by 6tting an ellipse to the repositioned edge points.
2. Tensor scale: theory and algorithms Tensor scale at an image point p is the parametric representation of the largest ellipse (ellipsoid in 3D) centered at p that is contained in the same object region. The process of tensor scale computation (in 2D) at any pixel p is schematically illustrated in Fig. 1. The method is based on tracing pairs of radially opposite sample lines emanating from p and locating the closest edge on every sample line. After locating edges on all sample lines, they are repositioned by making use of the ellipse’s axial symmetry. Finally, the tensor scale at p is determined by computing the best 6t ellipse to these repositioned edge points. The method involves the following four sequential steps—(1) intensity computation along each sample line, (2) location of the optimum edge on each sample line, (3) repositioning of the edge locations, and (4) computation of the best 6t ellipse to the repositioned edge locations. In the following, we describe the method in more detail. Although the method immediately extends to any 6nite dimension, we con6ne our description to 2D, the dimensionality of the present implementation. A digital object O=(Z 2 ; O ) is a fuzzy subset of Z 2 , where Z is the set of all integers and Z 2 represents a 2D digital
P.K. Saha, F.W. Wehrli / Pattern Recognition 37 (2004) 1935 – 1944
grid; and O : Z 2 → [0; 1] speci6es the membership value at each grid point, often referred to as a pixel. As pointed out earlier, the current application aims to compute orientation anisotropy of trabecular bone from -CT or in vivo -MR images. Following our target application, a digital object O will represent a trabecular bone network and at any pixel p, the membership value O will represent its bone volume fraction (BVF). An imaging system [27] acquires images containing information of a target object, often along with other co-objects. The process of extracting an object from an image, commonly referred to as image segmentation [28–31], has extensively been studied for decades. Following the main focus of the theory and methods addressed in this paper, we start with the de6nition of digital objects instead of digital images. However, in Methods, we will describe how BVF maps are generated for di3erent experiments from original image data. For now, let us use the above de6nition of digital objects. The support (O) of a digital object O is the set of all pixels with non-zero membership. Tensor scale is computed for all pixels p ∈ (O). (p) (p) (p) Let (1(p) ; 1 ); (2(p) ; 2 ); · · · ; (m(p) ; m ), where (p) (p) i ; i : [0; L] → R2 ; ∀i, denote m pairs of mutually opposite sample lines each having length L and emanating from the candidate pixel p ∈ Z 2 . The sample lines are intended to be uniformly distributed over the angular space around p to ensure that the edge points computed on these sample lines, eventually utilized to compute the tensor scale at p, are not skewed in any direction. A set of sample lines skewed in one direction may lead to a bias in the computed tensor scale in terms of its orientation. For example, if the sample lines were more densely distributed around the horizontal direction then the computed tensor scale would be less sensitive to vertical components of the local structure; instead, it could be biased to horizontal components of the structure. Each sample line is coupled with a radially opposite sample line; we thus refer to each such pair as a conjugate pair and to each sample line of the pair as a complementary sample line of the other. The membership pro6le on a sample line is computed by interpolating object membership values at points selected at regular intervals. During the second step, the optimum edge location is identi6ed on each sample line. For a raw (unsegmented) image, the edge point is located at the transition between two objects, i.e. at the zero crossing of the second-order derivative of intensity. Typically, voxel intensity is directly proportional to local object material density. However, at limited resolution, most voxels are only partially occupied by the object material, hence yielding only a fuzzy representation of the object. In these situations it is appropriate to presegment the images as a fuzzy object. For such fuzzily segmented objects, the edge point is not detected at the closest intersection between the boundary of the support of the object and the sample line. Rather, the fuzzy length of the locally connected object along a sample line is computed. Speci6cally, the object edge point i(p) on a sample line i(p) , emanating from p, is computed whose distance from p,
1937
denoted by i(p) −p, is equal to the integration of the membership values along i(p) until a non-object point (membership value = 0) is found. In other words, T (p;(p) ) (p) i di (t) (p) (p) i − p = O (i (t)) (1) dt; dt 0 where, T (p; i(p) ) is the parametric location of the 6rst (from p) zero membership point on i(p) . The edge location on the sample line i(p) is denoted by i(p) . The edge points, located as above, are intended to roughly describe the boundary of the ellipse representing the tensor scale at p. From the axial symmetry of an ellipse, two radially opposite points on its boundary are always equidistant from its center. However, the points along the boundary of a local structure do not always satisfy this condition of symmetry with respect to a candidate pixel. For example, in Fig. 1, the edge locations on the west- and north-bound sample lines are in general farther from the candidate pixel than the edge points on respective complementary sample lines. Therefore, to ensure axial symmetry of an ellipse, the edge points need to be repositioned. This step is accomplished by analyzing the edge points on every conjugate pair of sample lines. Speci6cally, the one closer to p between i(p) and
i(p) is selected and reMected on the complementary sample line. In Fig. 1, the edge points shown by pink dots are repositioned yielding the yellow dots; orange dots represent the edge points which remain unmoved after the repositioning process. The repositioned edge points will brieMy be referred to as r-edge points. The 6nal step in tensor scale computation is to 6t an ellipse to the r-edge points. This step is accomplished in two sub-steps—(i) determination of ellipse orientation, (ii) computation of the lengths of two semi-axes. The 6rst sub-step is completed using principal component analysis [32] of the r-edge points and assuming p as the center. Speci6cally, the orientation of the ellipse is computed as the angle of the eigenvectors associated to the larger eigenvalue of the covariance matrix de6ned as follows 1 1; 1 = (xi − xp )2 ; 2m i=1;2;:::;2m 2; 2 =
1 2m
(yi − yp )2 ;
(2)
i=1;2;:::;2m
and, 1; 2 = 2; 1 =
1 2m
(xi − xp )(yi − yp );
i=1;2;:::;2m
where (xi ; yi ) | i = 1; 2; : : : ; 2m are the coordinates of the r-edge points and (xp ; yp ) is the coordinate of p. Finally, the actual lengths a; b of the two semi axes of tensor scale ellipse are computed by minimizing the following error function u2 v2 ferror = 1 − i2 − i2 ; (3) a b i=1;2;:::;2m
1938
P.K. Saha, F.W. Wehrli / Pattern Recognition 37 (2004) 1935 – 1944
where, (ui ; vi ) | i = 1; 2; : : : ; 2m are the r-edge points after (i) applying a translation by −p (i.e., p is translated to origin), and then (ii) applying a rotation by angle − (i.e., the major semi axis of the ellipse is aligned to the horizontal coordinate axis). By setting the partial derivatives of ferror with respect to a and b to zero, it can be shown that the a and b which minimize ferror are given by the following expressions
was selected in this experiment to examine whether similar anisotropies are obtained in two di3erent sets of longitudinal slices within the same bone volume (which assumes that the trabecular structure is homogeneous within a given sample). In Experiment 4, designed to evaluate whether the technique was able to delineate local structural orientation in -MR images transverse MR images acquired in the distal
a=
xi2 yi2 i=1; 2; :::; 2m xi2 yi2 − i=1; 2; :::; 2m xi4 i=1; 2; :::; 2m yi4 ; 2 2 2 2 4 i=1; 2; :::; 2m xi yi i=1; 2; :::; 2m yi − i=1; 2; :::; 2m xi i=1; 2; :::; 2m yi i=1; 2; :::; 2m
b=
xi2 yi2 i=1; 2; :::; 2m xi2 yi2 − i=1; 2; :::; 2m xi4 i=1; 2; :::; 2m yi4 : 2 2 2 4 2 i=1; 2; :::; 2m xi i=1; 2; :::; 2m xi yi − i=1; 2; :::; 2m xi i=1; 2; :::; 2m yi
i=1; 2; :::; 2m
Finally, we shall use (p) to denote the tensor scale ellipse at any pixel p and Y1 (E) ¿ Y2 (E) to denote the two semi-axes of an ellipse E.
3. Methods Four di3erent experiments were conducted for preliminary evaluation of the proposed orientation anisotropy computation method and to demonstrate its potential in MR imaging applications of in vivo human trabecular bone: Experiments 1 and 2 were designed to evaluate the sensitivity of the computed orientation anisotropy to resolution and object rotation. The purpose of Experiment 3 was to demonstrate the characteristic di3erences in orientation anisotropy between transverse and longitudinal image slices. Finally, Experiment 4 focused on regional orientation anisotropy computed from 2D MR images of in vivo human trabecular bone. For all experiments, the following sequence of steps were applied—(1) computation of a BVF map from an acquired raw image, (2) tensor scale computation at every image location, and (3) computation of global or regional orientation anisotropy. 3.1. Material description and BVF map computation For Experiments 1–3, a -CT image was selected from a 3D stack of images that had been acquired previously (SCANCO Medical -CT 20, at 22 m isotropic voxel size) as part of a di3erent project [33] from samples cored from cadaveric bone of the distal radial metaphysis. The sample consisted of a cylindrical core of 9 mm height and 9 mm diameter with the cylinder axis parallel to the direction of the radius. The -CT image data set was processed to generate BVF maps. For Experiments 1 and 2, two central longitudinal slices, separated by 330 m (i.e., 15 slices), were selected. For Experiment 3, two additional longitudinal slices, shifted by 15 slices from the previous set, and two central transverse slices, separated by 15 slices from each other, were selected. A di3erent set of longitudinal slices
(4)
(5)
tibia (137 × 137 × 410 m3 voxel size) in two women, ages 53 and 49 years, using a 3D partial Mip-angle pulse sequence [34] were processed as described above. The objective of these experiments was to evaluate whether the method is able to quantify regional structural orientation from in vivo MR images. Because of high spatial resolution exceeding trabecular thickness by a factor of 5, the -CT images had a bimodal intensity distribution and thus could be segmented by selecting a threshold at the midpoint of the two modes representing the mean bone and mean marrow intensities. Subsequently, the largest 26-connected component [35] within the region thresholded for bone was computed. The intensities at the two modes were selected as the mean intensities for bone and marrow (denoted by B and M , respectively). At each pixel p, the BVF value f (p) was computed using the following equation: 0 if f(p) ¡ M ; f(p) − M if M 6 f(p) ¡ B ; (6) f (p) = B − M 1 otherwise: For Experiment 1, the BVF images of the selected slices were resampled at pixel sizes corresponding to integral multiples of the parent resolution, i.e., 44 × 44, 66 × 66, and 88 × 88 m2 . This strategy was chosen for generating 2D trabecular bone images at lower resolutions (rather than resampling the raw -CT slice data sets followed by computing the BVF images). The latter approach would cause signi6cant loss in structural detail in the segmented images using the simple threshold-based method used here. Although more sophisticated segmentation could be used, no other segmentation were attempted since this was deemed to be outside the scope of the present paper. For Experiment 2, the rotation of the raw 2D -CT images was performed before computing BVF since in this manner the maximum e3ect of rotation on the entire process can be assessed. Anticlockwise rotations of ∈ {15◦ ; 30◦ ; 45◦ ; 60◦ } were applied to each of the selected raw 2D -CT images around the center of each slice and the rotated images were
P.K. Saha, F.W. Wehrli / Pattern Recognition 37 (2004) 1935 – 1944
1939
redigitized into a square grid using linear interpolation. For both Experiments 2 and 3, BVF images were computed using Eq. (6). For the in vivo MR experiment, an actual region of interest (ROI), for which the tensor scale was computed, was generated by manually tracing boundaries on each of the two selected 2D images using an IDL (Interactive Data Language, Research Systems Inc.)-based graphical interface running on a Pentium III PC under MS Windows OS. The data encompassed by the ROIs were preprocessed by deshading and noise reduction using a histogram deconvolution method [36] to produce BVF images. 3.2. Tensor scale computation Tensor scale was computed at every image pixel as described in Section 2. Two parameters are required to compute tensor scale—(1) length L of each sample line, and (2) number of sample lines 2m. L should be suDciently large compared to trabecular thickness to capture local structural anisotropy. Here, we used L = 704 m or 32 pixels which is approximately 4–7 times larger than typical human trabecular thickness. The number of sample lines determines the precision of the computed tensor scales. Following the recommendation in [1], 30 sample lines were used for the experiments described. 3.3. Regional orientation and anisotropy Here, we describe the method to compute global as well as regional orientation of a trabecular bone network from its tensor scale map. Let p be any pixel inside the support (O) of the object O. The principal orientation of the trabecula at p may be de6ned by the unit vector Y1 (p)=|Y1 (p)|. Let H (i) | i ∈ {1; 2; : : : ; 180} represent the orientation distribution over the entire object where i is the value of an angle (anti-clockwise) in degrees. H (i) is computed as follows: p (i) H (i) = O (p) ; (7) |Y1 (p)| p∈(O)
where p (i) is the length of the “central chord” (a chord passing through the center) of the tensor scale ellipse (p) at p along the direction making i degrees anti-clockwise angle with a horizontal line. The length of p (i) can be expressed as follows: |Y1 (p)|2 × |Y2 (p)|2 ; (8) p (i) = |Y2 (p)|2 cos i + |Y1 (p)|2 sin i The value of p (i) as obtained by using the above equation is illustrated in Fig. 2. An ellipse, denoted by H , is computed from the rose plot of H (i) using a method similar to the one described in Section 2. Let Y1; H ¿ Y2; H denote the two semi-axes of H . Global orientation anisotropy of the TB network then is computed as 1 − |Y2; H |2 =|Y1; H |2 —the eccentricity of the ellipse H . Regional orientation anisotropy at a pixel p is obtained by calculating the orientation histogram Hp over
Fig. 2. Illustration of the length of the “central chord” (a chord passing through the center) of an ellipse along a given direction as de6ned in Eq. (8).
a neighborhood around p which, in this paper, has been chosen as a disk of radius 25 pixels. 4. Results In this section, we present results of the experiments described in the previous section. The motivation of the 6rst experiment was to assess the performance of tensor scale-based orientation anisotropy computation in various resolution regimes. Recall that four di3erent pixel sizes (22 × 22, 44 × 44, 66 × 66, and 88 × 88 m2 ) were chosen for this experiment. Figs. 3a,b show an original -CT image slice at 22 × 22 m2 voxel size and its BVF map computed as per Eq. (6). The BVF mapped image of the same slice resampled at 88 × 88 m2 resolution is presented in Fig. 3c. Tensor scale images, unlike other parametric images (e.g. BVF), are not scalar maps. Therefore, it is important to adopt a new display scheme that can adequately convey the visualization of the architectural make-up of the tissue. Since an ellipse centered at the origin is uniquely de6ned by three parameters—orientation of the major semi axis Y1 (p), anisotropy (= 1 − |Y2 (p)|2 =|Y1 (p)|2 ), and the length of the short axis (representing structure thickness) (=|Y2 (p)|) —we have adopted hue-saturation-intensity (HSI) color coding scheme to represent a 2D tensor scale image. Specifically, hue (H), saturation (S), and intensity (I) of the color vector at any pixel p encodes orientation, anisotropy, and thickness, respectively, of (p). Fig. 3d illustrates the color coding at full intensity (i.e., thickness is maximum). According to this coding convention, colors near the center of the disk represent locally isotropic structures while those near the periphery of the disk represent highly anisotropic structures. As illustrated in the disc, a color is repeated at 180◦ representing the axial symmetry of an ellipse.
1940
P.K. Saha, F.W. Wehrli / Pattern Recognition 37 (2004) 1935 – 1944
Fig. 3. E3ects of varying resolution on tensor scale-based TB orientation anisotropy computation: (a) an image slice from a 3D -CT dataset of a cadaveric TB specimen from the human distal radius acquired at 22 m isotropic voxel size; (b) a BVF mapped image computed from (a); (c) the same image as (b) but resampled at 88 × 88 m2 pixel size; (d) HSI color coding at full intensity used to represent tensor scale images; (e,f) color coded tensor scale images computed from (b,c), respectively; (g) rose plots and anisotropy ellipses of TB orientation distributions computed from tensor scale images at di3erent resolutions. The maximum variation in computed anisotropy for the di3erent resolutions examined was 0.70%.
Normalized thickness values are computed by dividing |Y2 (p)| with the maximum scale value L. HSI color-coded tensor scale images computed from Fig. 3b,c are illustrated in Fig. 3e,f, respectively. It may be noted that, although the tensor scale image at the lower resolution is somewhat blurred, orientation, anisotropy, and thickness of matching regions are well preserved. Rose plots of tensor scale-based orientation distributions computed from the sets of images at di3erent resolutions (22 × 22 m2 , 44 × 44 m2 , 66 × 66 m2 , and 88 × 88 m2 pixel size) are illustrated in Fig. 3g. The computed anisotropy ellipses, each 6tted to a rose plot, are also shown in the same 6gure. The maximum variation in anisotropy values at the four di3erent resolutions was 0.70%. The behavior of the proposed method under image rotation is illustrated by using the 2D image displayed in Fig. 4a, obtained by applying 30◦ anti-clockwise rotation on the raw the -CT TB slice image of Fig. 3a. The BVF map computed from Fig. 4a is shown in Fig. 4b while the tensor scale image of Fig. 4b is displayed in Fig. 4c using HSI color coding. It may be noted by comparing Fig. 4c with Fig. 3e that the color (orientation) of matching structures are signi6cantly di3erent in these two images. This mismatch is caused as the absolute orientation of a structure in an image is changed when the image itself is rotated and the amount
of rotation of the structure is exactly the same as that originally applied to the image. Therefore, it was expected that if the rotation applied to the original image were subtracted from tensor scale orientation of each pixel, the modi6ed tensor scales of matching structures in the original and in the rotated images would be similar. This is borne out by tensor scale image obtained by applying the orientation adjustment displayed in Fig. 4d. Rose plots of orientation distribution computed from the set of slices without rotation and those at 30◦ anticlockwise rotation are depicted in Fig. 4e; computed anisotropy ellipses are also shown in the same 6gure. The maximum variation in anisotropy values at di3erent rotations is 1.69%. The third experiment shows that the tensor scale-based orientation anisotropy analysis can reveal the known differences in structural anisotropy between longitudinal and transverse sections in the distal radius in images taken from the same 3D data. Trabeculae are preferentially oriented along the principal load direction and, since the cylindrical bone specimens were cored with their axes aligned with the radial shaft, greater orientation anisotropy is expected in longitudinal than in transverse sections [33,37]. Fig. 5a,b displays transverse and longitudinal slices selected from the same -CT data of a cadaveric TB specimen. Fig. 5c,d shows the BVF mapped images computed from Fig. 5a,b,
P.K. Saha, F.W. Wehrli / Pattern Recognition 37 (2004) 1935 – 1944
1941
Fig. 4. E3ects of image rotation on tensor scale-based TB orientation anisotropy computation: (a) the image of Fig. 3a after applying a 30◦ anticlockwise rotation; (b) a BVF mapped image computed from (a); (c) tensor scale image of (a); (d) the image of (c) after applying rotation-adjustment; (e) orientation rose-plots and anisotropy ellipses computed from respective tensor scale images demonstrating the method’s robustness to image rotation.
Fig. 5. Tensor scale-based TB orientation anisotropy shown for transverse and longitudinal image slices: (a,b) transverse and a longitudinal images from the same -CT image dataset; (c,d) BVF mapped images computed from (a,b); (d,e) tensor scale images derived from (c,d); (f) orientation rose-plots and anisotropy ellipses computed from respective tensor scale images.
1942
P.K. Saha, F.W. Wehrli / Pattern Recognition 37 (2004) 1935 – 1944
Fig. 6. Results of regional orientation analysis of in vivo TB -MR images using tensor scale: (a) cross-sectional image of the left ankle (one of 28 contiguous slices of 410 m thickness and 137 × 137 m2 pixel size from a 28-slice 3D dataset) acquired at the level of the malleolus (female, age 53 years). The foot is tilted approximately 30◦ ; (b) BVF image computed from (a) over a manually selected region; (c) color-coded tensor scale image of TB in (b); (d) regional orientation and anisotropy computed from (c); (e–h) same as (a–d) for a di3erent subject (female, age 49 years).
respectively, and Fig. 5e,f illustrate their color-coded tensor scale images. As visually apparent, the tensor scale image in Fig. 5e displays color almost randomly distributed over the entire color band while that in Fig. 5f shows greater frequency of colors between green (60◦ ) and blue (120◦ ) representing mostly vertical trabeculae (90◦ ). This trend is clearly visible in the rose-plots, shown in Fig. 5g, of the orientation distributions computed from the tensor scale images of selected sets of transverse and longitudinal slices. The computed orientation anisotropy values (0.49 and 0.88 for the sets of longitudinal and transverse slices, respectively) agreed well with these observations. Fig. 6 shows in vivo trabecular bone images from two subjects (female, age 53 and 49 years, respectively) taken in the distal tibia at the level of the malleolus (a lateral protrusion above the ankle), along with BVF maps and parametric images (tensor-scale and regional anisotropy maps). The images show distinct structural anisotropy with the trabecular orientation roughly parallel to the cortical rim. This structural arrangement reMects the multidirectional lateral forces acting on the ankle. 5. Conclusions A new tensor scale-based anisotropy computation method applied to images of trabecular bone has been presented. The method allows computation of tensor scale at any image location in a fuzzily segmented object representation and thus is applicable low resolution in vivo MR images. Since structural anisotropy has implications on the bone’s mechanical
properties, the method may have potential for studying regional changes in structure orientation during bone remodeling on the basis of in vivo MR images.
6. Summary Trabecular bone (TB) is a network of interconnected struts and plates that constantly remodels to adapt dynamically to the stresses to which it is subjected in such a manner that the trabeculae are oriented along the major stress lines (Wol3’s Law). Next to bone density, TB anisotropy has been found to be one of the most signi6cant determinants of the bone’s biomechanical competence. Typically, orientation anisotropy of TB networks is expressed in terms of the fabric tensor, obtained by measuring the mean intercept length between structure elements along test lines. This method, however, can provide only a global statistical average of TB orientation anisotropy and, in general, requires a large sampling volume. Here, we present a new method, based on the recently conceived notion of “tensor scale”, which provides regional information of TB orientation anisotropy. The method’s performance has been evaluated in terms of its sensitivity to resolution and image rotation and the characteristic di3erences between TB anisotropy computed from transverse and longitudinal sections have been examined. Finally, potential applications of the method to in vivo TB MR images are demonstrated. Although the method can be extended to any 6nite dimension, it is currently implemented only in two dimensions.
P.K. Saha, F.W. Wehrli / Pattern Recognition 37 (2004) 1935 – 1944
Tensor scale at an image point p is the parametric representation of the largest ellipse (ellipsoid in three dimensions) centered at p that is contained in the same object region. BrieMy, the method is based on tracing pairs of radially opposite sample lines emanating from p and locating the closest edge on every sample line. After locating edges on all sample lines, they are repositioned by making use of the ellipse’s axial symmetry. Finally, the tensor scale at p is determined by computing the best 6t ellipse to these repositioned edge points. Orientation anisotropy is computed by analyzing orientation of the major semi axis and the eccentricity of tensor scale ellipse at every image pixel location with non-zero fractional occupancy of TB. Experiments were conducted on the basis of microcomputed tomography images acquired at 22×22×22 m3 voxel size of a cadaveric TB specimen cored from the human distal radius with the axis of the cylinder aligned with that of the radius. The experimental results demonstrate the method to be remarkably robust over a wide range of resolution regimes (22 × 22, 44 × 44, 66 × 66, and 88 × 88 m2 ) and image rotations (anticlockwise rotations of 15◦ , 30◦ , 45◦ , 60◦ ). The maximum variation in computed tensor scale-based orientation anisotropy at di3erent resolutions examined was 0.70% while that for di3erent rotations was 1.69%. The method has been found to be successful to quantitatively characterize the di3erences in structural anisotropy for transverse and longitudinal sections of the same micro-CT data. Finally, the method’s applicability to in vivo MR imaging is demonstrated with data from the distal tibia. References [1] P.K. Saha, Novel theory and methods for tensor scale: a local morphometric parameter, in: Proceedings of SPIE: Medical Imaging, San Diego, CA, Vol. 5032, 2003, pp. 314–324. [2] D. Marr, Vision, W.H. Freeman and Company, San Francisco, CA, 1982. [3] A.P. Witkin, Scale-space 6ltering, in: Proceedings of the Eighth International Joint Conference Arti6cial Intelligence, Karlsruhe, West Germany, 1983, pp. 1019–1022. [4] J.J. Koenderink, The structure of images, Biol. Cybern. 50 (1984) 363–370. [5] T. Lindeberg, Scale-space for discrete signals, IEEE Trans. Pattern Recognition Machine Intell. 12 (1990) 234–254. [6] S.M. Pizer, D. Eberly, D.S. Fritsch, Zoom-invariant vision of 6gural shape: the mathematics of core, Comput. Vision Image Understand. 69 (1998) 55–71. [7] P.K. Saha, J.K. Udupa, D. Odhner, Scale-based fuzzy connected image segmentation: theory, algorithms, and validation, Computer Vision Image Understand. 77 (2000) 145–174. [8] P.K. Saha, J.K. Udupa, Scale based image 6ltering preserving boundary sharpness and 6ne structures, IEEE Trans. Med. Imaging 20 (2001) 1140–1155. [9] P.K. Saha, J.C. Gee, Z. Xie, J.K. Udupa, Tensor scale-based image registration, in: Proceedings of SPIE: Medical Imaging, San Diego, CA, Vol. 5032, 2003, pp. 743–753.
1943
[10] P.K. Saha, J.K. Udupa, Tensor scale-based fuzzy connectedness image segmentation, in: Proceedings of SPIE: Medical Imaging, San Diego, CA, Vol. 5032, 2003, pp. 1580–1590. [11] P. Carmeliet, R.K. Jain, Angiogenesis in cancer and other diseases, Nature 407 (2000) 249–257. [12] E.M. Brey, T.W. King, C. Johnston, L.V. McIntire, G.P. Reece, C.W. Patrick, A technique for quantitative three-dimensional analysis of microvascular structure, Microvasc. Res. 63 (2002) 279–294. [13] P.K. Saha, B.R. Gomberg, F.W. Wehrli, Three-dimensional digital topological characterization of cancellous bone architecture, Int. J. Imag. Syst. Tech. 11 (2000) 81–90. [14] F.W. Wehrli, P.K. Saha, B.R. Gomberg, H.K. Song, P.J. Snyder, M. Benito, A. Wright, R. Weening, Role of magnetic resonance for assessing structure and function of trabecular bone, Top. Magn. Reson. Imag. 13 (2002) 335–356. [15] M. Par6tt, Implications of architecture for the pathogenesis and prevention of vertebral fracture, Bone 13 (1992) S41–47. [16] C.L. Gordon, C.E. Webber, P.S. Nicholson, Relation between image-based assessment of distal radius trabecular structure and compressive strength, Bioengineering 49 (1998) 390–397. [17] S. Sarkar, B.H. Mitlak, M. Wong, J.L. Stock, D.M. Black, K.D. Harper, Relationships between bone mineral density and incident vertebral fracture risk with raloxifene therapy, J. Bone Miner. Res. 17 (2002) 1–10. [18] C.H. Chesnut, 3rd, S. Silverman, K. Andriano, H. Genant, A. Gimona, S. Harris, D. Kiel, M. LeBo3, M. Maricic, P. Miller, C. Moniz, M. Peacock, P. Richardson, N. Watts, D. Baylink, A randomized trial of nasal spray salmon calcitonin in postmenopausal women with established osteoporosis: the prevent recurrence of osteoporotic fractures study, PROOF study group., Am. J. Med. 109 (2000) 267–276. [19] M. Kleerekoper, A.R. Villanueva, J. Stanciu, D. Sudhaker Rao, A.M. Par6tt, The role of three-dimensional trabecular microstructure in the pathogenesis of vertebral compression fractures, Calcif. Tissue Int. 37 (1985) 594–597. [20] E. Legrand, D. Chappard, C. Pascaretti, M. Duquenne, S. Krebs, V. Rohmer, M.F. Basle, M. Audran, Trabecular bone microarchitecture, bone mineral density and vertebral fractures in male osteoporosis, J. Bone Miner. Res. 15 (2000) 13–19. [21] J.E. Aaron, P.A. Shore, R.C. Shore, M. Beneton, J.A. Kanis, Trabecular architecture in women and men of similar bone mass with and without vertebral fracture: II. Three-dimensional histology, Bone 27 (2000) 277–282. [22] M. Audran, D. Chappard, E. Legrand, H. Libouban, M.F. Basle, Bone microarchitecture and bone fragility in men: DXA and histomorphometry in humans and in the orchidectomized rat model, Calcif. Tissue Int. 69 (2001) 214–217. [23] J. Wol3, Das Gesetz der Transformation der Knochen, A. Hirschwald, Berlin, 1892. [24] A. Odgaard, J. Kabel, B. van Rietbergen, M. Dalstra, R. Huiskes, Fabric and elastic principal directions of cancellous bone are closely related, J. Biomech. 30 (1997) 487–495. [25] S.A. Goldstein, R. Goulet, D. McCubbrey, Measurement and signi6cance of three-dimensional architecture to the
1944
[26]
[27] [28] [29] [30] [31] [32]
P.K. Saha, F.W. Wehrli / Pattern Recognition 37 (2004) 1935 – 1944 mechanical integrity of trabecular bone, Calcif. Tissue Int. 53 (1) (1993) S127–32, discussion S132-3. H.W. Chung, F.W. Wehrli, J.L. Williams, S.D. Kugelmass, S.L. Wehrli, Quantitative analysis of trabecular microstructure by 400 MHz nuclear magnetic resonance imaging, J. Bone Miner. Res. 10 (1995) 803–811. Z.H. Chu, J.P. Jopes, M. Singh, Foundation of Medical Imaging, Wiley, New York, 1993. A. Rosenfeld, A.C. Kak, Digital Picture Processing, Vol. 1, Academic Press, Inc., Orlando, FL, 1982. A. Rosenfeld, A.C. Kak, Digital Picture Processing, Vol. 2, Academic Press, Inc., Orlando, FL, 1982. J.K. Udupa, G.T. Herman, 3D Imaging in Medicine, CRC Press, Boca Raton, FL, 1991. M. Sonka, V. Hlavac, R. Boyle, Image Processing, Analysis, and Machine Vision, PWS Publishing, Paci6c Grove, CA, 1999. E.L. Hall, Computer Image Processing and Recognition, Academic Press, New York, NY, 1979.
[33] S.N. Hwang, F.W. Wehrli, J.L. Williams, Probability-based structural parameters from 3D NMR images as predictors of trabecular bone strength, Med. Phys. 24 (1997) 1255–1261. [34] H.K. Song, F.W. Wehrli, In vivo micro-imaging using alternating navigator echoes with applications to cancellous bone structural analysis, Magn. Reson. Med. 41 (1999) 947–953. [35] P.K. Saha, B.B. Chaudhuri, 3D digital topology under binary transformation with applications, Computer Vision Image Understand. 63 (1996) 418–429. [36] S.N. Hwang, F.W. Wehrli, Estimating voxel volume fractions of trabecular bone on the basis of magnetic resonance images acquired in vivo, Int. J. Imaging Syst. Technol. 10 (1999) 186–198. [37] B.R. Gomberg, P.K. Saha, F.W. Wehrli, Topology-based orientation analysis of trabecular bone networks, Med. Phys. 30 (2002) 158–168.
About the Author—PUNAM KUMAR SAHA received the B.S. and M.S. degrees in computer science and engineering from Jadavpur University, Calcutta, India, in 1987 and 1989, respectively, and the Ph.D. degree from the Indian Statistical Institute, Calcutta, India, in 1997. He joined the Indian Statistical Institute as a Faculty Member in 1993. In 1997, he joined the Department of Radiology, Medical Imaging Section, University of Pennsylvania, Philadelphia, as a Postdoctoral Fellow. He is currently an Assistant Professor at the same university. His present research interests include biomedical imaging problems and the application of their solutions, scale, segmentation, 6ltering, image registration, inhomogeneity correction, digital topology and their applications, and estimation of trabecular bone strength from MR images. He has published over 40 papers in international journals. Dr. Saha is an Associate Editor of the Pattern Recognition Journal. He is a Member of the International Association for Pattern Recognition and a Member of the Governing Body of the Indian Unit for Pattern Recognition and Arti6cial Intelligence. He received a Young Scientist Award from the Indian Science Congress Association in 1996. About the Author—FELIX W. WEHRLI received the Ph.D. degree in physical organic chemistry from the Swiss Federal Institute of Technology, Zurich, Switzerland, in 1970. He has spent his entire research career in the 6eld of magnetic resonance in leading research and development positions at Varian Associates, Bruker Instruments, and General Electric Medical Systems. He is currently Professor of Radiologic Science, Biochemistry, and Biophysics and the Director of the Laboratory for Structural NMR Imaging in the Department of Radiology of the University of Pennsylvania, Philadelphia. Over the past 20 years, his research has been concerned with the development of magnetic resonance procedures for clinical diagnosis, with particular focus on the understanding of the relationship between structure and function of connective tissues. A major portion of his recent work has dealt with the development of image acquisition and processing tools for characterizing structure and strength of trabecular bone. He is the author or coauthor of over 120 scienti6c articles, three textbooks, and 30 book chapters and reviews. He has been the Editor-in-Chief of Magnetic Resonance in Medicine, the leading magnetic resonance research journal since 1991. Dr. Wehrli is a Fellow of the International Society for Magnetic Resonance in Medicine.
A robust method for measuring trabecular bone orientation anisotropy at in vivo resolution using tensor scale Punam Kumar Saha∗ , Felix W. Wehrli Department of Radiology, University of Pennsylvania, 4th Floor, Blockley Hall, 423 Guardian Drive, Philadelphia, PA 19104, USA Received 17 November 2003; accepted 19 December 2003
Abstract Trabecular bone (TB) is a network of interconnected struts and plates that constantly remodels to adapt dynamically to the stresses to which it is subjected in such a manner that the trabeculae are oriented along the major stress lines (Wol3’s Law). Next to bone density, TB anisotropy has been found to be one of the most signi6cant determinants of the bone’s biomechanical behavior. Typically, orientational anisotropy of TB networks is expressed in terms of the fabric tensor, obtained by measuring the mean intercept length between structure elements along test lines. This method, however, can provide only a global statistical average of TB orientation anisotropy and, in general, requires a large sampling volume. Here, we present a new method, based on the recently conceived notion of “tensor scale”, which provides regional information of TB orientation anisotropy. Regional structure is represented by local best 6t ellipsoid (ellipse in 2D) and the structural orientation is determined from the eigenvectors along the semi-axes. The method is found to be remarkably robust over a wide range of resolution regimes and image rotation as shown with micro-CT images from specimens of the human distal radius, the latter showing the characteristic di3erences in structural anisotropy for transverse and longitudinal sections. Finally, the method’s applicability to in vivo MR imaging is demonstrated with data from the distal tibia. ? 2004 Pattern Recognition Society. Published by Elsevier Ltd. All rights reserved. Keywords: Tensor scale; Fuzzy subset; Trabecular bone; Orientation; Anisotropy; -MRI-CT;
1. Introduction Analysis of regional architecture is of paramount interest in many medical imaging applications. This paper illustrates an application of the recently developed “tensor scale” concept [1] for quantitatively measuring tissue architectural parameters. Although tensor scale here is applied to the analysis of local structure morphology, the notion of scale has a long history in applications to computer vision and image processing. “Scale” may be thought of as the spatial resolution, or, more generally, a range of resolutions needed to ensure a suDcient yet compact object representation. Scale plays an important role in determining ∗ Corresponding author. Tel.: +1-215-662-6780; fax: +1-215898-9145. E-mail address: [email protected] (P.K. Saha).
the optimum trade-o3 between noise smoothing and perception of structures. Also, scale is helpful in breaking a computer vision and image-processing task into a hierarchy of tasks where those at higher levels deal with larger structures. Marr [2] pointed out many bene6ts of multi-scale image analysis and computer vision. Although the concept of scale originally evolved from scale-space theory [2–5], recently, the concept of “local scale” [6–8] has attracted signi6cant interest as it allows space-variant process control 6tting with local structures. The notion of local morphometric scale was presented in [7,8] using the spherical model. A major limitation of this model is its ignoring orientation and anisotropy, which makes it unsuited for many applications, especially in biomedical imaging where structures are inherently anisotropic. The local scale using a tensor model [1] (an ellipse in 2D and an ellipsoid in 3D) provides a uni6ed and parametric representation of size, orientation,
0031-3203/$30.00 ? 2004 Pattern Recognition Society. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.patcog.2003.12.019
1936
P.K. Saha, F.W. Wehrli / Pattern Recognition 37 (2004) 1935 – 1944
and anisotropy of local structures. Although tensor scale is a concept with wide-ranging applications in computer vision and image processing, here we focus on its application in quantitative analysis of regional morphology only; a few other applications of tensor scale have been reported in [9,10]. The focus of this paper is to explore the potential of tensor scale to analyze global and regional trabecular bone (TB) morphology, although the method may also be applicable to other tissues, e.g., the characterization of vascular micro architecture in during angiogenesis [11,12]. TB consists of a meshwork of interconnected struts and plates providing the bone high mechanical strength at minimum weight. Quantitative analysis [13] of trabecular bone structure by in vivo MR imaging [14] will be useful for fracture risk evaluation, early assessment of disease status, therapy response monitoring and screening of high risk subjects. The most common bone disorder is osteoporosis, which leads to enhanced bone fragility and increase the risk of fracture caused by both bone loss and micro architectural deterioration of bone tissue [15]. It is generally agreed [16] that 50–60% of the bone’s mechanical competence can be explained by variations in the apparent density (bone mass/tissue volume). However, clinical and experimental evidence for an independent contribution from trabecular architecture is compelling [17,18]. A few studies [19–22] involving histomorphometry in patients with and without vertebral fractures, matched for gender and bone mineral density (BMD), found the two groups of subjects to di3er in histomorphometric indices. Trabecular bone remodels, i.e., it adapts to the stresses to which it is subjected (Wol3’s Law [23]) in such a manner that the trabeculae follow the major stress lines. Thus, TB at most skeletal locations is highly anisotropic. Structural anisotropy (fabric) can be expressed in terms of a tensor and, although intuitively obvious, it has recently been shown that the stress tensor is aligned with the fabric tensor [24]. Next to material density (i.e., bone volume fraction) structural anisotropy most critically determines TB’s mechanical behavior [25]. Classically, the fabric tensor can be obtained from images by determining the mean intercept length (MIL) [26] of test lines along all spatial directions. This method, however, requires a suDciently large number of intercepts, a requirement that is only met for suDciently large sampling volume. Also, the MIL method provides only a statistical average of orientation anisotropy. This paper introduces a new method for computation of global and regional TB orientation anisotropy using the recently developed tensor scale theory and presents results of preliminary evaluative experiments. Speci6cally, the robustness of the method under varying resolutions and image rotations is studied and the characteristic differences between TB anisotropy computed from transverse and longitudinal sections are examined. Finally, potential applications of the method to in vivo TB MR images are demonstrated.
Fig. 1. Illustration of tensor scale computation at an image point (the red dot). Finitely many radially outward directed straight-line segments, called sample lines (blue), distributed at approximately uniform angular spacing, are traced and the optimum edge location (orange and pink dots) is identi6ed on each sample line. Exploiting the axial symmetry of the ellipse, the edge points on every pair of radially opposite sample lines are repositioned (orange and yellow dots) by selecting the one closer to the candidate point and then reMecting it onto the opposite sample line. Finally, tensor scale is computed by 6tting an ellipse to the repositioned edge points.
2. Tensor scale: theory and algorithms Tensor scale at an image point p is the parametric representation of the largest ellipse (ellipsoid in 3D) centered at p that is contained in the same object region. The process of tensor scale computation (in 2D) at any pixel p is schematically illustrated in Fig. 1. The method is based on tracing pairs of radially opposite sample lines emanating from p and locating the closest edge on every sample line. After locating edges on all sample lines, they are repositioned by making use of the ellipse’s axial symmetry. Finally, the tensor scale at p is determined by computing the best 6t ellipse to these repositioned edge points. The method involves the following four sequential steps—(1) intensity computation along each sample line, (2) location of the optimum edge on each sample line, (3) repositioning of the edge locations, and (4) computation of the best 6t ellipse to the repositioned edge locations. In the following, we describe the method in more detail. Although the method immediately extends to any 6nite dimension, we con6ne our description to 2D, the dimensionality of the present implementation. A digital object O=(Z 2 ; O ) is a fuzzy subset of Z 2 , where Z is the set of all integers and Z 2 represents a 2D digital
P.K. Saha, F.W. Wehrli / Pattern Recognition 37 (2004) 1935 – 1944
grid; and O : Z 2 → [0; 1] speci6es the membership value at each grid point, often referred to as a pixel. As pointed out earlier, the current application aims to compute orientation anisotropy of trabecular bone from -CT or in vivo -MR images. Following our target application, a digital object O will represent a trabecular bone network and at any pixel p, the membership value O will represent its bone volume fraction (BVF). An imaging system [27] acquires images containing information of a target object, often along with other co-objects. The process of extracting an object from an image, commonly referred to as image segmentation [28–31], has extensively been studied for decades. Following the main focus of the theory and methods addressed in this paper, we start with the de6nition of digital objects instead of digital images. However, in Methods, we will describe how BVF maps are generated for di3erent experiments from original image data. For now, let us use the above de6nition of digital objects. The support (O) of a digital object O is the set of all pixels with non-zero membership. Tensor scale is computed for all pixels p ∈ (O). (p) (p) (p) Let (1(p) ; 1 ); (2(p) ; 2 ); · · · ; (m(p) ; m ), where (p) (p) i ; i : [0; L] → R2 ; ∀i, denote m pairs of mutually opposite sample lines each having length L and emanating from the candidate pixel p ∈ Z 2 . The sample lines are intended to be uniformly distributed over the angular space around p to ensure that the edge points computed on these sample lines, eventually utilized to compute the tensor scale at p, are not skewed in any direction. A set of sample lines skewed in one direction may lead to a bias in the computed tensor scale in terms of its orientation. For example, if the sample lines were more densely distributed around the horizontal direction then the computed tensor scale would be less sensitive to vertical components of the local structure; instead, it could be biased to horizontal components of the structure. Each sample line is coupled with a radially opposite sample line; we thus refer to each such pair as a conjugate pair and to each sample line of the pair as a complementary sample line of the other. The membership pro6le on a sample line is computed by interpolating object membership values at points selected at regular intervals. During the second step, the optimum edge location is identi6ed on each sample line. For a raw (unsegmented) image, the edge point is located at the transition between two objects, i.e. at the zero crossing of the second-order derivative of intensity. Typically, voxel intensity is directly proportional to local object material density. However, at limited resolution, most voxels are only partially occupied by the object material, hence yielding only a fuzzy representation of the object. In these situations it is appropriate to presegment the images as a fuzzy object. For such fuzzily segmented objects, the edge point is not detected at the closest intersection between the boundary of the support of the object and the sample line. Rather, the fuzzy length of the locally connected object along a sample line is computed. Speci6cally, the object edge point i(p) on a sample line i(p) , emanating from p, is computed whose distance from p,
1937
denoted by i(p) −p, is equal to the integration of the membership values along i(p) until a non-object point (membership value = 0) is found. In other words, T (p;(p) ) (p) i di (t) (p) (p) i − p = O (i (t)) (1) dt; dt 0 where, T (p; i(p) ) is the parametric location of the 6rst (from p) zero membership point on i(p) . The edge location on the sample line i(p) is denoted by i(p) . The edge points, located as above, are intended to roughly describe the boundary of the ellipse representing the tensor scale at p. From the axial symmetry of an ellipse, two radially opposite points on its boundary are always equidistant from its center. However, the points along the boundary of a local structure do not always satisfy this condition of symmetry with respect to a candidate pixel. For example, in Fig. 1, the edge locations on the west- and north-bound sample lines are in general farther from the candidate pixel than the edge points on respective complementary sample lines. Therefore, to ensure axial symmetry of an ellipse, the edge points need to be repositioned. This step is accomplished by analyzing the edge points on every conjugate pair of sample lines. Speci6cally, the one closer to p between i(p) and
i(p) is selected and reMected on the complementary sample line. In Fig. 1, the edge points shown by pink dots are repositioned yielding the yellow dots; orange dots represent the edge points which remain unmoved after the repositioning process. The repositioned edge points will brieMy be referred to as r-edge points. The 6nal step in tensor scale computation is to 6t an ellipse to the r-edge points. This step is accomplished in two sub-steps—(i) determination of ellipse orientation, (ii) computation of the lengths of two semi-axes. The 6rst sub-step is completed using principal component analysis [32] of the r-edge points and assuming p as the center. Speci6cally, the orientation of the ellipse is computed as the angle of the eigenvectors associated to the larger eigenvalue of the covariance matrix de6ned as follows 1 1; 1 = (xi − xp )2 ; 2m i=1;2;:::;2m 2; 2 =
1 2m
(yi − yp )2 ;
(2)
i=1;2;:::;2m
and, 1; 2 = 2; 1 =
1 2m
(xi − xp )(yi − yp );
i=1;2;:::;2m
where (xi ; yi ) | i = 1; 2; : : : ; 2m are the coordinates of the r-edge points and (xp ; yp ) is the coordinate of p. Finally, the actual lengths a; b of the two semi axes of tensor scale ellipse are computed by minimizing the following error function u2 v2 ferror = 1 − i2 − i2 ; (3) a b i=1;2;:::;2m
1938
P.K. Saha, F.W. Wehrli / Pattern Recognition 37 (2004) 1935 – 1944
where, (ui ; vi ) | i = 1; 2; : : : ; 2m are the r-edge points after (i) applying a translation by −p (i.e., p is translated to origin), and then (ii) applying a rotation by angle − (i.e., the major semi axis of the ellipse is aligned to the horizontal coordinate axis). By setting the partial derivatives of ferror with respect to a and b to zero, it can be shown that the a and b which minimize ferror are given by the following expressions
was selected in this experiment to examine whether similar anisotropies are obtained in two di3erent sets of longitudinal slices within the same bone volume (which assumes that the trabecular structure is homogeneous within a given sample). In Experiment 4, designed to evaluate whether the technique was able to delineate local structural orientation in -MR images transverse MR images acquired in the distal
a=
xi2 yi2 i=1; 2; :::; 2m xi2 yi2 − i=1; 2; :::; 2m xi4 i=1; 2; :::; 2m yi4 ; 2 2 2 2 4 i=1; 2; :::; 2m xi yi i=1; 2; :::; 2m yi − i=1; 2; :::; 2m xi i=1; 2; :::; 2m yi i=1; 2; :::; 2m
b=
xi2 yi2 i=1; 2; :::; 2m xi2 yi2 − i=1; 2; :::; 2m xi4 i=1; 2; :::; 2m yi4 : 2 2 2 4 2 i=1; 2; :::; 2m xi i=1; 2; :::; 2m xi yi − i=1; 2; :::; 2m xi i=1; 2; :::; 2m yi
i=1; 2; :::; 2m
Finally, we shall use (p) to denote the tensor scale ellipse at any pixel p and Y1 (E) ¿ Y2 (E) to denote the two semi-axes of an ellipse E.
3. Methods Four di3erent experiments were conducted for preliminary evaluation of the proposed orientation anisotropy computation method and to demonstrate its potential in MR imaging applications of in vivo human trabecular bone: Experiments 1 and 2 were designed to evaluate the sensitivity of the computed orientation anisotropy to resolution and object rotation. The purpose of Experiment 3 was to demonstrate the characteristic di3erences in orientation anisotropy between transverse and longitudinal image slices. Finally, Experiment 4 focused on regional orientation anisotropy computed from 2D MR images of in vivo human trabecular bone. For all experiments, the following sequence of steps were applied—(1) computation of a BVF map from an acquired raw image, (2) tensor scale computation at every image location, and (3) computation of global or regional orientation anisotropy. 3.1. Material description and BVF map computation For Experiments 1–3, a -CT image was selected from a 3D stack of images that had been acquired previously (SCANCO Medical -CT 20, at 22 m isotropic voxel size) as part of a di3erent project [33] from samples cored from cadaveric bone of the distal radial metaphysis. The sample consisted of a cylindrical core of 9 mm height and 9 mm diameter with the cylinder axis parallel to the direction of the radius. The -CT image data set was processed to generate BVF maps. For Experiments 1 and 2, two central longitudinal slices, separated by 330 m (i.e., 15 slices), were selected. For Experiment 3, two additional longitudinal slices, shifted by 15 slices from the previous set, and two central transverse slices, separated by 15 slices from each other, were selected. A di3erent set of longitudinal slices
(4)
(5)
tibia (137 × 137 × 410 m3 voxel size) in two women, ages 53 and 49 years, using a 3D partial Mip-angle pulse sequence [34] were processed as described above. The objective of these experiments was to evaluate whether the method is able to quantify regional structural orientation from in vivo MR images. Because of high spatial resolution exceeding trabecular thickness by a factor of 5, the -CT images had a bimodal intensity distribution and thus could be segmented by selecting a threshold at the midpoint of the two modes representing the mean bone and mean marrow intensities. Subsequently, the largest 26-connected component [35] within the region thresholded for bone was computed. The intensities at the two modes were selected as the mean intensities for bone and marrow (denoted by B and M , respectively). At each pixel p, the BVF value f (p) was computed using the following equation: 0 if f(p) ¡ M ; f(p) − M if M 6 f(p) ¡ B ; (6) f (p) = B − M 1 otherwise: For Experiment 1, the BVF images of the selected slices were resampled at pixel sizes corresponding to integral multiples of the parent resolution, i.e., 44 × 44, 66 × 66, and 88 × 88 m2 . This strategy was chosen for generating 2D trabecular bone images at lower resolutions (rather than resampling the raw -CT slice data sets followed by computing the BVF images). The latter approach would cause signi6cant loss in structural detail in the segmented images using the simple threshold-based method used here. Although more sophisticated segmentation could be used, no other segmentation were attempted since this was deemed to be outside the scope of the present paper. For Experiment 2, the rotation of the raw 2D -CT images was performed before computing BVF since in this manner the maximum e3ect of rotation on the entire process can be assessed. Anticlockwise rotations of ∈ {15◦ ; 30◦ ; 45◦ ; 60◦ } were applied to each of the selected raw 2D -CT images around the center of each slice and the rotated images were
P.K. Saha, F.W. Wehrli / Pattern Recognition 37 (2004) 1935 – 1944
1939
redigitized into a square grid using linear interpolation. For both Experiments 2 and 3, BVF images were computed using Eq. (6). For the in vivo MR experiment, an actual region of interest (ROI), for which the tensor scale was computed, was generated by manually tracing boundaries on each of the two selected 2D images using an IDL (Interactive Data Language, Research Systems Inc.)-based graphical interface running on a Pentium III PC under MS Windows OS. The data encompassed by the ROIs were preprocessed by deshading and noise reduction using a histogram deconvolution method [36] to produce BVF images. 3.2. Tensor scale computation Tensor scale was computed at every image pixel as described in Section 2. Two parameters are required to compute tensor scale—(1) length L of each sample line, and (2) number of sample lines 2m. L should be suDciently large compared to trabecular thickness to capture local structural anisotropy. Here, we used L = 704 m or 32 pixels which is approximately 4–7 times larger than typical human trabecular thickness. The number of sample lines determines the precision of the computed tensor scales. Following the recommendation in [1], 30 sample lines were used for the experiments described. 3.3. Regional orientation and anisotropy Here, we describe the method to compute global as well as regional orientation of a trabecular bone network from its tensor scale map. Let p be any pixel inside the support (O) of the object O. The principal orientation of the trabecula at p may be de6ned by the unit vector Y1 (p)=|Y1 (p)|. Let H (i) | i ∈ {1; 2; : : : ; 180} represent the orientation distribution over the entire object where i is the value of an angle (anti-clockwise) in degrees. H (i) is computed as follows: p (i) H (i) = O (p) ; (7) |Y1 (p)| p∈(O)
where p (i) is the length of the “central chord” (a chord passing through the center) of the tensor scale ellipse (p) at p along the direction making i degrees anti-clockwise angle with a horizontal line. The length of p (i) can be expressed as follows: |Y1 (p)|2 × |Y2 (p)|2 ; (8) p (i) = |Y2 (p)|2 cos i + |Y1 (p)|2 sin i The value of p (i) as obtained by using the above equation is illustrated in Fig. 2. An ellipse, denoted by H , is computed from the rose plot of H (i) using a method similar to the one described in Section 2. Let Y1; H ¿ Y2; H denote the two semi-axes of H . Global orientation anisotropy of the TB network then is computed as 1 − |Y2; H |2 =|Y1; H |2 —the eccentricity of the ellipse H . Regional orientation anisotropy at a pixel p is obtained by calculating the orientation histogram Hp over
Fig. 2. Illustration of the length of the “central chord” (a chord passing through the center) of an ellipse along a given direction as de6ned in Eq. (8).
a neighborhood around p which, in this paper, has been chosen as a disk of radius 25 pixels. 4. Results In this section, we present results of the experiments described in the previous section. The motivation of the 6rst experiment was to assess the performance of tensor scale-based orientation anisotropy computation in various resolution regimes. Recall that four di3erent pixel sizes (22 × 22, 44 × 44, 66 × 66, and 88 × 88 m2 ) were chosen for this experiment. Figs. 3a,b show an original -CT image slice at 22 × 22 m2 voxel size and its BVF map computed as per Eq. (6). The BVF mapped image of the same slice resampled at 88 × 88 m2 resolution is presented in Fig. 3c. Tensor scale images, unlike other parametric images (e.g. BVF), are not scalar maps. Therefore, it is important to adopt a new display scheme that can adequately convey the visualization of the architectural make-up of the tissue. Since an ellipse centered at the origin is uniquely de6ned by three parameters—orientation of the major semi axis Y1 (p), anisotropy (= 1 − |Y2 (p)|2 =|Y1 (p)|2 ), and the length of the short axis (representing structure thickness) (=|Y2 (p)|) —we have adopted hue-saturation-intensity (HSI) color coding scheme to represent a 2D tensor scale image. Specifically, hue (H), saturation (S), and intensity (I) of the color vector at any pixel p encodes orientation, anisotropy, and thickness, respectively, of (p). Fig. 3d illustrates the color coding at full intensity (i.e., thickness is maximum). According to this coding convention, colors near the center of the disk represent locally isotropic structures while those near the periphery of the disk represent highly anisotropic structures. As illustrated in the disc, a color is repeated at 180◦ representing the axial symmetry of an ellipse.
1940
P.K. Saha, F.W. Wehrli / Pattern Recognition 37 (2004) 1935 – 1944
Fig. 3. E3ects of varying resolution on tensor scale-based TB orientation anisotropy computation: (a) an image slice from a 3D -CT dataset of a cadaveric TB specimen from the human distal radius acquired at 22 m isotropic voxel size; (b) a BVF mapped image computed from (a); (c) the same image as (b) but resampled at 88 × 88 m2 pixel size; (d) HSI color coding at full intensity used to represent tensor scale images; (e,f) color coded tensor scale images computed from (b,c), respectively; (g) rose plots and anisotropy ellipses of TB orientation distributions computed from tensor scale images at di3erent resolutions. The maximum variation in computed anisotropy for the di3erent resolutions examined was 0.70%.
Normalized thickness values are computed by dividing |Y2 (p)| with the maximum scale value L. HSI color-coded tensor scale images computed from Fig. 3b,c are illustrated in Fig. 3e,f, respectively. It may be noted that, although the tensor scale image at the lower resolution is somewhat blurred, orientation, anisotropy, and thickness of matching regions are well preserved. Rose plots of tensor scale-based orientation distributions computed from the sets of images at di3erent resolutions (22 × 22 m2 , 44 × 44 m2 , 66 × 66 m2 , and 88 × 88 m2 pixel size) are illustrated in Fig. 3g. The computed anisotropy ellipses, each 6tted to a rose plot, are also shown in the same 6gure. The maximum variation in anisotropy values at the four di3erent resolutions was 0.70%. The behavior of the proposed method under image rotation is illustrated by using the 2D image displayed in Fig. 4a, obtained by applying 30◦ anti-clockwise rotation on the raw the -CT TB slice image of Fig. 3a. The BVF map computed from Fig. 4a is shown in Fig. 4b while the tensor scale image of Fig. 4b is displayed in Fig. 4c using HSI color coding. It may be noted by comparing Fig. 4c with Fig. 3e that the color (orientation) of matching structures are signi6cantly di3erent in these two images. This mismatch is caused as the absolute orientation of a structure in an image is changed when the image itself is rotated and the amount
of rotation of the structure is exactly the same as that originally applied to the image. Therefore, it was expected that if the rotation applied to the original image were subtracted from tensor scale orientation of each pixel, the modi6ed tensor scales of matching structures in the original and in the rotated images would be similar. This is borne out by tensor scale image obtained by applying the orientation adjustment displayed in Fig. 4d. Rose plots of orientation distribution computed from the set of slices without rotation and those at 30◦ anticlockwise rotation are depicted in Fig. 4e; computed anisotropy ellipses are also shown in the same 6gure. The maximum variation in anisotropy values at di3erent rotations is 1.69%. The third experiment shows that the tensor scale-based orientation anisotropy analysis can reveal the known differences in structural anisotropy between longitudinal and transverse sections in the distal radius in images taken from the same 3D data. Trabeculae are preferentially oriented along the principal load direction and, since the cylindrical bone specimens were cored with their axes aligned with the radial shaft, greater orientation anisotropy is expected in longitudinal than in transverse sections [33,37]. Fig. 5a,b displays transverse and longitudinal slices selected from the same -CT data of a cadaveric TB specimen. Fig. 5c,d shows the BVF mapped images computed from Fig. 5a,b,
P.K. Saha, F.W. Wehrli / Pattern Recognition 37 (2004) 1935 – 1944
1941
Fig. 4. E3ects of image rotation on tensor scale-based TB orientation anisotropy computation: (a) the image of Fig. 3a after applying a 30◦ anticlockwise rotation; (b) a BVF mapped image computed from (a); (c) tensor scale image of (a); (d) the image of (c) after applying rotation-adjustment; (e) orientation rose-plots and anisotropy ellipses computed from respective tensor scale images demonstrating the method’s robustness to image rotation.
Fig. 5. Tensor scale-based TB orientation anisotropy shown for transverse and longitudinal image slices: (a,b) transverse and a longitudinal images from the same -CT image dataset; (c,d) BVF mapped images computed from (a,b); (d,e) tensor scale images derived from (c,d); (f) orientation rose-plots and anisotropy ellipses computed from respective tensor scale images.
1942
P.K. Saha, F.W. Wehrli / Pattern Recognition 37 (2004) 1935 – 1944
Fig. 6. Results of regional orientation analysis of in vivo TB -MR images using tensor scale: (a) cross-sectional image of the left ankle (one of 28 contiguous slices of 410 m thickness and 137 × 137 m2 pixel size from a 28-slice 3D dataset) acquired at the level of the malleolus (female, age 53 years). The foot is tilted approximately 30◦ ; (b) BVF image computed from (a) over a manually selected region; (c) color-coded tensor scale image of TB in (b); (d) regional orientation and anisotropy computed from (c); (e–h) same as (a–d) for a di3erent subject (female, age 49 years).
respectively, and Fig. 5e,f illustrate their color-coded tensor scale images. As visually apparent, the tensor scale image in Fig. 5e displays color almost randomly distributed over the entire color band while that in Fig. 5f shows greater frequency of colors between green (60◦ ) and blue (120◦ ) representing mostly vertical trabeculae (90◦ ). This trend is clearly visible in the rose-plots, shown in Fig. 5g, of the orientation distributions computed from the tensor scale images of selected sets of transverse and longitudinal slices. The computed orientation anisotropy values (0.49 and 0.88 for the sets of longitudinal and transverse slices, respectively) agreed well with these observations. Fig. 6 shows in vivo trabecular bone images from two subjects (female, age 53 and 49 years, respectively) taken in the distal tibia at the level of the malleolus (a lateral protrusion above the ankle), along with BVF maps and parametric images (tensor-scale and regional anisotropy maps). The images show distinct structural anisotropy with the trabecular orientation roughly parallel to the cortical rim. This structural arrangement reMects the multidirectional lateral forces acting on the ankle. 5. Conclusions A new tensor scale-based anisotropy computation method applied to images of trabecular bone has been presented. The method allows computation of tensor scale at any image location in a fuzzily segmented object representation and thus is applicable low resolution in vivo MR images. Since structural anisotropy has implications on the bone’s mechanical
properties, the method may have potential for studying regional changes in structure orientation during bone remodeling on the basis of in vivo MR images.
6. Summary Trabecular bone (TB) is a network of interconnected struts and plates that constantly remodels to adapt dynamically to the stresses to which it is subjected in such a manner that the trabeculae are oriented along the major stress lines (Wol3’s Law). Next to bone density, TB anisotropy has been found to be one of the most signi6cant determinants of the bone’s biomechanical competence. Typically, orientation anisotropy of TB networks is expressed in terms of the fabric tensor, obtained by measuring the mean intercept length between structure elements along test lines. This method, however, can provide only a global statistical average of TB orientation anisotropy and, in general, requires a large sampling volume. Here, we present a new method, based on the recently conceived notion of “tensor scale”, which provides regional information of TB orientation anisotropy. The method’s performance has been evaluated in terms of its sensitivity to resolution and image rotation and the characteristic di3erences between TB anisotropy computed from transverse and longitudinal sections have been examined. Finally, potential applications of the method to in vivo TB MR images are demonstrated. Although the method can be extended to any 6nite dimension, it is currently implemented only in two dimensions.
P.K. Saha, F.W. Wehrli / Pattern Recognition 37 (2004) 1935 – 1944
Tensor scale at an image point p is the parametric representation of the largest ellipse (ellipsoid in three dimensions) centered at p that is contained in the same object region. BrieMy, the method is based on tracing pairs of radially opposite sample lines emanating from p and locating the closest edge on every sample line. After locating edges on all sample lines, they are repositioned by making use of the ellipse’s axial symmetry. Finally, the tensor scale at p is determined by computing the best 6t ellipse to these repositioned edge points. Orientation anisotropy is computed by analyzing orientation of the major semi axis and the eccentricity of tensor scale ellipse at every image pixel location with non-zero fractional occupancy of TB. Experiments were conducted on the basis of microcomputed tomography images acquired at 22×22×22 m3 voxel size of a cadaveric TB specimen cored from the human distal radius with the axis of the cylinder aligned with that of the radius. The experimental results demonstrate the method to be remarkably robust over a wide range of resolution regimes (22 × 22, 44 × 44, 66 × 66, and 88 × 88 m2 ) and image rotations (anticlockwise rotations of 15◦ , 30◦ , 45◦ , 60◦ ). The maximum variation in computed tensor scale-based orientation anisotropy at di3erent resolutions examined was 0.70% while that for di3erent rotations was 1.69%. The method has been found to be successful to quantitatively characterize the di3erences in structural anisotropy for transverse and longitudinal sections of the same micro-CT data. Finally, the method’s applicability to in vivo MR imaging is demonstrated with data from the distal tibia. References [1] P.K. Saha, Novel theory and methods for tensor scale: a local morphometric parameter, in: Proceedings of SPIE: Medical Imaging, San Diego, CA, Vol. 5032, 2003, pp. 314–324. [2] D. Marr, Vision, W.H. Freeman and Company, San Francisco, CA, 1982. [3] A.P. Witkin, Scale-space 6ltering, in: Proceedings of the Eighth International Joint Conference Arti6cial Intelligence, Karlsruhe, West Germany, 1983, pp. 1019–1022. [4] J.J. Koenderink, The structure of images, Biol. Cybern. 50 (1984) 363–370. [5] T. Lindeberg, Scale-space for discrete signals, IEEE Trans. Pattern Recognition Machine Intell. 12 (1990) 234–254. [6] S.M. Pizer, D. Eberly, D.S. Fritsch, Zoom-invariant vision of 6gural shape: the mathematics of core, Comput. Vision Image Understand. 69 (1998) 55–71. [7] P.K. Saha, J.K. Udupa, D. Odhner, Scale-based fuzzy connected image segmentation: theory, algorithms, and validation, Computer Vision Image Understand. 77 (2000) 145–174. [8] P.K. Saha, J.K. Udupa, Scale based image 6ltering preserving boundary sharpness and 6ne structures, IEEE Trans. Med. Imaging 20 (2001) 1140–1155. [9] P.K. Saha, J.C. Gee, Z. Xie, J.K. Udupa, Tensor scale-based image registration, in: Proceedings of SPIE: Medical Imaging, San Diego, CA, Vol. 5032, 2003, pp. 743–753.
1943
[10] P.K. Saha, J.K. Udupa, Tensor scale-based fuzzy connectedness image segmentation, in: Proceedings of SPIE: Medical Imaging, San Diego, CA, Vol. 5032, 2003, pp. 1580–1590. [11] P. Carmeliet, R.K. Jain, Angiogenesis in cancer and other diseases, Nature 407 (2000) 249–257. [12] E.M. Brey, T.W. King, C. Johnston, L.V. McIntire, G.P. Reece, C.W. Patrick, A technique for quantitative three-dimensional analysis of microvascular structure, Microvasc. Res. 63 (2002) 279–294. [13] P.K. Saha, B.R. Gomberg, F.W. Wehrli, Three-dimensional digital topological characterization of cancellous bone architecture, Int. J. Imag. Syst. Tech. 11 (2000) 81–90. [14] F.W. Wehrli, P.K. Saha, B.R. Gomberg, H.K. Song, P.J. Snyder, M. Benito, A. Wright, R. Weening, Role of magnetic resonance for assessing structure and function of trabecular bone, Top. Magn. Reson. Imag. 13 (2002) 335–356. [15] M. Par6tt, Implications of architecture for the pathogenesis and prevention of vertebral fracture, Bone 13 (1992) S41–47. [16] C.L. Gordon, C.E. Webber, P.S. Nicholson, Relation between image-based assessment of distal radius trabecular structure and compressive strength, Bioengineering 49 (1998) 390–397. [17] S. Sarkar, B.H. Mitlak, M. Wong, J.L. Stock, D.M. Black, K.D. Harper, Relationships between bone mineral density and incident vertebral fracture risk with raloxifene therapy, J. Bone Miner. Res. 17 (2002) 1–10. [18] C.H. Chesnut, 3rd, S. Silverman, K. Andriano, H. Genant, A. Gimona, S. Harris, D. Kiel, M. LeBo3, M. Maricic, P. Miller, C. Moniz, M. Peacock, P. Richardson, N. Watts, D. Baylink, A randomized trial of nasal spray salmon calcitonin in postmenopausal women with established osteoporosis: the prevent recurrence of osteoporotic fractures study, PROOF study group., Am. J. Med. 109 (2000) 267–276. [19] M. Kleerekoper, A.R. Villanueva, J. Stanciu, D. Sudhaker Rao, A.M. Par6tt, The role of three-dimensional trabecular microstructure in the pathogenesis of vertebral compression fractures, Calcif. Tissue Int. 37 (1985) 594–597. [20] E. Legrand, D. Chappard, C. Pascaretti, M. Duquenne, S. Krebs, V. Rohmer, M.F. Basle, M. Audran, Trabecular bone microarchitecture, bone mineral density and vertebral fractures in male osteoporosis, J. Bone Miner. Res. 15 (2000) 13–19. [21] J.E. Aaron, P.A. Shore, R.C. Shore, M. Beneton, J.A. Kanis, Trabecular architecture in women and men of similar bone mass with and without vertebral fracture: II. Three-dimensional histology, Bone 27 (2000) 277–282. [22] M. Audran, D. Chappard, E. Legrand, H. Libouban, M.F. Basle, Bone microarchitecture and bone fragility in men: DXA and histomorphometry in humans and in the orchidectomized rat model, Calcif. Tissue Int. 69 (2001) 214–217. [23] J. Wol3, Das Gesetz der Transformation der Knochen, A. Hirschwald, Berlin, 1892. [24] A. Odgaard, J. Kabel, B. van Rietbergen, M. Dalstra, R. Huiskes, Fabric and elastic principal directions of cancellous bone are closely related, J. Biomech. 30 (1997) 487–495. [25] S.A. Goldstein, R. Goulet, D. McCubbrey, Measurement and signi6cance of three-dimensional architecture to the
1944
[26]
[27] [28] [29] [30] [31] [32]
P.K. Saha, F.W. Wehrli / Pattern Recognition 37 (2004) 1935 – 1944 mechanical integrity of trabecular bone, Calcif. Tissue Int. 53 (1) (1993) S127–32, discussion S132-3. H.W. Chung, F.W. Wehrli, J.L. Williams, S.D. Kugelmass, S.L. Wehrli, Quantitative analysis of trabecular microstructure by 400 MHz nuclear magnetic resonance imaging, J. Bone Miner. Res. 10 (1995) 803–811. Z.H. Chu, J.P. Jopes, M. Singh, Foundation of Medical Imaging, Wiley, New York, 1993. A. Rosenfeld, A.C. Kak, Digital Picture Processing, Vol. 1, Academic Press, Inc., Orlando, FL, 1982. A. Rosenfeld, A.C. Kak, Digital Picture Processing, Vol. 2, Academic Press, Inc., Orlando, FL, 1982. J.K. Udupa, G.T. Herman, 3D Imaging in Medicine, CRC Press, Boca Raton, FL, 1991. M. Sonka, V. Hlavac, R. Boyle, Image Processing, Analysis, and Machine Vision, PWS Publishing, Paci6c Grove, CA, 1999. E.L. Hall, Computer Image Processing and Recognition, Academic Press, New York, NY, 1979.
[33] S.N. Hwang, F.W. Wehrli, J.L. Williams, Probability-based structural parameters from 3D NMR images as predictors of trabecular bone strength, Med. Phys. 24 (1997) 1255–1261. [34] H.K. Song, F.W. Wehrli, In vivo micro-imaging using alternating navigator echoes with applications to cancellous bone structural analysis, Magn. Reson. Med. 41 (1999) 947–953. [35] P.K. Saha, B.B. Chaudhuri, 3D digital topology under binary transformation with applications, Computer Vision Image Understand. 63 (1996) 418–429. [36] S.N. Hwang, F.W. Wehrli, Estimating voxel volume fractions of trabecular bone on the basis of magnetic resonance images acquired in vivo, Int. J. Imaging Syst. Technol. 10 (1999) 186–198. [37] B.R. Gomberg, P.K. Saha, F.W. Wehrli, Topology-based orientation analysis of trabecular bone networks, Med. Phys. 30 (2002) 158–168.
About the Author—PUNAM KUMAR SAHA received the B.S. and M.S. degrees in computer science and engineering from Jadavpur University, Calcutta, India, in 1987 and 1989, respectively, and the Ph.D. degree from the Indian Statistical Institute, Calcutta, India, in 1997. He joined the Indian Statistical Institute as a Faculty Member in 1993. In 1997, he joined the Department of Radiology, Medical Imaging Section, University of Pennsylvania, Philadelphia, as a Postdoctoral Fellow. He is currently an Assistant Professor at the same university. His present research interests include biomedical imaging problems and the application of their solutions, scale, segmentation, 6ltering, image registration, inhomogeneity correction, digital topology and their applications, and estimation of trabecular bone strength from MR images. He has published over 40 papers in international journals. Dr. Saha is an Associate Editor of the Pattern Recognition Journal. He is a Member of the International Association for Pattern Recognition and a Member of the Governing Body of the Indian Unit for Pattern Recognition and Arti6cial Intelligence. He received a Young Scientist Award from the Indian Science Congress Association in 1996. About the Author—FELIX W. WEHRLI received the Ph.D. degree in physical organic chemistry from the Swiss Federal Institute of Technology, Zurich, Switzerland, in 1970. He has spent his entire research career in the 6eld of magnetic resonance in leading research and development positions at Varian Associates, Bruker Instruments, and General Electric Medical Systems. He is currently Professor of Radiologic Science, Biochemistry, and Biophysics and the Director of the Laboratory for Structural NMR Imaging in the Department of Radiology of the University of Pennsylvania, Philadelphia. Over the past 20 years, his research has been concerned with the development of magnetic resonance procedures for clinical diagnosis, with particular focus on the understanding of the relationship between structure and function of connective tissues. A major portion of his recent work has dealt with the development of image acquisition and processing tools for characterizing structure and strength of trabecular bone. He is the author or coauthor of over 120 scienti6c articles, three textbooks, and 30 book chapters and reviews. He has been the Editor-in-Chief of Magnetic Resonance in Medicine, the leading magnetic resonance research journal since 1991. Dr. Wehrli is a Fellow of the International Society for Magnetic Resonance in Medicine.