1. Introduction
Three-dimensional CT imaging technology is widely used in industrial detection and medical diagnosis due to its advantages of no damage and high resolution [1]. However, due to the influence of non-ideal factors (detector defects, noise, scattering, polychromatic rays, etc.) in industrial computed tomography (ICT) imaging system and imaging process, CT images generally contain scattering artifacts [2], beam hardening artifacts [3], noise [4], intensity inhomogeneity, boundary blurring, high data dimension and complex structure, which make it difficult for traditional segmentation methods to take into account the robustness and accuracy of segmentation.
The main challenge of image segmentation is that the segmented image usually contains noise, intensity inhomogeneity and artifacts. Most of the existing segmentation methods are two-dimensional, which can be divided into threshold-based segmentation method [5], region-based method [6], deep learning-based method [7] and partial differential equation-based method [8]-[10]. The segmentation method based on threshold and region is sensitive to noise. The deep learning model is particularly superior in processing volume data with complex structure and variable morphology. In order to improve the segmentation accuracy of brain tumors, Daimary [11] proposed a hybrid model combining multiple convolutional neural network architectures, such as U-Net and ResNet. In order to obtain a flexible and effective brain tumor segmentation system, Ranjbarzadeh [12] proposed a simple and efficient cascaded convolutional neural network (C-ConvNet/C-CNN). The model can extract local and global features on two different routes, but this method requires a large amount of dual-label data and a long time for pre-training. There are many types of industrial parts, and the structure is complex, so it is difficult to obtain an ideal training set. The active contour model (ACM) [13] does not require labeled data to guide segmentation, so it takes less time. With continuous exploration, ACM has become one of the most widely used and effective methods.
Osher and Sethian first proposed the level set (LS) method [14], and this method can be extended to higher dimensions, which provides a solid foundation for volume data segmentation. Volume data segmentation is to extract the target area with specific semantics from complex three-dimensional data. In the field of medical imaging, accurate organ or lesion segmentation is the basis of disease diagnosis and treatment planning. In the field of industrial inspection, component segmentation of 3D CT volume data is very important for product defect detection and quality evaluation. The numerical calculation of the level set method can be performed on a fixed Cartesian grid, which greatly promotes the further development of the level set method and its application in image segmentation [15]. For example, He et al. [16] proposed an improved three-dimensional edgeless active contour model (Chan-Vese model, referred to as CV model), which extended the level set method to three-dimensional space and applied it to the segmentation of CT volume data. This has brought great development to volume data segmentation based on level set. In order to achieve high-precision segmentation of ICT volume data, Huang et al. [13] proposed a level set segmentation method based on three-dimensional multiplicative additive model and achieved success.
The level set method can naturally handle topology changes and optimize boundaries, but it is sensitive to initialization. If the position of the initial contour is far away from the target area, it will increase the number of iterations and even lead to local optimum. The unreasonable form of the initial level set function will affect the stability and computational efficiency. The shape is greatly different from the target, and it is easy to miss details or over-segment. Therefore, designing a segmentation framework that integrates statistical modeling capabilities and boundary optimization capabilities is of great value for improving the segmentation performance of 3D volume data.
Markov random fields (MRFs) are widely used in image segmentation [17]. Its success is largely due to efficient algorithms, such as the combination of “data faithfulness” and “model smoothness” [18]. It is generally assumed that the intensity distribution of the region to be segmented in the image obeys the Gaussian distribution. Compared with the single Gaussian distribution, the Gaussian mixture model (GMM) method is more powerful in simulating the complex intensity distribution of the image. Wang et al. [19] proposed a GMM-based HMRF model (GMM-HMRF), which uses Gaussian mixture model to adapt to complex gray distribution, and realizes the alternating optimization of parameters and labels through EM algorithm, which improves the segmentation robustness of 3D volume data. However, the segmentation accuracy of this method at fine boundaries still needs to be improved, and it is difficult to meet the needs of high-precision applications.
In order to achieve accurate segmentation of volumetric data with intensity inhomogeneity, fuzzy boundaries, and noise interference, we propose a 3D CT image segmentation method based on HMRF-EM and level sets. First, the initial segmentation of the volumetric data is obtained using HMRF-EM, which is then used as the initial contour for the level set. Finally, a refined level set segmentation is performed, enabling the segmentation boundaries to fit the target region more precisely, thereby combining the rationality of the initial segmentation with the accuracy of the refined segmentation.
The organization of this paper is as follows. In Section 2, we introduce the related work of this paper and the proposed method. Section 3 presents the experimental results of our model. Section 4 is the conclusion.
2. The Proposed Model
2.1. HMRF-EM Framework
The HMRF-EM framework was first proposed for magnetic resonance brain image segmentation. Now, we use a segmentation method that combines the HMRF-EM framework with a Gaussian mixture model to segment three-dimensional volumetric data, and use it as the initialization for the level set method. First, we assume that the image is a two-dimensional grayscale image, and the intensity distribution of each region to be segmented follows a Gaussian distribution. Let the image
, where
is the number of pixels, and each
is the gray-level intensity of the pixel. we assume that each pixel corresponds to a specific label, i.e.
, where
,
is the set of all possible labels. Normally, in the binary segmentation problem, we label
. Henceforth, according to the MAP guideline, we are looking for a label that satisfies:
(1)
The prior probability
here is the Gibbs distribution, and the joint likelihood probability is:
(2)
where
is a Gaussian distribution with the parameter
.
In order to find the label
, we use the EM [20] algorithm to estimate the parameter set
, and the parameter set
and the label
are learned alternately.
Let
is a Gaussian distribution function with the parameter
:
(3)
We assume the prior probability is:
(4)
where
is a normalized constant, making the sum of
of all label configurations 1, where
is the set of all possible label configurations, and
is a priori energy function in the form of:
(5)
where
is the clique potential and
is the set of all possible cliques.
In the image domain, we assume that one pixel has at most 4 neighbors: the pixels in its 4-neighborhood. Then the clique potential is defined on pairs of neighboring pixels:
(6)
where
(7)
Thus assuming
(8)
where
represents the GMM likelihood probability of a single voxel,
is also a normalized constant. Given
and
, the likelihood energy is cast as:
(9)
Under these assumptions, the HMRF-EM algorithm is as follows:
1) First, use the kmeans method to initialize the parameter set
.
2) At the
-th iteration, calculate the likelihood distribution
.
3) The label is estimated using the current parameter set
by the above energy estimation method, that is, the maximized posterior probability is equivalent to the minimized total energy, which is achieved by the negative logarithmic transformation of the probability, combined with the definition of a priori probability and likelihood probability, we have:
(10)
4) Calculate the posterior distribution for all
and all pixels
using the Bayesian rule:
(11)
where
is the neighborhood configuration of
, and
(12)
Note that the posterior probability of label
given the neighborhood features
is defined as:
(13)
5) Use
to update the parameters:
(14)
(15)
In the previous assumption, the intensity distribution of each region to be segmented follows a Gaussian distribution with parameters
. However, this assumption is insufficient to capture the complexity of object intensity distributions in real-world scenarios. In contrast, a Gaussian Mixture Model (GMM) is more powerful than a single Gaussian distribution for modeling complex distributions. A Gaussian Mixture Model with
components can be represented by the parameters:
(16)
GMM now has weighted probabilities:
(17)
In our experiment, we suppose
, but we turn the m step of the em algorithm [20] into a Gaussian mixture model fitting problem. In fact, the GMM fitting problem itself can also be solved with the EM algorithm. In E-step, we determine which data should belong to which Gaussian component; In M-step, we recalculate the GMM parameters. By integrating the grayscale distribution modeling ability of the Gaussian hybrid model and the spatial constraints of the hidden Markov random field, the initial coarse segmentation of the three-dimensional volume data is realized.
2.2. Level Set Segmentation Model
Li et al. [21] proposed a variational level set framework for segmenting images and correcting intensity inhomogeneity. Now we consider a 3D grayscale image
on a given image domain
, for a voxel
in the image
, we define the local energy as:
(18)
is a kernel function with a standard deviation of
, truncated by a Gaussian function defined by:
(19)
where
is the radius, the spherical neighborhood centered on
.
is a normalization constant that makes
.
The local energy
can be equivalently expressed as the following calculation function in actual numerical calculations:
(20)
where
is the convolutional operation,
.
We define the following global energy as a data fitting term:
(21)
Let
be the level set function. We consider the image domain
being divided into two disjoint regions
and
. In this case, the regions
and
are represented by their membership functions
and
, respectively, where
is the Heaviside function, and
is the derivative of
. Defined as follows:
(22)
(23)
Then the energy in (21) can be expressed in the following level set formulation:
(24)
Two regularizations
and
are introduced:
(25)
(26)
where
.
In summary, the variational level set framework is:
(27)
where
and
are constants greater than zero. The above equation is minimized by the gradient descent method
(28)
Function
.
Fixing
and
, so that the energy
is the least optimal
is defined
(29)
where
.
Fixing
and
, so that the energy
is the smallest optimal
defined as:
(30)
where
and
.
2.3. Implementation
As we all know, the level set is sensitive to initialization. Now we directly map the coarse segmentation label matrix
to the initial implicit function
of the level set, avoiding the subjectivity of the traditional level set to manually set the initial contour. The level set function is initialized as follows :
(31)
where
is a constant,
represents the background region,
represents the target region, and
is an arbitrary voxel of the image domain
.
The main steps of the proposed model are as follows :
1) Find the label
.
2) Based on the label result, initialize the level set function
using the expression (31).
3) Evolve the level set using Equation (28), and then update
and
according to Equations (29) and (30).
3. Experimental Results
In this section, we will verify the ability of our method to segment ICT volume data. In this paper, let
, the value range of the length parameter is
, and the parameter of the Gaussian kernel function is
. Naturally, the length parameter
can be adjusted according to the types of segmented images. The general rule of setting
is: when you need to extract a smaller object, take a smaller value; if you need to segment a larger target, take a larger value. When the image is seriously uneven, a smaller
is used.
Tubular precision components, such as cylinders, are crucial in automobile manufacturing. Therefore, combined with the proposed segmentation method, accurate point cloud data can be obtained, making the error smaller. The car engine CT data and teapot CT data used in this paper are provided by the open scientific visualization data set (https://klacansky.com/open-scivis-datasets/category-ct.html).
Figure 1 shows the two-dimensional segmentation results of the proposed model for ICT volume data slice 27, slice 70, and slice 99 of engine cylinder samples (image size is 256 × 256 × 100). The segmentation details of these three slices are shown and compared with the segmentation results of HMRF-EM [19] (first column) and Li [21] (second column). In fact, our model can accurately identify the engine cylinder, depict the edge of the object, and detect tiny holes and thin walls to achieve superior segmentation performance. Figure 3(a) is the three-dimensional display of the engine segmentation results.
Figure 2 shows the two-dimensional segmentation results of the proposed model for the ICT volume data slice2 and slice10 of the teapot sample (the image size is 256 × 256 × 11), showing the segmentation details of the two slices. The first column represents the initial contour (HMRF-EM), the second and third
Figure 1. Segmentation results of slices 27, 70, and 99 of the engine cylinder ICT volume data (image size 256 × 256 × 100).
Figure 2. The segmentation process of the proposed model (image size is 256 × 256 × 11).
Figure 3. The 3D display of the segmentation results of simulated engine cylinder and teapot ICT volume data.
columns represent the segmentation in the intermediate process, and the fourth column shows the final evolution result of the LS function. It can be seen that the initial contour does not fully fit the edge of the teapot, but as the number of level set iterations increases, the segmentation results gradually become clear and tend to be stable. Figure 3(b) shows the three-dimensional display of the teapot segmentation results.
4. Conclusion
In this paper, a two-stage 3D volume data segmentation method combining GMM-HMRF and level set is proposed. The initial coarse segmentation with spatial continuity is realized by GMM-HMRF, which provides a reliable initial contour for the level set. Then, the optimized level set model is used to realize the fine optimization of the boundary, which effectively solves the problem of insufficient segmentation accuracy of the traditional method in the gray uneven and fuzzy boundary data. Experiments on real industrial CT data show that the proposed method is superior to the existing comparison methods in segmentation accuracy, robustness and efficiency. It can accurately extract the target area in three-dimensional volume data and provide high-quality segmentation results for subsequent tasks in the fields of medical image analysis and industrial detection.
Future work can further explore the multi-scale GMM-HMRF model to adapt to more complex grayscale distribution; at the same time, the parameter adaptive strategy of the level set energy function is optimized to improve the generalization ability of the algorithm in different types of volume data.
Acknowledgements
This work was supported by the high-level course “Algorithms and Data Structures” at Sichuan University of Science and Engineering.