Determination of the trabecular bone direction from digitised radiographs

Medical Engineering & Physics 25 (2003) 719–729 www.elsevier.com/locate/medengphy

Determination of the trabecular bone direction from digitised radiographs H. De´fossez a, R.M. Hall a,b,∗, P.G. Walker a, B.M. Wroblewski c, P.D. Siney c, B. Purbach c b

a School of Mechanical Engineering, University of Leeds, Leeds LS2 9JT, UK Academic Department of Orthopaedic Surgery, St. James’s Hospital, Leeds LS9 7TF, UK c The John Charnley Research Institute, Wrightington Hospital, Wigan WN6 9EP, UK

Received 26 July 2002; received in revised form 2 June 2003; accepted 2 July 2003

Abstract There is increasing evidence for monitoring the bone trabecular structure to explain, in part, the mechanical properties of bone. Despite the emergence of Computed Tomography, a radiograph is the standard format as it is cheap and used for assessing implant performance. Furthermore, various image-processing techniques developed to assess the trabecular structure from radiographs have regained interest owing to improvements in imaging equipment. This study assessed the precision and accuracy of the Co-occurrence and Run-length matrix, Spatial-frequency and Minkowski-fractal techniques to infer the trabecular direction from radiographs. Ten clinical images of femoral neck regions were obtained from digitised pelvic radiographs and subsequently analysed. These data were also used to generate synthetic images where the trabecular thickness, separation and directions were controlled in order to calculate the accuracy of the techniques. Additionally, a Laplacian noise was added in order to infer the precision of the techniques. All methods assessed the trabecular direction with a high degree of accuracy in these synthetic images including a single direction and no noise. However, only the Spatial-frequency and Co-occurrence matrix methods performed well on the clinical and heavily corrupted synthetic images. This demonstrated the possibility of inferring a linear trabecular direction in clinical conditions.  2003 IPEM. Published by Elsevier Ltd. All rights reserved. Keywords: Trabecular bone structure; Texture analysis; Anisotropy; Matrix; Spatial-frequency; Fractal; Radiographs

1. Introduction There are important clinical requirements for monitoring both bone modelling and remodelling such as providing osteoporosis diagnosis, predicting bone fracture or controlling the loosening of implants. Aspects of bone modelling and remodelling may be described in term of bone mineral density and bone trabecular structure. Bone mineral density is usually assessed using Dual Energy X-ray Absorptiometry (DEXA) [1], and more rarely Radiographic Absorptiometry (RA) [2,3] or Computed Tomography (CT) equipment [1]. The structure is less commonly evaluated than the density due to the dif∗ Corresponding author. Address: School of Mechanical Engineering, University of Leeds, Leeds LS2 9JT, UK. Tel.: +44-113-243-1751; fax: +44-113-243-4611. E-mail address: [email protected] (R.M. Hall).

ficulty of monitoring precisely the bone architecture. However, it is acknowledged that there is increasing evidence presenting important relationships between trabecular structure and bone toughness [4,5] and bone fracture [6]; and improvements in medical imaging equipment suggest an ameliorated assessment of the bone trabecular structure. The different techniques available to describe the bone trabecular structure were reviewed, with a particular attention towards those methods making use of radiographs, which still present many advantages over other recently introduced image formats. Indeed, even with the emergence of innovative image formats such as CT [1] or Magnetic Resonance Imaging (MRI) [7] offering the relatively easy assembly of three-dimensional data, radiographic images using film or more recently fully digital images provide the most common mode of assessment in most orthopaedic units. X-ray equipment is cap-

1350-4533/$30.00  2003 IPEM. Published by Elsevier Ltd. All rights reserved. doi:10.1016/S1350-4533(03)00123-1

720

H. De´ fossez et al. / Medical Engineering & Physics 25 (2003) 719–729

