2.1.

Bone Segmentation

Since bone structures exhibit high contrast on CT images, they can be easily identified with a conventional application of a threshold, which conveniently is also a fast operation. To automatically assign a threshold, we examined the histogram of the intensities of the CT volume data. When dealing with a normalized scale like the Hounsfield scale, it is straightforward to select a certain threshold that divides the image in bone and background. However, in this work we make no assumption on the CT data scaling so that the method can work seamlessly with different CT acquisition parameters and datasets, since in small animal imaging there exists less standardization between the data obtained, compared to clinical data. The analysis of many histograms of our CT data had shown that there are no significant features representing the intensity of bone tissue, like a local maxima or minima, which could be traced. The intensity of the bones is usually widely distributed throughout the histogram. Therefore we approximated our threshold $T_b$ by finding other distinct intensities, and assume that there is a linear relationship between those intensities and the threshold we need. Mathematically, this can be described by

When considering a typical histogram of a mouse CT, two distinct peaks can be noticed that could be used as reference points (Fig. 1 ). The highest peak can be found at the left side of the histogram and corresponds to voxels of very low density, in this case primarily the voxels corresponding to the air in the field of view surrounding the animal. A second significant peak, corresponding to water, can also be easily identified as a second maximum in the histogram. Soft tissue contains high amounts of water, and the area around that peak essentially indicates voxels corresponding to soft tissue. These two peaks can be employed for approximating the optimal threshold in Eq. 1. To determine the scaling factor $w$ in Eq. 1, we considered the Hounsfield scale, since it is a common standard for CT images. In the Hounsfield scale, air has an intensity of $-1000$ Hounsfield units (HU), water has an intensity of $0\mathrm{HU}$ , and bone structures start at $400\mathrm{HU}$ , i.e., 0.4 of the water-air difference. Thus in the Hounsfield scale, the threshold for bone structures $T_b$ at $400\mathrm{HU}$ can be determined by rewriting Eq. 1 as

Fig. 1

From the CT volume data, a histogram of the intensities was computed. In this example you can recognize the two distinct peaks (arrows) that represent soft tissue, consisting mostly of water and air.

2.2.

Ribcage Detection

For segmentation of lung and heart, we considered first the identification of orientation points to serve as initial points for subsequent segmentation steps. In this role, the ribcage serves as an easily identifiable structure that accurately delineates a big part of the outer surfaces of the lung and heart. To identify the ribcage, we analyzed the result of the bone segmentation by computing a histogram of the number of segmented bone voxels per axial slice. We treated the histogram as a signal $h_a\left(t\right)$ , where $t$ is the slice number. Within this signal the ribcage creates a distinct harmonic frequency [Fig. 2 ]. To detect the periodicity, we opted for the use of a Gabor filter; essentially a Gaussian function multiplied with a cosine, i.e.,

Fig. 2

(a) The histogram function $h_a\left(t\right)$ displays the number of segmented bone voxels per slice. Note the very harmonious frequency produced by the ribcage between the dotted lines. The lines mark the beginning and the end of the sternum. (b) The response of the Gabor filter has its highest peak close to the center of the ribcage. The solid line connecting all peaks is the interpolation of all maxima of the filter response.

This approach is similar to template matching, where the Gabor function describes the periodic oscillation of the ribcage. Since the frequency of the Gabor filter needs to match the frequency the ribcage produces in $h_a\left(t\right)$ , specific values for $\sigma$ and $\lambda$ had to be defined. To determine those values, we analyzed the frequency produced by the ribcage in three training datasets and adjusted $\sigma$ and $\lambda$ so that the Gabor filter fitted this frequency. When performing the ribcage detection on unknown test data, we used the differences of the voxel spacing and voxel size in the training and test data sets to compute a scaling factor for $\sigma$ and $\lambda$ . Thus the procedure is independent from scaled image data. Inherent in this procedure is the assumption that the size of the imaged mice does not vary significantly, and that the ribs have a distinct separation to each other. To apply the filter, we convoluted the histogram signal $h_a\left(t\right)$ with the Gabor filter $g_{\sigma ,\lambda }\left(t\right)$ . The result of the convolution $h_a{\left(t\right)}^{*}$ $g_{\sigma ,\lambda }\left(t\right)$ is a filter response that usually had its global maximum near the axial center of the ribcage [Fig. 2]. Furthermore, we used just the local maxima of the filter response to interpolate a new function [also Fig. 2]. In this function, the next minima to the left and right sides of the global maximum (the center of the ribcage) were selected as landmark points that defined a bounding box in the axial direction around the ribcage. Those landmarks do not necessarily mark specific anatomical points, but usually they appear near the top and bottom endings of the sternum.

