Volume-Based Fabric Tensors through
Lattice-Boltzmann Simulations
Rodrigo Moreno and ¨
Orjan Smedby
Center for Medical Image Science and Visualization (CMIV) Department of Medical and Health Sciences (IMH) Link¨oping University, Campus US, 58185 Link¨oping, Sweden
Email: {rodrigo.moreno,orjan.smedby@liu.se}
Abstract—This paper introduces a new methodology to com-pute fabric tensors from computational fluid dynamics simula-tions performed through the lattice-Boltzmann method. Trabecu-lar bone is modeled as a pipeline where a synthetic viscous fluid can flow from a single source located at the center of a spherical region of interest toward its boundaries. Two fabric tensors are computed from local velocities at the steady state estimated from the simulations, a tortuosity and a normalized tortuosity tensor. The main advantage of the proposed fabric tensors is that, unlike previous approaches, they intentionally disregard the trabecular termini in the computations, which do not play an important role in the estimation of trabecular bone quality. Thus, the proposed fabric tensors are less prone than previously proposed ones to unnecessary reductions of anisotropy related to the the presence of trabecular termini. The results of experiments conducted on synthetic and micro-computed tomography data in 2D and 3D show the artificial fluid flowing inside the trabecular bone has negligible velocities at trabecular termini, reducing in that way their influence in the estimation of the proposed fabric tensors.
Keywords—Fabric tensors, trabecular bone, computational fluid dynamics, lattice-Boltzmann method, tortuosity.
I. INTRODUCTION
Bone in humans is mainly composed by cortical and trabecular bone (cf. Fig. 1). As a matter of fact, skeletal sites that are more prone to osteoporosis-related fractures such as vertebrae the femoral neck or the radius have a relatively thin layer of cortical bone. Thus, the mechanical competence of bone in those sites mainly relies on the quality of trabecular bone. One of the ultimate goals of trabecular bone research in medicine is to determine the effect of pathological conditions, such as osteoporosis and osteoarthritis, and their treatments on the quality of trabecular bone. One of the parameters that can be used to evaluate the quality of trabecular bone is its anisotropy. For example, evidence supporting that changes in the anisotropy and orientation of trabecular bone are associated with osteoporosis has been reported [1], [2], [3].
Fabric tensors aim at modeling through tensors both anisotropy and orientation of a material of interest (usually referred to as phase in mechanics of materials) with respect to another one. In trabecular bone research, these two phases correspond to trabecular bone and bone marrow respectively. Fabric tensors are semi-global measurements in the sense that they are computed in relatively large neighborhoods, which are assumed quasi-homogeneous. It is worthwhile to mention that, from a clinical perspective, anisotropy seems more important than the orientation of the fabric tensors [4], [5], [6].
Fig. 1. 2D slices of micro computed tomography scans of a tibia (left) and a vertebra (right). Cortical bone is the compact tissue in the outer layer of the bone, while trabecular bone is the porous network inside the cortical bone.
Fabric tensors can be computed in many different ways (cf. [7] for a comprehensive review). The most used fabric tensor is the Mean Intercept Length (MIL) tensor, which aims at measuring anisotropy and orientation of the boundary between phases [8]. Related boundary-based tensors are the Generalized MIL tensor (GMIL) [8] and the global Gradient Structure Tensor (GST) [9]. However, boundary-based fabric tensors are unable to detect anisotropies and orientations generated by the shape of one of the phases. On the contrary, volume-based methods for computing fabric tensors are suitable for detecting such kinds of anisotropies. The proposed method in this paper belongs to this second family.
Most of the volume-based fabric tensors use a sampling procedure where statistics are computed from sampled points. For example, the tensor scale [10], the star volume and length distributions [11] and the tortuosity tensor proposed in [12] compute directional statistics from intercepting lines traced from the sampled points. Instead of lines, the sampling sphere orientation distribution [13] use spheres at the sampled voxels. The most important inconvenience of all these methods is that they treat trabecular termini (i.e., free-ended trabeculae) and struts (i.e., trabeculae that bear the main mechanical loads) in the same way. The difference between both is that, while struts bear the main responsibility for the mechanical competence of bones, trabecular termini have little or no contribution to the global mechanical competence of bone [14], [11], [15]. This can be seen in the simulation shown in Figure 21. It is
noteworthy to say that some studies have found that trabec-ular termini could have a potential power as a biomarker in osteoporosis, since it is related to the bone remodeling process [16]. However, our point here is that termini must be treated differently for computing fabric tensors, since such tensors
Fig. 2. von Mises stress computed on a synthetic 2D model with 2 termini. A horizontal displacement (1 m) was set at the top-left boundary of the strut, and the lower end boundary was fixed. As shown, the stress at the termini is very low compared to the stress at the strut.
must be correlated with mechanical properties of the bone [17], [18]. In addition to the aforementioned shortcoming, another common drawback of traditional volume-based methods is that their accuracy directly depends on the sampling scheme. Thus, usually a huge amount of testing is required to obtain a reasonable accuracy.
In this paper we explore a new methodology to com-pute volume-based fabric tensors from computational fluid dynamics (CFD) simulations. The proposed method treats a phase (e.g., trabecular bone) as a pipeline where a synthetic viscous fluid can flow from a single source located at the center of the region of interest toward the boundaries of a spherical region of interest. In particular, velocities and moment fluxes can be estimated locally from such simulations. In a second step two fabric tensors are computed: one from the orientation distribution of velocities and another from local moment fluxes. For the simulations, we take advantage of the lattice-Boltzmann method (LBM) [19], which has made CFD simulations affordable in many applications in the last few years. The LBM is more appropriate for fluids with low Reynolds number, which is actually the case for the laminar and creeping flows that are used in the current application.
Our approach has different advantages with respect to previously proposed methods. First and foremost, the influence of the termini in the estimation of fabric tensors is reduced by the fact that velocity of laminar flows in stagnant pockets is almost null. Second, unlike previous approaches, it does not require any sampling procedure, thus accurate results can be attained once the steady state of the CFD simulation is reached. The method is inspired in mechanics, so it is expected to yield better correlations with other mechanical properties, in particular with the stiffness tensor. CFD simulations have previously been used in bone research for modeling implanta-tion of bone cement in vertebrae [20], [21], [22]. However, to the best of our knowledge, our approach is the first attempt to compute fabric tensors by simulating synthetic fluid dynamics inside trabecular bone. Moreover, unlike the applications in bone cement fluid dynamics, a single CFD simulation is necessary, which makes it possible to compute fabric tensors in a reasonable amount of time.
The paper is organized as follows. Section II describes the proposed method for computing fabric tensors. Section III shows and discusses the results of experiments conducted on synthetic datasets and images acquired through Micro
Computed Tomography (µCT). Finally, Section IV makes some concluding remarks.
II. METHOD
A. Lattice-Boltzmann Simulation
The LBM has extensively been used in CFD simulations [23] and it is becoming more used in image processing applications, such as anisotropic filtering [24], dithering [25], inpainting [26] and active contour segmentation [27] among many others. The advantage of the LBM is that it can be coded for parallel machines and GPUs straightforwardly. Thus, implementations for dealing with complex hyperbolic partial differential equations (PDEs), such as the Navier-Stokes PDE, can be obtained efficiently.
The LBM simulates the interaction among synthetic par-ticles distributed in a lattice. At a microscopic level, such interactions are governed by two processes: namely, collision and streaming, where the particles are allowed to collide and propagate in a local neighborhood respectively.
Let f (~x, t) be the probability of finding a fluid particle at position ~x and time t. Such a probability can be discretized into several directional components, fα, each of them interacting
with the neighbor in the ~eα direction. By convention, f0
rep-resents the particles that are not propagated to the neighbors. The collision and streaming steps are governed by:
fα0(~x, t + ∆t) = fα(~x, t) + Ωα(~x, t), (1)
fα(~x + ~eα∆t, t + ∆t) = fα0(~x, t + ∆t) (2)
respectively, where ∆t is the time step and Ωαis the collision
operator. Ωαis of utmost importance in LBM, since different
PDEs can be solved by applying different operators. One of the most used operators in LBM is the Bhatnagar-Gross-Krook (BGK) function which is given by:
Ωα(~x, t) = τ (fαeq(~x, t) − fα(~x, t)), (3)
that is, the collisions are governed by a relaxation process with respect to an equilibrium function fαeq with a parameter
τ ∈ (0, 2) [19]. The equilibrium function can be derived from locally conserved quantities such as mass, momentum and energy for a specific lattice. As an example, the equilibrium function for the D2Q9 lattice (i.e., two dimensions and nine discrete velocity orientations) is given by:
fαeq = ωαρ 1 +3~e T α~v c2 + 9(~eT α~v)2 2c4 − 3(~vT~v)2 2c2 , (4)
with ωα being a set of lattice-dependent weights to achieve
Galilean invariance, ρ the density, ~v the estimated velocity and c the so-called lattice speed of sound, which is a lattice-dependent constant. With this equilibrium function, LBM becomes a solver for the Navier-Stokes equation at the macro-scopic limit for incompressible fluids [19].
Unlike the finite element method (FEM), LBM node creation is straightforward, as every voxel corresponds to a node. This has the additional advantage that boundary and equilibrium conditions can be set voxel-wise.
B. Estimation of Fabric Tensors
Once the probability functions fα are computed, different
parameters can be estimated, such as the local density ρ and velocity ~v. These parameters are given by:
ρ(~x, t) =X ∀α fα(~x, t), (5) v(~x, t) = 1 ρ(~x, t) X ∀α fα(~x, t)~eα, (6)
From these quantities, we propose to compute two fabric tensors: the tortuosity tensor (T) and the normalized tortuosity tensor eT. On the one hand, T is given by:
T = lim t→∞ Z ~ x∈V ~v(~x, t) ~v(~x, t)TdV, (7) where V is the spherical region of interest and the limit guarantees that the steady state has been attained. Following a similar approach as in [28], it can be shown that (7) estimates a tortuosity tensor for incompressible fluids and non-reentrant flows. Actually, (7) can be seen as a generalization of the tortuosity measure proposed in [28] where a tensor is computed instead of a scalar estimation of tortuosity in the main orientation of the flow.
One problem of T is that it gives more weight to voxels that are closer to the synthetic source of fluid, since the velocity decreases towards the boundaries of the region of interest. This effect is mainly due to the high viscosity of the simulated fluid. Unfortunately, the viscosity cannot be reduced in order to guarantee laminar flow. A possibility to tackle this problem is to use velocity directions instead of velocities in (7). Indeed, this approach is appealing, since all voxels inside the region of interest will have the same weight in the computations of the fabric tensor, as they should. However, low velocity areas (e.g., at the termini) will also be enhanced, which is contrary to the main purpose of the method. Assuming that the magnitude of the velocities follows a normal distribution, we propose to only consider the voxels with velocities larger than -2 standard deviations from the mean magnitude of the velocity. Thus, the proposed normalized tortuosity tensor eT is given by:
e T = lim t→∞ Z {~x∈V \||~v(~x,t)||≥κ} ~v(~x, t) ~v(~x, t)T ~v(~x, t)T~v(~x, t)dV, (8)
where κ is the mean magnitude of the velocities minus two standard deviations.
An additional problem can occur when the source of fluid is placed in a terminus. In such a case, that specific terminus will inevitably be considered in the fabric tensor. However, its influence in the total computation of the fabric tensor is usually negligible, due to the large amount of termini present in trabecular bone. Certainly, that problem can be solved by not placing the source of fluid in termini.
Notice that the proposed tensors can be computed indepen-dently for each phase. That means that in general four tensors can be computed with the proposed methodology: Tt, Tm, eTt
and eTm, where the subindices t and m stand for trabecular
bone and bone marrow respectively.
III. EXPERIMENTALRESULTS
Simulations have been coded in Python using the Sailfish CFD simulation platform2. This platform is especially
de-signed for performing lattice-Boltzmann simulations on GPUs. The simulations have been run in an NVIDIA Quadro 4000 with 2GB of RAM. Lattices D2Q9 (two dimensions, 9 velocity orientations in a 3 x 3 stencil) and D3Q19 (three dimensions, 19 velocity orientations in a 3 x 3 x 3 stencil) have been used in the 2D and 3D simulations respectively. Parameter τ in (3), ∆t and the Reynolds number of the simulated fluid has been set to 1, 1 × 10−4and 50 respectively. The following boundary conditions have been set: outflow nodes for the boundaries of the spherical (circular in 2D) region of interest and no-slip condition for nodes belonging to the complementary phase. In addition, a fluid source has been simulated by setting the feq
α = 0.95 for the nearest 18 (8 in 2D) neighbors of a central
voxel, where α indicates the direction from the center to these voxels. The center of the source has been put as far as possible from the complementary phase and as close as possible to the center of the sample. Special effort has been done for not to place the source of synthetic fluid in termini.
Anisotropy has been measured as ai= λi/λ1with λibeing
the i-th largest eigenvalue of the tensor. Experiments have been conducted on synthetic models and images of trabecular bone acquired through micro-CT. For all the experiments, segmentation of the micro-CT images has been performed as described in [29]. Notice that the anisotropy values computed with the proposed tensors are not comparable with any other fabric tensor. Actually, assessing the appropriateness of a fabric tensor is application-dependent. In particular, it is interesting to test the power of the proposed tensors for predicting the stiffness tensor. However, this aim is beyond the scope of this paper and is actually part of our ongoing research. Despite this, anisotropies of the MIL tensor (computed as in [8]) have been included in the results as a reference.
A. Experiments in 2D
Figure 3 shows the magnitude of the velocity vectors computed through the simulation and T and eT estimated on two synthetic images “cross” and “H”. The main orientation of both tensors was vertical for both images. Table I shows the anisotropies of these two tensors and the MIL tensor. Regarding the cross, T and eT are almost stick-shaped tensors, since a2 ≈ 0, compared to the non-depreciable 0.3 of the
MIL tensor. Thus, the aim of discarding the termini in the computation of the fabric tensors is well accomplished through the proposed tensors. The source of fluid has been located in different positions without significant changes in the obtained tensors. As for the “H” image, T was shown sensitive to the position of the source of fluid, while such a position did not have an important influence on eT. Notice that eT and the MIL tensor yield a similar anisotropy.
Figure 4 shows the magnitude of the velocity vectors computed through the simulation and Tm and eTm estimated
on a slice of an image of a radius acquired through micro-CT for the bone marrow. Notice that the simulation requires to have a single connected component for the phase of interest.
Fig. 3. Top: color encoded logarithm of the magnitude of the velocities computed through the lattice-Boltzmann simulation on two synthetic images. The steady state was attained at 5000 iterations. Each simulation took 5.3 s. Bottom: superimposed tensors Tm(left) andTem(right) for the two different
simulations. Plotted red and blue tensors correspond to the simulations shown at the top left and top right respectively.
Fig. 4. Top: color encoded logarithm of the magnitude of the velocities computed through the lattice-Boltzmann simulation on a 2D image of bone marrow using from two slightly different regions of interest. The steady state was attained at 25000 iterations. Each simulation took 17.3 s. Bottom: super-imposed tensors Tm(left) andeTm(right) for the two different simulations.
Plotted red and blue tensors correspond to the simulations shown at the top left and top right respectively.
In practice, that means that, in 2D, the proposed tensors can only be computed for bone marrow, since the connectedness of trabecular bone is only attained in 3D, as shown in the figure. As presented, the velocities are largely reduced in nearly closed areas.
An important issue for computing tensors for bone marrow is that it is less structured than trabecular bone. Thus, the estimated tensors are more dependent on the position of the synthetic source of fluid (cf. Figure 4 and Table I). A strategy to tackle this problem is to perform different simulations and perform an averaging of the computed tensors. It is necessary to consider the use of an appropriate averaging method in order to avoid the “swelling” effect described in [30].
B. Experiments in 3D
Figures 5 and 6 show renderings of the computed velocities for images of a radius and a vertebra acquired through
micro-TABLE I. ANISOTROPIESa2OF THE COMPUTED TENSORS FOR THE EXPERIMENTS IN2D. FORTANDTeTHE MEAN AND STANDARD
DEVIATION(IN PARENTHESIS)OFa2IS REPORTED FOR10DIFFERENT SIMULATIONS WHERE THE SOURCE OF THE SYNTHETIC FLUID IS PLACED
IN DIFFERENT POSITIONS.
Image T Te MIL Cross 0.0093 (0.001) 0.0140 (0.001) 0.2995 H 0.2587 (0.052) 0.1338 (0.007) 0.1399 Radius(bone marrow) 0.7169 (0.077) 0.8472 (0.072) 0.9541
Fig. 5. Rendering of the magnitude of the velocity of the simulated fluid flowing through the trabeculae of a section of a radius. Red and yellow represent high and low velocities respectively.
CT. As shown, the main struts of the trabeculae have higher values of velocity. Table II shows the anisotropies of the fabric tensors computed for these two images with different locations of the source of fluid. The problems described in the previous subsection for placing the source of fluid in bone marrow are also present in 3D. Fortunately, that is not the case for trabecular bone in 3D, since it is possible to find appropriate positions for the source of fluid. This point can be shown through the idealized rod-plate model of trabecular bone. This model assumes that trabecular bone is a network of rod-like and plate-like struts. It is clear that the flow at the rods will go either in one or the opposite direction, depending on the position of source of fluid. Thus, the contribution of the velocities at the rod to eT is the same, regardless of the position of source of fluid. Thus, it is ideal to place the source of fluid in a rod. It is important to remark that the rod-plate model is not accurate [31]. Despite this, our experiments have shown that eT has little variation with respect to the location of the source of fluid for the trabecular bone phase as shown in Table II. Thus, eT can be used for measuring anisotropies of trabecular bone.
Regarding the orientation, the angles between the main direction of the MIL tensor and the one of Tt and eTt were
below 5 degrees in all cases. As for bone marrow, the angular differences were larger, with a mean of 25.6 degrees and a standard deviation of 10.2 degrees.
IV. DISCUSSION
A new methodology to compute fabric tensors from CFD simulations performed through LBM has been presented. The phase of interest is modeled as a pipeline where a synthetic viscous fluid can flow from a single source toward the bound-aries of a spherical region of interest. Fabric tensors are computed from local velocities estimated from the simulations. The main advantage of the proposed method is that it reduces
TABLE II. ANISOTROPIESa2ANDa3OF THE COMPUTED TENSORS FOR THE EXPERIMENTS IN3DFOR TRABECULAR BONE(TB)AND BONE MARROW (BM). FORTANDeTTHE MEAN AND STANDARD DEVIATION(IN PARENTHESIS)OFa2ANDa3IS REPORTED FOR10DIFFERENT SIMULATIONS WHERE THE
SOURCE OF THE SYNTHETIC FLUID IS PLACED IN DIFFERENT POSITIONS. THEMILTENSOR IS INCLUDED AS A REFERENCE. THEMILTENSOR FOR TRABECULAR BONE AND BONE MARROW ARE EQUIVALENT.
Image T eT MIL
a2 a3 a2 a3 a2 a3
Radius (tb) 0.8233 (0.101) 0.5129 (0.121) 0.4796 (0.023) 0.3374 (0.031) 0.9091 0.6404 Radius (bm) 0.6166 (0.120) 0.3410 (0.150) 0.7825 (0.031) 0.6239 (0.070) id. id. Vertebra (tb) 0.3879 (0.121) 0.3293 (0.077) 0.6638 (0.031) 0.6274 (0.033) 0.941 0.1519 Vertebra (bm) 0.8039 (0.321) 0.4340 (0.220) 0.8599 (0.041) 0.7215 (0.050) id. id.
Fig. 6. Rendering of the magnitude of the velocity of the simulated fluid flowing through the trabeculae of a section of a vertebra. Red and yellow represent high and low velocities respectively.
the influence of the termini in the estimation of fabric tensors. Discarding the termini might result in better estimations of mechanical competence of trabecular bone through fabric tensors.
Results of experiments conducted on synthetic and micro-computed tomography data show that eT is promising for com-puting anisotropy and orientation of trabecular bone. Results also show that the tensors computed on bone marrow are less reliable with a large dependence on the position of the simulated source of the synthetic fluid. The reason for this is that trabecular bone is more structured than bone marrow, so the fluid movement is largely restricted, making it less dependent on the source of the fluid.
GPU implementations of the LBM have made CFD simu-lations affordable for computing fabric tensors. However, it is necessary to improve its efficiency, in order to close the current computational effort gap between volume- and boundary-based methods. As an example, processing a volume of 200 x 200 x 200 voxels with the proposed method requires 4.18 minutes using D3Q19 or 3.35 minutes using D3Q15 compared to 1.4 seconds taken by the CPU implementation of the GMIL tensor [8]. As already mentioned, other volume-based fabric tenors are also much more expensive than their boundary-based counterparts.
As expected, anisotropies obtained from eT are in general different from those extracted from the MIL tensor. As already mentioned, this result is also expected with any other type of fabric tensor. In general, the proposed tensors tend to be more anisotropic since they intentionally disregard the termini in the computations.
It is worthwhile to mention that fabric tensors have two po-tential uses in trabecular bone research. On the one hand, com-bined with histomorphometric parameters, anisotropy
measure-ments can potentially be used for diagnosis and following up of treatments in clinical practice [1], [2], [3]. Under this perspective, clinical studies must be performed in order to assess the power of different fabric tensors for diagnostic purposes. On the other hand, fabric tensors have been used to predict mechanical properties of trabecular bone [17], [18]. In this case, studies must be performed in order to assess the power of the proposed fabric tensors for predicting the stiffness tensor with respect to previous approaches. Both type of studies are beyond the scope of this paper and are part of our ongoing research in fabric tensors.
In particular, we are interested in testing with mechanical data the main hypothesis we made in this paper, i.e., termini are not relevant for predicting mechanical competence of trabecular bone. In [16], a relatively small amount of termini per volume was reported for healthy bones. This result suggests that termini are generated as a product of the normal bone remodeling process without any related mechanical function. Since the termini are largely increased in osteoporotic bones [16], we expect eT to be more appropriate for predicting the stiffness tensor than other fabric tensors, especially in osteoporotic bones.
Furthermore, our ongoing research also include proposing a variant of eT that can be directly applied to gray-scale images. Such an extension would allow us to avoid the segmentation step, which can introduce inaccuracies in images acquired in vivo.
ACKNOWLEDGMENT
We thank Prof. Osman Ratib from the Service of Nuclear Medicine of the Geneva University Hospitals for providing the micro-CT data of the vertebra and Andres Laib and Torkel Brismar for providing the micro-CT data of the radius. This research has been supported by the Swedish Research Council (VR), grants no. 2011-5197 and 2012-3512, and the Swedish Heart-Lung Foundation (HLF), grant no. 2011-0376.
REFERENCES
[1] T. E. Ciarelli, D. P. Fyhrie, M. B. Schaffler, and S. A. Goldstein, “Variations in three-dimensional cancellous bone architecture of the proximal femur in female hip fractures and in controls,” J Bone Miner Res, vol. 15, no. 1, pp. 32–40, 2000.
[2] J. Homminga, B. Van-Rietbergen, E. M. Lochm¨uller, H. Weinans, F. Eckstein, and R. Huiskes, “The osteoporotic vertebral structure is well adapted to the loads of daily life, but not to infrequent “error” loads,” Bone, vol. 34, no. 3, pp. 510–516, 2004.
[3] J. M. Kreider and S. A. Goldstein, “Trabecular bone mechanical properties in patients with fragility fractures,” Clin Orthop Relat Res, vol. 467, no. 8, pp. 1955–1963, 2009.
[4] R. Huiskes, “If bone is the answer, then what is the question?” J Anatomy, vol. 197, pp. 145–156, 2000.
[5] R. Moreno, M. Borga, E. Klintstr¨om, T. B. Brismar, and ¨O. Smedby, “Correlations between fabric tensors computed on cone beam and micro computed tomography images,” in Proc Comput Vis Med Image Process (VIPIMAGE). CRC Press, 2013, pp. 393–398.
[6] ——, “Anisotropy estimation of trabecular bone in gray-scale: com-parison between cone beam and micro computed tomography data,” in Developments in Medical Image Processing and Computational Vision, J. Tavares and R. N. Jorge, Eds. Springer, 2014, in press.
[7] R. Moreno, M. Borga, and ¨O. Smedby, “Techniques for computing fabric tensors: A review,” in Visualization and Processing of Tensors and Higher Order Descriptors for Multi-Valued Data, C. F. Westin, A. Vilanova, and B. Burgeth, Eds. Springer, 2014, pp. 271–292. [8] ——, “Generalizing the mean intercept length tensor for gray-level
images,” Med Phys, vol. 39, no. 7, pp. 4599–4612, 2012.
[9] Z. Tabor and E. Rokita, “Quantifying anisotropy of trabecular bone from gray-level images,” Bone, vol. 40, no. 4, pp. 966–972, 2007. [10] Y. Liu, P. K. Saha, and Z. Xu, “Quantitative characterization of
trabecular bone micro-architecture using tensor scale and multi-detector CT imaging,” in Med Image Comput Comput-Assist Interv (MICCAI), ser. LNCS, vol. 7510, 2012, pp. 124–131.
[11] A. Odgaard, J. Kabel, B. van Rietbergen, M. Dalstra, and R. Huiskes, “Fabric and elastic principal directions of cancellous bone are closely related,” J Biomech, vol. 30, no. 5, pp. 487–495, 1997.
[12] J. Bear and Y. Bachmat, Introduction to modeling of transport phenom-ena in porous media. Kluwer, 1990.
[13] P. Varga and P. K. Zysset, “Assessment of volume fraction and fabric in the distal radius using HR-pQCT,” Bone, vol. 45, no. 5, pp. 909–917, 2009.
[14] D. Stauffer and A. Aharony, Introduction to percolation theory, 2nd ed. Taylor & Francis, 1994.
[15] Z. Tabor, “Novel algorithm detecting trabecular termini in µCT and MRI images,” Bone, vol. 37, no. 3, pp. 395 – 403, 2005.
[16] J. E. Aaron, P. A. Shore, R. C. Shore, M. Beneton, and J. A. Kanis, “Trabecular architecture in women and men of similar bone mass with and without vertebral fracture: II. three-dimensional histology,” Bone, vol. 27, no. 2, pp. 277–282, 2000.
[17] P. K. Zysset, “A review of morphology-elasticity relationships in human trabecular bone: theories and experiments,” J Biomech, vol. 36, no. 10, pp. 1469–1485, 2003.
[18] T. Gross, D. Pahr, and P. Zysset, “Morphology-elasticity relation-ships using decreasing fabric information of human trabecular bone
from three major anatomical locations,” Biomech Model Mechanobiol, vol. 12, no. 4, pp. 793–800, 2013.
[19] Y. H. Qian, D. D’Humi`eres, and P. Lallemand, “Lattice BGK models for Navier-Stokes equation,” Europhys Lett, vol. 17, no. 6, pp. 479–484, 1992.
[20] T. Zeiser, M. Bashoor-Zadeh, A. Darabi, and G. Baroud, “Pore-scale analysis of newtonian flow in the explicit geometry of vertebral trabec-ular bones using lattice Boltzmann simulation,” J Eng Med, vol. 222, pp. 185–194, 2008.
[21] R. Singh, P. D. Lee, T. C. Lindley, R. J. Dashwood, E. Ferrie, and T. Imwinkelried, “Characterization of the structure and permeability of titanium foams for spinal fusion devices,” Acta Biomater, vol. 5, no. 1, pp. 477–487, 2009.
[22] A. Syahrom, M. R. Abdul Kadir, J. Abdullah, and A. chsner, “Perme-ability studies of artificial and natural cancellous bone structures,” Med Eng Phys, vol. 35, no. 6, pp. 792–799, 2013.
[23] S. Succi, The lattice-Boltzmann equation for fluid dynamics and beyond. Oxford University Press, 2001.
[24] B. D. Jawerth, P. Lin, and E. Sinzinger, “Lattice Boltzmann models for anisotropic diffusion of images,” J Math Imaging Vis, vol. 11, no. 3, pp. 231–237, 1999.
[25] K. U. Hagenburg, M. Breuß, O. Vogel, J. Weickert, and M. Welk, “A lattice Boltzmann model for rotationally invariant dithering,” in Int Symp Visual Comput (ISVC), ser. LNCS, vol. 5876, 2009, pp. II:949–959. [26] Y. Chen, “A lattice-Boltzmann method for image inpainting,” in Int
Congr Image Signal Process (CISP), vol. 3, 2010, pp. 1222–1225. [27] X. Sun, Z. Wang, and G. Chen, “Parallel active contour with lattice
Boltzmann scheme on modern GPU,” in Int Conf Image Process (ICIP), 2012, pp. 1709–1712.
[28] A. Duda, Z. Koza, and M. Matyka, “Hydraulic tortuosity in arbitrary porous media flow,” Phys Rev E Stat Nonlin Soft Matter Phys, vol. 84, no. 3, pp. 036 319–(1–8), 2011.
[29] E. Klintstr¨om, ¨O. Smedby, R. Moreno, and T. B. Brismar, “Trabecular bone structure parameters from 3D image processing of clinical multi-slice and cone-beam computed tomography data,” Skeletal Radiol, vol. 43, no. 2, pp. 197–204, 2014.
[30] M. Moakher, “On the averaging of symmetric positive-definite tensors,” J Elast, vol. 82, pp. 273–296, 2006.
[31] R. Moreno, M. Borga, and ¨O. Smedby, “Evaluation of the plate-rod model assumption of trabecular bone,” in Proc Int Symp Biomed Imaging (ISBI), 2012, pp. 470–473.
