<figure><img src="image_1.png" /></figure>

HAL
open science

# Segmentation of 3D cellular networks from SR-micro-CT images

Alexandra Pacureanu, Jerome Rollet, Chantal Revol-Muller, Vasile Buzuloiu, Max Langer, Françoise Peyrin

► To cite this version:

Alexandra Pacureanu, Jerome Rollet, Chantal Revol-Muller, Vasile Buzuloiu, Max Langer, et al. Segmentation of 3D cellular networks from SR-micro-CT images. 2011 IEEE International Symposium on Biomedical Imaging: From Nano to Macro, Mar 2011, Chicago, United States. pp.1970-1973, 10.1109/ISBI.2011.5872796 . hal-01825257

HAL Id: hal-01825257

https://hal.science/hal-01825257

Submitted on 10 Apr 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
SEGMENTATION OF 3D CELLULAR NETWORKS FROM SR-MICRO-CT IMAGES

Alexandra Pacureanu¹,², Jérôme Rollet¹, Chantal Revol-Muller¹, Vasile Buzuloiu¹, Max Langer¹,²,
Françoise Peyrin¹,²

¹ CREATIS-LRMN, INSA-Lyon; Inserm, U630; CNRS, UMR5220;
Université de Lyon; Université Lyon 1 ; F-69621, Lyon, France
² ESRF, Imaging Group, 38043 Grenoble Cedex

ABSTRACT

Bone fragility involved in diseases such as osteoporosis implicates many mechanisms at the cellular level. It was recently shown that the lacunar-canalicular network interconnecting osteocytes has a major role in mechanosensitivity. So far, this system has only been studied from 2D microscopic images. In a previous work, we demonstrated the feasibility of synchrotron radiation micro-CT with a voxel size of 0.28µm, to image the lacunar-canalicular porosity in 3D. Nevertheless, the segmentation of this dense network of slender channels with average diameters of ~300-600 nanometers, at the limit of the spatial resolution, is difficult. In this work we propose a level set based method to automatically segment this complex system. To this aim, we designed an automatic initialization process and we apply a post-processing filter. Quantitative results on a ground truth prepared image are presented. On real data sets, expert evaluation showed good results.

***Index Terms*** — 3D image segmentation, level sets, bone, lacunar-canalicular network, X-ray -µCT.

**1. INTRODUCTION**

In appearance rigid and inactive, bone is actually a dynamic tissue
with remarkable properties of self repairing and permanent
remodeling throughout the life. The mechanisms managing these
processes are complex and despite the numerous research studies
and advances realized during the last century, bone remodeling and
failure is not completely elucidated. Understanding how to control
bone resorption and formation at cellular level opens the way
towards a new stage in therapy and diagnosis of bone diseases.
One of the most stringent maladies of our times is osteoporosis,
currently causing fractures in 1 of 3 women and 1 of 5 men over 50
years old, while in the next 40 years the incidence of hip fracture is
expected to increase by 310% in men and 240% in women [IOF].

First considered a passive cell embedded in the bone matrix,
osteocyte is today largely accepted like being the key element in
bone remodeling [1]. By sensing the external mechanical stimuli,
osteocytes trigger bone renewal. These bone cells are hosted in
spaces called lacunae which are regularly spread throughout the
bone matrix. Osteocytes are interconnected via a dense network of
processes hosted in slender (100-600nm in diameter) channels
called canaliculi. Canaliculi enable signal transmission and
recently it was reported [2] that mechanosensitivity is significantly
higher at the level of osteocyte processes comparing to the cell
body, hence the quantity of canaliculi could be correlated with the

<figure><img src="image_3.png" /><figcaption>**Figure 1:** a) A slice in the reconstructed volume from SR-µCT showing Haversian canals, osteocyte lacunae and canaliculi in cortical human bone. b) Detail in the slice from (a); c) Maximum intensity projection over 100 slices, in the region from (b).</figcaption></figure>

tissue capacity to sense mechanical stimuli. Therefore, quantitative
investigation of the lacunar-canalicular network is necessary for
bone quality assessment and for both development and evaluation
of new therapy methods.

Achievement of quantitative results relies on appropriate
imaging techniques and on accurate and robust segmentation
methods. Imaging bone tissue at a sub-micrometric scale keeping
in the same time a sufficiently large field of view to allow a
consistent evaluation is the first difficulty to surmount. Very few
quantitative studies on lacunar-canalicular structure exist in
literature and they were obtained by manual histomorphometry
from optical or confocal microscopy. These imaging methods are
limited to very thin sample slices (~3-30µm), covering only one
layer of cells in depth. In addition, the resolution is not
homogeneous, preventing a valid 3D evaluation.
We demonstrated recently the feasibility of imaging the lacunar-canalicular network with Synchrotron Radiation micro-tomography (SR-µCT) [3][4] at ESRF (ID19) which enables 3D imaging of bone tissue with a homogenous resolution of 0.28µm with a FOV of (573µm)³.

