ABSTRACT

Brain ideal midline estimation is vital in medical image processing, especially in analyzing the severity of a brain injury in clinical environments. We propose an automated computer-aided ideal midline estimation system with a two-step process. First, a CT Slice Selection Algorithm (SSA) can automatically select an appropriate subset of slices from a large number of raw CT images using the skull’s anatomical features. Next, an ideal midline detection is implemented on the selected subset of slices. An exhaustive symmetric position search is performed based on the anatomical features in the detection. In order to enhance the accuracy of the detection, a global rotation assumption is applied to determine the ideal midline by fully considering the connection between slices. Experimental results of the multi-stage algorithm were assessed on 3313 CT slices of 70 patients. The accuracy of the proposed system is 96.9%, which makes it viable for use under clinical settings.

1. Introduction

Human brain has two hemispheres with an approximate bilateral symmetry distinguished by the ideal midline (IML), which is the longitudinal fissure marked by the falx cerebri in the mid-sagittal plane [1,2]. The computer-aided estimation of IML has attracted a great deal of attention in the recent two decades [1,3,4]. The interhemispheric fissure line segments have been widely used to detect the ideal midline on MRI images which usually has a high visibility on the fissure line [5,6]. In the case of brain CT slices, longitudinal fissure cannot be used as a primary index for detection due to the low-to-zero visibility of the fissure which can seriously affect the accuracy of detecting the Mid-Sagittal Plane (MSP) or IML. Moreover, some pathological symptoms, such as a tumor, may curve the fissure and completely change the direction of the fissure. To avoid the above limitation, skull symmetry has been included as another important anatomical feature in MSP/IML detection [6,7]. G. Ruppert et al. extracted the MSP based on bilateral symmetry maximization [6]. W. Chen et al. combined bone symmetry and direct detection of the anatomical features in CT images in IML detection [8]. This method works effectively and accurately on a single CT slice but lacks connection or comparison with the detection results from other CT slices of the same patient. For a computer aided medical image processing system to detect IML, the capability of automatically selecting relevant CT slices is essential. Although dozens of CT images can be acquired from a patient’s brain scan, only those slices depicting clear anatomical features and limited inherent noise are used in the IML quantification. Currently in clinical practice, the process of selecting appropriate images for IML diagnosis is performed manually by physicians [7,8]. From our research, we did not find any existing automated method that is in implementation to perform this task.

In this work, we propose a two-step algorithm for the automated IML detection. As the first step, a CT Slice Selection Algorithm (SSA) is proposed, wherein the algorithm finds an appropriate subset of slices from a large number of raw CT images. In the second step, brain IML is detected accurately and effectively by considering both anatomical features and the connections among CT slices. A database of 3133 CT slices of 70 patients with TBI cases yields highly desirable accuracy and efficiency when tested with the proposed method.

2. Methodology

The flowchart of the ideal midline (IML) detection system is shown in Figure 1. The two steps, CT Slice Selection Algorithm (SSA) and Ideal midline detection (IML detection) comprise the core of the system. SSA is used to greatly reduce the slice number while the IML detection is aimed to accurately detect the IML position and rotation angle. The following subsection 2.1 and 2.2 will describe the details of each step.

2.1. CT Slice Selection Algorithm

With head CT scan in the clinical environment, dozens of raw images can easily be acquired for one patient. However, not all images are ideal for IML detection. As shown in Figure 2, there are 42 raw CT images obtained from one patient’s CT scan. Some images taken from the lower section of the head contains too much noise from other organs, such as the eye and nose in slice No.15. Some images capture a small intracranial area because the scan position is too close to the calvaria, as seen in slice No.36. Some images capture integrated skull contours and large intracranial area but lack good convexity there rendering them improper for IML detection, such as slice No.19. From the viewpoint of anatomical features, the ideal CT slices usually contain larger intracranial area, integrated skull bone contour and good convexity of the skull, such as No.22 through No.30 in Figure 2. Therefore, CT slice selection should ideally be based on the above mentioned features.

