Research on a bifurcation location algorithm of a drainage tube based on 3D medical images.

Based on patient computerized tomography data, we segmented a region containing an intracranial hematoma using the threshold method and reconstructed the 3D hematoma model. To improve the efficiency and accuracy of identifying puncture points, a point-cloud search arithmetic method for modified adaptive weighted particle swarm optimization is proposed and used for optimal external axis extraction. According to the characteristics of the multitube drainage tube and the clinical needs of puncture for intracranial hematoma removal, the proposed algorithm can provide an optimal route for a drainage tube for the hematoma, the precise position of the puncture point, and preoperative planning information, which have considerable instructional significance for clinicians.

Introduction

Trauma can cause blood vessels to burst in the brain or between the skull and brain tissue. Subsequently, an intracranial hematoma may form from blood pooling up in the brain or between the skull and brain, which compresses the brain tissue. An intracranial hematoma is a common but serious secondary damage mode of craniocerebral injury. The incidence of closed craniocerebral injury is approximately 10%, while the incidence of severe craniocerebral injury is approximately 40%–50%. The complications of intracranial hematomas include their effect on cerebral blood flow, cerebral hernia, cerebral edema, and Cushing’s reaction, which can seriously damage brain tissue, and some of these damage modes are irreversible and life-threatening [1-4]. Puncture and ablation of the intracranial hematoma is a widely used medical method. With the increasingly busy modern life, heavy work stress, mental pressure, and lack of physical labor, the incidence rates show that the age of onset of this condition is decreasing. Once it occurs, the disease places a burden on individuals and families, so all countries and the international community have attached considerable importance to it [2]. The advantages of puncture operation include effective application, high security, quick recovery, short operating time, and relatively few postoperative complications [5-7]. At present, the methods of positioning such as the stereotactic technique [8] and neuronavigation technique [9] mainly involve a freehand technique for insertion of a drainage tube, which is based on fixed anatomical landmarks, does not consider individual variations, and often exhibits insufficient precision. Hence, a time-efficient and low-cost technique to localize the hematoma puncture point and to provide path planning will be beneficial, especially when highly sophisticated and expensive navigation systems cannot be made available in developing regions. In the present study, according to a dataset of brain computerized tomography (CT) images, we reconstructed the 3D model of a patient’s brain by using a 3D Slicer software [10-14] and extracted the area of intracranial hematoma as well. To quickly find the best puncture point in theory, an algorithm based on k-means clustering [15] was proposed to optimize the search space, eliminate redundant computation, effectively sort out the search space, and reduce the possibility of particles falling into a local extremum [16]. Then, the algorithm is improved from the viewpoint of the search mode, and the point cloud search algorithm of adaptive weighted particle swarm optimization [17-19] is proposed, which greatly reduces the calculation time and the number of iterations of the algorithm. In the experimental environment, the algorithm proposed in this study can greatly improve the search efficiency of this method in obtaining the global optimal solution. Compared with the direct optimization algorithm [20], the average maximum acceleration ratio is 1288.67%. According to the characteristics of the multitube drainage tube [21] and the clinical needs of the puncture process for intracranial hematoma removal, the bifurcation location algorithm provides the precise position of the puncture point, the optimum location of probes, and the best route under ideal conditions. In this study, the 3D model is effectively combined with the multitube drainage tube; furthermore, a preoperative simulation is proposed, which can provide significant guidance and value in formulating the puncture operation plan and decreasing the risk of blindness.

Segmentation of the intracranial hematoma

In clinical medicine, for patients with an intracranial hematoma, most hospitals conduct CT, hematoma location, puncture marking, followed by other processes [5, 6]. We used the 3D Slicer software to implement three-dimensional reconstruction based on brain CT data. This approach can greatly reduce the treatment time and improve efficiency. 3D Slicer software can reconstruct anatomical parts such as the hematoma, blood vessels, skull, and nerve bundles. In this way, doctors can avoid important parts such as blood vessels and nerve bundles when they mark the puncture points and accurately and efficiently puncture the hematoma, thus reducing treatment time [12-14]. The 3D Slicer program is a software platform for the analysis (including registration and interactive segmentation) and visualization (including volume rendering) of medical images and for research in image guided therapy. The platform originated in an MSc thesis program in 1998 and is jointly sponsored by the Surgery Plan Laboratory at Brigham and Women’s Hospital and the Artificial Intelligence Laboratory at the Massachusetts Institute of Technology [10-12].