To define the bounding box also in sagittal and coronal directions, we computed the histograms $h_s\left(t\right)$ and $h_c\left(t\right)$ indicative of the number of segmented bone voxels in these directions. Note that we only used the slices between the axial landmark points to compute those histograms. In the histograms, we searched for the global maxima and the first slices left and right to them, where $h_s\left(t\right)$ and $h_c\left(t\right)$ respectively equal zero, that is, were the ribcage ends. We confined our bounding box only around the ribcage and excluded artifacts outside the mouse that sometimes occur. If this detection scheme experiences difficulties due to noise and artifacts, it is possible to employ searches that define areas where $h_s\left(t\right)$ and $h_c\left(t\right)$ become smaller than a predetermined value greater than zero to get a bounding box tight around the ribcage. Overall, ribcage determination is an essential step for further detection of anatomic structures, as described in the following.

2.3.

Lung Segmentation

To enable lung segmentation, we utilized a seed growing algorithm within the confines of the detected ribcage. For this purpose, possible seed points inside the lungs had to be found automatically. The whole respiratory system can be naturally recognized on XCT images by means of its low density and the corresponding high contrast to surrounding tissue. However, image intensity and contrast alone did not suffice for accurate detection. This was because the bounding box of the seed growing algorithm was usually still wide enough to occasionally contain regions of air outside the mouse or parts of the digestive system that showed also a low intensity and contrast. Other challenges involved the blurring of borders due to possible moving artifacts. To refine the region of interest and achieve a correct segmentation result, we created a spherical region of interest (SROI) inside the bounding box. The SROI was initialized as a sphere with a radius of zero at the center of the bounding box [see Fig. 3 ], and was allowed to grow so that the sphere radius measured 90% of the distance between the center of the sphere and the boundary of the mouse. In Fig. 3 a CT image is displayed with the according bounding box (solid blue line) and the SROI (bright circle).

Fig. 3

(a) The image shows a CT slice including the SROI (brighter spherical region), the computed bounding box (solid line), and possible seed points (white spots). The dotted line marks the middle of the bounding box, dividing right and left lobes of the lung. (b) The graph shows the intensity distribution within the region of interest. The peak represents water in the soft tissue.

To find seed points inside the SROI, we computed an intensity histogram from all the voxels inside the SROI [Fig. 3]. Here, the voxels with the lowest intensity $I_{\text{low}}$ mark the dark bronchial tubes, and the high peak $I_{\text{peak}}$ marks soft tissue. We took these easy to detect points as references to compute an interval $[I_1,I_2]$ , where

For the seed growing algorithm itself, a mean intensity ${\overline{I}}_a$ was defined using the chosen seed point $s$ and its neighboring voxels. Also, a confidence interval was defined by

Because we also wanted to be sure that both the right and left lobes of lung are segmented, we chose an equal number of seed points from both. We distinguished between the right and left lobe by simply dividing our bounding box in the middle.

2.4.

Heart Position Approximation and Segmentation

The last procedure of our framework is the approximation of the heart position and its segmentation. For this purpose we propose the use of a shape model of the heart generated from manually segmented training data. In this work we used one manually segmented volume dataset to gain this model. The model is a closed mesh that consists of a number of vertices connected through edges. To yield a rough initial position and scaling factor for the model, we used the bounding box from the ribcage detection as reference. Using the training dataset, we examined the heart position relative to the borders of the bounding box by considering the box as a normalized cube with a side length of 1. We initialized the heart model at the same position in bounding boxes of other CT volume images. A scaling factor for controlling the size of the heart model was approximated using the sizes of the bounding boxes around the ribcages of training and test datasets as references. Scaling factors were computed for all three directions and averaged to get the main scaling factor. This averaging results in a more robust scaling, for our experiments had shown that the sizes of the bounding boxes varied enough to receive unusual heart shapes.

This operation generally placed the heart model close to its supposed position. However, it also attained regions where the heart model was overlapping the other segmentations of the lung and bone structures. This is because of the rough initial position and scale approximation. To adjust the model position, we searched for all of the overlapping voxels and created for each one a unit vector that points to the center of gravity of the heart model. Thus a vector field was created. The field represents forces that push the model away from overlapping sections. After the heart model was translated by the vector field, the procedure was repeated iteratively. A decreasing weighting factor thereby ensures the convergence of the procedure. The iterative process was stopped when either no more voxels with segmentation overlap were detected, or the translational improvement was beneath a specified threshold. The latter usually occurs when there is a balance between forces from opposite sides, i.e., lung and ribcage/sternum, which means that the heart model is too large to fit. We then scaled down the heart model to 95% of its size and restarted the iterative position adjustment until finally no more regions with overlapping segmentations remained.

We note that this algorithm does not provide a segmentation of the heart that fully incorporates wide shape variations. Since the heart model is static, it cannot fully fit the actual image data. Nevertheless, it still can be used as an approximation of a segmentation result and as a new initial position for further, more advanced segmentation algorithms like active contour models that have yet to be implemented.

