3D CT Image Segmentation Based on HMRF-EM and Level Set

Abstract

Aiming at the problems of intensity inhomogeneity, boundary blurring and noise interference in the segmentation of three-dimensional volume data (such as medical images and industrial CT data). In this paper, the hidden Markov random field based on Gaussian mixture model is used for the initial segmentation of volume data, and the statistical and spatial information of the image is effectively used to obtain a reasonable initial contour and region division. Then, the level set method is used for fine segmentation to accurately capture the complex boundary shape in the image, which has good adaptability to the change of the target boundary and makes the segmentation boundary more accurately fit the target area. This improvement takes into account the rationality of the initial segmentation and the accuracy of the fine segmentation. The experimental results verify the effectiveness of the proposed method in image segmentation.

Share and Cite:

Liang, D. L. and Liu, C.J. (2026) 3D CT Image Segmentation Based on HMRF-EM and Level Set. Open Access Library Journal, 13, 1-12. doi: 10.4236/oalib.1114708.

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 Y=( y 1 ,, y N ) , where N is the number of pixels, and each y i is the gray-level intensity of the pixel. we assume that each pixel corresponds to a specific label, i.e. X=( x 1 ,, x N ) , where x i L , L is the set of all possible labels. Normally, in the binary segmentation problem, we label L=0,1 . Henceforth, according to the MAP guideline, we are looking for a label that satisfies:

X * = argmax X { P( Y|X,Θ )P( X ) } (1)

The prior probability P( X ) here is the Gibbs distribution, and the joint likelihood probability is:

P( Y|X,Θ )= i P( y i |X,Θ ) = i P( y i | x i , θ x i ) (2)

where P( y i | x i , θ x i ) is a Gaussian distribution with the parameter θ x i =( μ x i , σ x i ) .

In order to find the label X * , we use the EM [20] algorithm to estimate the parameter set Θ={ θ l |lL } , and the parameter set Θ and the label X are learned alternately.

Let G( z; θ l ) is a Gaussian distribution function with the parameter θ l =( μ l , σ l ) :

G( z; θ l )= 1 2π σ l 2 exp( ( z μ l ) 2 2 σ l 2 ) (3)

We assume the prior probability is:

P( X )= 1 Z exp( U( X ) ) (4)

where Z= Xχ exp( U( X ) ) is a normalized constant, making the sum of P( X ) of all label configurations 1, where χ is the set of all possible label configurations, and U( x ) is a priori energy function in the form of:

U( X )= cC V c ( X ) (5)

where V c ( X ) is the clique potential and C 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:

V c ( x i , x j )= 1 2 ( 1 I x i , x j ) (6)

where

I x i , x j ={ 0, if x i x j 1, if x i = x j (7)

Thus assuming

P( Y|X,Θ )= i P( y i | x i , θ x i ) = i G( y i ; θ x i ) = 1 Z exp( U( Y|X ) ) (8)

where G( y i ; θ x i ) represents the GMM likelihood probability of a single voxel, Z = Y exp( U( Y|X ) ) is also a normalized constant. Given Y and Θ , the likelihood energy is cast as:

U( Y|X,Θ )= i U( y i | x i ,Θ ) = i [ ( y i μ x i ) 2 2 σ x i 2 +ln σ x i ] (9)

Under these assumptions, the HMRF-EM algorithm is as follows:

1) First, use the kmeans method to initialize the parameter set Θ ( 0 ) .

2) At the t -th iteration, calculate the likelihood distribution P ( t ) ( y i | x i , θ x i ) .

3) The label is estimated using the current parameter set Θ ( 0 ) 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:

X ( t ) = argmax Xχ { P( Y|X, Θ ( t ) )P( X ) } = argmin Xχ { ln[ P( Y|X, Θ ( t ) )P( X ) ] } = argmin Xχ { lnP( Y|X, Θ ( t ) )lnP( X ) } = argmin Xχ { U( Y|X, Θ ( t ) )+U( X ) } (10)

4) Calculate the posterior distribution for all lL and all pixels y i using the Bayesian rule:

P ( t ) ( l| y i )= G( y i ; θ l )P( l| x N i ( t ) ) P ( t ) ( y i ) (11)

where x N i ( t ) is the neighborhood configuration of x i ( t ) , and

P ( t ) ( y i )= lL G( y i ; θ l )P( l| x N i ( t ) ) (12)

Note that the posterior probability of label l given the neighborhood features x N i is defined as:

P( l| x N i ( t ) )= 1 Z exp( j N i V c ( l, x j ( t ) ) ) (13)

5) Use P ( t ) ( l| y i ) to update the parameters:

μ l ( t+1 ) = i P ( t ) ( l| y i ) y i i P ( t ) ( l| y i ) (14)

( σ l ( t+1 ) ) 2 = i P ( t ) ( l| y i ) ( y i μ l ( t+1 ) ) 2 i P ( t ) ( l| y i ) (15)

In the previous assumption, the intensity distribution of each region to be segmented follows a Gaussian distribution with parameters θ x i =( μ x i , σ x i ) . 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 g components can be represented by the parameters:

θ l ={ ( μ l,1 , σ l,1 , w l,1 ),,( μ l,g , σ l,g , w l,g ) } (16)

GMM now has weighted probabilities:

G min ( z; θ l )= c=1 g w l,c G( z; μ l,c , σ l,c ). (17)

In our experiment, we suppose g=1 , 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 I:Ω R 3 on a given image domain Ω , for a voxel x( x,y,z )Ω in the image I , we define the local energy as:

e i ( x )= K( yx ) | I( x )b( y ) c i | 2 dy (18)

K is a kernel function with a standard deviation of σ , truncated by a Gaussian function defined by:

K( u )={ 1 a e | u | 2 / 2 σ 2 , | u |ρ, 0, | u |>ρ. (19)

where ρ is the radius, the spherical neighborhood centered on y . a is a normalization constant that makes K( u )du =1 .

The local energy e i ( x ) can be equivalently expressed as the following calculation function in actual numerical calculations:

e i ( x )= I 2 1 K 2 c i I( bK )+ c i 2 ( b 2 K ) (20)

where is the convolutional operation, 1 K ( x )= K( yx )dy .

We define the following global energy as a data fitting term:

E= i=1 2 Ω i ( K( yx ) | I( x )b( y ) c i | 2 dx )dy ,x,yΩ (21)

Let ϕ be the level set function. We consider the image domain Ω being divided into two disjoint regions Ω 1 and Ω 2 . In this case, the regions Ω 1 and Ω 2 are represented by their membership functions M 1 ( ϕ )=H( ϕ ) and M 2 ( ϕ )=1H( ϕ ) , respectively, where H is the Heaviside function, and δ ϵ ( x ) is the derivative of H ϵ ( x ) . Defined as follows:

H ϵ ( x )= 1 2 [ 1+ 2 π arctan( x ϵ ) ] (22)

δ ϵ ( x )= H ϵ ( x )= 1 π ϵ ϵ 2 + x 2 (23)

Then the energy in (21) can be expressed in the following level set formulation:

E= i=1 2 Ω i ( K( yx ) | I( x )b( y ) c i | 2 dy ) M i ( ϕ( x ) )dx ,x,yΩ (24)

Two regularizations ( ϕ ) and p ( ϕ ) are introduced:

( ϕ )= | H( ϕ ) |dx (25)

p ( ϕ )= p( | ϕ | )dx (26)

where p( s )= 1 2 ( s1 ) 2 .

In summary, the variational level set framework is:

E= i=1 2 Ω i ( K( yx ) | I( x )b( y ) c i | 2 dy ) M i ( ϕ( x ) )dx +ν Ω δ ϵ ( ϕ( x ) )| ϕ( x ) |dx +μ Ω p( | ϕ | )dx ,x,yΩ (27)

where ν and μ are constants greater than zero. The above equation is minimized by the gradient descent method

ϕ t = δ ϵ ( ϕ )( e 1 e 2 )+ν δ ϵ ( ϕ )diν( ϕ | ϕ | )+μdiν( d p ( | ϕ | )ϕ ) (28)

Function d p ( s ) p ( s ) s .

Fixing ϕ and b , so that the energy E is the least optimal c is defined

c ^ i = ( bK )I u i dydx ( b 2 K ) u i dydx ,i=1,2 (29)

where u i ( x )= M i ( ϕ( x ) ) .

Fixing ϕ and c , so that the energy E is the smallest optimal b defined as:

b ^ = ( I J ( 1 ) )K J ( 2 ) K (30)

where J ( 1 ) = i=1 N c i u i and J ( 2 ) = i=1 N c i 2 u i .

2.3. Implementation

As we all know, the level set is sensitive to initialization. Now we directly map the coarse segmentation label matrix X * to the initial implicit function ϕ 0 ( x,y,z ) 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 :

ϕ 0 ( x,y,z )={ c 0 , X * ( x,y,z )=1, c 0 , X * ( x,y,z )=2. (31)

where c 0 =1 is a constant, X * ( x,y,z )=1 represents the background region, X * ( x,y,z )=2 represents the target region, and ( x,y,z ) is an arbitrary voxel of the image domain Ω .

The main steps of the proposed model are as follows :

1) Find the label X * .

2) Based on the label result, initialize the level set function ϕ 0 ( x,y,z ) using the expression (31).