Hematoma modeling

After the CT data have been preprocessed, three-dimensional reconstruction of the brain can be realized to obtain a three-dimensional image. In 3D Slicer, the region of interest (ROI) function is used to extract the hematoma area in the image, and the Segmentations module is selected to model the hematoma, providing models in the STL and OBJ formats to extract the “optimal external central axis” needed in this study.

Hematoma location locking

First, the Paint function in Slicer is applied to draw the hematoma, and then, the Fast Marching Method is used to draw the models of the intracerebral hematoma within a three-dimensional sphere, which is of great help in the follow-up work. Finally, the three-dimensional sphere is initialized by the Fast Marching function. Because the mask area around the hematoma is very large after initialization, the mask is narrowed by adjusting the slider for Segment volume, and the position of the hematoma is locked (Fig. 1).

Fig. 1

Process of distinguishing the hematoma in the yellow circle. a Thresholding; b Using smoothing to optimize staining; c Using Islands to separate the hematoma; and d The fusion of the hematoma and the three-dimensional map

Applying the threshold method to segment the hematoma

We select the Threshold function, in which the threshold is set to 50–100, set the selected intensity range to an editable intensity for the mask, and then, convert the masked area into a 3D model. For multiple small masks that appear in CT images, we use the Smoothing function to make the hematoma easier to segment by removing protrusions and filling holes to make the segmentation boundaries smoother. After smoothing, we use the Islands function to segment the hematoma from the numerous masks. We select the Keep selected island and click on the hematoma to segment it (Fig. 2). For specific modeling steps and methods, please refer to ref. [12].

Fig. 2

Schematic diagram of hematoma extraction. a Mask effect; b Smoothing; and c hematoma segmentation

Methods

Basic concept and definition

As an effective tool to describe the morphological characteristics of the geometry, the center axis has many advantages; for example, it can describe the geometry information inside the model and at any position of the boundary [22-24]. Furthermore, it provides a method to intuitively describe the geometric characteristics of a complex model and represent the topological relation of the model. The definition of the external central axis is as follows: Let Ω be a model in n-dimensional Euclidean space R, and assume it is homeomorphic to a 3D closed sphere; we define Ω as an n − D model. The external center axis of Ω is the set of the centers of the maximum circumscribed spheres [25]. Inspired by the concept of the external center axis, we found the optimal external center axis suitable for this study. The maximum Euclidean distance in the hematoma point cloud is determined, and the line segment joining the two points is defined as the optimal external center axis.

Algorithm of optimal external central axis extraction

The optimal external central axis can be obtained using a direct optimization algorithm [20]. However, a large amount of computation is required because the time complexity increases at the N2 level with an increasing quantity of point clouds. In the case of a large-scale geometric model, the approach to obtain the central axis is excessively time consuming. Therefore, to reduce the computation redundancy, a point cloud clustering algorithm based on k-means is proposed to optimize the search space [15, 16]. The algorithm presents strong adaptive ability, robustness, and low computation requirements. However, it does not consider the search method, and thus, it is not a suitable method for extracting the axis. Based on the above analysis, we propose the point cloud search arithmetic method for modified adaptive weighted particle swarm optimization [17–19, 26–30]. The algorithm reduces computational redundancy and greatly enhances the efficiency of searching axes in complex geometries. Our proposed algorithm is written in MATLAB (The MathWorks, Inc., MA).

Algorithm flow

The flow of the point cloud search arithmetic method for modified adaptive weighted particle swarm optimization based on k-means clustering is as follows:

Step 1. First, define the data set D in 3D Euclidean space and randomly select n items (n is the number of selected clusters) as centroids C(i = 1, 2, ⋯, n), where each centroid is the clustering center of a category. Then, calculate the Euclidean distance between C and D(1, 2, ⋯, k) in D. For example, if D is the nearest to C, it is classified as cluster i.

Step 2. Through the first step, the data set D is initially divided into k classes. Calculate the mean value of each dimension of all data items in each cluster. A new centroid is formed and updated to the new centroid.

Step 3. Repeat steps 1 and 2 to obtain the new centroid until the centroid of each class no longer changes.

