Watertight 2-manifold 3D bone surface model reconstruction from CT images based on visual hyper-spherical mapping

1.   Introduction

Bone medical images generated from Computed Tomography (CT) scanning of patients are widely used as a powerful tool to provide diagnosis and analysis of bone mineral density, bone injury, bone tumor, etc. ^[1,2,3].Since CT medical imaging data is composed of multi-layer 2D images and cannot provide enough 3D spatial geometry and topology information, watertight 2-manifold 3D bone surface model reconstruction from 2D CT images is crucial to make follow-up surgical plans, finite element modeling & analysis and so on ^[4,5,6]. Because exposure to radiation at high doses may cause harm to the human body, low dose and thick slice (greater than 0.625mm) are often used in clinical CT scan. Consequently, CT artifacts ^[7,8] during CT scanning inevitably result in problems such as noise, blurry tissue boundaries and locally non-uniform distribution of HU (Hounsfield Unit) values in CT images as shown in Figure 1. HU values corresponding to the femoral head and femoral condyle regions have obvious non-uniform distribution (Figure 1(A), (D) and (E)), and the 3D reconstructed surface models using such CT images thus have precision problems, such as redundant and abnormal shapes, holes, and stair-step shapes (Figure 1(H) and (I)). The 3D bone surface model reconstruction from CT image includes two key steps: image segmentation ^[9] (Figure 1(B)) and surface model generation ^[10] (Figure 1(F) and (G)). CT image segmentation methods can be classified into three typical kinds ^[11,12,13] including manual, semi-automatic ^[11], and automatic ^[12,13]. Considering the complexity of CT images, many segmentation algorithms are designed for certain kind of bones or organs ^[14,15], and algorithms for full body skeleton still have its limitations ^[16], which are not yet applied widely. 3D bone surface model reconstruction from point cloud data ^[10] is accomplished by interpolating triangular meshes ^[17,18] or implicit function fitting methods ^{[19,20,21,22,23]}, and the model precision mainly depends on the segmentation results of CT images. On the one hand, bone surface reconstruction is extremely sensitive to the resolution of CT images and parameters selected for CT image segmentation. On the other hand, there are problems such as non-obvious bone boundaries, overlapping of threshold ranges between regions corresponding to the bone and its nearby cartilage or between adjacent bones. The factors related above bring many difficulties in CT image segmentation and 3D bone surface model reconstruction from CT images, and so that research interest is attracted.

Researchers have never stopped studying medical image segmentation algorithms ^[11,12,13] since the invention of CT scan. The common used methods are in the top-down form ^{[24,25,26,27]}, such as region-growing ^[24], active contour model ^[25], statistical shape ^[26], and graph cuts methods ^[27]. The top-down methods have excellent performance on segmentation of certain organs, but they are not suitable for accurate segmenting bones of the full body skeleton, due to complicated distribution of bone HU values in CT images, especially at the presence of fractures or bone tumors.

Recently, some researchers introduced bottom-up methods based on statistical machine learning on medical image segmentation ^{[13,16,28,29,30,31,32,33]}, which are mainly classified into two categories: clustering based methods ^[28] and classification based methods ^[29,30,31]. Clustering based methods, such as Mean-shift and K-means, assign pixels with similar properties to the same clustering block which is corrected iteratively until the clustering result is convergence. Clustering based methods belong to unsupervised learning algorithm where no predefined training samples are needed, and thus the segmentation result is usually coarse. On the contrary, the classification based methods, CNN ^[29,30] and SVM ^[31] for example, which are used in supervised learning algorithm where predefined training samples are needed, have good segmentation result. Study on combination of the above two categories' advantages to segment medical images is one of the hot topics in recent years ^[16,32,33]. The irregular shape such as bone tumor and complicated HU values distribution will bring great inaccuracy to the bottom-up methods such as CNN, SVM, and there are few general-purpose mature algorithms that can automatically segment all types of bones from CT images. That means proper CT image segmentation methods should be chosen when segmenting different types of bones or different regions of the same bone.

After segmentation of CT images, bone's voxel or pixel data and its corresponding volume sampling points can be obtained. The next step is to extract bone's boundary information or surface points from the volume sampling points used to reconstruct bone surface model. Since healthy bone surfaces have the characteristics of continuity and smoothness, implicit surface reconstruction algorithms such as the tangent planes ^[19], RBF ^[20], MLS ^[21], MPU ^[22], and Poisson equation ^[23], are usually used to reconstruct the watertight 2-manifold bone surface model. Implicit surface reconstruction algorithms can effectively reconstruct surface models from point cloud data that have complicated topologies. Nevertheless, due to the existence of noise, stair-step effect, data loss, etc., accurate bone surface models cannot be achieved by directly applying surface reconstruction operation on the extracted surface points. Therefore, studies on methods of accurate bone surface model reconstruction are meaningful.

We draw inspiration from methods including thresholding ^[34], point cloud processing ^[23,35], shape modeling ^[36,37] and computer vision ^[38]. This paper proposes a method to reconstruct accurate watertight 3D bone surface model from CT images. The volume sampling points are used in the surface point extraction and accurate surface generation stages of the proposed method, which are also named as initial volume sampling points with narrower threshold range and target volume sampling points with wider threshold range respectively. Overview of the algorithm is described as follows:

Firstly, resample the CT images according to the requirement of reconstruction precision; segment bone's ROI (Region of Interest) by thresholding with combination of 3D morphology operation; extract the initial and target volume sampling points respectively. Secondly, calculate visible surface points from the initial volume sampling points by method based on visual hyper-spherical mapping proposed in this paper. Thirdly, generate the watertight 2-manifold surface model from the visible surface points by implicit surface reconstruction. Finally, take the subset of the target volume sampling points that are near the actual bone surface as constraint point set; deform the initial surface under given constraints to generate the accurate bone surface model.

The remainder of this paper is organized as follows: section 2 describes details of the proposed 3D bone surface reconstruction method including volume sampling point extraction, visible surface point calculation, initial surface reconstruction, and accurate surface generation; section 3 presents experiments; section 4 discusses result; section 5 draws conclusions.

2.   Methods

2.1.   Volume sampling point extraction