To effectively select a few appropriate CT slices from a large number of CT scan images acquired for each patient, the CT Slice Selection Algorithm (SSA) was designed. As the flowchart shows in Figure 3, this algorithm analyzes every slice by examining multiple anatomic features. Each as aspect of the flowchart has been described in the following sections.

2.1.1. Skull Detection

As the first step in SSA algorithm, the skull detection is firstly implemented on every raw CT slice as shown in Figure 4. In this step the raw CT images are treated as a raw matrix I with m rows and n columns (Equation (1)).

(1)

where I_i,j represents the intensity of the pixel at the ith row and jth column. Using a threshold method, the possible bone pixels can be extracted from the raw matrix to build up a new binary matrix B as shown in Figure 4(b).

(2)

where the pixels with their original intensity I_i,j is larger than the threshold of T. In this study, based on experimentation, the value for T is set to 250, which lies within the common rage for bone intensity within CT images.) (See Equation (3)).

(3)

where B_i,j is the element at ith row and jth column in the new binary matrix B. Then those possible bone pixels constitute a certain number of connected regions C₁, C₂, ∙∙∙, C_P by means of the connected component algorithm (CCA) [9]. We choose the ath connected region C_a (1 ≤ a ≤ p) which contains the largest number of elements as the candidate skull as shown in Equation (4).

(4)

