**6. Automatic image processing, topology analysis and measurement of statistical features**

This chapter describes an approach to automatic micro-CT image processing using computer vision techniques and the Mercox vascular corrosion casts of the intestinal mucosa described in the previous chapter. The input data were acquired with Xradia XCT 400 (Xradia, Pleasanton, CA, USA) by scanning Mercox vascular corrosion casts. The resulting data set consisted of a DICOM stack of approximately 1000 slices, each 1000 times 1000 pixels, such that the scanned volume contained 109 volume elements (voxels), each represented by a 16-bit signed number. The depicted volume in reality was a cube with an edge of approximately 35 mm, whereas the voxel edge was 36.53 µm. For image processing, we used the Insight Toolkit library for C++ (www.itk.org) and MATLAB with its Image Processing Toolbox (The MathWorks, Inc., Natick, MA, USA). The first step of the image processing was an analysis of the density distribution of the data, using a density histogram; see Fig. 4.

Segmentation is a process in computer vision that divides the voxels of the source image into two subsets – foreground and background – depending on their affinity to the objects of the real world (Sonka et al., 1998). In our case, the foreground was formed by the vascular corrosion cast, whereas all of the other data were considered to be background.

Thresholding is a common and widely used image segmentation method in machine vision. This technique outputs a binary image *g(i,j,k)* that classifies the voxels of the source image *f(i,j,k),* in which *i, j, k* represent the spatial indices of the data, according to equation 2:

$$\begin{aligned} \text{sg}(\mathbf{i}, \mathbf{j}, \mathbf{k}) &= \mathbf{1} \quad \text{for } f(\mathbf{i}, \mathbf{j}, \mathbf{k}) \ge T \\ &= \mathbf{0} \quad \text{for } f(\mathbf{i}, \mathbf{j}, \mathbf{k}) < T \end{aligned} \tag{2}$$

In this equation, *T* is the threshold value; *g(i,j,k)*=1 for image elements of objects and *g(i,j,k)*=0 for image elements of the background (Sonka et al., 1998). The threshold value *32,768*, used for segmentation in our case, is highlighted in red in Fig 4.

Correlating Micro-CT Imaging with Quantitative Histology 183

**Figure 5.** A label image histogram. The object labels each stand for voxels belonging to a single

**Figure 6.** A three-dimensional reconstruction of the entire vessel system (A), and of the largest

suitable for large data sets and produces a 1-voxel-thin 3-D skeleton. The topology analysis of the skeleton involves identifying the node and terminal points of the tree. The points are located based on their 26-connectivity to neighboring voxels. Based on the properties of the skeletonization algorithm, we can find nodes of valence 3 (i.e., nodes at which three microvascular segments join). Existing nodes of valence 4 are decomposed into two nodes of valence 3. Fig. 7 shows the resulting skeleton with located node points. In this sample, 149 nodes were detected. Using the number of branching points with a known valence, the numerical density of the microvessels *NV(cap/ref)* can be estimated, as described by

under study.

continuous object (B).

Lokkegaard et al. (2001).

contiguous object, i.e., the count of each label denotes the size of the object it represents, and the number of unique labels corresponds to the number of contiguous objects in the volume. Object labels higher than 100 were omitted for a better overview. There were 1,967 unique objects within the image data set

**Figure 4.** A density histogram of the source image data. In this histogram, which has been smoothened for visualization, we can observe a local minimum prior to a peak representing the density of desired vessel voxels. According to the histogram, the segmentation methods based on density values appear to be suitable for the current image data set.

After segmentation was performed, it was necessary to further pre-process the binary image to eliminate artifacts, e.g., holes and bays caused by irregular spreading of the casting material, and measurement errors. These errors were resolved by applying the morphological operation of hole filling, which converts cavities inside of objects into object voxels. The next step of the data analysis used the labeling algorithm, which assigned a unique label to each contiguous region of a binary image (Sonka et al., 1998). This procedure gave us a large amount of useful information about object counts and sizes in the volume. Part of the label image histogram is shown in Fig. 5.

The label histogram (Fig. 5) demonstrates that there are many objects of insignificant size that are produced simply by measurement noise. We can convert these objects into background and eliminate them from further processing. The size of the objects at which the segmentation of the blood vessel regions becomes unreliable, reveals also the resolution limits of the current micro-CT scan. Acquiring a reliable representation of all of the microvessels would therefore require another scan with a different objective and micro-CT settings. We can also deduct the count of dominant objects in the volume. For the purposes of visualization, in the subsequent processing steps, we will only work with the largest object. Fig. 6 demonstrates the difference between a 3-D model of the whole vessel system (Fig. 6A), and the largest continuous object only (Fig. 6B).