Segmentation of this network is particularly delicate due to the structure complexity and density (Fig.1c), with objects of different geometrical shapes. Contrast and signal to noise ratio are limited by the detector and dose effects on the sample and the slenderness of canaliculi determines partial volume effect as they are only one to two voxels thick. Furthermore the inhomogeneous texture of the bone matrix increases the difficulty to separate the cellular network.

Here we propose a level sets based segmentation method adapted for cellular networks in bone. To surmount the difficulties raised by the segmentation of highly dense, one voxel thick, faintly contrasted, curvilinear structures, we use a shape orientated enhancement filter. Residual noise is removed from the segmented image by filtering of small size connected components. The method is tested on a ground truth phantom constructed from experimental images.

# 2. SEGMENTATION METHOD

## 2.1. Level sets

Let $f(\mathbf{x})$ with $\mathbf{x} = (x, y, z)$ be a 3D image. Level set methods are based on the evolution of a level set function $\phi$ initialized by $\phi(x, y, z, 0) = \phi_0(x, y, z)$:

$$ \frac{\partial \phi}{\partial t} = |\nabla \phi| \cdot V \quad (1) $$

where $V$ is the speed factor controlling the propagation. Several formulations of the speed factor have been proposed and different implementations of the level set methods [5] [6]. In this work we chose the standard approach where the speed factor contains a propagation term accelerating the evolution in homogeneous regions, a regularization term based on the mean curvature and an advection term. The propagation equation can then be written as:

$$ \frac{\partial \phi}{\partial t} = |\nabla \phi| \left( (\lambda_c k - \lambda_p) d_f - \lambda_a \nabla d_f \cdot \frac{\nabla \phi}{|\nabla \phi|} \right) \quad (2) $$

where

$$ d_f = s \circ g_f \quad (3) $$

$$ g_f(x, y, z) = g(|\nabla(G_\sigma * f)(x, y, z)|) \quad (4) $$

$$ s_{\alpha, \beta}(u) = \frac{1}{1 + \exp(-\frac{u - \beta}{\alpha})} \quad (5) $$

where $s$ is a sigmoid function and $G_\sigma(x, y, z)$ is a Gaussian function of standard deviation $\sigma$, and $k$ is the mean curvature defined by:

$$ k = \mathrm{div} \left( \frac{\nabla \phi}{|\nabla \phi|} \right) \quad (6) $$

and $\lambda_a, \lambda_c, \lambda_p$ are parameters controlling the balance between the terms. In order to set the parameters values, we quantitatively evaluated (DICE coefficient, recall, connectivity) the results on the prepared phantom image (Table 1).

<figure><img src="image_4.png" /><figcaption>**Figure 2:** Initialization details in a slice. Arrows indicate points of interest: a) Original image; b) Max entropy threshold on original image; c) Max entropy threshold on filtered image d) Final initialization image</figcaption></figure>

## 2.2 Initialization

Initialization showed to be a crucial part of the segmentation in this case and we designed it as a simple and fast automatic pre-segmentation. Initializing the algorithm with a germs map as starting point for level sets propagation proved to be the only solution to attain good results.

First step consists in finding an optimal threshold to create a base image for propagation. Different methods have been evaluated and we obtained significantly better results by maximizing the inter-class entropy to find the threshold [7]. However, the thresholded image conserves some background noise and small parts from the bone matrix (Fig. 2b, Fig. 3b). In order to remove these residues we erode the image using a cross structuring element of radius 1.

To maximize the likelihood to segment all the canaliculi we seek to put a germ in each of them, especially in the sections where contrast is very weak, thus not appearing in the thresholded image. For this aim we apply a shape orientated filter to enhance curvilinear structures. This filter is based on eigenvalue analysis of the local Hessian matrix and several formulations have been proposed [8] [9]. We applied the criterion proposed in [8] which amplifies the gray level in tubular structures, especially in regions with faint contrast, enabling to reconnect parts of disrupted canaliculi after thresholding (Fig. 2c). On the other hand, the filter removes the lacunae. To obtain the final initialization map we combine the two thresholded images by logical disjunction operation.

## 2.3 Connected component analysis

Ideally, the lacunar-canalicular network is highly connected, so analysis of structure connectivity is an important factor in the evaluation of the segmentation quality. Furthermore, this analysis enables filtering of individual components of chosen sizes, allowing removal of residues in the segmented image, without affecting the size and shape of the other components.