where C_k (k = 1, ∙∙∙, p) is the kth connected region and f(C_k) is the number of elements in region C_k. Next, all the other connected regions from the image are removed except for the candidate skull. However, some small holes still possibly exist in the candidate skull Ca as shown in Figure 4(c). To remove those small holes in bone, the binary matrix is copied and inverted to form a new matrix, called the inverted matrix M (Equations (5) and (6),

(5)

with

(6)

where M_i,j is the converted intensity of the pixel at the ith row and the jth column in the inverted matrix M. Using the connected component algorithm (CCA) again, q pieces connected regions (D₁, D₂, ∙∙∙, D_q) are obtained from the inverted matrix M. Using the identified region seach of the component of D_k. L_k (k = 1, 2, ∙∙∙, q) is used to represent the pixels within the matrix for each component of D_k.

(7)

with

(8)

After finding each of the components which does notbelong to the skull, the area of these components is computed. Using the computed areas, only those connected regions with an area less than a set threshold S is considered as a hole within the bone structure of the original scan. For this study based on the relative sizes of the objects found in brain CT scans, S has been set to 200 pixels which is a fair estimate of possible hole size. Once these holes have been identified inside the candidate skull structure, they are filled with bone intensity (equal to 1 in binary the matrix). This helps unify the overall identified bone structure by covering all the holes. A subset of the inverted matrices which are identified as holes within the bone structure is given as H_k (k = 1, 2, ∙∙∙, q) to express the bone hole regions.

(9)

where f(D_k) is the number of elements in region D_k and O is a zero m × n matrix. Then, we can obtain the final detected skull G by combining the candidate skull matrix (J − M) with all whitened small holes matrices L_k as shown in Equations (7), (8) and Figure 4(d).

(10)

where M is the converted matrix defined in Equation (6), J is the matrix with every element equal to one. Thus, (J − M) represents the candidate skull corresponding to the connected region C_a.

2.1.2. Closed Skull Inspection

Followed by the skull detection, the second step in SSA algorithm is the closed skull inspection. This process aims to remove the slices with either unclosed skull or with the skull containing too many separated regions. The non-integrated skull affects the following IML identification since symmetry value calculation through the exhaustive symmetric position search process is sensitive to the skull contour.

We define a new measure, called skull closing level F, using the number of zero matrices among all hole matrices H_k (k = 1, 2, ∙∙∙, q).

(11)

with

(12)

where H_k is the kth hole matrix which was defined in Equation (9), O is a zero m × n matrix, and thus g(H_k) is a binary state-variable used to express whether H_k is a zero matrix. When H_k is a zero matrix, it means that the kth connected-converted region D_k belongs neither to the candidate skull nor to the small bone hole. Therefore, those zero matrices H_k should be the regions separated by the detected skull G as the black regions either inside or outside the detected skull shown in the bottom three figures of Figure 5.

If the computed skull closing level F is equal to 2, it implies that the skull is integrated and ideal for the following steps of detection, such as slice No.27 shown in Figure 5. However, if F is not equal to 2, it means that the image cannot be used in the detection of the ideal midline due to an inappropriate scan position, such as the slices No.41 and No.15 in Figure 5. The skull closing level measure can quantitatively evaluate the integrity of the skull in head CT images. After closed skull inspecttion, all images with F ¹ 2 are removed from the slice subset.

2.1.3. Intracranial Area Detection

Based on clinical experience, the ideal CT slice for IML detection generally has larger intracranial area, such as the slices No.22 - 30 in Figure 2. Hence in this step, the area of the inner region surrounded by the detected skull, namely the intracranial pixels, is calculated and sorted for all remaining slices in the subset.

In the subset of images acquired after the closed skull detection step, every CT image should contain only two black regions which are separated by the detected skull. One of them is the intracranial region which contains the mass center of the skull and the other is the region outside of the skull. In order to calculate the intracranial area, the intracranial region has to be distinguished from region outside of the skull. This can be achieved using the coordinate of the skull’s mass center.

Generally, the image moment m_pq of the order p + q of

the digital image G can be defined as below,

(13).

where G_ij with the value of 1 or 0 is the intensity of the element at the ith row and jth column in the detected skull matrix G. Then, the coordinate of the mass center (x₀, y₀) of the detected skull can be obtained by

(14)

Thus the connected-converted region D_in containing the coordinate of the skull mass center (x₀, y₀) is the in tracranial region. The intracranial area S_in of this image is given by

(15).

The intracranial area of every slice in the subset is calculated and sorted in descending order. The first l slices with larger intracranial area are selected out for the following inspection. This number of l is a variable that depends on the number of slices for one patient or physician’s requirement.

(16)

where γ is a default number of selected slices and η is the selected percentage of the whole number of slices. After experimenting with various values, γ = 10 and η = 25% were finally chosen for this study.

2.1.4. Convexity Inspection

According to practical experience, another important characteristic of an ideal head CT image for IML detection is generally a good convexity for the intracranial region. For instance, the intracranial regions of the slices No.19 - 20 which are not ideal for IML detection both have integrated skull and larger intracranial area but show partially concavity. In contrast, intracranial regions of slices 26 through 28 have good convexity (see Figure 2). In addition, the concave shape of the intracranial region could affect the accuracy of the exhaustive symmetric position search, which is performed in the subsequent IML detection.

To qualitatively measure and evaluate the convexity of the intracranial region, we define a new measure Λ, called the intracranial convex measure. As shown in Figure 6(c), we extract the contour of the intracranial region. Then we can scan those pixels row by row. We define the far left and far right junctions (the blue points) of the ith row line (the upper red dash line) and the intracranial contour (the black curve) as points LP_i and RP_i in Figure 6(c). We use the function ξ(i,j) to describe the out-of-intracranial-region pixels between LP_i and RP_i on the ith row as the green bold line shown in Figure 6(c).

(17)

where R represents the intracranial region. Then the number of the out-of-intracranial-region pixels on the ith row φ(i) is given as below

(18)

We define the intracranial convex measure Λ using the total number of the out-of-intracranial-region pixels on all m rows of the image.

(19)

Then, we can rotate the image by φ degree as shown in Figure 6(d). The intracranial convex measure at angle φ can be calculated and noted as Λ_φ.

(20)

where Ψ_φ(i) the number of the out-of-intracranial-region pixels on the ith row in the image rotated by angle φ. With the sum of the all Λ_φ at all rotating angles, the total intracranial convex measure Λ_Total is given as below.

(21)

Larger values of the total intracranial convex measure Λ_Total represent increasingly worse convexity of the in tracranial region in a CT image. Using the sorted Λ_Total, we can keep the better l CT slices in the slice-subset by removing the ones with worse convexity. Value of l can be decided by Equation (16). For the example under study, using γ = 6 and η = 15% around 6 slices were obtained for the slice subset after convexity inspection.

With the implementation of the SSA algorithm on the raw slices, the number of slices greatly decreases. For

instance, in the case of Figure 1, the slice number of this patient’s scans decreases from 42 to only 6. The reduction in the number of slices effectively enhances the efficiency and saves computational time in the steps that follow.

2.2. Ideal Midline Detection

After the slice selection is performed using SSA algorithm, nearly all the slices of the obtained SSA subset are considered to be appropriate for the ideal midline detection (IML detection). Ideal midline detection is a twostep procedure, which includes the detection of the approximate ideal midline and subsequently refining the detected ideal midline, as shown in Figure 7. In the first step, the ideal midline can be approximated based on the skull symmetry, however including other features of the skull and the brain can help improve the accuracy of theapproximation. Thus, in the refined IML detection step, the bone protrusion on the upper part of the skull and the falx cerebri in the lower part are used to accurately detect the position of the ideal midline. To fully consider the connection among the slices in the subset, we utilize a global rotation assumption in both steps to determine the rotation angle of the skull. This method can further reduce the detection error due to individual non-ideal image.

2.2.1. Approximate IML Detection

In the calculation of the SSA algorithm, the detected skull G and its mass center (x₀, y₀) have been determined by Equations (10) and (14), respectively. To find the approximate ideal midline, we use the exhaustive symmetric position search algorithm which was developed for a prior work by our research group [8].

The row symmetry is defined as the difference in distance between each side of the skull edge and the current approximate midline. The CT image is rotated around the mass center of the skull. The symmetry cost S_θ of the image at the rotation angle θ is calculated as the sum of all row symmetry in the resulting image as shown in Equation (22).

(22)

where m is the number of rows in the image with the rotation angle θ (−45˚ < θ < +45˚ as used in this study) and measure l_i and r_i are the distance between the edge of the skull on the left/right side and the current approximate midline at the ith row. More details can be found in [8]. Finally, the rotation angle θ with the minimum symmetry cost S_θ determines the rotation direction of the midline of the brain for each particular CT slice.

(23)

where θ_p is the rotation angle of the midline on the pth slice and  is the symmetry cost of the pth slice at the rotation angle θ_pj.

All 6 CT slices in the SSA subset are processed one at a time using exhaustive symmetric position search. However, due to the non-uniform-thickness of the skull to serious deformation of the skull on one side after injury, it is hard to get an accurate position of the midline by processing only one slice. In this work, a global rotation assumption is used to decide the approximate ideal midline of all the CT images from one patient with full consideration of the connection among all the slices.

In the global rotation assumption, we assume that all CT images of one patient have the same rotation angle of the ideal midline due to the fixed posture of the patient during scanning. The rotation direction of the approxi-

mate ideal midline is determined by the median value of the rotation angles of all 6 slices in the SSA subset as shown below.

(24)

where the angles θ₁, θ₂, ∙∙∙, θ_λ are the rotation angles of the midlines determined by the exhaustive symmetric position search, l is the number of slices in the SSA subset, and θ_a is the approximate ideal midline of the whole set of slices. At the end of the approximate IML detection, the approximate ideal midline on each slice is calibrated to the vertical direction by rotating the skull by −θ_a angle.

2.2.2. Refined Ideal Midline Detection

Once the approximate midline is estimated and calibrated, brain anatomical features, such as the position of the falx cerebri and protrusion of skull bone, are used to refine the detection. In the detection of the falx cerebri and protrusion, we use the same algorithm from our previous work [8].

The falx cerebri is a strong arched fold of dura mater that descends vertically in the longitudinal fissure between the left and right cerebral hemispheres (Figure 8). In this work, edge detection method and Hough transform are used to detect this anatomical feature quickly and accurately. On the other hand, a bone protrusion is located in the anterior section of the skull. As shown in Figure 8, the bone protrusion curves down to a minimum point which is considered to be the upper starting point of the falx cerebri. To locate the lowest point of the protrusion curve, the derivative of the curve is calculated in a limited neighborhood area, which has been chosen to be 10 - 15 pixels in this work. The local minima point a is determined by the following equation.

(25)

where the function is the extracted curve of the interior bone edge and w is the neighborhood width. In fact, several small local minimal points may exist around the neighbor area of the protrusion due to the noise of the image or the irregularities of the skull. Using the maximal second derivative of the curve as the point a, Equation (25) is used to successfully extract the true protrusion minimal point by avoiding the influence of noise. More details of the detection of the falx cerebri and the bone protrusion can be found in [8].

Using the detected falx cerebri and the bone protrusion, we can obtain the refined rotation angle θ_q of the midline on each slice. Again, the global rotation assumption is used to determine the refined ideal midline of the whole set of slices. Rather than using the median method in the approximate midline detection, the weighted average method is used in this refine detection step. The rotation angle θ_f of the refined ideal midline of all the slices is given by

(26)

where θ_q is the refined rotation angle of the midline on the qth slice and μ_q is the weight of θ_q in the refined IML detection calculation.

(27).

where the values of v₁ and v₂ are both in the range of 0 - 1. We set v₁ = 0.2 and v₁ = 0.3 in this work.

At the end of the refined IML detection, the ideal midline on each slice is calibrated again to the vertical direction by rotating the skull by −θ_f angle. Therefore, in the two-step ideal midline detection, the ideal midline is centered by the mass center of the skull and rotated by an angle of −(θ_a + θ_f) from the original position in the slice.

3. Results and Discussion

3.1. Data

This database contains 3133 axial CT scan slices acquired across 70 patients with cases of both mild and severe Traumatic Brain Injuries (TBI). All the available CT scans have been utilized in testing the system’s processes and in the estimation of ideal midline.

3.2. Experimental Results and Discussion

This algorithm is primarily based on the anatomical characteristics of the skull and closely simulates the process of manual CT slice selection and decision making in IML by physicians.

In our dataset, all CT slices selected out by the SSA algorithm have been found to be acceptable for use in IML detection with the physician’s confirming. For instance, the patient in Figure 1 has 42 raw CT slices after CT scan. With the implementation of the SSA algorithm, only 6 slices were selected to be in the SSA subset as shown in the upper figures in Figure 9. All 6 slices are found to be ideal for the IML detection step due to the features such as integrated skull, large intracranial area, and good convexity of the skull. The result of the IML detection is displayed in the lower images of Figure 9. We can see that the detected ideal midline is accurately located in the middle of the skull and that the skulls in each scan are calibrated correctly.

In order to quantitatively measure the performance of the proposed system, the collaborating physician manually estimated IML is used as the ground truth. With a strict definition of accuracy, which is an allowed error of three pixels in horizontal distance δ between the estimated IML and the ground truth, the accuracy of IML estimation in this system is 96.9%, which is higher than 95% reported in [8]. In order to evaluate the result of the

mid-sagittal plane estimation, Ruppert et al. used an average z-distance measure to indicate the displacement between the resulting plane and the ground truth plane inside one image [1]. Therefore, the average z-distance measure has the similar physical meaning as the error δin our method. The lesser the mean value of error δ, the closer the estimated IML is to the ground truth. As the comparison in Table 1, the mean value of the error δ in our method is only 2.1 pixels, which is much less than the ones reported in [1,8,10]. The above experimental results demonstrate the high reliability of the proposed system.

4. Summary

In this paper, we proposed a system to identify the ideal midline (IML) using CT scans of patients with head injuries. The proposed SSA algorithm is used to closely simulate the process of manual selection of CT slice by physicians. With the implementation of SSA, a vast reduction can be achieved of the number of slices that is used for the computation of IML detection steps. Having fewer and more appropriate slices effectively increases the efficiency of the algorithm and also has the potential to save the cost and time required in practice. The global rotation assumption fully considers the connection among CT slices and thereby compensates the error generated by using a single CT slice. The obtained results show the potential for such a system to be implemented in clinical settings.

REFERENCES

NOTES

^*Corresponding author.