Let XYZ-O denote the 3D Cartesian coordinate system of the CT images with origin at the center of the first image slice. Z axis parallels to the scanning direction. X and Y are parallels to U and V of the image's 2D pixel coordinate system UV-O respectively. The Cartesian coordinates (xyz) of a sample point ${\boldsymbol{p}}_{\boldsymbol{m}, \boldsymbol{n}, \boldsymbol{l}}\in {\mathbb{R}}^{1\mathrm{x}3}$ are denoted as ${\left\{{x}_{m}\right\}}_{0}^{{N}_{x-1}}$, ${\left\{{y}_{n}\right\}}_{0}^{{N}_{y-1}}$ and ${\left\{{z}_{l}\right\}}_{0}^{{N}_{z-1}}$, where ${x}_{m} = mΔx$, ${y}_{n} = nΔy$, ${z}_{l} = lΔz$. $Δx, Δy$ are spatial distances between two adjacent pixels along X axis and Y axis respectively, and $Δz$ is the slice thickness. Generally, $Δx$ and $Δy$ are equal to each other, and $Δz$ varies greatly for different CT scans and clinical requirements. ${N}_{x}, {N}_{y}$, and ${N}_{z}$ denote CT image slice numbers along X, Y, and Z axis, respectively. Thus, the number of volume sampling points is ${N = N}_{x}\times {N}_{y}\times {N}_{z}$. All sample points form a point set $\boldsymbol{P} = \left\{{\boldsymbol{p}}_{\boldsymbol{m}, \boldsymbol{n}, \boldsymbol{l}}\right\}$, $\boldsymbol{P}\in {\mathbb{R}}^{\mathrm{N}\times 3}$. The voxel corresponding to ${\boldsymbol{p}}_{\boldsymbol{m}, \boldsymbol{n}, \boldsymbol{l}}$ can be denoted as ${\boldsymbol{w}}_{\boldsymbol{m}, \boldsymbol{n}, \boldsymbol{l}}\in {\mathbb{R}}^{1\times 3}$, and the voxel set corresponding to $\boldsymbol{P}$ is $\boldsymbol{W} = \left\{{\boldsymbol{w}}_{\boldsymbol{m}, \boldsymbol{n}, \boldsymbol{l}}\right\}$, $\boldsymbol{W}\in {\mathbb{R}}^{\mathrm{N}\times 3}$. The HU value of ${\boldsymbol{w}}_{\boldsymbol{m}, \boldsymbol{n}, \boldsymbol{l}}$ is represented by ${H}_{m, n, l} = F\left({x}_{m}, {y}_{n}, {z}_{l}\right)$.

if $Δx = Δy\ne Δz$, the extracted voxel data exhibits anisotropy, which will cause surface of the generated 3D bone model being unnatural. In order to obtain isotropic voxel data ($Δx = Δy = Δz = d$), numerical interpolation and resampling on CT images is needed. The commonly used methods ^[39] include Nearest-Neighbor, linear, quadratic, B-spline, cubic, Lagrange, and Gaussian interpolations. Because B-spline interpolation method has the advantages of good robustness, high precision and efficiency, which is easy to implement, this paper adopts it to obtain evenly spaced resampling points data with better resolution. Figure 2 shows the 3D bone surface models reconstructed before and after CT images being resampled. High resolution resampling with small thickness simultaneously on all three coordinate directions can effectively reduce the stair-step effects both in CT images and reconstruction results, and enhances the clarity of CT images.

After interpolation and resampling, the following step is to segment bone's ROI and extract its corresponding voxel data or volume sampling points from CT images. The two-step thresholding method is adopted in this paper because of its high efficiency and robustness. Let ${\boldsymbol{W}}_{H} = \{{H}_{m, n, l}\in \left[{\theta }_{min}, {\theta }_{max}\right], {\boldsymbol{w}}_{m, n, l}\in \boldsymbol{W}\}$, where ${\theta }_{min}$ and ${\theta }_{max}$ are the lower and upper limit value of the threshold range. ${\theta }_{max}$ for a bone tissue is generally fixed, so the quality of segmentation results mainly depends on ${\theta }_{min}$, and the segmentation results have the behavior of outward expansion or inward shrinking with the change of threshold range.

As ${\theta }_{min}$ increase gradually, the connection regions between adjacent bones in the binary images will be shrunk, but the holes and regions with data loss will become large in the corresponding binary images. When ${\theta }_{min}$ is a value bigger enough to segment the bone separately from adjacent bones while still making the segmentation retain the main surface profile of the bone as much as possible, such as ${\theta }_{min} = {H}_{0}$, the initial volume sampling points ${\boldsymbol{P}}_{H0}$ corresponding to voxel data ${\boldsymbol{W}}_{H0}$ are obtained.

As the threshold range become wider with the decrease of ${\theta }_{min}$, the holes and those regions with data loss of the segmentation results will be gradually filled and recovered, and also noise data from cartilage and adjacent bones will be introduced which cause regions of adjacent bones being merged together. When ${\theta }_{min}$ is a value smaller enough to make the segmentation result fully contain the region of the target bone while introducing noise points as few as possible, such as ${\theta }_{min} = {H}_{T}$, the target volume sampling points ${\boldsymbol{P}}_{HT}$ corresponding to voxel data ${\boldsymbol{W}}_{HT}$ are obtained.

The segmentation results can be visually presented as adjusting ${\theta }_{min}$ dynamically, so ${H}_{0}$ and ${H}_{T}$ are easy to be determined interactively. The narrower threshold range $\left[{H}_{0}, {\theta }_{max}\right]$ is used to extract volume sampling points for reconstructing the initial surface model, and the wider threshold range$\left[{H}_{T}, {\theta }_{max}\right]$ is used in the iterative deformation stage to further restore local shape of the bone surface, especially for shape near the joint position. Figure 3 shows the binary image segmentation results corresponding to different threshold ranges.

In order to reduce the negative influence on the volume sampling point extraction such as noises, holes and data loss caused by directly using thresholding method, this paper applies the morphology operation on the segmentation result to fill holes, smooth boundaries and break small connections between adjacent bones.

Let ${\boldsymbol{W}}_{H}^{k}$ denote the voxel set of bone with index $k$ ($k\in \left[\mathrm{0,206}\right]$, 206 is the total bone number of Asian), and its corresponding volume sampling data point set is ${\boldsymbol{P}}_{H}^{k}$. Let ${Erode}^{m}\left(\boldsymbol{W}\right)$, ${{Dilate}}^{m}\left(\boldsymbol{W}\right)$, ${{Open}}^{m}\left(\boldsymbol{W}\right)$, and ${{Close}}^{m}\left(\boldsymbol{W}\right)$ represent the corresponding $\mathrm{E}\mathrm{r}\mathrm{o}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}$, $\mathrm{d}\mathrm{i}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}$, $\mathrm{o}\mathrm{p}\mathrm{e}\mathrm{n}\mathrm{n}\mathrm{i}\mathrm{n}\mathrm{g}$ and $\text{Closing}$ morphology operations on voxel set $\boldsymbol{W}$ with $m(m\ge 1)$ times iteration respectively. The following describes the process of volume sampling point extraction.

Firstly, a hybrid filter composed of Gaussian filter and $\rho$-percentile gradient filter is employed on resampled CT images in order to enhance boundaries of bone tissue and fill potential holes.