able of an excellent spatial resolution of 50 µm when compared with the 30 µm spatial resolution of CT equipment and 300 µm spatial resolution of MRI equipment. Additionally, a significant part of the information contained in 3D images is also contained in the corresponding radiograph [8]. Furthermore, only radiographs can be used to image bone with metallic implants, since the metal blocks the X-rays impeding the image reconstruction on CT equipment, and interferes with the magnetic field in the MRI system. A number of image-processing techniques have been used to describe the trabecular structure. Three sets of approaches applicable to radiographic images were distinguished. The first set relies on the matrix analyses that can differentiate fractured bones, and are derived from the Mean Intercept Length Fabric Tensor method [8,9]. The Run-length matrix method enables the inference of the trabecular direction [10–12], while the Co-occurrence matrix provides texture descriptors as well as the trabecular direction [10–12]. The second set of approaches is based on Spatial-frequency analyses that have been correlated to bone mechanical properties. The Intensity Orientation Method allows the inference of the dominant trabecular direction [4,6,13], while other power-spectrum methods allow the inference of the roughness of the trabeculae [10,14]. The final set of methodologies is derived from fractal analyses where a specific feature of the image is found at different scales [15–20]. The Minkowski-fractal approach has the advantage of enabling the inference of the trabecular direction as well as providing a texture descriptor correlated to bone mechanical properties [16]. The aim of the present study was to assess the potential of the matrix, Spatial-frequency and fractal techniques to infer the trabecular direction from radiographic images. It was decided to assess the accuracy and the precision of these techniques towards the analysis of the trabecular structure found in human femoral necks. Ten images of femoral necks were derived from digitised clinical pelvic radiographs. The radiographs were obtained from patients at Wrightington hospital who were selected for a clinical trial involving a unilateral hip replacement. Only the nonoperated femur, that was diagnosed as healthy by the surgeons, was used to create the radiographic images. These data were employed to generate synthetic test images in which a perfectly linear trabecular texture was created. The orientation of the texture was controlled at a desired angle in order to allow the calculation of the accuracy of the techniques, which quantifies how close a measure is from the real value. Additionally, a random Laplacian white noise was introduced in the synthetic images at various controlled levels in order to mimic a degradation of the apparent trabecular structure [21]. This noise also served to calculate the precision of the techniques, which quantifies the repeatability of the measurement in the same image, at

the same noise level but for several random noise patterns. Eventually, the clinical images were analysed to statistically measure the validity of the image-processing methods in clinical conditions.

2. Materials and methods 2.1. Images Ten clinical pelvic radiographs were digitised at 512 dpi on a 8 bits grey scale using a ScanMaster DX scanner (Howtek, New Hampshire, USA). These clinical data were used to develop a series of synthetic test images. In each radiograph, a clear trabecular bone region of 1 cm2 was randomly selected in the femoral neck as it allowed the visual observation of the longitudinal trabeculae that follow a linear directional pattern in the femoral neck (Fig. 1a). This area was subsequently processed using a histogram equalisation in order to enhance the contrast of the image (Fig. 1b) [10,11]. A single trabecula was then selected at random, and its measured width and separation from the other trabeculae were employed as a model to create three synthetic test images including a trabecular thickness of 10 pixels and separation of 10 pixels. Image A included a perfect single and linear trabecular bundle orientated at 60° (Fig. 2a). Image B was created similar to A, for the exception of the trabecular structure being fragmented (Fig. 2b). Image C included two perfect trabecular bundles at 90° or 110° (Fig. 2c). This trabecular structure, representative of ‘perfect’ clinical conditions, allowed the calculation of the accuracy of the image-processing techniques since the trabecular direction was known. Subsequently, a random Laplacian white noise [21] was introduced into the synthetic images, which enabled the precision of the techniques to be assessed in addition to being able to elicit further information on the accuracy of the methodologies. Six different noise levels were introduced into each synthetic image. Fig. 3 illustrates this pattern on the test image A. These levels, which correspond to the Laplacian coefficients 5, 50, 150, 500, 1000 and 10,000, were arbitrarily chosen to make the visual appearance of the images ranging from unaffected to fully noisy. Eventually, the clinical images were analysed to demonstrate the possibility to infer the trabecular direction in clinical conditions. 2.2. Image-processing techniques 2.2.1. Run-length matrix The Run-length matrix of an image I of Np × Np pixels of Ng grey levels, which has been rotated at the angle q, is defined by the term L(i,j) representing the occurrence of j pixels in succession on the matrix rows having the grey intensity i [10]. The Run-length matrix of a

H. De´ fossez et al. / Medical Engineering & Physics 25 (2003) 719–729

721

Fig. 2. Synthetic images inspired from the radiographic image CL1: (a) image A: perfect includes a singular linear trabecular structure at 60°, a width of 10 pixels and a separation of 10 pixels, (b) image B: fragmented is similar to image A, expect for the addition of voids and trabecular extensions, and (c) image C: dual direction is similar to image A, but it includes two linear trabecular directions at 90° and 110°.

Fig. 1. (a) Pelvic radiograph which was digitised at 512 dpi and in which a bone portion was selected in the right femur (white box). (b) Close-up of the selected trabecular bone portion after histogramequalisation. In this clinical image CL1, the black dashed helps to visualise the trabecula that was selected at random to estimate the trabecular width, while the white dots help to visualise the trabecula that was selected at random to estimate the trabecular separation.

radiographic image was computed for all angles from 0° to 180° in 1° steps. It was not necessary to complete the computations from 180° to 360° due to the symmetry of the matrix system. The Long Run Emphasis (LRE) term was subsequently computed as