3) Evolve the level set using Equation (28), and then update c and b 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 μ=1 , the value range of the length parameter is ν=[ 0.005× 255 2 ,0.02× 255 2 ] , and the parameter of the Gaussian kernel function is σ=[ 2,4 ] . 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.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1] Peng, J., Li, W., Liao, S. and Zhu, Y. (2024) Microcrack Segmentation of 3D CT Images Based on Multi-Task Learning. IEEE Access, 12, 138192-138200.[CrossRef]
[2] Yang, F., Zhang, D., Huang, K., Shi, W. and Wang, X. (2018) Scattering Estimation for Cone-Beam CT Using Local Measurement Based on Compressed Sensing. IEEE Transactions on Nuclear Science, 65, 941-949.[CrossRef]
[3] Tang, S., Huang, K., Cheng, Y., Mou, X. and Tang, X. (2018) Optimization Based Beam-Hardening Correction in CT under Data Integral Invariant Constraint. Physics in Medicine & Biology, 63, Article ID: 135015.[CrossRef] [PubMed]
[4] Liu, C., Liu, W. and Xing, W. (2019) A Weighted Edge-Based Level Set Method Based on Multi-Local Statistical Information for Noisy Image Segmentation. Journal of Visual Communication and Image Representation, 59, 89-107.[CrossRef]
[5] Zheng, J., Zhang, D., Huang, K. and Sun, Y. (2018) A CBCT Series Slice Image Segmentation Method. Journal of X-Ray Science and Technology: Clinical Applications of Diagnosis and Therapeutics, 26, 815-832.[CrossRef] [PubMed]
[6] Qian Zhao, Y., Hong Wang, X., Fang Wang, X. and Shih, F.Y. (2014) Retinal Vessels Segmentation Based on Level Set and Region Growing. Pattern Recognition, 47, 2437-2446.[CrossRef]
[7] Zheng, J., Tang, C. and Sun, Y. (2024) Thresholding-Accelerated Convolutional Neural Network for Aeroengine Turbine Blade Segmentation. Expert Systems with Applications, 238, Article ID: 122387.[CrossRef]
[8] Zhang, W., Wang, X., You, W., Chen, J., Dai, P. and Zhang, P. (2020) RESLS: Region and Edge Synergetic Level Set Framework for Image Segmentation. IEEE Transactions on Image Processing, 29, 57-71.[CrossRef] [PubMed]
[9] Liu, G., Jiang, Y., Chang, B. and Liu, D. (2022) Superpixel-Based Active Contour Model via a Local Similarity Factor and Saliency. Measurement, 188, Article ID: 110442.[CrossRef]
[10] Braga, A.M., Marques, R.C.P., Medeiros, F.N.S., Neto, J.F.S.R., Bianchi, A.G.C., Carneiro, C.M., et al. (2021) Hierarchical Median Narrow Band for Level Set Segmentation of Cervical Cell Nuclei. Measurement, 176, Article ID: 109232.[CrossRef]
[11] Daimary, D., Bora, M.B., Amitab, K. and Kandar, D. (2020) Brain Tumor Segmentation from MRI Images Using Hybrid Convolutional Neural Networks. Procedia Computer Science, 167, 2419-2428.[CrossRef]
[12] Ranjbarzadeh, R., Bagherian Kasgari, A., Jafarzadeh Ghoushchi, S., Anari, S., Naseri, M. and Bendechache, M. (2021) Brain Tumor Segmentation Based on Deep Learning and an Attention Mechanism Using MRI Multi-Modalities Brain Images. Scientific Reports, 11, Article No. 10930.[CrossRef] [PubMed]
[13] Huang, K., Li, Z., Tang, S., Zeng, Y., Ye, W. and Yang, F. (2024) A Level-Set Method with the 3D Multiplicative-Additive Model for CT Volume Data Segmentation. Measurement, 229, Article ID: 114442.[CrossRef]
[14] Osher, S. and Sethian, J.A. (1988) Fronts Propagating with Curvature-Dependent Speed: Algorithms Based on Hamilton-Jacobi Formulations. Journal of Computational Physics, 79, 12-49.[CrossRef]
[15] Paragios, N. and Deriche, R. (2000) Geodesic Active Contours and Level Sets for the Detection and Tracking of Moving Objects. IEEE Transactions on Pattern Analysis and Machine Intelligence, 22, 266-280.[CrossRef]
[16] Fubin He, and Yi Sun, (2015) Segmentation of Noisy CT Volume Data Using Improved 3D Chan-Vese Model. 2015 IEEE 7th International Conference on Awareness Science and Technology (iCAST), Qinhuangdao, 22-24 September 2015, 31-36.[CrossRef]
[17] Lei Zhang, and Qiang Ji, (2010) Image Segmentation with a Unified Graphical Model. IEEE Transactions on Pattern Analysis and Machine Intelligence, 32, 1406-1425.[CrossRef] [PubMed]
[18] Wang, Q. and Boyer, K.L. (2012) The Active Geometric Shape Model: A New Robust Deformable Shape Model and Its Applications. Computer Vision and Image Understanding, 116, 1178-1194.[CrossRef]
[19] Wang, Q. (2012) GMM-Based Hidden Markov Random Field for Color Image and 3D Volume Segmentation. arXiv: 1212.4527.
[20] Wang, Q. (2012) HMRF-EM-Image: Implementation of the Hidden Markov Random Field Model and Its Expectation-Maximization Algorithm. Computer Science, 94, 222-233.
[21] Li, C.M., Huang, R., Ding, Z.H., Gatenby, J.C., Metaxas, D.N. and Gore, J.C. (2011) A Level Set Method for Image Segmentation in the Presence of Intensity Inhomogeneities with Application to MRI. IEEE Transactions on Image Processing, 20, 2007-2016.[CrossRef] [PubMed]

Copyright © 2026 by authors and Scientific Research Publishing Inc.

Creative Commons License

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.