Secondly, the initial voxel set ${\boldsymbol{W}}_{H0}^{k} = \{{H}_{m, n, l}\in \left[{H}_{0}, {\theta }_{max}\right], {v}_{m, n, l}\in \boldsymbol{W}\}$ is obtained, and ${\boldsymbol{H}}_{0}$ is chosen on the basis of retaining as many bone tissue regions as possible while minimizing the connection regions. If ${\boldsymbol{W}}_{H0}^{k}$ still includes voxel data from adjacent bones, erosion operation will be iteratively executed on ${\boldsymbol{W}}_{H0}^{k}$ until achieving the sub-region ${Erode}^{m}\left({\boldsymbol{W}}_{H0}^{k}\right)$ which is disconnected with adjacent bones, ${{ Erode}^{m}\left({\boldsymbol{W}}_{H0}^{k}\right)\subset \boldsymbol{W}}_{H0}^{k}$. However, because the erosion operation also erode correct regions while removing redundant connection regions between adjacent bones, the sub-region${Erode}^{m}\left({\boldsymbol{W}}_{H0}^{k}\right)$ should be applied dilation operation in order to get the temporary region ${\boldsymbol{W}}_{DE}$ which contain the eroded correct regions, ${{\boldsymbol{W}}_{DE} = Dilation}^{m}\cdot {Erode}^{m}\left({\boldsymbol{W}}_{H0}^{k}\right)$. Then, the bone's voxel data is obtained by operation ${\boldsymbol{W}}_{H0}^{k} = {\boldsymbol{W}}_{H0}^{k}\cap {\boldsymbol{W}}_{DE}$. Voxel set ${\boldsymbol{W}}_{Hi}^{k}$ corresponding to threshold value ${\theta }_{min} = {H}_{i}$ is also obtained by the same way, where ${H}_{i}\in \left[0, {H}_{0}\right]$, and${H}_{i} < {H}_{i-1}$.

Finally, opening and closing operations are applied on ${\boldsymbol{W}}_{Hi}^{k}$ in order to further improve the segmentation results: ${\boldsymbol{W}}_{Hi}^{k} = {Open}^{l}{\cdot close}^{n}\left({\boldsymbol{W}}_{Hi}^{k}\right)$. Map the voxel set ${\boldsymbol{W}}_{Hi}^{k}$ to the corresponding volume sampling points ${\boldsymbol{P}}_{Hi}^{k}$, and then, the initial and target volume sampling points ${\boldsymbol{P}}_{H0}^{k}$ and ${\boldsymbol{P}}_{HT}^{k}$ are obtained by ${\theta }_{min}$ = ${H}_{0}$ and ${\theta }_{min}$ = ${H}_{T}$ respectively.

2.2.   Bone surface points extraction

Because the volume sampling point ${\boldsymbol{P}}_{H}^{k}$ is a point set which contains not only the surface data, but also inner structure data of the bone, as well as noise data. It is difficult to reconstruct an accurate watertight 2-manifold bone surface model from ${\boldsymbol{P}}_{H}^{k}$ directly through geometric or non-geometric reconstruction algorithms. Further operations on ${\boldsymbol{P}}_{H}^{k}$ are needed in order to extract visible surface points ${S(\boldsymbol{P}}_{H}^{k})$ from ${\boldsymbol{P}}_{H}^{k}$ used to construct the bone surface model, ${S(\boldsymbol{P}}_{H}^{k})\subset {\boldsymbol{P}}_{H}^{k}$. Based on the principle of retinal spherical imaging, this paper proposed the visible surface point extraction method from volume sampling points.

When human observe objects, objects with shorter distance to eyes will block those with longer distance, and the visible regions are inverted and projected onto the spherical retina. According to the visual imaging principle of human eyes, the eye's imaging system can be simplified to that of the pinhole camera model when the imaging plane is replaced by spherical imaging surface of retina. This means that for object points on the same ray line emitted from the eye, the one closest to the eye is visible, and will be projected onto the retinal spherical imaging surface. Projections of the other invisible vertices fall inside the retinal spherical imaging surface.

When applying this mechanism to a point cloud with different observation locations, we can develop an algorithm of mapping visible surface points to the inside of a super-thin spherical shell with thickness $\delta \ll R$ while mapping invisible data points to the inside of spherical shell's inner surface, where $R$ is the outer radius of the shell.

2.2.1.   Hyper-spherical surface projection mapping

Imaging screen of human eye retina is similar to a concave spherical surface and the image formed through the eye vision system is concave. Spherical imaging surface of retina has better adaptability and wider view field compared with planar imaging plane. Based on the combination of the eye vision system and pinhole imaging model of camera, we propose the hyper-spherical projection mapping method to extract visible surface points, which is described as follows in Figure 4.

Firstly, the imaging plane of the pinhole imaging system is replaced by a spherical imaging surface corresponding to the concave retina in Figure 4(B) and (E).

Secondly, visible points of the object is inverted and projected near to the spherical imaging surface; invisible points inside the object are inverted and projected to the inside of the spherical imaging surface in Figure 4(B) and (E).

Thirdly, for the convenience of computation and analysis, the object and its corresponding flipped and magnified hyper-spherical projection mapping model are placed on the same side of the projection center (observation center) $\mathrm{O}$ in Figure 4(A) and (C), Figure 4(D) and (F). Thus, the hyper-spherical projection mapping system is constructed.

The demonstration of the hyper-spherical projection mapping system is shown in Figure 4. As the radius $R$ of the hyper-sphere increases in Figure 4(E), more visible surface points will be projected near to the hyper-spherical surface, and invisible points off the hyper-spherical surface with a relative large distance. The invisible points can be filtered automatically in the following convex hull reconstruction stage.

Based on pinhole imaging principle and the hyper spherical projection mapping system described above, the mapping between a point ${\boldsymbol{p}}_{i}\in \boldsymbol{P}$ and its projection to the hyper-spherical surface ${\boldsymbol{p}}_{i}^{\mathrm{*}}\in {\boldsymbol{P}}_{}^{\mathrm{*}}$ is linear, namely ${{\boldsymbol{P}}_{}^{\mathrm{*}} = \{{p}_{i}^{\mathrm{*}} = {f}_{O}\left(R, {\boldsymbol{p}}_{i}\right)|{\boldsymbol{p}}_{i}\in \boldsymbol{P}\}}_{}$, where $R$ is radius of the hyper sphere and $\mathrm{O}$ is the observation position. The function has characteristics that visible points are projected near the hyper spherical surface. Therefore, the mapping function can be written as

where ${p}_{i}^{\mathrm{*}}$ is the mapped point of ${p}_{i}$.

2.2.2.   Surface points extraction

Let $\widehat{\boldsymbol{P}} = {\boldsymbol{P}}_{}^{\mathrm{*}}\cup \mathrm{O}$. Then, the convex hull of $\widehat{\boldsymbol{P}}$ is calculated by using Qhull algorithm ^[40] in Figure 5(A) and (B). A point is visible if its mapping point is corresponding to a vertex of the convex hull, and then all vertices of the convex hull except $\mathrm{O}$ are corresponding to the visible surface points $S\left(\boldsymbol{P}\right)$ of $\boldsymbol{P}$ extracted from the given view direction in Figure 5(B) and (C). The convex hull varies with the hyper sphere's radius $R$ changing, and more visible points will appear as $R$ increases in Figure 4. When the value range of $R$ is $R\in [{10}^{3.5}, {10}^{3.8}]$ ^[35], the surface points can be extracted sufficiently.

When we observe a point cloud with simple shape, taking a vase for an example as shown in Figure 5, we can get a full view of the object from the left and right view directions, and most of the vase model's surface data can be extracted. For object with complex shapes such as bones, six or more view directions are necessary in Figure 6. After the visible points from different view directions are extracted and added together, the surface points $S\left(\boldsymbol{P}\right)$ are obtained.

Figure 6(C) shows the sub-surface patches corresponding to the visible points extracted from different view directions of the thighbone's volume sampling data. As can be seen from Figure 6(D), most of the surface points can be extracted from eight view directions. The optimal view direction with most visible points is calculated by PCA analysis method.