Fig. 3. Concise illustration, the Laplacian noise being simultaneously incorporated from left to right in image A at six increasing levels. In practice, the noise was separately introduced in all the synthetic images at each noise level.

冋冘 冘

(1)

lar direction, but also of multiple orientations when more than one is present.

The LRE term was shown to obtain its highest value when the matrix angle was taken following the trabecular orientation [10,11]. It was also shown that this parameter enables not only the inference of a unique trabecu-

2.2.2. Co-occurrence matrix The Co-occurrence matrix of an image I of Np × Np pixels of Ng grey levels, which has been rotated at the angle q, is defined by a step of Dc pixels not greater

Ng⫺1 Np

i ⫽ 0j ⫽ 2

册冒冋 冘 冘 Ng⫺1 Np

(L(i,j) ⫻ j2)

册

(L(i,j)) .

i ⫽ 0j ⫽ 2

H. De´ fossez et al. / Medical Engineering & Physics 25 (2003) 719–729

722

than Np, and the term C(i,j) being the occurrence of the pair of pixel Mi and Nj of grey intensity of respectively i and j being present in the matrix rows and separated by the distance Dc [10]. The normalised associated matrix is calculated using P(i,j) ⫽ C(i,j) / C

(2)

where C is the sum of all occurrences. The Co-occurrence matrix method requires a prior knowledge concerning the trabecular thickness and separation that is subsequently used to calculate an optimal value for the matrix step dimension Dc. This is important, since a correct value of this step is required for the correct inference of the trabecular direction. Given that the thickness and separation have to be known a priori, test images including such a singular linear trabecular structure were used. The matrix of this image is subsequently computed at the angle perpendicular to the trabeculae and for various values of the step ranging from 2 pixels to the image length given in pixels. The computed contrast value, which is representative of the difference of intensity between pixels, is calculated as:

冘冘

Ng⫺1Ng⫺1

Contrast ⫽

(P(i,j) ⫻ (i⫺j)2).

(3)

i ⫽ 0j ⫽ 0

The contrast was shown to obtain a maximum value for a particular step dimension which is used to define the optimal step. Once this optimal value has been obtained for an image whose trabecular thickness and separation are known, the trabecular direction of any other image including a trabecular structure of same dimensions is possible, whatever the trabecular direction may be. The finding of the direction is achieved by computing the cooccurrence matrix at all angles from 0° to 180° in 1° steps. The associated contrast obtains a minimum value for the matrix taken following the texture direction [10,11]. 2.2.3. Spatial-frequency This technique studies the spatial frequencies of the image by computing its Fast Fourier Transform (FFT) and subsequently analysing the associated power-spectrum. The centre of the power-spectrum groups together the slowly varying components representing the main features of the image, while the outer parts of the spectrum gather together the high frequency elements representing the fine details and noise being present in the image. The trabecular structure of the bone is organised to sustain the applied physiological stresses. The longitudinal trabeculae are mainly grouped together in a parallel manner, and this arrangement produces a linear pattern. Therefore their associated frequency distributions are represented by a wave of high amplitude which lies in the centre of the power-spectrum, and which is

aligned perpendicularly to the longitudinal trabecular direction. The transverse trabeculae that are perpendicular to the longitudinal ones, are significantly smaller, and abruptly break the continuity of the longitudinal trabeculae. This structural organisation leads to a higher spatial variance, which then becomes apparent in the middle of the power-spectrum in the form of a second wave of amplitude smaller than the first wave, and aligned perpendicularly to the transverse trabeculae. Overall, the directional variation in the frequency components of the trabecular texture leads to an ellipsoidal pattern, whose minor axis lies in the direction of the trabeculae. The Intensity Orientation Method has been shown to be a robust tool to assess the direction of this axis [4,6,13]. This method was implemented by evaluating the intensities of the pixels situated on straight rays defined from the centre of the power-spectrum up to an arbitrary chosen distance DFFT, in a radial manner from 0° to 180° in 1° steps. The ray of highest intensity represented the major axis of the ellipse, indicating the transverse direction of the trabeculae. Furthermore, the theory behind the Fourier transform enables the adoption of a correct value for DFFT, and explains at which minimum resolution the radiograph needs to be scanned at [22]. So, a pixel out of the middle of the power-spectrum by i pixels (i苸[0⫺n / 2] for an image of n × n pixels) is representative of the frequency f(Digitalisation—Rate) ⫻ i / n.

(4)

