Long range node-strut analysis of trabecular bone microarchitecture
LLong Range Node-Strut Analysis of Trabecular Bone Micro-architecture
Tanya Schmah
Department of Computer Science, University of Toronto,Toronto, Ontario M5S 3G4, Canada
Norbert Marwan
Potsdam Institute for Climate Impact Research (PIK), D-14412 Potsdam, Germany, andInterdisciplinary Center for Dynamics of Complex Systems, University of Potsdam, D-14415 Potsdam,Germany
Jesper Skovhus Thomsen
Department of Connective Tissue Biology, Institute of Anatomy, University of Aarhus, DK-8000, ˚Arhus C,Denmark
Peter Saparin
Perceptive Informatics, a PAREXEL technology company, Am Bahnhof Westend 15, D-14059 Berlin,Germany
Purpose
We present a new morphometric measure of trabecular bone micro-architecture, called mean nodestrength (NdStr), which is part of a newly-developed approach called long range node-strut analysis . Ourgeneral aim is to describe and quantify the apparent “lattice-like” micro-architecture of the trabecular bonenetwork.
Methods
Similar in some ways to the topological node-strut analysis introduced by Garrahan et al., ourmethod is distinguished by an emphasis on long-range trabecular connectivity. Thus, while the topologicalclassification of a pixel (after skeletonisation) as a node, strut, or terminus, can be determined from the 3 × node strength . The node strength is averaged overa region of interest to produce the mean node strength (NdStr) of the region. Results
We have applied our long range node-strut analysis to a set of 26 high-resolution peripheralquantitative computed tomography (pQCT) axial images of human proximal tibiae acquired 17 mm belowthe tibial plateau. We found that NdStr has a strong positive correlation with volumetric trabecular bonemineral density (BMD). After an exponential transformation, we obtain a Pearson’s correlation coefficient of r = 0 .
97. Qualitative comparison of images with similar BMD but with very different NdStr values suggeststhat the latter measure has successfully quantified the prevalence of the “lattice-like” micro-architectureapparent in the image.Moreover, we found a strong correlation ( r = 0 . Conclusions
The newly introduced morphometric measure allows a quantitative assessment of the long-range connectivity of trabecular bone. One advantage of this method is that it is based on pQCT images thatcan be obtained noninvasively from patients, i.e. without having to obtain a bone biopsy from the patient.PACS numbers: 87.59.bd, 05.45.-a, 07.05.PjKeywords: trabecular bone, osteoporosis, structure analysis, histomorphometry, pQCT
I. INTRODUCTION
Several studies have shown that measures of trabec-ular bone micro-architecture and bone strength arecorrelated.
Together with loss of bone mass, changesin the trabecular bone micro-architecture occur duringageing, during development of osteopenia and osteoporo-sis as well as in connection with immobilisation or spaceflight, and can lead to an increased risk of bone fracture.The vertebral bodies and the epiphyses and metaphysesof the long bones consist mainly of trabecular bone sur-rounded by a thin cortical shell.
A dramatic changein the state of the trabecular bone leads to an increased fracture risk. Bone mineral density (BMD) is the most commonly usedpredictor of bone strength and fracture risk, and alsothe most commonly used general descriptor of the stateof the bone. Non-linear relationships have been estab-lished between volumetric BMD and compressive bonestrength and elastic modulus.
However, for trabecu-lar bone, it has been established that a part of the vari-ation in the strength of the bone cannot be explained byBMD alone, but is instead due to the micro-architectureof the trabecular network. For example, a relationshiphas been established between the mechanical propertiesof the bone and the shape, orientation, bone trabecular a r X i v : . [ phy s i c s . m e d - ph ] A ug volume fraction, and thickness of the trabeculae. A series of new methodologies based on techniques fromnonlinear data analysis has also been introduced in orderto study the relationship between the complexity of thetrabecular bone network and bone strength.
Saparinet al. established such a relationship by use of structuralmeasures of complexity based on symbolic encoding.
Furthermore, several studies using numerical modellingof the trabecular bone and finite element analysis havealso confirmed the importance of the trabecular bonemicro-architecture to bone strength.
In the present study we propose a new method for anal-ysis of trabecular bone micro-architecture from high-resolution quantitative computed tomography (QCT)images, as at present it is still not possible to performmicro CT ( µ CT) on humans in vivo and in situ. In con-trast to classical histomorphometry, CT images can beobtained in a nondestructive and noninvasive manner,which is preferable in a clinical setting. Our methoduses a new approach called long range node-strut analy-sis that quantifies the apparent nodes and struts in thetrabecular network. In contrast to the node-strut anal-ysis of Garrahan et al., our method emphasises longrange connectivity of the trabecular network, over a dis-tance controlled by a parameter in the algorithm. Theanalysis has the dual aims of describing the shape of thetrabecular network and predicting bone strength.The trabecular bone compartment consists basically oftwo components: bone and marrow. Due to the limitedresolution of present day CT and peripheral quantita-tive CT (pQCT) scanners it is not possible to completelyresolve the trabecular bone micro-architecture. This re-sults in variations of the CT values of the trabecular vox-els, even if the intrinsic bone density is constant through-out the trabecular network. Consequently, our methodtakes this apparent variation in bone density into consid-eration rather than segment or binarise the image, whichwould imply a loss of this additional information. Weassume that higher CT values of a given trabecula will,all else being equal, result in a higher compressive bonestrength. We note that our method does not requireskeletonisation of the image.We apply our method to 2-dimensional pQCT imagesof human proximal tibiae and quantify the trabecularbone micro-architecture at different levels of bone in-tegrity ranging from normal healthy bone to osteoporoticbone (as assessed by their BMD). We compare our re-sults with the node/terminus ratio (Nd/Tm), com-puted using traditional histomorphometry performed onbone biopsies obtained at the same location as the pQCTscans. II. MATERIALS
The study population comprised 18 women aged 75–98years and 8 men aged 57–88 years. At autopsy, the tibialbone specimens were placed in formalin for fixation. Foreach specimen a pQCT image and a bone biopsy wereobtained from the same location.For each proximal tibia, an axial QCT slice was acquired17 mm below the tibial plateau with a Stratec XCT-2000pQCT scanner (Stratec GmbH, Pforzheim, Germany),with an in-plane pixel size of 200 µ m × µ m and aslice thickness of 1 mm. In some cases, the scans wereperformed after the biopsies were taken. Therefore, theholes left from the bone biopsy appear in some of thepQCT images. A standardised image pre-processing pro-cedure was applied to exclude the cortical shell from theanalysis. One of the resulting images is shown inFig. 1.
FIG. 1. An axial pQCT slice of human proximal tibia ac-quired 17 mm below the tibial plateau. The cortical shellhas been removed from the image. The trabecular BMD ofthis sample is 107 mg/cm . The hole visible on the left of theimage is the result of a cylindrical biopsy. Cylindrical bone samples with a diameter of 7 mmwere obtained 17 mm distal from the centre of the me-dial facet of the superior articular surface by drillingwith a compressed-air-driven drill with a diamond-tippedtrephine at either the right or the left proximal tibia.These bone biopsies were embedded undecalcified inmethyl methacrylate, cut in 10- µ m-thick sections on aJung model K microtome (R. Jung GmbH, Heidelberg,Germany), and stained with aniline blue (modified Mas-son trichrome). The mounted sections were placed in anflat-bed image scanner and 2540 dpi digital 1 bit imagesof the sections were obtained as previously described indetail. The resulting pixel size is 10 µ m × µ m.The trabecular BMD of each pQCT slice was calculatedusing a linear relationship derived on the basis of experi-mental calibration with the European Forearm Phantom,as described by Saparin et al. The trabecular BMDs ofthe slices in this study range from 30 to 150 mg/cm ,with a median of 97 mg/cm . III. ANALYTICAL METHOD
Most of this section describes the new image analysismethod, long range node-strut analysis , which includesthe new measure, mean node strength (NdStr). Thismethod was applied to the pQCT sections described inSection II, after removal of the cortical shell. We alsocomputed a standard measures from the same regions ofinterest of the same pQCT images: trabecular volumet-ric bone mineral density (BMD), calculated as describedby Saparin et al. For further comparison, topological 2-dimensional node-strut analysis was performed on the histological sec-tions described in Section II, using a custom-made com-puter program. The trabecular bone profile was iter-atively eroded until it was only one pixel thick using aHilditch skeletonisation procedure, and nodes and ter-mini were automatically detected by inspecting the local3 × × Consequently, only the node-to-terminus ratio (Nd/Tm) from this analysis of the his-tological sections is used in the present study.
A. Basic definitions
We start with the algorithm to find strands in each image.A strand is a connected trabecular path, i.e., a chain ofone or more struts, with the whole path going in approxi-mately the same direction. In our algorithm, we use eightdirections, labelled by the points of the compass: north(N), northeast (NE), east (E), southeast (SE), south (S),southwest (SW), west (W), or northwest (NW). “North”is the anterior direction, which corresponds to “up” inthe images shown in this article. A node is a pixel thatis joined to strands in at least three of these eight di-rections, any two of which are at least 90 degrees apart(e.g. N, E, and SW; but not N, NE, and E). The nodestrength of a pixel is 0 if it is not a node at all, but oth-erwise depends on the lengths of the strands that meetat the node and the pQCT values of the pixels in thesestrands.
B. CT values of bone
These basic definitions could be implemented in a varietyof ways, depending for example on the definitions of “con-nected” and “same direction”, and on how the pQCTvalues of the pixels are used. In our algorithm, the firststep is to remove marrow from further consideration, byapplying a bone threshold filter. The thresholded pQCTvalues will be referred to in this section as
CT values ofbone , denoted b . Thus if a is the pQCT value of a pixel,then b = (cid:40) a − a threshold if a > a threshold . (1)We choose a threshold = 275, corresponding to a BMDvalue of 24 mg/cm . This is the soft tissue threshold usedin symbol encoding for complexity measures by Saparinet al. The threshold was determined as the pQCT valueof the densest non-bone tissue tested, plus the standarddeviation of the noise (a calibration protocol for findingthe threshold based on a phantom will be developed inthe future).
C. Strands
In order to explain the central algorithm, we begin byconsidering a fictitious image consisting of just one row.The CT values of the pixels in the row form a sequence b , b , b , b , · · · . We will define a corresponding sequence of strandstrengths , s , s , s , s , · · · , where each s j describes the pattern of densities to theleft of pixel j in the row. We begin by setting s = 0,and then we recursively define s j = ( b j − + T s j − ) (cid:18) min { b j , b j − } b j − (cid:19) , (2)where T is a transmission constant between 0 and 1.Note that if all of the CT values have a common value b then s j = ( b + T s j − ) = ( b + T ( b + T s j − ))= b (1 + T ) + T s j − = · · · = b (cid:0) T + · · · + T j (cid:1) , since s = 0. As j increases, s j approaches an upperbound of b − T . Thus, for the simple case of a uni-formly dense row, we have defined the strand strengthto be 11 − T times the CT values of bone. This constant11 − T can be interpreted as a characteristic length ofthe method, which we can vary depending on the lengthof strand that we believe to be most important to thestrength of the bone. For the present study, we have cho-sen T = 0 .
95, corresponding to a characteristic length of20 pixels, i.e. 4 mm. For the general case of a variableCT values row, the presence of the transmission constant T in the recursive formula ensures that the CT value ofpixels closest to the j th pixel have the greatest effect on s j . Finally, the factor min { b j , b j − } b j − ensures that s j de-pends strongly on b j , so that any weak link in a chain ofhigh-density pixels lowers the strand strength dramati-cally.The difficulty in generalising this definition to a 2-dimensional image is of course that there are infinitelymany directions but only four that are parallel with thepixel grid. We have resolved this in a practical but adhoc fashion. As already stated, we consider eight direc-tions. Consider first the definition of strand strength inthe leftwards (W) direction. ¿From pixel ( i, j ) (at row i and column j ) we consider that there are five length-3paths leading approximately leftwards:(i) ( i + 1 , j − ← ( i + 1 , j − ← ( i, j ) , (ii) ( i + 1 , j − ← ( i, j − ← ( i, j ) , (iii) ( i, j − ← ( i, j − ← ( i, j ) , (iv) ( i − , j − ← ( i, j − ← ( i, j ) , (v) ( i − , j − ← ( i − , j − ← ( i, j ) , as illustrated in Fig. 2. j –2 j –1 j j –2 j –1 j j –2 j –1 ji +1 i –1 ii +1 i –1 i (i) (ii) (iii)(iv) (v) FIG. 2. The five possible length-3 paths in a leftwards(“west”) direction from pixel ( i, j ). The start pixel ( i, j ) ismarked in black.
D. Strand strength
For each of these paths, we define b j , b j − and b j − as theCT values of the pixels in the path in columns j, j −
1, and j −
2. We then apply the recursive formula s j = (cid:18) b j − + T ( b j − + T s j − ) (cid:18) min { b j − , b j − } b j − (cid:19)(cid:19) (3) × (cid:18) min { b j , b j − } b j − (cid:19) , which is the formula in Eq. (2), applied twice. (We mustnow set s = s = 0.) We now have five s j values,one for each of the five paths, i.e., s ( i ) j , s ( ii ) j , . . . , s ( v ) j .We multiply all but the third (Fig. 2(iii)) of thesefive values by a bending coefficient γ between 0 and1, to penalise strands that bend away from the hor-izontal direction. Then the maximum of these fivenumbers is our leftwards (W) strand strength S Wij =max { γs ( i ) j , γs ( ii ) j , s ( iii ) j , γs ( iv ) j , γs ( v ) j } .Strand strengths in the other seven directions S NWij , S
Nij , S
NEij , . . . , S
SWij are determined in a similarmanner. For diagonal directions, the geometry of thefive length-3 paths used in the computation is slightlydifferent, so a different bending coefficient γ is used. Wechose γ = 0 . γ = 0 . E. Node strength
Finally, we calculate the node strength σ ij at each pixel.There are k = 16 possible ways of choosing 3 out ofthe 8 directions in such a way that each pair of chosendirections makes an angle of at least 90 degrees. Forexample, E, N, and SW comprise an allowable choice,but E, N, and NE do not. At each pixel, for each al-lowable choice of three directions, we calculate the min-imum of the three strand strengths in these directions,e.g., ˆ S kij = min { S Eij , S
Nij , S
SWij } . The node strength σ ij is the maximum of the 16 minima ˆ S kij subtracted by aminimum strength constant S , σ ij = max k =1 ... (cid:8) ˆ S kij (cid:9) − S . (4)Pixels with a positive node strength are called nodes . Themean node strength over the region of interest (ROI) iscalled the node strength of the region,NdStr = 1ROI (cid:88) ij ∈ ROI σ ij . (5)The purpose of subtracting the minimum strength con-stant S is to allow us to ignore short “strands”, which (a) (b) (c)(d) (e) (f) FIG. 3. Node strength algorithm illustrated on a region of trabecular bone taken from the full section shown in Fig. 1. (a)original image; (b) horizontal strands; (c) vertical strands; (d) diagonal strands (“northwest - southeast”); (e) diagonal strands(“northeast - southwest”); (f) node strength. are often just transverse sections of trabeculae with anapparent width of more than one pixel. The trabeculaein our images have an apparent width of approximatelytwo pixels, or 0.4 mm, which is higher than the true av-erage trabecular width, due to the 1 mm thickness of theCT image and the 0.2 mm pixel size. Thus, we wish toignore strands of bone that are only two pixels long. Thestrength of such a strand depends on the pQCT values ofthe two pixels. Following Eq. (2), the strength s of sucha strand with CT value b = b = b is s = b . Based onexamination of the resulting node strength images (suchas Fig. 3 (f)), we find ad hoc a minimum strength con-stant of 225 (corresponding to a CT value of 500, Eq. (1),corresponding to a BMD value of 331 mg/cm , which ishigher than the mean CT value for trabecular bone). Sostrands of lower strength than this will be ignored by thealgorithm.We note that our complete algorithm involves two thresh-olding steps: one at the beginning, when pQCT val-ues are converted to CT values of bone; and one at theend, when the minimum strength constant is subtracted.These are not equivalent to one thresholding step witha higher threshold. The purpose of the first thresholdfiltration is to ignore marrow. The purpose of the secondthreshold filtration is to ignore short strands, and it canbe seen as an alternative to skeletonisation.Numerical experiments have shown that relative (cross-subject) values of NdStr are stable with respect to thechoice of the bone threshold value ( a threshold ) and min-imum strength constant. Specifically, if either of theseconstants is varied by ±
10% from the values used in thisarticle, and NdStr is recomputed for all subjects, the re-sulting NdStr values are very strongly linearly correlatedwith the original NdStr values (Pearson correlation co- efficient r > .
995 in all cases). The absolute values ofNdStr depend on the choice of these constants: increasingthe bone threshold by 10% leads to a 40% mean decreasein NdStr (percent decrease averaged over all subjects),while increasing the minimum strength constant by 10%leads to an average decrease in NdStr of 11%. Never-theless, and as already mentioned, the relative (cross-subject) variation of NdStr remains constant.A key property of NdStr is that it depends on both thegeometry of the trabecular network and the CT valuesof the trabeculae. NdStr depends linearly on thresholdedCT values, if the CT values of all pixels are varied by thesame factor. It follows that images with the same geom-etry but different BMD will have different NdStr values.On the other hand, two images with the same BMD butdifferent geometry can have different NdStr values. Thisis apparent from the description of the method, and con-firmed by Figs. 4 and 5, as discussed below, and also bythe results discussed in Section IV. F. Illustration
To illustrate the analytical method, we now present theresults in visual form for an enlarged region near thebottom (posterior) of the slice in Fig. 1. The enlargedregion of the original image is shown in Fig. 3(a). Parts(b) to (e) of Fig. 3 show directional strand strengths,and part (f) shows the final node strength plot. Eachdirectional strand strength plot shows the sum of twostrand strengths at every pixel, in opposite directions:east/west; north/south; northwest/southeast; and north-east/southwest. In each of the directional strand strength (a) (b)
FIG. 4. Trabecular micro-architecture of a section with BMD 107 mg/cm . (a) original image; (b) the node strength of eachpixel. The mean node strength is 71.2. plots, the strands in the given direction are shown withthe highest intensity, but most of the trabeculae are stillvisible, even if faintly. In contrast, in the node strengthplot (part (f)), most of the trabeculae are invisible. Thisis because of the subtraction of the minimum strengthconstant. In this example, there are almost no nodes inthe right half of the image. This correctly describes themicro-architecture of the original image in that region,which contains many trabeculae but few that cross eachother to make a lattice-like micro-architecture. The lefthalf of the image contains many nodes. Notice that, inthe node strength plot, the nodes seem to be thicker thanin the original image. This is because the trabeculae inthe original image are actually slightly thicker than theyappear, the outer pixels being dimmer (i.e. lower CT val-ues) and thus not easily registered by the eye. Since theouter pixels near the apparent nodes in the original imageare almost as well-connected as pixels in the centres ofthe nodes, they have large node strengths, and are veryvisible in the node strength plot.The two specimens depicted in Figs. 4 and 5 have compa-rable bone mineral densities, but their trabecular micro-architecture are visibly different. In Fig. 4, the specimenhas a trabecular BMD of 107 mg/cm , which is near themedian value (97 mg/cm ) of the specimens in this study.On the left is the original image, and on the right is thenode strength plot. Notice that there are a lot of nodesin most of the outer areas, with the notable exception ofa region near the bottom left. The mean node strengthis 71.2.The specimen in Fig. 5 has a trabecular BMD of 94mg/cm , which is only 12% lower than that of the speci-men shown in Fig. 4, but it has substantially fewer nodes.The mean node strength is only 42.2, which is 40% lowerthan that for the specimen shown in Fig. 4. This reflectsthe lack of a strong lattice-like micro-architecture in theoriginal image. IV. RESULTS
The long range node strength of each pixel of each imagewas computed as described in the previous section, and amean node strength, NdStr, was calculated for each im-age. We also computed the bone mineral density (BMD)from the pQCT slices, as explained in Section III. A scat-ter plot showing NdStr versus BMD for all 26 specimensis shown in Fig. 6. There is a strong positive correlation,which we quantified in three ways. Firstly, Pearson’s cor-relation coefficient is r = 0 .
87, indicating a very stronglinear correlation. Secondly, since the scatter plot clearlysuggests a nonlinear relationship, we fitted an exponen-tial curve to the data and then found Pearson’s correla-tion coefficient to be r = 0 .
97. Thirdly, the Spearmanrank correlation is ρ = 0 .
98, indicating a very strong cor-relation. The Spearman rank correlation coefficient is arobust nonparametric correlation measure that is appro-priate when little is known about the distributions andnature of the correlation between the variables. We also compared NdStr with the node-terminus ra-tio, Nd/Tm, in the topological node-strut analysis intro-duced by Garrahan et al. The two measures are simi-lar in philosophy, because they both quantify the nodesin the trabecular network. However, the definition ofNd/Tm is highly localised: after the skeletonisation pro-cess has eroded the trabecular network to a thickness ofone pixel, each pixel is classified as node, strut, or termi-nus depending on its 3 × µ m,so the classification is made on the basis of a 30 µ m × µ m region. In contrast, node strength is semi-global,taking into account longer strands to a degree controlledby the transmission strength constant. In the presentstudy, with this constant set to 0 .
95, the method has acharacteristic length of 4 mm, in the sense described inSection III.The node-terminus ratio, Nd/Tm, was calculated usinghistomorphometry performed on bone biopsies as de-scribed in Section III. Recall that these biopsies were (a) (b)
FIG. 5. Trabecular micro-architecture of a section with BMD is 94 mg/cm . (a) original image; (b) the node strength of eachpixel. The mean node strength is 42.2. from the same regions of the same donors as the pQCTslices from which NdStr has been computed. Pearson’scorrelation coefficient for the relationship between NdStrand Nd/Tm is r = 0 .
62, and the Spearman rank correla-tion coefficient is ρ = 0 .
65. We also measured the corre-lation between Nd/Tm and trabecular BMD: Pearson’scorrelation coefficient is r = 0 .
64, and the Spearman rankcorrelation coefficient is ρ = 0 .
61. Using either measureof correlation, mean node strength is more strongly corre-lated with trabecular BMD than the node-terminus ratiois correlated with either of these variables.
20 40 60 80 100 120 140 BMD50100150NdStr
FIG. 6. Scatter plot of mean node strength vs. trabecularbone mineral density (mg/cm ). V. DISCUSSION
We have introduced a new morphometric measure forcharacterising the micro-architecture of trabecular bone, long range node strength , which measures the degree towhich a pixel in a 2-dimensional bone image has long-range connectivity in three or more directions, each atleast 90 degrees from the others. In addition, we havecalculated the mean node strength, NdStr, for each ofthe 26 bone samples considered in the study. We havefound that NdStr has a strong positive correlation with trabecular BMD ( r = 0 .
97, after exponential transforma-tion). Furthermore, we have ascertained a strong corre-lation ( r = 0 .
62) between NdStr and the established his-tomorphometric measure, node-terminus ratio (Nd/Tm).Moreover, qualitative comparison of images with similartrabecular BMD but different mean NdStr (see Figs. 4and 5) suggest that NdStr successfully quantifies how“lattice-like” the micro-architecture is.Further studies, including either clinical data or mea-sured bone strength, are needed in order to determinethe utility of this measure relative to other existing mor-phometric measures. Such studies could also determinethe most useful choice of the parameters that appearin the algorithm: the transmission constant, minimumstrength constant, and bending coefficients. For example,these constants could be chosen to maximise correlationswith bone strength. In the future, the sensitivity of themethod on these parameters and on different CT settings(e.g., CT pixel size, slice thickness, mean CT value ofbone) should be systematically investigated. Moreover,a future study on synthetic data could be used to verifythat NdStr measures structure and not just bone volumeor mass, and specifically that it finds nodes. Also differ-ent skeletal sites should be analysed with our approachin order to test its potential to describe structural differ-ences.We would like to note that our algorithm for computingnode strength is dependent on pixel size, and is thus un-suitable for absolute comparisons between studies. How-ever, this limitation is not unique to the node strengthmeasure. For example, Guggenbuhl et al. have showedusing CT images with different thickness (1 mm, 3 mm,5 mm, and 8 mm) that the outcome of texture analy-sis depends substantially on the slice thickness . In thepresent study we have used pQCT equipment with anin-plane pixel size of 200 µ m × µ m and a slice thick-ness of 1 mm. We have not conducted a formal inves-tigation into the influence of slice thickness on the longrange node strength, but it is fair to assume that thelong range node strength is similarly affected by the slicethickness. If further studies were to confirm the practicalvalue of the measure, an implementation could be devel-oped to produce node strengths that are broadly compa-rable between images with different pixel size. However,the technological development of high-resolution pQCTequipment is progressing very quickly, and already high-resolution pQCT scanners exist that can image a tibia atan isotropic pixel size of approximately 100 µ m, whichis comparable to the trabecular thickness in the humanproximal tibia , thus making the slice thickness less ofan issue.Some previous studies have investigated the trabecularbone micro-architecture using texture analysis appliedon X-ray images of bone . However, these tech-niques are not comparable to the method presented inthe present study as they are based on projections ofa 3-dimensional trabecular network on a 2-dimensionalplane, whereas our method is applied to 2-dimensionalsections obtained through the 3-dimensional trabecularnetwork. Nevertheless, Ranjanomennahary et al. veryrecently showed some significant correlation did exist be-tween radiograph based texture analysis and µ CT basedunbiased 3-dimensional measures of trabecular micro-architecture . Apostol et al. compared 3-dimensionalmeasures of micro-architecture based on 3-dimensionalsynchrotron radiation CT with more than 350 textureparameters obtained from simulated radiographs thatwere created from the synchrotron radiation CT datasets . They found using multiple regression analysisthat a combination of a subset of texture analysis param-eters correlated to the 3-dimensional measures of micro-architecture. However, they had to use at least three tex-ture analysis parameters in the correlation with each ofthe 3-dimensional measures of micro-architecture. Cortet et al. investigated CT images of the distal radius usingtexture analysis . The CT image used by Cortet et al. used similar pixel size to those used in the present study.They analysed the CT images using the traditional node-strut analysis and texture analysis including the greylevel run length method . As illustrated in the presentwork there is a moderate correlation ( r = 0 .
64) betweenthe node-strut analysis and the NdStr measure, whichwe believe illustrates that the NdStr measure capturessomewhat different information from that obtained withthe node-strut analysis. Grey level run length is basedon runs that have exactly the same grey level, and istherefore very sensitive to the choice of discretisation ofthe grey levels. In contrast, the NdStr measure treatsthe grey level as a continuous parameter, and thus it isable to detect runs with variable grey levels, and is con-sequently less sensitive to discretisation choices.Another advantage of the method is that the proposedmethod utilises all of the CT value information in thepQCT images, as no binarization of the images is per-formed.In the present study we have not applied the developedmethod to the histological sections directly. The rea-son for this is that the histological sections only cover a limited area of interest and thereby a limited trabecularlength, which renders the NdStr less meaningful. How-ever, this is a limitation of the histological examinationprocedure where it is only feasible to investigate a smallersample (biopsy) of a larger structure (the proximal tibia)and not a limitation of the proposed method. In ad-dition, the histological images are segmented into boneand marrow by their nature and it is thus not possibleto assign an intensity value to pixels, analogous to a CTvalue, which is needed by the algorithm in its presentform. However, it is possible to apply the concept oflong range node-strut analysis to such binary images butthis is outside the scope of the present investigation.In the present study the long range node-strut analy-sis was applied to pQCT images of the proximal tibia.However, we would like to stress that the method is notlimited to this skeletal site and thus can be applied to2-dimensional CT images obtained from any skeletal sitelike e.g. the vertebral body or the calcaneous.Finally, we note that the general method of long rangenode-strut analysis provides more than just the meannode strength. In the present study, we have focusedon mean node strength for simplicity, but the intermedi-ate measures of directional strand strength , used in thecomputation of node strength, may be useful in them-selves as a measure of directional strength. In the presentstudy the long range node-strut analysis has been appliedto 2-dimensional pQCT images obtained in the horizon-tal plane. However, the trabecular micro-architecture ofthe proximal tibia is mostly isotropic in the horizontalplane, whereas the micro-architecture in vertical direc-tion is highly anisotropic compared with the horizon-tal plane . As µ CT scanners become more and moreprevalent and as the pixel size and imaging capacity ofpQCT scanners steadily improve it will be an importantfuture task to generalise the long range node-strut anal-ysis into three dimensions and to apply the technique to3-dimensional data sets obtained from such equipment.In particular, the directional strand strength could beused to investigate anisotropic differences of the trabec-ular micro-architecture of such 3-dimensional data sets.
ACKNOWLEDGEMENTS
The data acquisition parts of this project were made pos-sible by grants from the Microgravity Application Pro-gram/Biotechnology from the Manned Spaceflight Pro-gram of the European Space Agency (ESA) (ESA project S. A. Goldstein, R. Goulet, and D. McCubbrey, “Measure-ment and Significance of Three-Dimensional Architecture to theMechanical Integrity of Trabecular Bone,” Calcified Tissue ,S127–S133 (1993). Li. Mosekilde, A. Viidik, and Le. Mosekilde, “Correlation Be-tween the Compressive Strength of Iliac and Vertebral TrabecularBone in Normal Individuals,” Bone , 291–295 (1985). G. Delling and M. Amling, “Biomechanical stability of theskeleton – it is not only bone mass, but also bone structurethat counts,” Nephrology Dialysis Transplantation , 601–606(1995). J.-Y. Rho, L. Kuhn-Spearing, and P. Zioupos, “Mechanical prop-erties and the hierarchical structure of bone,” Medical Engineer-ing & Physics , 92–102 (1998). T. Hildebrand, A. Laib, R. M¨uller, J. Dequeker, andP. R¨uegsegger, “Direct Three-Dimensional Morphometric Analy-sis of Human Cancellous Bone: Microstructural Data from Spine,Femur, Iliac Crest, and Calcaneus,” Journal of Bone and MineralResearch , 1167–1174 (1999). C. S. Rajapakse, J. S. Thomsen, J. S. Espinoza Ortiz, S. J.Wimalawansa, E. N. Ebbesen, Li. Mosekilde, and G. H.Gunaratne, “An expression relating breaking stress and den-sity of trabecular bone,” Journal of Biomechanics (2004),10.1016/j.jbiomech.2003.12.001. J. H. Kinney, J. S. St¨olken, T. S. Smith, J. T. Ryaby, and N. E.Lane, “An orientation distribution function for trabecular bone,”Bone , 193–201 (2005). R. Eastell, Li. Mosekilde, S. F. Hodgson, and B. L. Riggs, “Pro-portion of human vertebral body bone that is cancellous,” Jour-nal of Bone and Mineral Research , 1237–1241 (1990). H. Ritzel, M. Amling, M. P¨osel, M. Hahn, and G. Delling, “TheThickness of Human Vertebral Cortical Bone and Its Changesin Aging and Osteoporosis: A Histomorphometric Analysis ofthe Complete Spinal Column from Thirty-Seven Autopsy Speci-mens,” Journal of Bone and Mineral Research , 89–95 (1997). P. McDonnell, P. E. McHugh, and D. O’Mahoney, “Vertebralosteoporosis and trabecular bone quality,” Annals of BiomedicalEngineering , 170–189 (2007). D. R. Carter and W. C. Hayes, “The Compressive Behaviour ofBone as a Two-Phase Porous Structure,” The Journal of Boneand Joint Surgery , 954–962 (1977). E. N. Ebbesen, J. S. Thomsen, H. Beck-Nielsen, H. J. Nepper-Rasmussen, and Li. Mosekilde, “Lumbar Vertebral Body Com-pressive Strength Evaluated by Dual-Energy X-ray Absorptiom-etry, Quantitative Computed Tomography, and Ashing,” Bone , 713–724 (1999). Y. Jiang, J. Zhao, P. Augat, X. Ouyang, Y. Lu, S. Majum-dar, and H. K. Genant, “Trabecular bone mineral and calcu-lated structure of human bone specimens scanned by peripheralquantitative computed tomography: Relation to biomechanicalproperties,” Journal of Bone and Mineral Research , 1783–1790(1998). S. Majumdar, M. Kothari, P. Augat, D. C. Newitt, T. M. Link,J. C. Lin, T. Lang, Y. Lu, and H. K. Genant, “High-resolutionmagnetic resonance imaging: Three-dimensional trabecular bonearchitecture and biomechanical properties,” Bone , 445–454(1998). A. Odgaard, J. Kabel, B. van Rietbergen, M. Dalstra, andR. Huiskes, “Fabric and elastic principal directions of cancellousbone are closely related,” Journal of Biomechanics , 487–495(1997). D. Ulrich, B. van Rietbergen, A. Laib, and P. R¨uegsegger, “Theability of three-dimensional structural indices to reflect mechan-ical aspects of trabecular bone,” Bone , 55–60 (1999). T. Baum, J. Carballido-Gamio, M. B. Huber, D. M¨uller, R. Mon-etti, C. R¨ath, F. Eckstein, E. M. Lochm¨uller, S. Majumdar, E. J.Rummeny, T. M. Link, and J. S. Bauer, “Automated 3d trabec-ular bone structure analysis of the proximal femur—predictionof biomechanical strength by ct and dxa,” Osteoporosis Interna-tional , 1553–1564 (2010). J. S. Thomsen, E. N. Ebbesen, and Li. Mosekilde, “RelationshipsBetween Static Histomorphometry and Bone Strength Measure-ments in Human Iliac Crest Bone Biopsies,” Bone , 153–163(1998). J. S. Thomsen, E. N. Ebbesen, and Li. Mosekilde, “PredictingHuman Vertebral Bone Strength by Vertebral Static Histomor-phometry,” Bone , 502–508 (2002). P. I. Saparin, W. Gowin, J. Kurths, and D. Felsenberg, “Quan-tification of cancellous bone structure using symbolic dynamicsand measures of complexity,” Physical Review E , 6449 (1998). G. Dougherty, “A comparison of the texture of computed tomog-raphy and projection radiography images of vertebral trabecularbone using fractal signature and lacunarity,” Med Eng Phys ,313–321 (2001). S. Prouteau, G. Ducher, P. Nanyan, G. Lemineur, L. Benhamou,and D. Courteix, “Fractal analysis of bone texture: a screen-ing tool for stress fracture risk?” European Journal of ClinicalInvestigation , 137–142 (2004). P. I. Saparin, J. S. Thomsen, S. Prohaska, A. Zaikin, J. Kurths,H.-C. Hege, and W. Gowin, “Quantification of spatial structureof human proximal tibial bone biopsies using 3D measures ofcomplexity,” Acta Astronautica , 820–830 (2005). P. I. Saparin, J. S. Thomsen, J. Kurths, G. Beller, andW. Gowin, “Segmentation of bone ct images and assessment ofbone structure using measures of complexity,” Med. Phys. ,3857–3873 (2006). N. Marwan, J. Kurths, and P. Saparin, “Generalised RecurrencePlot Analysis for Spatial Data,” Physics Letters A , 545–551(2007). N. Marwan, P. Saparin, and J. Kurths, “Measures of complexityfor 3D image analysis of trabecular bone,” The European Phys-ical Journal – Special Topics , 109–116 (2007). N. Marwan, J. Kurths, J. S. Thomsen, D. Felsenberg, andP. Saparin, “Three dimensional quantification of structures intrabecular bone using measures of complexity,” Physical ReviewE , 021903 (2009). B. van Rietbergen, S. Majumdar, W. Pistoia, D. C. Newitt,M. Kothari, A. Laib, and P. R¨uegsegger, “Assessment of can-cellous bone mechanical properties from micro-FE models basedon micro-CT, pQCT and MR images,” Technology and HealthCare , 413–420 (1998). X. E. Guo and C. H. Kim, “Mechanical Consequence of Trabec-ular Bone Loss and Its Treatment: A Three-dimensional ModelSimulation,” Bone , 404–411 (2002). L. Pothuaud, B. Van Rietbergen, L. Mosekilde, O. Beuf,P. Levitz, C. L. Benhamou, and S. Majumdar, “Combinationof topological parameters and bone volume fraction better pre-dicts the mechanical properties of trabecular bone,” Journal ofBiomechanics , 1091–1099 (2002). U. Wolfram, L. O. Schwen, U. Simon, and M. Rumpf, “Statisticalosteoporosis models using composite finite elements: A parame-ter study,” Journal of Biomechanics , 2205–2209 (2009). N. J. Garrahan, R. W. E. Mellish, and J. E. A. Compston, “Anew method for the two-dimensional analysis of bone structure inhuman iliac crest biopsies,” Journal of Microscopy , 341–349(1986). P. Saparin, W. Gowin, and D. Felsenberg, “Comparison of boneloss with the changes of bone architecture at six different skele-tal sites using measures of complexity,” Journal of GravitationalPhysiology , 177–178 (2002). J. S. Thomsen, A. Laib, B. Koller, S. Prohaska, Li. Mosekilde,and W. Gowin, “Stereological measures of trabecular bone struc-ture: comparison of 3D micro computed tomography with 2Dhistological sections in human proximal tibial bone biopsies,”Journal of Microscopy , 171–179 (2005). J. S. Thomsen, E. N. Ebbesen, and Li. Mosekilde, “A NewMethod of Comprehensive Static Histomorphometry Applied onHuman Lumbar Vertebral Cancellous Bone,” Bone , 129–138(2000). L. Lam and C. Y. Suen, “Thinning methodologies – a comprehen-sive survey,” IEEE Transaction on Pattern Analysis and MachineIntelligence , 869–885 (1992). J. E. Compston, K. Yamaguchi, P. I. Croucher, N. J. Garrahan,P. C. Lindsay, and R. W. Shaw, “The effects of gonadotrophin-releasing hormone agonists on iliac crest cancellous bone struc-ture in women with endometriosis,” Bone , 261–267 (1995). W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flan-nery,
Numerical Recipes in C — the Art of Scientific Computing ,2nd ed. (Cambridge University Press, Cambridge, 1992). M. Ding, A. Odgaard, F. Linde, and I. Hvid, “Age-related vari-ations in the microstructure of human tibial cancellous bone,”Journal of Orthopaedic Research , 615–621 (2002). P. Guggenbuhl, D. Chappard, M. Garreau, J.-Y. Bansard,G. Chales, and Y. Rolland, “Reproducibility of ct-based bonetexture parameters of cancellous calf bone samples: Influenceof slice thickness,” European Journal of Radiology , 514–520(2008). M. Ding and I. Hvid, “Quantification of age-related changes inthe structure model type and trabecular thickness of human tibialcancellous bone,” Bone , 291–295 (2000). W. G. M. Geraets, P. F. van der Stelt, C. J. Netelenbos, andP. J. M. Elders, “A new method for automatic recognition of theradiographic trabecular pattern,” Journal of Bone and MineralResearch , 227–233 (1990). D. Chappard, P. Guggenbuhl, E. Legrand, M. F. Basl´e, andM. Audran, “Texture analysis of x-ray radiographs is correlatedwith bone histomorphometry,” Journal of Bone and Mineral Metabolism , 24–29 (2005). D. Chappard, F. Pascaretti-Grizon, Y. Gallois, P. Mercier, M. F.Basl´e, and M. Audran, “Medullar fat influences texture analysisof trabecular microarchitecture on x-ray radiographs,” EuropeanJournal of Radiology , 404–410 (2006). M. B. Huber, J. Carballido-Gamio, K. Fritscher, R. Schubert,M. Haenni, C. Hengg, S. Majumdar, and T. M. Link, “Devel-opment and testing of texture discriminators for the analysis oftrabecular bone in proximal femur radiographs,” Medical Physics , 5089–5098 (2009). C. M. Korstjens, Li. Mosekilde, R. J. Spruijt, W. G. M. Ger-aets, and P. F. van der Stelt, “Relations between radiographictrabecular pattern and biomechanical characteristics of humanvertebrae,” Acta Radiologica , 618–624 (1996). P. Ranjanomennahary, S. S. Ghalila, D. Malouche, A. Mar-chadier, M. Rachidi, C. Benhamou, and C. Chappard, “Com-parison of radiograph-based texture analysis and bone mineraldensity with three-dimensional microarchitecture of trabecularbone,” Medical Physics , 420–428 (2011). L. Apostol, V. Boudousq, O. Basset, C. Odet, S. Yot, J. Tabary,J.-M. Dinten, E. Boiler, P.-O. Kotzki, and F. Peyrin, “Relevanceof 2d radiographic texture analysis for the assessment of 3d bonemicro-architecture,” Medical Physics , 3546–3556 (2006). B. Cortet, N. Boutry, P. Dubois, P. Bourel, A. Cotten, andX. Marchandise, “In vivo comparison between computed tomog-raphy and magnetic resonance image analysis of the distal radiusin the assessment of osteoporosis,” Journal of Clinical Densito-metry , 15–26 (2000). B. Cortet, P. Dubois, N. Boutry, P. Bourel, A. Cotten, andX. Marchandise, “Image analysis of the distal radius trabecu-lar network using computed tomography,” Osteoporosis Interna-tional , 410–419 (1999). A. Chu, C. M. Sehgal, and J. F. Greenleaf, “Use of gray valuedistribution of run lengths for texture analysis,” Pattern Recog-nition Letters , 415–420 (1990). M. M. Galloway, “Texture analysis using gray level run lengths,”Computer Graphics and Image Processing4