Step 4. Randomly initialize the particles in the solution space.

Step 5. Each particle updates its velocity and position by tracking its own and the group’s optimal solution. The velocity updating formula of particle i in the dth dimension is as follows:

The position updating formula is

Repeat step 5 and incrementally iterate until the optimal position of the given threshold is found.

Note: is the d − th dimensional component of the k − th iterative particle i flight velocity vector; is the d − th dimensional component of the position vector; c1c2 is the maximum learning step parameter; r1r2 is the random function in the range of [0 − 1]; and ω is the inertia weight of the search, which has a considerable impact on the search scope and method.

Step 6. Compute the mean distance between d1, d2, ⋯, d: , where denotes the distance between two clustering centers.

Step 7. Calculate and obtain the difference between D and D: S = D − D, where ω denotes the ratio of S to ∑S as the generating weight of the small block ω: .

Step 8. Dynamically adjust the inertial weight of particles according to the NAIW method proposed by Li et al. [30]. where ωmax and ωmin are 0.9 and 0.4, respectively [20], Δx is the distance between the ith particles and the extreme position of the group, and k is the iteration coefficient. The NAIW method balances the weight of the particles under different circumstances, and it has advantages in convergence and robustness. Figure 3 is the flow diagram of the algorithm.

Fig. 3

Flow diagram of the algorithm to extract the optimal external middle axis

Optimal path planning

According to clinical demand and geometry symmetry, the straight line, where the optimal external axis is, is the path of the drainage tube into the brain. The multitube drainage tube can be divided into two subtubes and three subtube drainage tubes. The purpose of this study is to provide the corresponding optimal path planning for a specific multitube drainage tube [21]. The Table 1 is the explanation of symbols.

Table 1

Notation

Symbol	Description
M1(x1, y1, z1), M2(x2, y2, z2)	The coordinate of the joint point of the optimal external axis and the boundary of the hematoma
d12 =  ∣ M1M2∣	The distance from M1 to M2
d34 =  ∣ M3M4∣	The maximum diameter of the vertical surface through M1/2 in the hematoma
M3(x3, y3, z3), M4(x4, y4, z4)	The coordinates of the endpoints of d34 in the hematoma
D, d	Ddenotes the diameter of the duct in the main tube; d denotes the diameter of the hole in the subtube
dO	The distance between the two central points of the holes
dB	The distance between the two bifurcations
d_B_i	The distance between the main probe tip and the i − th bifurcation-tube probe tip (we define the first subtube to be nearest to the main probe tip)

Note: Only some common symbols are listed here; the rest are explained later

The known conditions can be determined by M1 and M2

According to the intersection coordinates M1(x1, y1, z1) and M2(x2, y2, z2), we obtain the parameter equation of the center axis as follows:

(6)

The midpoint coordinates of M1M2 are .

The distance from M1 to M2 is .

The plane equation of M3M4 is

The distance from M3 to M4 is .

A comparison of the distance between d12 and can offer directional evidence for puncture surgery.

The flow diagram of the bifurcation location algorithm of the drainage tube is as shown in Fig. 4.

Fig. 4

Flow diagram of the bifurcation location algorithm of the drainage tube

Optimal path planning of two subtube drainage tubes

According to the clinical demand, the best position of the drainage tube in the hematoma is the midpoint between the two holes, coinciding with the midpoint of the central axis, that is, , and the path of the main drainage tube coincides with the optimal external center axis.

Calculating the coordinates of B1 and B2

Substituting t1, t2 into Eq. (6), we obtain the coordinates of and .

Calculating the direction vector of the subtube

To find two line segments through M3 with a fixed angle between this axis and the center axis, let the intersection points between the line segment and the axis be M31, M32; these points should meet the following limits:

Substituting t31, t32 into Eq. (6), we obtain the coordinates of M31, M32; the direction vector of subtubes are and , and the parameter equation of the two subtubes is as follows:

Scenario one:

(12)

Scenario two:

(13)

The intersection points between the four rays and the boundary of the hematoma are B11, B21, B12, and B22.

In the first scenario, the coordinates of the probe tip of the subtube are as follows:

(15)

In the second scenario, the coordinates of the probe tip of the subtube are.

(16)

The optimal path planning of three subtube drainage tubes