2.3.   Initial watertight bone surface model reconstruction

Let $M$ represent a watertight 2-manifold surface mesh model.$\boldsymbol{V} = ({\boldsymbol{v}}_{0}, \cdots {, \boldsymbol{v}}_{n})$ is the vertex set of $M$, $\boldsymbol{V}\in {\mathbb{R}}^{\mathrm{n}\times 3}$, ${\boldsymbol{v}}_{i}\in {\mathbb{R}}^{1\times 3}$.

The initial bone surface model ${M}_{0}$ reconstructed from surface points ${S(\boldsymbol{P}}_{H0}^{k})$ should be closed and watertight. Because ${S(\boldsymbol{P}}_{H0}^{k})$ extracted by method proposed in section 2.2 generally have problems such as local data loss, holes and noises, geometric reconstruction algorithms are not suitable for reconstruction of this kind data. In this paper, implicit surface reconstruction method is adopted to generate ${M}_{0}$ which can fill holes, restore local shapes with data loss and suppress noise existed in ${S(\boldsymbol{P}}_{H0}^{k})$.

Implicit surface reconstruction method uses spatial indicator function to describe the inner, outer, and boundary information of the model. It can reconstruct surface from complicated point clouds by simple function. Denote ${\chi }_{M}$ as the spatial indicator function of surface model $M$. ${\chi }_{M}>0$ indicates points inside the model $M$, ${\chi }_{M} = 0$ means points on boundary or surface of $M$, and ${\mathrm{\chi }}_{\mathrm{M}} < 0$ indicates points outside $M$. In other words, implicit surface reconstruction is to solve the indicator function ${\mathrm{\chi }}_{\mathrm{M}}$.

The implicit surface reconstruction algorithm based on Poisson equation combines both the advantages of global and local fitting methods, and the reconstructed surface shows more feature details and less noise. Let vector field $\stackrel{⃑}{\boldsymbol{U}}$ denote the oriented points of the sampling points with computed normal and ${\nabla \chi }_{M}$ represent the gradient of ${\chi }_{M}$. There is a corresponding integral relationship between $\stackrel{⃑}{\boldsymbol{U}}$ and ${\chi }_{M}$. Thus, the problem of surface reconstruction from point cloud $\boldsymbol{P}$ reduces to computing the indicator function ${\chi }_{M}$ of which the gradient ${\nabla \chi }_{M}$ best approximates $\overrightarrow{U}$.

After applying divergence operator on Eq (2), this variational problem becomes a standard Poisson problem:

Solution of Eq (3) can finally reduce to a well-conditioned sparse linear system. The goal of this section is to reconstruct a watertight and smooth initial surface model ${M}_{0}$, and therefore an octree of depth 7 or smaller is chosen in order to restore the contour shape efficiently at regions with data missing while approximating the visible surface points as much as possible.

Figure 7(A) and the enlarged view in Figure 7(C) show that there are large holes in the femur head and condyle regions. Figure 7(B) shows a watertight 2-manifold bone surface model is generated by Poisson surface reconstruction algorithm with holes filled and noise smoothed, but the details such as linea aspera are missing.

2.4.   Accurate surface reconstruction

Although the initial surface model ${M}_{0}$ reconstructed in section 2.3 is watertight and smooth, it is still have large shape deviation from the actual surface model at certain regions, and is not accurate enough to be used for medical analysis.

Denote a point set on the actual bone surface ${M}_{T}$ as ${\boldsymbol{B}}_{P}$, ${\boldsymbol{B}}_{P}\subset {\boldsymbol{P}}_{H}^{k}$. In order to reconstruct an accurate surface model $M$, ${M}_{0}$ need to be deformed in order to best approximate or interpolate ${\boldsymbol{B}}_{P}$ in the form of non-rigid deformation that follows the shape characteristics of bone, and the position deviation between ${M}_{0}$ and $M$ should be within a reasonable range. Generally, extra constraints are imposed during solving for $M$ in order to satisfy the spatial relationship between adjacent bones and also reduce the influence of noises existing in ${\boldsymbol{P}}_{H}^{k}$. The constraints usually cause the solution process to be non-linear and non-convex, which is difficult to achieve a global optimization result. To solve this kind of problem, local-global iterative surface deformation algorithms are used to make $M$ best approximate ${\boldsymbol{B}}_{P}$

where $i$ is the number of iterations and ${M}_{i}$ is the corresponding deformation model.

As described in section 2.1, when adopting thresholding method to extracted bone volume sampling points from CT images, the volume sampling points gradually cover the entire bone tissue space as expansion of the threshold range, and ${\boldsymbol{P}}_{HT}^{k}$ is used as the reference point set to generate the accurate surface model. But due to the characteristics of CT image data, it is inevitable that the sampling points ${\boldsymbol{P}}_{HT}^{k}$ are mixed with data from its own cartilages, inner structures and adjacent bones. As a result, it is difficult to extract such a surface point set ${\boldsymbol{B}}_{P}$ of which every point is on the actual bone surface, while the actual extracted surface point set ${\boldsymbol{B}}_{P}^{*}$ also includes extra noise ${\boldsymbol{B}}_{N}$, ${\boldsymbol{B}}_{P}^{*}\approx {\boldsymbol{B}}_{P}\cup {\boldsymbol{B}}_{N}$, where ${\boldsymbol{B}}_{N}$ is the noise point data. In this paper, we use the subset points ${\boldsymbol{B}}_{P}^{*}$ of ${\boldsymbol{P}}_{HT}^{k}$ near the initial surface as the deformation constraint point set, ${\boldsymbol{B}}_{P}^{*} = S\left({\boldsymbol{P}}_{H0}^{k}\right)\bigcup ({\boldsymbol{P}}_{HT}^{k}-{\boldsymbol{P}}_{H0}^{k})$.

Therefore, the accurate surface reconstruction transform into solving $M$ through iterative deformation operation from ${M}_{0}$, so that $M$ can best approximates the sub set ${\boldsymbol{B}}_{P}$ of ${\boldsymbol{B}}_{P}^{*}$, while the negative influences of noise ${\boldsymbol{B}}_{N}$ are suppressed.

A constraint energy function is designed to implement the iterative deformation: ${M}_{i} = f\left({M}_{i-1}, {\boldsymbol{B}}_{P}^{*}\right)$. The constraints include proximity matching constraint to implement the match between vertices of deformation surface model and the constraint point set ${\boldsymbol{B}}_{P}^{*}$, shape deformation constraint to deform and optimize the deformation surface model while preserving its local feature details, and shape sharpness constraint to penalize relative large displacement from the previous iteration step. The constraint function is expressed as follows

where $E\left(\boldsymbol{V}\right)$ is the total constraint energy, $\alpha$ and $\beta$ are coefficients, and ${E}_{Proximity}$, ${E}_{deformation}$ and ${E}_{sharpness }$ are constraint energies corresponding to proximity matching, shape deformation, and shape sharpness, respectively, which are described as follows.

2.4.1.   Proximity matching constraint energy

The proximity matching of vertex set $\boldsymbol{V}$ with the constraint point set ${\boldsymbol{B}}_{\mathrm{P}}^{\mathrm{*}}$ is obtained by minimizing the following constraint energy function