Also, in a radiograph digitised at 512 dpi, the highest inferable frequency, which is called the Nyquist limit, is found at i = n / 2, where f = (512) × ((n / 2) / n) = 256 dpi, corresponding to a resolution of ⬵100 µm. Since the width of the human trabeculae varies from 1100 to 100 µm, the Nyquist limit suggests digitising radiographs at a resolution as high as 512 dpi, and to monitor the rays for a length DFFT as long as half the image size. Several additional processes were used to enhance the standard Intensity Orientation Method. Prior to determining the FFT, the image was processed using a Hanning window, which eliminates the artefacts introduced in the power-spectrum [13]. These artefacts arise in finite digital systems where the FFT reproduces the infinite space by tilting the image. The Hanning windowing process smoothens the edges of the image to a null intensity, allows continuity and removes the detrimental artefacts. In addition, the values of the power-spectrum, which are defined by P(u,v) where u and v are the frequency variables, were used after the following log transformation PL(u,v) ⫽ log(1 ⫹ P(u,v)).

(5)

This image-processing ploy compensates for the difficulty in observing the high frequencies in images where the power-spectrum decreases rapidly. Finally, in order

H. De´ fossez et al. / Medical Engineering & Physics 25 (2003) 719–729

to avoid aliasing problems, a bilinear interpolation was used to define the rays [23]. 2.2.4. Minkowski-fractal Fractal methods consist of identifying an image feature reproduced at different scales [15]. Various features have been studied in the past, using a box element counting method [20], a power-spectrum analysis [18] or a ‘dilatation/erosion’ operational scheme [16]. The latter method was adopted since the Minkowski dimension derived from these morphologic operations has been correlated to bone strength [16]. The Minkowski dimension of an image I of Np × Np pixels of Ng grey levels, which has been rotated at the angle q, is defined as log(PS(e) / e3) MD[I,q] ⫽ lim , log(1 / e) e→0

(6)

where S is a structuring element at scale e and PS(e) is the pixel intensity between the two processed versions of the image [16]. The term PS(e) is obtained as

冘冘 Np

PS(e) ⫽

Np

[(I丣eS)⫺(I丢eS)],

(7)

i ⫽ 0j ⫽ 0

where (I丢eS) and (I丣eS) are the eroded and dilated version of the image, using the structuring element S at the scale e. Therefore, if a nonsymmetric structural element such as a line array is adopted, the Minkowski dimension is function of the image angle. For the analysis of the radiographic images, a line array, ranging in scales from 2 to 10 pixels, was used as a structural element and employed at all angles from 0° to 180° in 1° steps. The Minkowski dimension was expected to obtain a maximum when the angle followed the trabecular direction, permitting the inference of the texture orientation. Additionally an in-house modification of this method was introduced by using only the intensity between the two processed versions of the image. 2.3. Statistical performance of the techniques The performance of a technique is usually quantified by computing its accuracy and precision [24]. In the present case, the accuracy quantifies the ability of an image-processing technique to measure the trabecular orientation as close as its real value. The synthetic images were used to infer the accuracy since the trabecular orientation was defined at a known angle. The clinical images, in which the trabecular direction was only visually estimated, were therefore only used to demonstrate the validity of the methods in clinical conditions. The accuracy error was calculated in the synthetic images as the absolute value of the difference between the ‘true’ trabecular direction and the measured trabecular direction. In the case of repeatable measurements of a same image, it was calculated as the absolute

723

value of the difference between the true direction and the mean of the repeated measurements. The precision of a technique is the quantification of the variation of repeated measurements of a same source [24]. This is calculated as the standard deviation defined as SDj ⫽