According to the clinical demand, the best position of the drainage tube in the hematoma is the midpoint of the second hole, coinciding with the midpoint of the central axis, i.e., , with the path of the main drainage tube coinciding with the optimal external center axis.

Calculating the coordinate of B1, B2, and B3

Substituting t11, t12 and t13 into Eq. (6), we obtain the coordinates , , and .

To find a line segment through M3 with a fixed angle between the axis M1M2 and the center axis, let the intersection points between the line segment and the axis be M31(x31, y31, z31); they should satisfy the following limits:

We obtain .

Because the angle between two subtubes and the angle between the subtube and the main tube are fixed values, we can calculate the other two direction vectors of the subtubes based on the geometrical relationship.

The parametric equation of the three subtubes is as follows.

There are three cases for the direction vector of each subtube; let the coordinates of the subtube and the boundary of the hematoma be B11, B12, B13; B21, B22, B23; B31, B32, B33.

When we find the maximum value , the coordinates of the bifurcation subtube probe in hematoma can be obtained as , , and .

Results and discussion

According to the CT data, the hematoma area was segmented, and an independent three-dimensional diagram was obtained, as shown in Fig. 5.

Fig. 5

Three-dimensional intracranial hematoma model

The STL format of a three-dimensional point cloud is obtained by hematoma modeling. Figure 6 shows the three-dimensional point cloud diagram in the hematoma area, and Fig. 7 shows the triangulation diagram of this area [24].

Fig. 6

The point cloud model

Fig. 7

Schematic diagram of triangulation of the hematoma in STL format

In this case, the point cloud search arithmetic method for modified adaptive weighted particle swarm optimization to search the hematoma region uses 81,226 point-chains. The coordinates of the intersection between the central axis and the boundary of the hematoma are M1(x1, y1, z1) and M2(x2, y2, z2). The optimal external central axis of the hematoma is shown in Fig. 8. The intersection points of the central axis and the boundary of the hematoma are M1(28.555, 9.623, 7.123) and M2(43.384, 21.295, 48.887).

Fig. 8

Optimal external central axis of the hematoma. a Two-dimensional projection; b Three-dimensional projection; and c The exact position in skull

Therefore, the midpoint coordinates of the center axis M1M2 are M1/2(35.970, 5.836, 28.010), and the distance from M1 to M2 is d12 =  ∣ M1M2 ∣  ≈ 54.0377. Figure 8 is optimal external central axis of the hematoma.

The route planning for the puncture operation is as follows:

Because , it is suggested that doctors choose three subtube drainage tubes when performing puncture operations. From the coordinate of the center axis, the plane equation through M3M4 can be obtained as

According to the point cloud search arithmetic method for modified adaptive weighted particle swarm optimization, we obtain the coordinates as M3(26.201.451, 18.457) and M4(36.353, −7.797, 37.816), and the distance from M3 to M4 is d34 = 38.134751.

Through the above algorithm process, the coordinates of the three bifurcation points can be calculated

The three direction vectors of the subtube are as follows:

There are three cases for the direction vector of each subtube; let the coordinates of the subtube and the boundary of the hematoma be B1, 1, B1, 2, B1, 3; B2, 1, B2, 2, B2, 3; B3, 1, B3, 2, B3, 3. The maximum value from the bifurcation point to the boundary of the hematoma is

The coordinates of the bifurcation probe to reach the hematoma area are

Figure 9 is a schematic diagram of the path planning process.

Fig. 9

Route planning. M2 denotes the accurate position of the puncture point

In this study, 3D Slicer software was used to reconstruct a three-dimensional model of an intracranial hematoma, and an improved algorithm based on clustering and adaptive particle swarm optimization to extract the optimal external central axis of the hematoma was proposed. According to the characteristics of the drainage tube and the specific morphology of the intracranial hematoma, optimal path planning under ideal conditions was performed, and the preoperative simulation was provided, which will be useful in clinical medicine.

The bifurcation localization algorithm of the drainage tube proposed in this study assumes an ideal state. It is necessary to consider the distribution of peripheral nerves and blood vessels and the specific location of the hematoma in the brain to apply this method to clinical medicine. The focus of subsequent research work will be to optimize the algorithm by considering the above factors, to implement the diffusion model and simulation of the drug injection, to plan the optimal surgical project, and to provide a dynamic simulation before operation.