We implemented the technique proposed in [10] providing the number and size of connected components. Isolated elements smaller than 30 voxels (0.6 µm³) were removed. Below we outline the proposed segmentation method in a flowchart diagram and specify the chosen parameters for different steps:
<figure><img src="image_5.png" /></figure>

where

i – Input 3D image

ft – Image after thresholding on shape oriented filter result

it – Image resulted after erosion of threshold(i)

g – Initialization image for level sets

s – Segmented image after level sets

sc – Final segmentation result, after elimination of small connected components

Operators:

F – Shape orientated filter (σ=1, α₁=0.5, α₂=2)

T – Maximum entropy thresholding

E – Binary morphological erosion with a cross element of radius 1

OR – logical operator OR

LS – level sets algorithm based on [6] with

$$ \lambda_a = 0, \lambda_c = 40, \lambda_p = 70, \alpha = -1, \beta = 20 $$

CC – connected component analysis, filtering individual components smaller than 30 voxels.

# 3. SEGMENTATION EVALUATION

## 3.1. Construction of a realistic phantom

For evaluation we prepared a phantom based on a volume extracted from a real image (256³ voxels) in order to conserve the structure characteristics. The binary reference image was created by manual thresholding followed by removal of residues not belonging to the structure and filling up of gaps in canaliculi. Subsequently we set the background gray level at 100, the structure level at 160 and added a Gaussian noise (σ=25), the chosen values being estimated from real images.

## 3.2. Figures of merit

In order to evaluate the segmentation we made use of the DICE coefficient, the recall metric and connectivity error defined as below.

$$ DICE = \frac{2|R \cap S|}{|R| + |S|} \quad \text{recall} = \frac{tp}{tp + fn} \quad \text{Err}_{cc} = \frac{cc_S - cc_R}{cc_R} \quad (7) $$

where R, S are the sets of white voxels in respectively ground truth volume and segmented volume; tp stands for true positive detections, fn is false negative detections; cc_R

and cc_S are the No. of connected components in the reference image and respectively in the segmented image.

# 4. RESULTS

In Figure 2 is illustrated a detail in the experimental image and the different initialization steps, showing that by using the shape orientated filter, gaps in canaliculi are filled up and some low contrasted canaliculi are detected after the structure enhancement.

Table 1 contains quantitative evaluation results of the initialization and the segmentation for different values of the β parameter in sigmoid function (expression (5)). The same quantitative evaluation was applied to set all the level sets parameters presented in the caption of the flowchart diagram. We can observe that after maximum entropy thresholding we keep too

<figure><img src="image_6.png" /><figcaption>**Figure 3:** Results on the phantom (MIP over 10 slices): a) Input image; b) Max entropy threshold on (a); c) Max entropy threshold on (a) filtered image; d) Segmentation result (β=20) overlaid on the input image.</figcaption></figure>

many voxels comparing to the reference ground truth image and there are too many unconnected components (633% error). The results improve considerably in the final initialization image (22% error) and after optimal the level sets segmentation we reduce the error to 3.76%.

The recall metric shows very good results, indicating that we detect almost the integral structure. The slightly lower values of the DICE coefficient reflect the fact that level sets tends to thicken with one voxel on each side the canaliculi, introducing a relative important error in the number of voxels, considering the actual thickness of canaliculi which corresponds to one or two voxels. Computation time for the entire pipeline of the method applied on a (256)³ volume is about 20 minutes. Segmentation result over 10 projected slices in the phantom image are shown in Figure 3d.

## 4.2. Application to experimental data

The method was applied to experimental data obtained with high resolution SR-μCT at ESRF, representing 3D images of the lacunar-canalicular structure in human femur bone tissue. For the sample presented in this paper, the acquisition resolution was

<table><caption>TABLE I
SEGMENTATION EVALUATION</caption><thead><tr><th>Image</th><th>No. CC</th><th>No. Vox</th><th>Err (%) No. CC</th><th>DICE (%)</th><th>Recall (%)</th></tr></thead><tbody><tr><td>Reference</td><td>2634</td><td>287490</td><td></td><td></td><td></td></tr><tr><td>MaxEntrT</td><td>19310</td><td>366307</td><td>633.11</td><td>85.15</td><td>96.83</td></tr><tr><td>Initialization</td><td>3231</td><td>207467</td><td>22.67</td><td>80.75</td><td>69.51</td></tr><tr><td>Segmentation</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>β=1</td><td>2774</td><td>448867</td><td>5.32</td><td>74.55</td><td>95.48</td></tr><tr><td>β=5</td><td>2780</td><td>447852</td><td>5.54</td><td>74.64</td><td>95.46</td></tr><tr><td>β=10</td><td>2809</td><td>444137</td><td>6.64</td><td>74.96</td><td>95.39</td></tr><tr><td>β=15</td><td>2778</td><td>437441</td><td>5.47</td><td>75.48</td><td>95.17</td></tr><tr><td>β=20</td><td>2535</td><td>429529</td><td>-3.76</td><td>75.96</td><td>94.72</td></tr><tr><td>β=50</td><td>1986</td><td>1117321</td><td>-24.60</td><td>38.51</td><td>94.09</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>**Figure 4:** Volume of interest showing the lacunar-canalicular network around Haversian canals, in a human femur cortical bone sample (image width ~ 448 µm). Isosurface of the segmentation achieved with the proposed method.</figcaption></figure>