The volumetric image (Fig. 6) is useful for visualization. However, it is not desirable for topological investigation. A shape simplification is needed that would preserve the tree connectivity and geometric conditions. In this study, we used the algorithm introduced by Lee et al. (1994). This algorithm is based on parallel 3-D binary image thinning, which is

182 Injury and Skeletal Biomechanics

be suitable for the current image data set.

Part of the label image histogram is shown in Fig. 5.

(Fig. 6A), and the largest continuous object only (Fig. 6B).

In this equation, *T* is the threshold value; *g(i,j,k)*=1 for image elements of objects and *g(i,j,k)*=0 for image elements of the background (Sonka et al., 1998). The threshold value *32,768*, used

**Figure 4.** A density histogram of the source image data. In this histogram, which has been smoothened for visualization, we can observe a local minimum prior to a peak representing the density of desired vessel voxels. According to the histogram, the segmentation methods based on density values appear to

After segmentation was performed, it was necessary to further pre-process the binary image to eliminate artifacts, e.g., holes and bays caused by irregular spreading of the casting material, and measurement errors. These errors were resolved by applying the morphological operation of hole filling, which converts cavities inside of objects into object voxels. The next step of the data analysis used the labeling algorithm, which assigned a unique label to each contiguous region of a binary image (Sonka et al., 1998). This procedure gave us a large amount of useful information about object counts and sizes in the volume.

The label histogram (Fig. 5) demonstrates that there are many objects of insignificant size that are produced simply by measurement noise. We can convert these objects into background and eliminate them from further processing. The size of the objects at which the segmentation of the blood vessel regions becomes unreliable, reveals also the resolution limits of the current micro-CT scan. Acquiring a reliable representation of all of the microvessels would therefore require another scan with a different objective and micro-CT settings. We can also deduct the count of dominant objects in the volume. For the purposes of visualization, in the subsequent processing steps, we will only work with the largest object. Fig. 6 demonstrates the difference between a 3-D model of the whole vessel system

The volumetric image (Fig. 6) is useful for visualization. However, it is not desirable for topological investigation. A shape simplification is needed that would preserve the tree connectivity and geometric conditions. In this study, we used the algorithm introduced by Lee et al. (1994). This algorithm is based on parallel 3-D binary image thinning, which is

for segmentation in our case, is highlighted in red in Fig 4.

**Figure 5.** A label image histogram. The object labels each stand for voxels belonging to a single contiguous object, i.e., the count of each label denotes the size of the object it represents, and the number of unique labels corresponds to the number of contiguous objects in the volume. Object labels higher than 100 were omitted for a better overview. There were 1,967 unique objects within the image data set under study.

**Figure 6.** A three-dimensional reconstruction of the entire vessel system (A), and of the largest continuous object (B).

suitable for large data sets and produces a 1-voxel-thin 3-D skeleton. The topology analysis of the skeleton involves identifying the node and terminal points of the tree. The points are located based on their 26-connectivity to neighboring voxels. Based on the properties of the skeletonization algorithm, we can find nodes of valence 3 (i.e., nodes at which three microvascular segments join). Existing nodes of valence 4 are decomposed into two nodes of valence 3. Fig. 7 shows the resulting skeleton with located node points. In this sample, 149 nodes were detected. Using the number of branching points with a known valence, the numerical density of the microvessels *NV(cap/ref)* can be estimated, as described by Lokkegaard et al. (2001).

Correlating Micro-CT Imaging with Quantitative Histology 185

Let us summarize the advantages and disadvantages of our automatic image processing approach. Compared with stereological methods currently in use (see the next chapter), automatic image processing requires virtually no user interaction. This approach is able to compute the distributions of the volumes, surfaces and lengths of vessels of given diameters with a high precision. Moreover, the stereological methods may be applied to our volumetric model to achieve results comparable to those of stereological measurements performed by a human, but with no interaction at all. A disadvantage is the necessity of the choice of the threshold used for the segmentation, due to the sensitivity of certain results to the threshold settings. As stated before, the threshold value does not greatly influence the topology of the vessel tree but strongly correlates with the diameters of the vessels. A proposed solution to this issue is to compute the skeleton with a higher threshold and the