2.5.

Validation of Segmentation Results

As a reference for evaluation of segmentation results, we used gold-standard manual segmentation revised by an expert specialized in mouse anatomy. We segmented the whole skeleton, both lobes of the lung, and the heart of a CT volume image with a size of $267\times 242\times 452$ voxels on a $64\text{-bit}$ PC with a quad core CPU $\left(2.67\mathrm{GHz}\right)$ and $4\mathrm{GB}$ of RAM. Notice that results of the bone segmentation will always be constant. The lung segmentation algorithm, on the other hand, picks randomly only a few of many possible seed points, thus producing different results. Since the heart position approximation depends on the lung segmentation result, these results vary too. To compensate for this fact, we performed the segmentation process several times to yield the mean performance.

As a main criterion for the evaluation, we used the Dice coefficient $0⩽s⩽1$ with

4.1.

Segmentation

Figure 4 shows the empiric results of the bone segmentation. Notice that very thin bone structures like the blade bones exhibit holes. In our CT images, these structures show lower intensities than bone usually does due to blurring artifacts. Overall, the results yielded Dice coefficients of 0.8721. FRR (0.1062) and FAR (0.1485), which show that these operations resulted in oversegmentation. When we visually evaluated the result, we recognized that nearly all segmentation errors occurred along the borders. This is mainly due to blurring artifacts at the borders between different tissues. Thus we considered these errors to be within normal uncertainty bounds. The segmentation of the bones took $3.3\sec$ , and the recognition of the ribcage took $1.2\sec$ , which is very fast for data volumes of such large size.

Fig. 4

The result of typical bone segmentation: (a) the original CT slice, (b) the corresponding slice of the segmented data, and (c) surface model of the skeleton computed from the segmentation result.

For the lung segmentation, we analyzed 30 segmentations of our reference image data. We experienced that five seed points per lobe of the lung usually were enough to achieve an accurate and robust segmentation result while still being time efficient. The results are displayed in Fig. 5 . The framework achieved a mean Dice coefficient of 0.766 with a variance of 0.007. Nonsegmented voxels (FRR 0.3096) had the greatest influence on this result, while oversegmentation was much smaller (FAR 0.1091). Falsely accepted voxels were usually part of the bronchial tubes outside the lung. The falsely rejected voxels were mostly voxels with a considerably higher intensity, where the lung tissue showed pathologies. The speed of the lung segmentation differs, since the number of iterations of the seed growing algorithm depends on the initial seed point. Usually the segmentation was done in less than $30\sec$ , including the search for appropriate seed points.

Fig. 5

Result of the lung segmentation: (a) the original CT slice, (b) the corresponding slice of the segmented data, and (c) surface model of the skeleton and the lung computed from the segmentation results.

The heart segmentation was also done 30 times. In Fig. 6 you can see the adapted heart model inside the ribcage. The mean Dice coefficient was 0.7647 and had a variance of only 0.0004. Considering that we only used a static model build from one single training dataset, we consider this a very good result. Most notably, this result was due to the quite high FRR of 0.3378, while only 0.0936 of the segmented voxels were falsely accepted. The time needed for the approximation of the heart also depends on the initial position. On our test data it took less than $45\sec$ to adjust the heart model.

Fig. 6

Result of the heart segmentation. The image shows surface models of all three segmented anatomical structures.

4.2.

Reconstruction

Figure 7 shows the results of the utilization of XCT anatomical information as priors in an FMT inversion scheme. The FMT images are laid over the corresponding CT slice. To simplify matters, only one slice out of a reconstructed volume is presented for each approach. For the evaluation of the reconstruction improvement using anatomical priors, we simulated a situation of inflamed lungs [Fig. 7], modeled after previous studies of lung inflammation,^{2, 34} and used three different reconstruction procedures, i.e., 1. no regularization, 2. inversion using Thikonov regularization, and 3. the Laplace regularization. However, segmentation of the in-vivo CT imaging data was done using our framework, and no simulated segmentation was used. We note that the first two approaches do not utilize the segmented information image priors, and that no noise was added in the simulation.

Fig. 7

(a) Simulated fluorescence signal in the lung. (b) Result without regularization. (c) Tikhonov regularization. (d) Laplace regularization showing the best imaging performance in this case.

Figure 7 shows the inversion obtained without regularization. In this case the inversion generates significant artifacts, especially on the borders leading to a highly inaccurate reconstruction. Figure 7 depicts high blurring of the fluorescent signal. The intensity of the signal is also too low, and a prominent spot can be recognized in one lobe of the lung while the intensity should be homogeneous. Finally, Fig. 7 shows the best reconstruction results due to the priors. The fluorescence intensity was reconstructed accurately; it is distributed homogeneously in the lung and only small blurring artifacts occur along the borders.