where ${\omega }_{i}$ is constraint weight of ${v}_{i}$ and ${\omega }_{i}>0$, $d\left({\boldsymbol{v}}_{i}\right)$ is calculated as the distance of ${\boldsymbol{v}}_{i}$ and its projection $\psi \left({\boldsymbol{v}}_{i}\right)$ on the bone surface implied in${\boldsymbol{B}}_{P}^{*}$. Formally, $\psi \left({\boldsymbol{v}}_{i}\right)$ is taken as moving ${\boldsymbol{v}}_{i}$ in the minimal way in order to achieve proximity matching.

Then Eq (6) becomes

where $\boldsymbol{V} = ({\boldsymbol{v}}_{1}, \dots, {\boldsymbol{v}}_{n})$, $\boldsymbol{\psi }\left(\boldsymbol{V}\right) = \left(\psi \right({\boldsymbol{v}}_{n}), \dots, \psi ({\boldsymbol{v}}_{\boldsymbol{n}}\left)\right)$, and $\boldsymbol{W}$ is a diagonal matrix. When ${\omega }_{i} = 1$, Equation (8) becomes

By minimization of Eq (9), the proximity projection of $\boldsymbol{V}$ can be solved, which match the constraint points ${\boldsymbol{B}}_{P}^{*}$ as closely as possible without deviating from $\boldsymbol{V}$ too far away.

2.4.2.   Shape deformation constraint energy

The purpose of the shape deformation is to make the shape of $M$ match with that of the actual bone surface as accuracy as possible. The selection of deforming strategy affects significantly on the final result. Because Laplacian surface editing can retain most of the original shape features and avoids introducing new overlaps when the surface is deformed, this paper uses Laplacian surface editing to design the deformation energy function.

Denote one-ring neighbors of ${\boldsymbol{v}}_{i}$ as ${N}_{i} = \left\{{\boldsymbol{v}}_{ij}\right|{\boldsymbol{v}}_{ij}ϵ\boldsymbol{V}, j\le {k}_{i}\}$, and ${k}_{i}$ is the vertex number of ${N}_{i}$. The coordinate vector from the weighted average coordinates of ${\boldsymbol{v}}_{i}$'s one-ring neighbors to ${\boldsymbol{v}}_{i}$ can be expressed as

where ${w}_{ij}$ is the weight corresponding to ${\boldsymbol{v}}_{ij}$ which is the $j$th one-ring neighbor of ${\boldsymbol{v}}_{i}$. ${\boldsymbol{\delta }}_{i}$ is also called the differential coordinate of absolute coordinate ${\boldsymbol{v}}_{i}$.Write Eq (10) in vector form while taking the weight value as one for convenience:

where ${\bf{\Delta }} = ({\boldsymbol{\delta }}_{1}, \dots, {\boldsymbol{\delta }}_{n})$, ${\boldsymbol{\delta }}_{i} = {\boldsymbol{v}}_{i}-(1/k)\sum {\boldsymbol{v}}_{ij}$, $\boldsymbol{A} = {\{0, 1\}}^{n\times n}$ is the adjacency matrix of $\boldsymbol{V}$ on surface $M$, and $\boldsymbol{D} = dig({k}_{1}, \dots, {k}_{n})$ is a diagonal matrix. $\boldsymbol{L} = \left(\boldsymbol{I}-{\boldsymbol{D}}^{-}\boldsymbol{A}\right)$ is called Laplacian operator which is translation-invariant, and the rank of matrix $\boldsymbol{L}$ is $\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\left(\boldsymbol{L}\right) = n-1$ for 2-manifold mesh model. For meshes which have fixed adjacency matrix, $\boldsymbol{V}$ can be obtained by solving linear system as long as one vertex is determined.

When deforming surface mesh ${M}^{\mathrm{*}}$ to obtain a new surface mesh $M$ in condition of keeping its adjacency matrix unchanged, the deformation constraint is described by the displacement vector ${\boldsymbol{c}}_{i}$ of vertex ${\boldsymbol{v}}_{i}^{\mathrm{*}}$ from ${M}^{\mathrm{*}}$ to $M$, ${\boldsymbol{c}}_{i} = {\boldsymbol{v}}_{i}-{\boldsymbol{v}}_{i}^{\mathrm{*}}$, with a deforming weight ${\omega }_{i}\in [0, 1]$. If ${\omega }_{i} = 0$, the vertex ${\boldsymbol{v}}_{i}^{\mathrm{*}}$ is called free-deform vertex, and its displacement completely depends on the adjacency matrix and those vertices whose weights are not zero. The vertices with non-zero weights control the deformation, which are named control vertices. As a result, the weight ${\omega }_{i}$ and displacement ${\boldsymbol{c}}_{i}$ determine the shape of the new surface. When modeling with Laplacian operator, the first thing to do is to determine the absolute coordinates of control vertices, i.e., ${\boldsymbol{v}}_{i} = {\boldsymbol{v}}_{i}^{\mathrm{*}}+{\boldsymbol{c}}_{i}, i\in \{m, n\}$, where $m < n$. Then calculate new coordinates of free-deform vertices. Solving for free-deform vertices is transform into solve the following deformation energy function

By minimizing the deferential coordinate deviations and the compliance with the positional constraints, the quadratic minimization problem becomes solving sparse linear equations as below

where ${\bf{\Omega }} = \mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\omega }_{1}, \cdots, {\omega }_{n})$ is a diagonal matrix composed of weights, and ${\boldsymbol{V}}^{\mathrm{*}}$ is initialized with ${\boldsymbol{V}}_{0}$ of ${M}_{0}$.

2.4.3.   Shape sharpness constraint energy

Because excessive large displacement will introduce unnatural and sharp shape deformation, the shape sharpness constraint energy is used to control the amount of vertex's relative displacement from previous step. The displacement can be penalized by using the following energy constraint function

where ${\boldsymbol{V}}_{pre}$ represents the position state in last iteration.

2.4.4.   Solution of deformation function

When the partial derivation of $E\left(\boldsymbol{V}\right)$ satisfies $\partial E\left(\boldsymbol{V}\right)/\partial \boldsymbol{V} = 0$, the total constraint energy $E\left(\boldsymbol{V}\right)$ in Eq (5) gets its extreme value, and then a linear system is obtained as follows

The direct solver for sparse matrix is used to solve Eq (15). During the computation, at earlier iteration stage, some of the displacements are relatively large, which will increase the unreliability of the corresponding projections on the implied surface of the constraint point set. In order to suppress the negative effects caused by ${E}_{Proximity}\left(\boldsymbol{V}\right)$, initial values for $\alpha$ and $\beta$ are set bigger to increase the influence strength of ${E}_{deformation}\left(\boldsymbol{V}\right)$ and ${E}_{sharpness}\left(\boldsymbol{V}\right)$, e.g., $\alpha = 1, \beta = 1$.

As the iteration steps increase, $\boldsymbol{V}$ gradually approaches to the target position, and the proximity match becomes more reliable. When solutions of two consecutive iterations are very close and result in little progress, values for $\alpha$ and $\beta$ are decreased to reduce the influence strength of ${E}_{deformation}\left(\boldsymbol{V}\right)$ and ${E}_{sharpness}\left(\boldsymbol{V}\right)$ on the solution and indirectly increase the weight of ${E}_{Proximity}\left(\boldsymbol{V}\right)$, in order to make the deformation result approximate the local feature details of the target surface more precisely.