Current micro-CT devices are typically bundled with sophisticated software packages that offer a number of automated quantification procedures. However, correlating the micro-CT results with quantitative histology favors the use of unbiased stereological methods, which are highly standardized and widely accepted in biomedical microscopy research (Howard and Reed, 1998; Mouton, 2002). This chapter illustrates the stereological assessment of micro-CT scans of bone scaffolds and microvascular corrosion casts, including the quantification of the volume fraction (*VV*, dimensionless), surface density (*SV*, m-1), length

In bone tissue samples, the micro-CT resolution is currently capable of providing images that can be used for both analysis of bone vascular canals, and counting individual osteocyte lacunae. Quantification of bone microporosities is used for testing their effect on the viscoelastic properties of bone tissue (Tonar et al., 2011). The microporosity has at least two functional levels, the vascular porosity (related to the vascular canals; the order of magnitude is 10–1000 µm) and the lacunar-canalicular porosity (surrounding the osteocytes; the order of magnitude is 0.1–10 µm). Mechanical experiments clearly demonstrated that the hierarchical organization of bone architecture is crucial and that bone structural heterogeneity varies with the scale of magnification. Whereas Fig. 8 demonstrates twodimensional sections of compact bone produced by histology and micro-CT, Fig. 9 shows

The vascular corrosion casts described in section 4 and used in section 5 also can be assessed with spatial stereological methods. Image data acquired by micro-CT are demonstrated in

After manually tracing the microvessel profiles within a series of consecutive twodimensional micro-CT sections (software Ellipse, ViDiTo, Košice, Slovak Republic), a three-

density (*LV*, m-2), orientation and anisotropy of microvessels (Kochová et al., 2011).

diameters with a lower one.

Fig. 10.

**7. Quantitative micro-CT, histology and stereology** 

the results of 3-D micro-CT reconstruction of cancellous bone.

dimensional system of oriented lines can be acquired (Fig. 11).

**Figure 7.** Vessel skeleton with branching node points.

Currently, we are able to segment a vessel tree formed by a corrosive cast. In the example analyzed in this study, continuous regions of the volumetric model were counted to detect the number of individual objects within the volume. The volume of each individual object was estimated. Using skeletonization, the number of branching nodes was counted, and the number of vascular segments between the nodes was estimated, as was the length of each segment. Knowing the vessel length and volume, it is also possible to compute the average diameter of the vessel. However, there are still challenges for future work, such as estimating the vessel diameter in each voxel of the skeleton using a distance transform. With this information, the diameter distribution with respect to the vessel length or the vessel volume could be acquired. With a known diameter, the surface area of the vessels can also be computed. Tracking the skeleton leading its binary graph construction determines the spatial distribution of the branching nodes from the proximal vascular segments to the periphery of the vascular tree.

The approach presented in this chapter is based on several assumptions that deserve to be discussed, as they also represent the limitations of this study. We assume, for our computation, that the vessels are in the form of generalized cylinders, which means that the cross-section orthogonal to the medial axis is a circle. The voxels in our data set were cubical; however, a data set with unequal voxel edges may be resampled into cubical voxels.

An important question is the choice of the threshold value. The segmented vessel diameter is partially linearly dependent on the chosen threshold value. Setting an excessively low threshold causes too many artifacts to appear and objects to merge into each other; contrariwise, an excessively high threshold makes vessels very thin, such that small vessels disappear. This effect implies that the threshold value influences the absolute diameter of the vessels and biases statistical markers. However, the construction of the skeleton and the topological analysis do not appear to be affected by the selection of the threshold.

Let us summarize the advantages and disadvantages of our automatic image processing approach. Compared with stereological methods currently in use (see the next chapter), automatic image processing requires virtually no user interaction. This approach is able to compute the distributions of the volumes, surfaces and lengths of vessels of given diameters with a high precision. Moreover, the stereological methods may be applied to our volumetric model to achieve results comparable to those of stereological measurements performed by a human, but with no interaction at all. A disadvantage is the necessity of the choice of the threshold used for the segmentation, due to the sensitivity of certain results to the threshold settings. As stated before, the threshold value does not greatly influence the topology of the vessel tree but strongly correlates with the diameters of the vessels. A proposed solution to this issue is to compute the skeleton with a higher threshold and the diameters with a lower one.