0.28µm, the energy was set at 20.5keV and 2000 projections were taken with a counting time of 0.8 seconds.

We show below the results obtained for a 1600x1000x256 sub-volume extracted from a (2048)³ image. Figure 3 illustrates the segmentation result of the lacunar-canalicular system, from 3D SR-µCT. For this volume, the initialization image contained 71255 connected components while the segmented image contained only 37225 individual components, showing that the level sets propagation filled gaps between canaliculi. In the result image we could calculate the ratio between the volume occupied by the lacunar-canalicular network and the bone matrix volume, excluding the Haversian canals. We estimated this value at 2.65%.

## 5. CONCLUSION

SR-micro-CT enables now imaging of bone tissue at sub-cellular scale, making possible to visualize and evaluate in 3D the lacunar-canalicular network. The interest to study this structure is increasing rapidly, given the recent results on the importance of osteocytes in bone remodeling and bone diseases. Considering the

number of lacunae and canaliculi, the density and complexity of this network, an automatic method of segmentation and quantification is necessary. Previous studies have been done manually, therefore they were highly time consuming and applied on very small regions of tissue, limiting the consistency of the results.

Here we propose an automatic 3D segmentation method adapted for cellular networks. The correctness of the method was evaluated with several metrics and proved good results, allowing some immediate simple quantification. One of the key points of the method stands in the initialization part based on maximum entropy threshold selection and geometry orientated filtering. Level sets propagation increases the connectivity of the network, reducing the number of individual components by a factor 2. Finally, connected component analysis permits filtering of small residual components belonging to the bone matrix texture. Isosurfaces of the segmentation results are presented and some preliminary quantification values are reported.

## 6. REFERENCES

[1] E.H. Burger and J. Klein-Nulend, "Mechanotransduction in bone - role of the lacuno-canalicular network," *FASEB J.*, vol. 13, May. 1999, pp. 101-112.

[2] T. Adachi, Y. Aonuma, M. Tanaka, M. Hojo, T. Takano-Yamamoto, and H. Kamioka, "Calcium response in single osteocytes to locally applied mechanical stimulus: Differences in cell process and cell body," *Journal of Biomechanics*, vol. 42, Aug. 2009, pp. 1989-1995.

[3] A. Pacureanu, A. Larrue, Z. Peter, and F. Peyrin, "3D non-linear enhancement of tubular microscopic bone porosities," *Biomedical Imaging: From Nano to Macro*, 2009. ISBI '09. IEEE International Symposium on*, 2009, pp. 602-605.

[4] M. Salomé, F. Peyrin, P. Cloetens, C. Odet, A.M. Laval-Jeantet, J. Baruchel, and P. Spanne, "A synchrotron radiation microtomography system for the analysis of trabecular bone samples," *Medical Physics*, vol. 26, Oct. 1999, pp. 2194-2204.

[5] J.A. Sethian, *Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science*, Cambridge University Press, 2000.

[6] V. Caselles, R. Kimmel, and G. Sapiro, "Geodesic active contours," *Computer Vision, IEEE International Conference on*, Los Alamitos, CA, USA: IEEE Computer Society, 1995, p. 694.

[7] P.K. Sahoo, S. Soltani, and A.K.C. Wong, "A survey of thresholding techniques," *Computer Vision, Graphics, and Image Processing*, vol. 41, Feb. 1988, pp. 233-260.

[8] Y. Sato, S. Nakajima, N. Shiraga, H. Atsumi, S. Yoshida, T. Koller, G. Gerig, and R. Kikinis, "Three-dimensional multi-scale line filter for segmentation and visualization of curvilinear structures in medical images," *Medical Image Analysis*, vol. 2, Jun. 1998, pp. 143-168.

[9] A. Frangi, W. Niessen, K. Vincken, and M. Viergever, "Multiscale Vessel Enhancement Filtering," *Medical Image Computing and Computer-Assisted Intervention — MICCAI'98*, 1998, p. 130.

[10] J. Hoshen and R. Kopelman, "Percolation and cluster distribution. I. Cluster multiple labeling technique and critical concentration algorithm," *Physical Review B*, vol. 14, Oct. 1976, p. 3438.