Figure 8 shows an example of deforming an initial surface model in Figure 8(A) to the target model in Figure 8(B). As can be seen, local feature details can still be kept in regions where the displacement deviations between initial surface model and the reference model are remarkable, and the constraint of proximity matching plays a leading role during in the later deforming stage. Figure 8(E) and (H) show that the final deformation result approximates the target model accurately.

3.   Experiments

3.1.   Resampling and segmentation

Because of the narrow interosseous space between the hip and the caudal vertebrae, and the complicated HU values distribution in regions such as acetabular fossa, femoral head, acetabula notch and trochanter fossa, accurate segmentation and reconstruction of the femur and hip are very difficult and representative. Therefore, we take both the healthy and diseased CT images containing femur and hip as examples to analyze the proposed algorithm.

The slice thickness of CT image is usually between 0.625 mm and 10 mm. It has great influences on CT image segmentation and bone surface reconstruction. Big thickness not only tends to seriously blur interosseous regions that make it hard to segment adjacent bones separately with each other, but also results in obvious stair-step shape on the reconstructed surface in Figure 2(A). Therefore, in order to avoid the negative effects mentioned above, it is necessary to resample the spatial space of CT images by smaller thickness, which can achieve image enhancement, improve image clarity especially for regions at joint location, and alleviate the stair-step effect in Figure 2(B). Experiments show that when the resampling thickness is taken as $0.2~0.3\mathrm{m}\mathrm{m}$ for long bones and $0.1~0.2\mathrm{m}\mathrm{m}$ for vertebrae and short bones, better segmentation and reconstruction results are achieved.

According to biomedical characteristics of the bone to be segmented and application scenarios of the segmentation results, the HU threshold range used for bone's ROI segmentation from CT images is determined. The upper limit value of the threshold range usually remains unchanged. For the lower limit value, a bigger value is chosen to extract initial volume sampling points for reconstructing the initial surface model, and a smaller value for the target volume sampling points which is used as the reference point set for generating accurate surface model. Taking femur as an example in Figure 10, the threshold ranges [300-2500] and [150-2500] are chosen for surface points extraction and surface deformation respectively.

After resampling operation, for most cases, the adjacent bones can be effectively disconnected by applying once or twice consecutive erosion morphological operation on the interosseous regions. For some special cases, such as the threshold range of interosseous region overlaps too much with that of the bone's ROI, or the interosseous region is too small and narrow, we adjust the threshold range to that of the corresponding cartilage in order to identify and segment out the region containing cartilages, and some interaction may be required. After subtracting the cartilage region from CT images, the bone's ROI can be segmented correctly.

3.2.   Bone surface points extraction

Since the mapping points and the source points are matched one-to-one, inverse mapping of the convex hull reconstructed from the mapping points is corresponding to the local surface patch of the visible region on the source points. After projection operation, convex hull reconstruction and inverse mapping operation are sequentially applied on the source points from multi-view directions, and the visible surface points can be extracted efficiently. Generally, the object can be fully observed from six view directions including top, bottom, front, back, left and right. For model with complex shapes in Figure 6, additional view directions are added until the whole picture of the model is observed. The surface mesh patches reconstructed from the visible surface points usually contain noise especially at boundaries, which can be removed according to the smoothness and connection characteristics of surface model.

When the visible surface points extracted from multi-view directions of the bone's initial volume sampling points are fitted by implicit function in Figure 7, the depth of the octree used in Poisson surface reconstruction determines the details of the reconstructed model. The bigger the depth is, the more geometric details will be reconstructed. On the contrary, the smaller the depth is, the more geometric details will be smoothed out. Smaller depth makes the Poisson surface reconstruction procedure more effective while removing noise and filling holes, and then a smooth, closed and continuous surface model is achieved.

3.3.   Watertight surface reconstruction

The initial watertight surface model ${M}_{0}$ is reconstructed from the visible surface points, so the geometry size of ${M}_{0}$ is usually smaller than that of the actual bone surface model. The subsequent deformation modeling behavior is similar to expanding ${M}_{0}$ outward in Figure 11(K).

The deformation result mainly depends on the sampling points that are distributed near the bone's actual surface, and the sampling points inside the surface have little contribution to the deformation result. The constraint points ${\boldsymbol{ }\boldsymbol{B}}_{P}^{*}$ used in the iterative deformation stage are achieved through set operations ${\boldsymbol{B}}_{P}^{*} = S\left({\boldsymbol{P}}_{H0}^{k}\right)\bigcup ({\boldsymbol{P}}_{HT}^{k}-{\boldsymbol{P}}_{H0}^{k})$.The reference points ${\boldsymbol{P}}_{HT}^{k}$ with lower accuracy requirements also can be obtained by using binary mask to 3D mesh reconstruction algorithm in ITK.

As shown in Figure 8, when the tooth model in Figure 8(A) is deformed from the initial to the final state, the vertex's movement direction of the deformation model is determined by its relative position to the target model, e.g., moving inward (corresponding to blue region), moving outward (corresponding to red region) and fixed (on the surface of target model) respectively.

Directional search is performed to obtain the matching pairs $({\boldsymbol{v}}_{i}, \psi ({\boldsymbol{v}}_{i}\left)\right)$ with given search radius $R$. Vertices in the matching pairs are used as control points, and those beyond the search radius are free deformation points. The directional search center ${\boldsymbol{v}}_{c}$ corresponding to ${\boldsymbol{v}}_{i}$ is determined by offsetting a distance$\delta$ along the normal direction, ${\boldsymbol{v}}_{\boldsymbol{c}} = {\boldsymbol{v}}_{i}+\stackrel{⃑}{\boldsymbol{n}}*\delta$, when moving outward $\delta \in \left[0.5{d}^{*}, {d}^{*}\right]$, otherwise$\delta \in \left[-{d}^{*}, {-0.5d}^{*}\right]$.The range of the search radius $R$ is usually taken as$R\in [{d}^{*}, \sqrt{3}{d}^{*}]$, which can ensure that the deformation amplitude is controlled in reasonable range. ${d}^{*}$ is the average edge length of the deformation surface model.

During the iterative deformation process, the weights are assigned with ${\omega }_{i} =$1 and ${\omega }_{i} =$0 for vertices in control region and free deformation region, respectively. After each iteration, ${\omega }_{i}$ will be updated with 1 for vertex with little or no movement. When the deformation result reaches a stable stage, all weights are updated with 1. Let $N$ denote the number of total iterative steps. Good results are achieved when the value of $N$ is taken between ^[6,10]. The deformation vector ${\boldsymbol{c}}_{i}$ is defined as follow:

where $j$ is the current iteration step.

