Automatic 3D pelvimetry framework in CT images and its validation

The proposed framework consists of four modules: (1) 3D pelvic model reconstruction; (2) detection of spatial femoral head centers; (3) recognition of the superior sacral endplate; (4) calculation of pelvic parameters. The second and third modules can be executed in parallel, and structure of the proposed framework is demonstrated in Fig. 1.Fig. 1Schematic diagram of the proposed framework. Four modules are included: 3D pelvic model reconstruction (red box), detection of spatial femoral head centers (black box), recognition of the superior sacral endplate (purple box), and calculation of pelvic parameters (blue box).3D pelvic model reconstructionPelvic CT images, including skeleton, soft tissue, organs and other components, are usually acquired from different scanners, thus extraction of bone regions by a classic thresholding method is not accurate enough to obtain a satisfactory pelvic model. Therefore, the K-means algorithm15 was chosen to segment the CT images into background and bone parts. By incorporating the weighted quality evaluation function, this module divided the image pixels into two clusters based on the gray scale feature of CT images. Segmentation quality \(E\) was measured by the weighted distance between the pixels and their cluster centers, which can be formulated as follows:$$\begin{array}{c}\begin{array}{c}E={\sum }_{i=1}^{\text{K}}\frac{{n}_{i}}{N}{\sigma }_{\text{i}}\end{array}\end{array}$$
(1)
where \({\sigma }_{i}\) represents the standard deviation of the ith cluster, and \(N\) is the total number of pixels that need to be segmented,Firstly, two different gray values were randomly selected as the initial clustering centers on the image to be segmented, and the weighted distance \(L\left(P, {O}_{i}\right)\) between each pixel and each clustering center was calculated as follows:$$\begin{array}{c}\begin{array}{c}L\left(P, {O}_{i}\right)=\frac{1}{2\text{N}}{\left(\frac{1}{{\sigma }_{i}}d\left(P,{O}_{i}\right)\right)}^{2}\end{array}\end{array}$$
(2)
where \(d(P,{O}_{i})\) is the Euclidean distance between object \(P\) and object \({O}_{i}\). Then, the pixel was put into the category with the smallest distance, and the clustering centers of each category were updated by using the average value. Next, the cluster quality \(E\) was calculated by Eq. (1). If \(E\) reached the expected value or the number of iterations reached the maximum number, the process would be stopped; otherwise, the clustering process was re-iterated. Finally, we got the segmented pelvic images. Introduction of these two equations ensures that points at the boundary are more likely to be divided into large clusters, which would make the boundary of image blocks clear and improve the quality of the reconstruction model.As an open-source software for manipulating and displaying scientific data, the Visualization Toolkit (VTK) comes with state-of-the-art tools for 3D rendering and a suite of widgets for 3D interaction. In addition, VTK is also convenient to use across platforms, thus VTK was chosen as the primary tool of the proposed framework to realize visual operations and data operations. Here the preprocessed CT images were imported into the data pipeline of VTK, and the built-in marching cubes (MC) algorithm was used for 3D reconstruction of pelvic models.The reconstructed 3D model can be saved into a vtk format file for direct read later, and can also be displayed directly in the system using the visualization pipeline of VTK. An original image, the extracted image of skeleton, and the 3D reconstruction model of pelvis are shown in Fig. 2.Fig. 2Process of the 3D reconstruction. (a) Original image; (b) Bone part image extracted from the original image; (c) 3D reconstruction model of pelvis.Detection of spatial femoral head centersThe femoral head is round, and its top half is similar to a sphere. Central points of the two femoral heads are essential reference points in the pelvic parameter measurement. To automatically locate the femoral heads on the established pelvic model, the femoral heads was first segmented on the 2D images, and centers of the femoral heads were then determined and mapped to the 3D model. Core of this module was the segmentation of the femoral heads on 2D images by deep learning technology.Structure of the DRINetTo recognize the circular regions (i.e., femoral heads) on CT images, the DRINet was used as image segmentation network in this module13. The DRINet is based on the ideas of DenseNet16 and Inception-Resnet17, and consists of three parts: Dense Connection Block (DC_Block), Residual Inception Block (RI_Block) and Unpooling Block (Unpooling_Block). DC_Block is the encoder of the DRINet, and its main function is to extract and aggregate image features. RI_Block and Unpooling_Block are the decoders of the DRINet, and mainly upsample the image and output the final image. RI_Block and Unpooling_Block both adopt multi-branch deconvolution operations, and finally concatenate the features to realize image sampling. BN layers are also used to prevent the gradient from disappearing.Three DC_Blocks, six RI_Blocks and three Unpooling_Blocks are combined to get the DRINet. The structure of each block and the whole network are shown in Fig. 3.Fig. 3DRINet trainingImage data from local hospital and DeepLesion18 were used to train the DRINet, and there were 2500 CT scans containing femoral heads. Pelvic image sequences of 15 patients were locally acquired using a Siemens Definition AS 40 CT system, and scanning parameters are listed as follows: tube voltage, 120 kV; tube current, 154 mA; scanning field of view (DFOV), 407.0 mm; axial slice thickness, 5 mm; inter-slice gap, 5 mm; reconstruction thickness, 1 mm; reconstruction interval, 1 mm; reconstruction matrix, 512 \(\times\) 512. This study was approved by the Institutional Review Board of Anhui Medical University and in accordance with the Declaration of Helsinki, and informed consent was obtained from all patients.Before training, the original CT scans were converted into gray images in BMP format to improve training speed. Each image has a size of \(512\times 512\), and its corresponding label is an image that contains only the circular area of the femoral head with the same size. In the training stage, the learning rate was set to 0.001 and training epochs to 250, the optimizer was set as Adam, the loss function was set as Dice_loss, and both images and labels were normalized. The image data were divided into a training set (80%) and a validation set (20%), and K-fold cross-validation19 was used in training. Definition of Dice_Loss is formulated as follows:$$\begin{array}{c}\begin{array}{c}Dic{\text{e}}_{\text{loss}}= 1- \frac{2\left|X\cap Y\right|}{\left|X\right|+\left|Y\right|}\end{array}\end{array}$$
(3)
where \(X\) represents the label image, \(Y\) represents the output image, \(\left|X\cap Y\right|\) represents the intersection of the label and output images, and \(\left|X\right|\) and \(\left|Y\right|\) represent the number of elements in \(X\) and \(Y\) respectively.For an image containing femoral heads, its corresponding output of DRINet should contain circular regions; and for an image in the absence of femoral heads, its corresponding output of DRINet should be blank. Femoral head can be approximated as circle, and the Hough transform algorithm20 was used to detect the centers and radii of the circles in the output images. Figure 4 shows examples of an original image, output image of the DRINet, and image after circle detection.Fig. 4Examples of an original image, output image of the DRINet, and image after circle detection. (a) Original image; (b) Output image; (c) Circle detection.3D sphere fitting of femoral headFor images containing femoral heads, center and radius of each femoral head region have been obtained in “DRINet training”. For the left and right femoral heads of each patient, centers and radii of the femoral head regions in all CT images were traversed, and therefore the maximum radius and its corresponding center were respectively taken as the radius and hemispherical center of each femoral head.According to the pixel spacing of the continuous images, the center and radius of each femoral head were mapped into 3D space, and a sphere at the corresponding position was drawn using VTK to fit the femoral head. Because the centers and radii determined on the 2D images are pixel coordinates, we use \({pixel}_{x}\) to represent the coordinates on the x-axis of the image, \({x}_{spacing}\) to represent the pixel spacing on the x-axis of the image, and \({space}_{x}\) to represent the space coordinate obtained after mapping. We set the serial number of the image which contains a central point of the femoral head as \({current}_{z}\), which represents the position of the central point on the z-axis in the pelvic image sequence. The space point mapping is formulated as follows:$$\begin{array}{c}\begin{array}{c}\left\{\begin{array}{c}\begin{array}{c}spac{e}_{x}=pixe{l}_{x}\times {x}_{spacing}\\ spac{e}_{y}=pixe{l}_{y}\times {y}_{spacing}\\ spac{e}_{z}=curren{t}_{z}\times {z}_{spacing}\end{array}\end{array}\right.\end{array}\end{array}$$
(4)
Thereafter, center of the femoral head is denoted as \({P}_{f}\) with coordinates of \(({space}_{x}, {space}_{y}, {space}_{z})\).Recognition of the superior sacral endplateFor the superior sacral endplate, normal vector \({V}_{s}\) and central point \({P}_{s}\) are also important variables for calculating pelvic parameters. However, it is time-consuming and laborious to manually determine the position of endplate by conventional methods. In this module, automatic recognition of the superior sacral endplate was carried out by deep neural network, and there were three steps: (1) recognizing a transition image between L5 (i.e., the fifth lumbar vertebra) and S1 (i.e., the first sacral vertebra) by neural network; (2) determining the point set \({S}_{p}\) of the superior sacral endplate in the image by analysis of image connected region; (3) mapping the 2D points in \({S}_{p}\) into 3D space, then all points and their normal vectors of the endplate were obtained by the plane growth algorithm to calculate \({P}_{s}\) and \({V}_{s}\).Image classificationRecognition of vertebral types in images by neural networks is a classification problem. The VGG16 network14 was empirically chosen as the classification network, and several convolution kernels with size of \(3\times 3\) were used to replace the large convolution kernel in the standard CNN. Compared with the standard CNN, the VGG16 not only reduces the number of parameters, but also learns more complex features. However, the standard VGG16 has a large number of neurons in the last fully-connected layer, which may lead to too many parameters and reduce the efficiency of the network, so the number of neurons in the fully-connected layer was reduced appropriately.A number of 2000 images were selected from DeepLesion (mentioned in “DRINet training”) as the training set of the VGG16 network, and the labels were the classification number of each image rather than images. All pelvic slice images were classified into three categories: lumbar vertebrae (ID: 0), S1 (ID: 1) and femoral heads (ID: 2). Because of the particularity of pelvic spatial structure, an image may contain two objectives, so we stipulated that if there were both vertebral L5 and S1 in the image, its category was lumbar vertebrae; if there were both S1 and femoral heads in the image, its category was S1. Figure 5 shows typical images of these three classification categories.Fig. 5Typical images of the three classification categories. (a) An image with L5 (ID: 0); (b) An image with Sl (ID: 1); (c) An image with femoral heads (ID: 2).In the training stage, we set the number of neurons in the last three fully-connected layers of the VGG16 network as 128, 256 and 3 respectively. The activation function in the network was the Relu function, the learning rate was 0.001, the optimizer was Adam, the loss function was MSE, and the training epochs was 50. Finally, the training accuracy of the VGG16 network was over 95%, and the network converged with its parameters saved in an h5 file. Output of the network was a vector of \(1\times 3\), and each component in the vector represented the probability of classification. Subscript of the component with the highest probability was chosen as the final classification result.Image connected region analysisSuperior sacral endplate locates below the lumbar L5, and it is a slope, so it spans several 2D pelvic slices. But there is a rule: there must be some points belonging to this plane in several images of the junction of L5 and S1. We chose the images belonging to the junction of L5 and S1 according to the classification results of the VGG16 (generally select one to two images before the image classified as S1), observed them and found that the largest areas near the centerline of these images belong to S1, and the points on the edge of the top part of S1 was on the superior sacral endplate. Therefore, as long as some points of the superior sacral endplate in the 2D image were mapped to the 3D space as anchors, plane-growth can be carried out in the 3D model, and center points and normal vectors of the plane can be obtained.As for how to get these points, our initial idea was to use deep learning to get the results directly, but the number of images belonging to the junction of L5 and S1 in the dataset was too small to allow the network to be fully trained. So we put forward another method according to the observation: we first used the 4-connected region algorithm to find the largest connected region \({R}_{\text{l}}\) in an image of the junction of L5 and S1, then we calculated the central point \({O}_{R}\) of \({R}_{\text{l}}\), and took out 20 adjacent points which locates above \({O}_{R}\) and at the edge of \({R}_{\text{l}}\). Most of these 20 points were points on the superior sacral endplate.Finally, the points on these images were mapped into 3D space by Eq. (4), and these points would be saved in point set \({S}_{a}\) to be used in the next step. Figure 6 shows the process of finding points of the top endplate of the S1 on a 2D image.Fig. 6Process of finding points on the superior sacral endplate in a 2D image. (a) The image at the junction between L5 and S1; (b) The largest connected area; (c) The points on the upper plane of the sacrum.Plane growth on the 3D modelThe point set \({S}_{a}\) was obtained in the previous step, and in this step, \({S}_{a}\) would be used to find the spatial points on the superior sacral endplate in the 3D pelvic model. The plane growth algorithm is described in Algorithm 1.Calculation of pelvic parametersBased on the center point \({P}_{f}\) and the radius of femoral head (“Detection of spatial femoral head centers”), and the center point \({P}_{s}\) and the normal vector \({V}_{s}\) of superior sacral endplate (“Recognition of the superior sacral endplate”), pelvic parameters were calculated. It should be noted that there were two central points of the femoral heads, which were denoted as \({P}_{f1}\) and \({P}_{f2}\) respectively.Here we specified that the normal unit vector of the horizontal plane of space was \(Z\), then the formulas for calculating the three 3D pelvic parameters were as follows:$$\begin{array}{c}\begin{array}{c}{P}_{mid}= \left({P}_{f1}+{P}_{f2}\right)/2\end{array}\end{array}$$
(5)
$$\begin{array}{c}{\text{PI}}_{3\text{D}}=arccos\frac{{(P}_{\text{s}}-{P}_{mid})\cdot {\text{V}}_{\text{s}}}{\Vert {P}_{s}-{P}_{mid}\Vert \Vert {\text{V}}_{\text{s}}\Vert }\end{array}$$
(6)
$$\begin{array}{c}{\text{PT}}_{3\text{D}}= arccos\frac{{(P}_{\text{s}}-{P}_{mid})\cdot \text{Z}}{\Vert {P}_{s}-{P}_{mid}\Vert }\end{array}$$
(7)
$$\begin{array}{c}\begin{array}{c}{\text{SS}}_{3\text{D}}= \text{arcsin}\frac{{V}_{s}\cdot Z}{\Vert {V}_{s}\Vert }\end{array}\end{array}$$
(8)
The calculation of the 3D pelvic parameters is an extension of the 2D pelvic parameters12.Specifically, when automated measurements fail to meet requirements, the system offers a manual adjustment module. The manual adjustment can utilize the functionality provided by VTK to manually move the predicted femoral head sphere to the actual femoral head position. The manual adjustment module operates on the premise of automated measurements, meaning that before using the manual adjustment module, the corresponding automated measurement module must be called. When using the femoral head center position adjustment module, the 3D sphere representing the femoral head already drawn in the 3D scene was first cleared. Then, the femoral head coordinates obtained from the automated measurement module were modified according to the direction of movement. Finally, based on the new femoral head coordinates, new femoral head was drawn in the 3D scene.To verify the accuracy of the pelvic parameters automatically generated by the proposed framework, pelvic image sequences of 15 patients (mentioned in “DRINet training”) were processed by an orthopedic surgeon. Firstly, the 3D pelvic slice images were mapped to the sagittal plane, and a clinician was asked to measure the pelvic parameters manually on the 2D image. The clinician was then asked to identify the superior sacral endplate and the femoral head in a completely manual manner in the initial framework and then obtained the 3D pelvic parameters. Next, these images were imported into the proposed framework, and the 2D and 3D pelvic parameters were calculated automatically and compared to verify the validity of the pelvic parameters automatically generated by the proposed framework.