¯) , 冪 冘 冉(xn⫺x ⫺1 冊 nj

ij

i⫽1

j

2

(8)

j

for a source j (an image), where nj is the number of repeatable measurements (the angle of the trabecular direction), xij the result of the ith measurement, and x¯ j the mean of all xij measurements (the angle given in degrees). The general definition tacitly implies that the source to be measured, or the measuring equipment, is subjected to variations making the possibility of the repeatable measurements to be different from each other. In the present case, an image (the source) which is analysed several times by the same image-processing algorithm (the measuring equipment) will lead to the same repeated measurements since no variation affects neither the source nor the equipment; the precision will therefore appear to be maximal (SD of 0°). Therefore, it was decided to introduce some variations in the source data by adding some random noise in the image, as previously described. The precision of a technique was computed at each of the six noise levels for each synthetic image. For a desired noise level and image, the random noise was introduced, the image analysed and the process repeated 25 times. Note that the precision error was only introduced in this study via the random incorporation of noise. In the analyses of a radiographic image, other factors will contribute to the precision error. These can include film type, radiographic setting, film development process, operator dependence when selecting the ROI or the adopted thresholding algorithms. The performance of a technique in clinical conditions is of critical importance. This provides the possibility to define if a new method is appropriate to supersede an old one, to quantify how better a technique performs over another one in clinical conditions. In these conditions, the true measurements remain unknown and therefore the typical use of correlation is insufficient. Instead, the standard deviation of the mean of the difference between two techniques should be used [25]. For each image-processing technique, the various clinical images were analysed and the differences between the visually estimated trabecular orientation values and the measured values were employed to calculate the standard deviation. This indicates the probability to measure a value contained between ± 1.96 × SD with 95% confidence.

724

H. De´ fossez et al. / Medical Engineering & Physics 25 (2003) 719–729

2.4. Software and computing resources Custom written software, developed in the IDL programming environment (Research Systems International, Crowthorne, UK), was used to implement the various image-processing techniques. All processing was undertaken on a PC (AMD Athlon clocked at 1.2 GHz with 256 MB ram). 3. Results For clarity, it was decided to present the results in two main sections. The first one describes the typical analyses of the synthetic and clinical images by each imageprocessing method. The second section introduces the complete analyses of all images which served to statistically calculate the performance of the techniques. 3.1. Typical analyses of synthetic and clinical images 3.1.1. Run-length matrix analyses The typical results of clinical and synthetic images are displayed in Fig. 4, where the LRE value is given as a

Fig. 4. Typical Run-length matrix analyses of synthetic and clinical images: (a) synthetic images A, B and C, (b) clinical image CL1.

function of the angle at which the image has been analysed. In the case of the synthetic image A, which includes a single perfect linear trabecular direction, the LRE term obtained a sole peak allowing the accurate inference of the trabecular direction (Fig. 4a). In the image B, which includes a fragmented version of the trabecular structure of image A, the LRE term also obtained a sole peak allowing the accurate inference of the trabecular direction (Fig. 4a). In the image C, two peaks of similar intensity were distinguishable and accurately indicative of both trabecular directions (Fig. 4a). In the clinical image CL1, the LRE value still obtained a peak allowing the detection of the texture orientation, even though this maximum value was not as distinct as in the synthetic test images which used an ideal trabecular texture (Fig. 4b). 3.1.2. Co-occurrence matrix analyses It was possible to infer the optimal step dimension used in the Co-occurrence matrix of the image where the ideal trabeculae with the expected width and separation were introduced. The results presented in the Fig. 5a showed the typical pattern of the contrast values as a function of the step dimension Dc. So, for the synthetic image A, which includes a perfect unique linear trabecular direction, the contrast increases as the step dimension increases, until reaching a maximal value defining the optimal step dimension. Subsequently, the contrast values drop until reaching a minimum, and this replicates to form a cyclic pattern. This phenomenon was possible due to the perfect symmetry of the trabeculae of the test image where the distance required to ‘jump’ one or more trabeculae is a multiple of the fundamental step dimension. Furthermore, the contrast was null when the matrix was defined using the fundamental step dimension. This phenomenon is due to the step value which leads to the distance between the points forming the pairs of the matrix to be such that the two points lie inside a same trabeculae (or separation locations) having a homogenous intensity. Note that the optimal step value needed to be calculated for a matrix taken perpendicular to the trabecular direction, in order to normalise this scheme. Otherwise, if the image was taken as such, another optimal value was found leading to an incorrect inference of trabecular direction. The definition of the matrix being perpendicular to the trabeculae was possible since a perfect image was generated. In the case of an image where the trabecular orientation is unknown, two solutions are possible. The first one consists in using an already known optimal step dimension corresponding to the expected texture dimensions; such knowledge would be gained from experience. The other solution consists of inferring the optimal step dimension using the textural orientation of the image that

H. De´ fossez et al. / Medical Engineering & Physics 25 (2003) 719–729

725

Fig. 5. Typical Co-occurrence matrix analyses of clinical and synthetic images: the optimal step was determined to be 9 pixels for the synthetic image A (a) and 15 pixels for the clinical image CL1 (b); the trabecular orientation was determined to be 60° for image A (c) and 64° for image CL1 (d).

has previously be determined by another method, such as the Spatial-frequency method. This latter scheme was adopted for analysing the step dimension of the clinical image CL1. The contrast in function of the step dimension followed a pattern similar to the one of the ideal synthetic image (Fig. 5b). However, the cyclic pattern was not present due to the imperfect symmetry of the image trabecular pattern. Furthermore, the contrast was not null for the fundamental step since the trabeculae were not of a homogeneous intensity. In this case the first peak was used for defining the optimal step value. Finally, the analysis of the trabecular direction was straightforwardly performed using the optimal step for the synthetic and clinical images (Fig. 5c and d). 3.1.3. Spatial-frequency analyses The Spatial-frequency analysis enabled the inference of the trabecular direction since the intensity of the rays obtains a distinct maximum at the angle corresponding to the direction perpendicular to the trabecular direction. Fig. 6 illustrates this phenomenon for the perfect linear trabecular structure present in the image A, and the fragmented linear structure of the image B. In the clinical image CL1, a similar pattern was also observed (Fig. 6). The peak was still distinguishable even though the signal became slightly noisy due to the imperfect trabeculae. 3.1.4. Minkowski-fractal analyses For the synthetic image A, which includes a single perfect trabecular direction, the computed Minkowski dimension obtained its maximum at the angle following

Fig. 6. Typical Spatial-frequency analyses of the synthetic images A and B and the clinical image CL1.

726

H. De´ fossez et al. / Medical Engineering & Physics 25 (2003) 719–729

the trabecular direction (Fig. 7a). Furthermore, it was observed that the sum of the processed surfaces at the different scales appeared sufficient to infer the trabecular direction as it obtained a minimum when the image angle was taken following the trabecular direction (Fig. 7b). The same patterns were present for the synthetic image B that includes a single but fragmented trabecular direction. The analysis of the clinical image CL1 by the Minkowski and modified fractal analysis showed similar patterns for the inference of the trabecular direction; the signal was slightly noisy due to the imperfect trabecular texture (Fig. 7a and b).

3.2. Analyses of all images providing statistical performance of the techniques in artificial and clinical conditions 3.2.1. Analyses of the uncorrupted synthetic images The synthetic image A, without the addition of noise, was perfectly analysed by all the image-processing techniques. The same excellent accuracy was obtained for the image B, even though its texture was fragmented. However, only the Run-length matrix method was able to accurately determine the two trabecular directions that are present in the image C. 3.2.2. Analyses of the synthetic images corrupted with noise The accuracy and precision of the techniques when employed in the synthetic images which include a single linear trabecular structure are presented in Fig. 8a for image A and Fig. 8b for image B. – The Run-length matrix method started to fail to infer the direction even with the addition of the smallest noise (accuracy error of 12.20°, 9.68° and standard deviation of 15.02°, 15.23°, respectively, for images A and B). – The Co-occurrence matrix method was the most accurate and precise method, up to a noise of 1000 (accuracy error of 0.68°, 0.80° and standard deviation of 5.64°, 5.18°, respectively, for images A and B). – The Spatial-frequency method was excellent at inferring the trabecular orientation in the noisy synthetic images up to a noise of level 500 (accuracy error of 0.04°, 0.08° and standard deviation of 0.20°, 0.28°, respectively, for images A and B). – The Minkowski and modified fractal approaches allowed the finding of the trabecular direction up to a noise of 150 (accuracy error better than 1° and standard deviation better than 2° for the images A and B). The analyses of the image C, which includes two trabecular bundles, showed that no method was able to infer the two trabecular directions when noise was introduced. Instead, only one direction was picked up, this being a value between the two ‘true’ directions. This failure occurred even for the smallest noise level, while the Range-length method successfully inferred the two trabecular directions. An explanation for this failure is presented in Discussion.

Fig. 7. Typical fractal analyses of synthetic and clinical images: (a) Minkowski analyses of the synthetic images A and B and the clinical image CL1, (b) modified fractal analyses of the images A, B and CL1.

3.2.3. Analyses of the clinical images The analyses of the clinical images including portions of the femoral neck are presented in Fig. 9. The accuracy error, here being the mean difference, was computed using the trabecular direction that was visually estimated. This estimation was conducted by selecting five trabeculae at random and visually estimating the five

H. De´ fossez et al. / Medical Engineering & Physics 25 (2003) 719–729

Fig. 8.

727

Results of the analyses of the synthetic images without noise and once corrupted at six various noise levels: (a) image A, (b) image B.

Fig. 9. Results of the analyses of the clinical images.

associated linear directions. The mean of these five directions defined the trabecular direction in the image. The results indicated that: – The Run-length matrix method appeared fairly accurate at evaluating the trabecular direction. The mean difference was found to be 22.1° and the standard deviation 14.47°. The method performs better in clinical images than expected by the theoretical analysis of the noisy synthetic images. In the clinical images, there is however an imperfect trabecular structure, and there is still a significant amount of trabeculae rep-

resented by pixels of same intensity, while this is not the case with the noisy synthetic images. – The Co-occurrence matrix method was very efficient with a mean difference of 4.0° and a standard deviation of 3.0°. This confirmed the assumption given by the theoretical analysis forecasting that this method should be valid in clinical conditions. However, the analysis is quite difficult since the definition of the optimal step dimension is needed to search for the trabecular direction. – The Spatial-frequency method appears to be also very efficient with a mean difference of 7.2° and a standard deviation of 5.7°. This confirmed the assumption given by the theoretical analysis forecasting that this method should be valid in clinical conditions. – The Minkowski-fractal appears to be quite efficient with a mean difference of 8.7° and a standard deviation of 5.3°. The modified fractal analysis is slightly less efficient (mean difference of 10.6° and a standard deviation of 10.1°). The theoretical analysis on the synthetic images allowed the correct forecasting of the efficiency of this method.

728

H. De´ fossez et al. / Medical Engineering & Physics 25 (2003) 719–729

4. Discussion This study provided new insights into the relative performance of three different broad strategies used to determine trabecular structure from radiographic images. It was demonstrated that whilst each of the methods work well when applied to idealised synthetic images, they had a different performance with the addition of varying levels of noise to these images and in clinical conditions. – Concerning the Run-length matrix method, previous studies had shown a consistent correlation between this analysis and the bone quality in the calcaneus [12], and that this method was adequate at quantifying the trabecular directions in the femoral neck [10]. The present study showed the poor accuracy and precision of the technique in noisy synthetic images. However, it also showed that this method, although being the worst of the various adopted approaches, performs better than expected in clinical conditions studying the femoral neck. This is due to the Laplacian noise introduced, which considerably reduces the number of successive pixels of same intensity, on which the method relies. A strong reduction in these pixels is achieved even at the smallest noise level due to its artificial nature. However, although the trabecular texture found in clinical conditions is imperfect, there is still a significant number of these relevant pixels to give the results observed. – For the Co-occurrence matrix, Ferrand showed that the method could be used to infer the trabecular direction, but the information concerning the precision and accuracy were not provided [10]. Kaufman demonstrated that textural analysis of radiographs could be employed to quantify architectural changes under clinical circumstances in disuse osteopenia [12]. The present study presented the very high accuracy and precision of the technique applied not only to noisy synthetic images but also to femoral neck clinical images. Furthermore, the difficulty and time required to employ this method was highlighted. – The Spatial-frequency method has been recognised as giving good results [4,6,13]. Only the current study provided information concerning the accuracy and precision of this method. Furthermore, this method was relatively easy to implement, and the most rapid to analyse the images on our imaging equipment; this took less than 30 s compared to circa 10 min for the matrix analyses. – Various fractal methods were studied [15–20]. The correlation between the fractal dimension and the bone quality was demonstrated; however, there was no mention of the trabecular direction, except for the study of Jiang [16]. His research showed that the Minkowski-fractal method provides a trabecular direction

in addition to the Minkowski dimension [16]. These values, added to mineral density, were highly relevant to improve the explanation of the bone strength in the femoral neck. The study also highlighted a significant discrepancy between the measured and estimated trabecular direction that was measured in radiographs of bone cubes cut from femoral necks. Interestingly, the present analysis, which studied the femoral neck obtained from pelvic radiographs, showed only a small discrepancy between the measured and estimated trabecular direction. The present radiographs were digitised at a higher resolution that the one used by Jiang, and this might be an explanation for the difference between the results.

5. Conclusion This study not only reiterated the possibility to infer the trabecular direction in the femoral neck in clinical conditions using pelvic radiographs, but also extended the correct body of knowledge in two ways. The first way consisted in comparing the various image-processing techniques to infer the trabecular direction by calculating their accuracy and precision in synthetic images as well as their statistical performance in clinical conditions. The second way was to analyse pelvic radiographs which were digitised at 512 dpi; this resolution being high enough to assure that even the thinnest trabeculae were visible in the image. The best method to employ was the Spatial-frequency method. The Co-occurrence matrix approach performed only marginally better, but it was difficult and time consuming to implement, while the use of the Spatial-frequency method was straightforward. There was also a necessity to scan the radiographs at a high resolution of 512 dpi if one wants to assure the monitoring of the smallest trabeculae in order to conduct robust trabecular analyses. The theoretical analyses performed on synthetic images was very useful to predict the performance of the methods in clinical conditions. However, the study of the clinical radiographs remained indispensable to ensure the validity of the predictions since the theoretical situation never fully encompasses the clinical conditions. Finally, the image-processing methods presented in this study only allow the inference of a linear trabecular structure. The trabecular pattern can however appear linear, fan shaped, curved or simply random at various bone locations. Nonetheless, this trabecular arrangement is always optimised to make the bone suitable to sustain the applied stresses. In the femur, the neck is mostly studied, since this location include longitudinal trabecular that follow a clear linear direction. The clarity of the trabeculae present in this linear pattern region can be visually appreciated, and is an indicator of bone quality.

H. De´ fossez et al. / Medical Engineering & Physics 25 (2003) 719–729

Bone mineral density studies focusing on the humeral neck are indeed used to assess risk fracture in the hip. A more complex trabecular structure would be analysed using an extension of the presented algorithms as well as other methods that need to be developed. Acknowledgements This work was partly supported by the John Monk Hip Research Fund, the Peter Kershaw and John Charnley Trusts’. References [1] Ebbesen E, Thomsen J, Beck-Nielsen H, Nepper-Rasmussen H, Mosekilde L. Lumbar vertebral body compressive strength evaluated by Dual energy X-ray Absorptiometry, quantitative computed tomography, and ashing. Bone 1997;25(6):713–24. [2] Strid K, Kalebo P. Bone mass determination from microradiographs by computer-assissted videodensitometry. I Methodology. Acta Radiol 1998;29(Fasc. 4):611–7. [3] Strid K, Kalebo P. Bone mass determination from microradiographs by computer-assissted videodensitometry. II Aluminium as a reference substance. Acta Radiol 1998;29(Fasc. 4):465–72. [4] Libergall M, Simkin A, Neeman V, Mosheiff R, Jabschinsky C, Leichter I. Bone radiograph derived variables—can they replace DEXA? In: 46th Annual Meeting. Orthopaedic Research Society; 2000. [5] Wigderpwitz CA, Paterson CR, Dashti H, McGurty D, Rowley DI. Prediction of bone strength from cancellous structure of the disal radius: can we improve on DEXA? Osteoporosis Int 2000;11:840–6. [6] Hagiwara H, Inoue N, Matsuzaki H, Cohen D, Kostuik J, Chao E. Relationship between structural anisotropy and bone anisotropy and bone mineral density in the human vertebral body. In: 46th Annual Meeting. Orthopaedic Research Society; 2000. [7] Majumdar S, Kothari M, Augat P, Newitt DC, Link TM, Lin JC et al. High-resolution magnetic resonance imaging: three-dimensional trabecular bone architecture and biomechanical properties. Bone 1998;22(5):445–54. [8] Luo G, Kinney J, Kaufman J, Haupt D, Chiabrera A, Siffert R. Relationship between plain radiographic patterns and threedimensional trabecular architecture in the human calcaneus. Osteoporosis Int 1999;(9):339.

729

[9] Odgaard A. Three-dimensional methods for quantification of cancellous bone architecture. Bone 1997;20(4):315–28. [10] Ferrand D. Traitement bas niveau sur les images medicales: application a la detection des trabecules dans une radiographie de la hanche. France: Laboratoire d’Etudes et de Recherches sur les Signaux et Images, Faculte des Sciences et Techniques de SaintJerome, Universite d’Aix-Marseile III, 1999. [11] Dubus JP. Mesure par analyse d’image: analyse statistique et texturelle. Techniques de l’Ingenieur 1998;R630. [12] Kaufman J, Mont A, Ohley N, Lundhal T, Soifer T, Siffert S. Texture analysis of radiographic trabecular patterns in disuse osteopenia. Trans Orthopaed Res Soc 1987;12:265. [13] Hamilton R. Texture mapping of trabecular structure using digital image processing. British Research Orthopaedic Society, 1999. [14] Caligiuri P, Giger M, Favus M, Jia H, Doi K, Dixon L. Computerized radiographic analysis of osteoporosis: preliminary evaluation. Radiology 1993;186(2):471–4. [15] Berry J, Webber R, Horton R, Santago P, Pope T. Fractal dimension as a measure of change in trabecular bone structure. BED Adv Bioengng ASME 1993;26. [16] Jiang C, Giger M, Chinander M, Martell J, Kwak S, Favus M. Characterization of bone quality using computer extracted radiographic features. Med Phys 1999;26(6):872–9. [17] Webber R, Underhill T, Horton R, Dixon R, Pope T. Predicting osseous changes in ankle fractures. IEEE Engng Medicine Biol 1993;103–10. [18] Ruttimann U, Webber R, Hazelrig J. Fractal dimensions from radiographs of peridental alveolar bone. Oral Surg Oral Medicine Oral Pathol 1992. [19] Benhamou C, Lespessailles E, Jacquet G, Harba R, Jennane R, Loussot T. Fractal organization of trabecular bone images on calcaneus radiographs. J Bone Mineral Res 1994;9(2):1909–18. [20] Majumdar S, Weinstein R, Prasad R. Application of fractal geometry techniques to the study of trabecular bone. Med Phys 1993;20(6):1611–9. [21] Pitas I. Digital image processing algorithms. Prentice-Hall, 1993 p. 445. [22] Gonzalez R, Wintz P. Digital image processing, 2nd ed. Reading, MA: Addison-Wesley, 1987. [23] Hamilton R. Design and testing of techniques to measure anisotropy in trabecular bone structure. 4H MSc Project Report, Trevelyan College, University of Durham, 1999. [24] Gluer C, Blake G, Lu Y, Blunt B, Jergas M, Genant H. Accurate assessment of precision errors: how to measure the reproducibility of bone densitometry techniques. Osteoporosis Int 1995;5:262–70. [25] Bland J, Altman D. Statistical methods for assessing agreement between two methods of clinical agreement. Lancet 1986;8:307–11.