As shown in Figure 9, the spherical point cloud is obtained by union of sampling points of spherical surfaces with radiuses of 100,120,140 and 160 mm. The spherical point cloud are copied and then translated by 300 mm and 190 mm along axis X to obtain the local overlapped point cloud as shown in Figure 9(A) and Figure 9(B) respectively. The points in the overlapped region are distributed irregularly. The sphere with radius of 100 mm in Figure 9(A0) and (B0) are iteratively deformed outwards by the proposed method. After performing once, twice and thrice iterative directional deformation operation, the corresponding deformation results are obtained in Figure 9(A1)-(A3); Figure 9(B1)-(B3). The radii of spheres in Figure 9(A1)-(A3) are approximately equal to 120,140 and 160 respectively, and the same for Figure 9(B1)-(B3).

Figure 8 and Figure 9 show that the deformation algorithm used in this paper has the characteristics of directional deformation, and the deformation amplitude, direction and result are controllable.

As shown in Figure 8(I), the shape deviation between the deformation result and the target model with a maximum of 0.024 mm is smaller than the average edge distance 0.0256mm of the target model.

Figure 9 shows that when the sphere surface is deformed and approaches to the surface data of the point cloud, the shape of the sphere surface is well preserved under action of the constraint energy function. The proposed reconstruction algorithm has the ability to control the deformation scale within reasonable range, and will not produce extra noise or excessive deformation. The maximum shape deviation 9.214 mm in Figure 9(D) between the deformation result and its corresponding standard fitting sphere surface is smaller than the local average sampling density 11.305 mm of the point cloud's irregular region.

Because the proposed algorithm is iterative and the deformation process are controllable, the iterative steps, direction, scale and area of the deformation can be specified or changed during the reconstruction process as needed. As can be seen from Figures 8-11 that the deformation algorithm proposed in this paper can detect the accurate surface shape from the constraint point set. The mechanism of the reconstruction algorithm can be explained theoretically and experimentally from the constraint energy function in Figures 8-11.

4.   Results and discussion

4.1.   Shape analysis of reconstruction results

Figure 10 shows the surface reconstruction procedure of a femur, and Figure 11 for a hip bone. The CT images are resampled with thickness value of 0.3mm. The segmentation threshold ranges used to extract the initial volume sampling points and the target volume sampling points are [300,2500] and [150,2500] respectively. The initial surface models in Figure 10(C) and Figure 11(C), and the cross section models in Figure 10(J) and Figure 11(H) are all generated by Poisson surface reconstruction algorithm. The octree depth used to reconstruct the surface models in Figure 10(J) and Figure 11(H) is 10. The surface of the reconstructed model from the target volume sampling points indicates where the actual bone surface is located in Figure 10(L) and Figure 11(J).

The enlarged view in Figure 10(F) shows that the stair-step surface model is obtained without considering the shape deformation constraint energy ($\mathrm{\alpha } = 0$) in the last iteration step. As can be seen from the stair-step shapes at the femoral head, which are consistent with the 3D slice contours of the femoral head directly extracted from CT images in all three directions, the proposed proximity matching algorithm can detect the data of bone surface correctly from the volume data. The accurate femur surface model in Figure 10(G) is obtained when the iterative deformation reaches a stable state. Enlarged view in Figure 10(G) shows that the linea aspera of the femur is restored clearly.

From the cross section contours in Figure 10(J) and Figure 11(H), it can be seen that there are local shape missing, abnormal and redundant shapes in the model reconstructed directly from the corresponding volume sampling points by Poisson surface reconstruction.

Cross section contours in Figure 10(K) and Figure 11(I) show that the initial bone surface model has large differences in shape compared with the actual shape contained in the corresponding Poisson reconstructed surface model, especially at the region circled with blue border rectangle.

Figure 10(L) and Figure 11(J) are the generated accurate femur and hip bone surface models respectively. As can be seen from the enlarged views in Figure 10(L) and Figure 11(J), the cross-section contour curve of the accurate bone surface model correctly approximates the actual contour curve contained in the Poisson reconstructed surface model while avoiding the influences of noise data effectively. The local features of the femur such as trochanteric fossa, linea aspera are also restored correctly in Figure 10(G) and (L).

The cross section contour comparisons in Figure 10(L), (M) and Figure 11(J), (K) show that the proposed surface reconstruction algorithm in this paper can effectively reduce the influences of noise data, and it can also restore the local shape with data loss and reconstruct the feature details accurately. The premise is only to roughly segment out bone's ROI from CT images, in order to generate a watertight 2-manifold initial bone surface model.

In order to get an accurate bone's ROI disconnected from its adjacent bones, the resampling thickness should be smaller than the smallest interosseous space or thickness of cartilages between adjacent bones. Because the interosseous spaces corresponding to long bones and hip bones are relatively obvious, a bigger resampling thickness can be chosen. The interosseous spaces between adjacent short bones or vertebrae are much smaller compared to long bones or hip bones, so it needs a smaller resampling thickness to effectively segment the bone's ROI from CT images.

However, because there are non-obvious segmentation boundaries between adjacent sub-bones of the skull, the segmentation of sub-bone's ROI is very difficult. The segmentation and reconstruction of the skull's sub-bones is the goal of future research.

4.2.   Computational complexity

The proposed surface reconstruction algorithm mainly consists of three parts: volume sampling point extraction, visible surface point extraction, initial surface model reconstruction, and accurate surface generation.

Both CT images resampling and segmentation by thresholding have linear time $O(N$) complexity. The time complexities corresponding to hyper-spherical projection and convex hull reconstruction included in surface data point extraction algorithm are $O\left(N\right)$ and $O\left(NlogN\right)$ respectively.

Poisson surface reconstruction algorithm has $O\left(Nlogn\right)$ time complexity with $\le N/{2}^{lnN-{\rm H}+1}$, where $N$ is the points number, ${\rm H}$ is the depth of the octree. The depth ${\rm H}$ used to reconstruct the initial surface model is chosen as 7 in this paper, so it greatly reduces the computational time taken for the surface reconstruction.

For the accurate surface generation algorithm, the most time-consuming items are the proximity matching pairs' searching and solving of the sparse linear system of the deformation algorithm with a $O\left(N\right)$ time Cholesky factorization. The $O\left(Nlogk\right)$ time KNN (k-nearest neighbors) algorithm ^[41] is used to search for matching pairs in this paper.

Therefore, the corresponding time complexity of the proposed reconstruction algorithm can be measured by combinations of $O\left(NlogN\right)$ and $O\left(N\right)$. The number of iteration steps required for deformation algorithm mainly depends on the density and search radius of the constraint point set. Generally, stable results can be achieved after 5-10 iterations. The reconstruction results show that the proposed method is accurate and efficient.

4.3.   Comparison and quantitative analysis

The 3D bone surface model reconstruction algorithm of CT images based on thresholding and implicit surface fitting has characteristics of strong stability and robustness compared with other algorithms, and it is the mainstream algorithm used by 3D medical image processing software such as the most popular Materialise Mimics.

Figure 12(A) shows a patient's local CT image with femoral tumor, which includes a healthy right hip bone, a healthy right femur, a diseased left hip bone, and a left femur with bone tumor and fracture. The slice thickness, slice increment and pixel size of CT image in Figure 12(A) are 1mm, 0.4mm and 0.39mm respectively. The diseased region and shape size of bone tumor are usually uncertain and random. Bone tumor can lead to abnormal distribution of HU values in the diseased region, which will cause the segmentation boundary between bone tumor and its surrounding soft tissue being blurred and unclear (see the enlarged view of Figure 12(A)). The irregular shape and complicated HU values distribution of bone tumor will bring great inaccuracy to the algorithm based on differential geometric properties, which are calculated according to HU values. So, we use CT images indicated in Figure 12 (A) to compare and analyze the reconstruction effects between the proposed algorithm and Materialize Mimics software.

The statistical threshold ranges of compact bone, spongial bone and muscle tissue are [580,2500] and [145,2500], [-25,140] respectively. Because the lower limit of spongial bone and upper limit of muscle tissue are very close. It is hard to classify pixels or voxels with HU value between [130,160]. Therefore, lots of holes and segmentation boundary missing problems will appear when segmenting CT images by thresholding, especially in regions distributing with more spongial bone tissue. Figure 12(B)-(D) are segmentation results corresponding to threshold range of [150,2500], [226,2500] and [330,2500] respectively. Figure 12(E)-(G) show the cross section contour curves of the 3D bone models reconstructed directly by Mimics materialize (version 20.0) after thresholding segmentation. As can be seen from Figure 12, while the contour curves of the reconstructed model are faithful to the segmentation data, there are lots of abnormal shapes, holes and internal redundant structure in the reconstructed bone models, which can't be used directly for subsequent artificial prosthesis design, internal honeycomb structure design and FEM (Finite Element Method) analysis.

In Materialize Mimics (version 20.0), the watertight 3D surface model reconstruction of CT images is realized by the combination method of "thresholding" (Figure 13(B)), "Calculate part" (Figure 13(C)), and "Wrap" (Figure 13(D)-(F)), and the watertight surface reconstruction procedure and purpose is similar to that of this paper.

Figure 13 shows the reconstruction procedure of the watertight 3D bone surface model by Materialize Mimics software (Version 20.0). Figure 13(B) shows the segmentation results with threshold range of [226,2500], and the colorful region corresponding to each bone are separated with adjacent bones manually. Figure 13(C) shows the 3D bone models reconstructed directly from the segmentation results in Figure 13(B) by Mimics. Figure 13(D)-(F) are corresponding to the wrapped results of Figure 13(C) with parameters "smallest detail = 1mm" and "Gap closing distance = 2, 8, 15mm" respectively. Figure 13(G)-(I) are the cross section contour curves corresponding to Figure 13(D)-(F), and the intersection plane marked is the plane of CT slice in Figure 13(A).

Enlarged views of Figure 13(G)-(I) at the same region of the CT image show the comparison between the red cross section contour curves and the blue manually selected actual contour curves. As can be seen from Figure 13(J), with the increase of "Gap closing distance", the contour curves are evolved from multiple rings to a single closed contour curve with holes being filled and gaps being closed, and a watertight 3D bone surface model can be achieved when the contour curve corresponding to each CT slice is a single closed curve in Figure 13(F) and (I).

Figures 15-18 show the watertight surface wrapping procedure of the right femur, the right hip bone, the left femur, and the left hip bone respectively. The red curves are cross section contour curves, and the blue curves represent the manually selected actual contour curves. As can be seen from the enlarged views in Figures 15-18, while the contour curves are evolved as the "Gap closing distance" increases, the holes, internal redundant structures and gaps are gradually disappeared, but local surface shape in regions such as trochanteric fossa (Enlarge view in Figures 14 and 16) and acetabula notch (Enlarge view in Figures 15 and 17) are also over deformed in the direction of increasing deviation error from the correct position. When "Gap closing distance" = 15mm, the wrapped result for each bone are watertight.

Because the wrap algorithm relies on the geometric information of the bone model to be wrapped, rather than the actual CT segmentation boundary information, it will introduce new shape errors which have the same quantity level with the value of "Gap closing distance"' while closing the shape. For bone models in Figures 15-18, Table 1 show that the maximum shape errors of the wrapped watertight surface model are usually between 5mm and 15mm, which will cause assembly accuracy problems between the remaining healthy bones after tumor resection and the corresponding artificial prosthesis.

Figure 18 shows the watertight 3D bone surface model reconstruction by the method proposed in this paper. As can be seen from Figure 18(G) and (H), the watertight and smooth bone surface models are obtained while local details are well preserved.

Figure 19 shows the red cross section contour curves of Figure 19(I)-(L) for each bone in Figure 18(F) and the corresponding manual selected actual contour curves in Figure 19(M)-(P). Figure 19(A)-(D) indicate the location of the intersection plane in Figure 19(U) for each bone respectively. The enlarged views in Figure 19(Q)-(T) show comparisons between the cross section contour curve and the actual contour curve.

Figure 20 shows the quantitative analysis of the deviation error between the cross section contour curve and the manually selected actual contour curve. Figure 20(A) and Figure 20(B) show that the maximum deviation of both the healthy right femur and the healthy right hip bone is less than 1 mm, and the deviation of 95% sampling points is less than 0.5mm ($\approx 1.28$pixel size). Figure 20(C) show that for the left diseased hip bone and the left femoral tumor, there are about 1~2% sampling points with deviation error lager than 1 mm, deviation error of 85% sampling points is less than 0.5mm. Figure 20(B) shows the deviation curve comparison between the diseased left and healthy right hip bone. Because of irregular changes of local shapes and HU values corresponding to the diseased bone, deviation curve of the diseased hip bone changes much sharply than that of the healthy hip bone.

For the cases used in this paper, statistical results demonstrate that the deviation error of 98% sampling points for all bones is less than 1mm (2.5pixels), and the reconstruction accuracy is coincide with CT slice thickness (1.0 mm). The maximum deviation is only 1.3 mm, which is far less than that of the reconstruction results by Mimics in Table 1.

In fact, the volume sampling data corresponding to CT images are regular point set with equidistance between adjacent points. The search radius of the deformation algorithm is easy to be determined according to the CT slice thickness. Therefore, the deformation amplitude of each iteration can be controlled within the CT slice thickness.

4.4.   Limitations

The premise of the bone surface reconstruction algorithm proposed in this paper is to firstly segment the bone's ROI from CT images in order to achieve the initial volume sampling points, which are used to reconstruct the initial watertight surface model. When dealing with bones that have non obvious segmentation boundary in CT images such as sub-bones of the skull or sternum, it is difficult to segment sub-bone's ROI from CT images directly by thresholding method, and a more efficient and accurate segmentation method needs to be designed.

The proposed algorithm is designed for surface model reconstruction. The reconstructed model doesn't include the bone model's inner honeycomb structure. Therefore, to reconstruct a bone model with detailed inner structure, the algorithm needed to be redesigned based on volumetric mesh model

5.   Conclusions

Aiming at the problems that exist in the 3D bone surface model reconstructed from CT images, such as noises, holes, stair-step shapes, abnormal and redundant shapes, this paper proposes a watertight 2-manifold 3D bone surface reconstruction method. The proposed method includes three main steps: two-step thresholding, initial watertight surface reconstruction and shape optimization. From the experimental results of the cases used in this paper, the proposed algorithm can achieve good reconstruction result for bone with obvious interosseous space. For special cases, some manual interaction is still needed to achieve the desired deformation results.

Acknowledgments

Funding: This work was supported by the National Natural Science Foundation of China [grant number 51705183].

Conflict of interest

The authors declare there is no conflict of interest.