A fast and efficient numerical algorithm for image segmentation and denoising

1.   Introduction

Image segmentation is the process of dividing an image into homogenous regions based on a characterization of the given image, and represents one of the most fundamental procedures of image processing ^{[1,2,3,4,5,6,7,8,9,10,11,12]}. Currently, the prevalent applications of pattern recognition technologies such as medical image recognition, robotics, and visual field monitoring are achieved through the implementation of image segmentation algorithms ^[1]. Furthermore, Ahmad et al. ^[2] presented an augmented Lagrangian technique for the image deblurring problem in its primal form based on mean curvature. Song et al. ^[3] proposed tensor Conjugate-Gradient schemes, and applied the proposed scheme to restore images and videos. Lee et al. ^[4] and Li et al. ^[5] studied the restoration of fingerprint images, and they utilized the phase-field model for image processing. Pan and Feng ^[6] conducted research on parameter estimation necessary for image restoration. They proposed algorithms to improve the accuracy of parameter estimation, and verified the advantages when using the estimated parameters for image restoration.

Numerous phase-field methods have been developed for image segmentation ^{[7,8,9,10,11,12]}. Chen et al. ^[7] proposed a novel image segmentation model designed for complex two-dimensional images with intensity inhomogeneity. The proposed phase-field model uses local information. Furthermore, their proposed scheme excels at minimizing background interference compared to previous methods. Fang et al. ^[8] introduced an integrated vector-valued active contour model for image segmentation. They combined the advantages of edge-based and region-based active contour models to propose an integrated method that uses adaptively adjustable weights to effectively obtain boundary information. The proposed scheme efficiently segmented vector-valued images. Jeong et al. ^[9] solved a modified Allen-Cahn (AC) equation with a fidelity term. They proposed a fast and accurate hybrid method that possesses the advantages of simplicity and versatility in arbitrary computational domains. Liu et al. ^[10] presented an Allen-Cahn Cahn-Vese (ACCV) model that deals with multi-phase image segmentation. They validated the effectiveness of the proposed method by comparing it to previous models. Thasneem et al. ^[11] proposed an algorithm that consists of three steps: Slice selection, segmentation, and interpolation. In the segmentation step, they presented a segmentation method using a modified AC equation. Yang et al. ^[12] proposed a modified Cahn-Hilliard (CH) equation for image segmentation. A significant feature of the proposed model is that it can interpolate omitted contours over wide distances. They improved the tailored-finite-point method for the proposed image segmentation model. Image denoising is a process that aims to remove unwanted noise from an image and enhances its visual quality. Noise, often caused by factors like low-light conditions or sensor limitations, can distort image details. Denoising algorithms can remove this noise while preserving essential image features. Common techniques include wavelet transforms, non-local means, and deep learning-based approaches. In deep learning, convolutional neural networks (CNNs) are frequently used for their ability to learn complex noise patterns and restore clean images ^[13]. Image denoising plays a crucial role in various applications, such as medical imaging, surveillance, and photography to have clearer and more accurate visual representations ^[14,15].

The main purpose of this study is to present a novel, highly efficient, and fast numerical scheme specifically developed for a phase-field mathematical model of image segmentation. The proposed computational scheme is based on an operator splitting method ^[16]. The three split operators are solved by using closed-form analytic solutions and a finite difference method (FDM) with an alternating direction explicit (ADE) method ^[17].

The outline of this paper is as follows: In Section 2, we present the proposed numerical solution algorithm. In Section 3, we conduct a series of numerical experiments to validate the notable efficiency and rapid computational performance of the proposed numerical algorithm. Finally, conclusions are given in Section 4.

2.   Proposed numerical solution algorithm

We consider the following phase-field model for image segmentation and denoising:

where $u({\bf x}, t)$ is an order parameter at space ${\bf x}$ and time $t$, $\epsilon$ is an interfacial thickness related parameter, $\lambda$ is a fidelity parameter, and $f({\bf x})$ is the given image. In addition, $c_1(t)$ and $c_2(t)$ are the averages of the given image $f$ in the region ($u\ge0$) and ($u < 0$), respectively, and are defined as

The phase-field model for image segmentation and denoising can be derived from the phase-field approximation of the Mumford-Shah functional

Given our focus on a stable numerical algorithm in this study, for a more detailed derivation and definitions, please refer to ^[18].

Equation (2.1) is discretized in a domain $\Omega = (L_x, R_x) \times (L_y, R_y)$, which is discretized as follows: $\Omega_{h} = \{ (x_i, y_j) : x_i = L_x+ih, \; y_j = L_y + jh, \; 0 \le i \le N_{x}, \; 0 \le j \le N_y \}$, where $h = (R_x-L_x)/N_x = (R_y-L_y)/N_y$ for some integers $N_x$ and $N_y$. For simplicity of notion, let $u_{ij}^n = u(x_i, y_j, n \Delta t)$, where $\Delta t$ is the uniform time step and $n$ is the number of temporal iteration steps.

Using the operator splitting method, we first solve the following term in the governing Eq (2.1):

which is analytically solved as follows and the solution after $\Delta t$ time is given as for $i = 1, 2, \ldots, N_x-1$ and $j = 1, 2, \ldots, N_y-1$:

Here, $c_{1}^n$ and $c_{2}^n$ are defined as

Next, we solve the diffusion term in the governing Eq (2.1)

which is numerically solved by a semi-implicit Euler method

where $\Delta_h u_{ij}^{n+1, 2} = (u_{i-1, j}^{n+1, 2} +u_{i+1, j}^{n+1, 2} -4u_{ij}^{n+1, 2}+u_{i, j-1}^{n+1, 2}+u_{i, j+1}^{n+1, 2})/h^2$. To efficiently solve the discrete Eq (2.6), we apply the ADE method ^[17] as follows: If $n = 4m+1$ for some integer $m$, then for $j = 1, 2, \ldots, N_y-1$ and $i = 1, 2, \ldots, N_x-1,$

Equation (2.7) can be rewritten as

If $n = 4m+2$ for some integer $m$, then for $j = 1, 2, \ldots, N_y-1$ and $i = N_x-1, N_x-2, \ldots, 1$,

which can be rewritten as

If $n = 4m+3$ for some integer $m$, then for $j = N_y-1, N_y-2, \ldots, 1$ and $i = 1, 2, \ldots, N_x-1$,

which can be rewritten as

If $n = 4m$ for some integer $m$, then for $j = N_y-1, N_y-2, \ldots, 1$ and $i = N_x-1, N_x-2, \ldots, 1$,

which can be rewritten as

For simplicity of exposition, we use the Dirichlet boundary condition

We should note that there are many other numerical methods for solving the diffusion equation such as Fourier spectral method and multigrid methods ^[19]. In this study, we choose the above scheme because it is simple and practically stable.

Finally, we solve the following nonlinear term in the governing Eq (2.1):

which is solved by using a separation of variables and its closed-form solution is given as

We redefine $f$ as $f = (f-f_{min})/(f_{max} -f_{min})$, where $f_{max}$ is the maximum and $f_{min}$ is the minimum of $f$. In this study, we use $u^0 = 2f-1$, which implies that $-1 \le u^0 \le 1$. Because we use $(u^2-1)^2 / (4\epsilon^2)$ in the total cost functional ${\mathcal{E}}(u)$, the values of $u$ in the bulk phases are approximately $\pm 1$.

In summary, the proposed numerical algorithm follows these steps:

We note that the main motivation of the proposed algorithm is as follows: It is an explicit scheme having practical stability because of using an ADE method for solving the diffusion equation, which was previously solved by a fully explicit scheme.

3.   Computational tests

Now, we conduct several computational tests using the proposed numerical method on synthetic and real images.

3.1.   Comparison test with a previous method

In the previous study, Jeong et al. ^[9] proposed the explicit hybrid method for image segmentation. We compare the previous and our proposed methods for the segmentation of a real blood vessel image. The image consists of $64$ pixels. Figure 1 displays the image segmentation results for a real blood vessel image with the previous method (top row) and our proposed method (bottom row). The model parameters used are $\epsilon = \epsilon_{4}$ and $\lambda = 10$. We use the time step size $\Delta t = 0.1$ and $\Delta t = 0.5$ for the previous and proposed methods, respectively. Here, $\epsilon_m = hm/[2\sqrt{2} \tanh^{-1}(0.9)]$; see ^[9] for more detailed information. The previous explicit hybrid method is fast and accurate, however it has the time step restriction because it is a fully implicit scheme. Our proposed method achieved successful segmentation through fewer iterations with more stability for the time step size.

Figure 1.  Temporal evolutions of segmentation results of the left anterior descending vessel with the explicit hybrid method ^[9] (top row) and the proposed method (bottom row), the number of iterations and time are shown below each figure.

DownLoad: Full-Size Img PowerPoint

3.2.   Image segmentation

In this section, we investigate the segmentation performance of the proposed numerical algorithm. To validate segmentation results, we consider images of $256 \times 256$ pixels. First, we consider a brain MR image as shown in Figure 2(a). Our proposed algorithm can be used for the analysis of medical images to provide essential information for medical treatment. Figure 2 shows the temporal evolution of segmentation results for a solid brain tumor. Here, we use $\Delta t = 0.5$, $h = 1$, $\epsilon = \epsilon_4$, and $\lambda = 4$. From the result in Figure 2(c), we can observe that the proposed method segments solid brain tumors in the brain MR image.

Figure 2.  Temporal evolution of segmentation results of a solid brain tumor image, the number of iterations and time are shown below each figure.

DownLoad: Full-Size Img PowerPoint

Next, we consider a satellite image of Europe's night-lights, which is complex scattered data, to perform segmentation using the proposed method. The given image of Europe's night-lights is shown in Figure 3(a). Figure 3(b), (c) shows the segmentation results at $t = 10$ and $t = 50$, respectively. We note that the homogeneous Neumann boundary condition is applied in this test. Here, interfacial parameter $\epsilon = \epsilon_{80}$, $\lambda = 2$, time step $\Delta t = 0.5$, and $h = 1$ are used. The proposed algorithm generates a clear segmentation result for the European terrain.

Figure 3.  (a) Given image of Europe night-lights; temporal evolution of contour at (b) $t = 10$ and (c) $t = 50$.

DownLoad: Full-Size Img PowerPoint

3.3.   Image denoising

Next, we consider image denoising by using the proposed numerical algorithm for the phase-field model. The governing equation contains a diffusion term, where diffusion is effective in noise removal from images due to its inherent smoothing property. In image processing, diffusion processes, such as the diffusion equation, redistribute pixel values by averaging neighboring intensities over time. This averaging removes high-frequency variations associated with noise and promotes a more uniform and visually coherent appearance. Diffusion operates as a low-pass filter and reduces high-frequency noise components while preserving essential image features. The gradual blending of pixel values minimizes abrupt intensity changes, resulting in a smoother image appearance and effectively reducing the impact of noise for enhanced visual clarity.

Figure 4 shows an image denoising process on a text image corrupted with $10\%$ salt-and-pepper noise level using the proposed algorithm. Figures 4(a), (b), and (c) are the initial image, computational results after $1$ iteration, and $4$ iterations, respectively. Here, the denoising is performed in a domain $\Omega = (0, 4)\times(0, 1)$ with grid sizes $N_x = 256$ and $N_y = 128$, $\epsilon_{2}$, $\Delta t = 2 \times 10^{-5}$, and $\lambda = 10^4$ are used. From this test result, we can confirm that the proposed algorithm is highly efficient because it can obtain a good result with a small number of iterations.

Figure 4.  (a) Initial noisy image, (b) and (c) are the computational results after $1$ iteration and $4$ iterations, respectively.

DownLoad: Full-Size Img PowerPoint

3.4.   Fingerprint image

We consider two fingerprint images with $10\%$ salt-and-pepper noise. The proposed scheme is applied to the noisy fingerprint images, and the results are shown in Figure 5. Figure 5 displays image processing for corrupted fingerprint images on a domain $\Omega = (0, N_x) \times (0, N_y)$. In the first row from Figure 5, interface parameter $\epsilon_{2}$, time step $\Delta t = 0.2$, $N_x = 256$, $N_y = 256$, $h = 1$, and $\lambda = 1$ are used. In the second row from Figure 5, interface parameter $\epsilon_{2}$, time step $\Delta t = 0.3$, $N_x = 300$, $N_y = 300$, $h = 1$, and $\lambda = 20$ are used. We note that image noise was removed with only a small number of iterations because we can use a practically larger time step.

Figure 5.  (a) Initial images with noise; (b) $t = \Delta t$ ($1$ iteration); (c) $t = 3\Delta t$ ($3$ iterations).

DownLoad: Full-Size Img PowerPoint

3.5.   Bar code image

We conduct denoising test with a seriously damaged image. Figure 6 shows the results of image denoising. Figure 6(a) is the original image. Figure 6(b) is the corruption of Figure 6(a) with scratch and $10\%$ salt-and-pepper noise level. Using the proposed algorithm, we obtain a denoised result as shown in Figure 6(c). Here, the used parameters are the interface parameter $\epsilon_{2}$, time step $\Delta t = 0.2$, $N_x = 380$, $N_y = 250$, $h = 1$, and $\lambda = 0.1$. As a result, we obtain a denoised image that can be scanned with a bar code scanner.

Figure 6.  (a) Original bar code image; (b) damaged bar code image; (c) $t = 9\Delta t$ ($9$ iterations).

DownLoad: Full-Size Img PowerPoint

4.   Conclusions

In this article, we proposed a novel, highly efficient, and fast computational algorithm, which is based on an operator splitting method, specifically developed for a modified AC equation of image segmentation. The split operators were solved by using two closed-form analytic solutions and the FDM with an ADE method. To validate the superior efficiency and rapid computational performance of the proposed computational method, we conducted several numerical tests, such as noise removal in a blood bessel image, a text image, and fingerprint images. The computational tests confirmed a high efficiency of the proposed algorithm.

Use of AI tools declaration

The authors have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The first author (Y. Z. Jin) was supported by the Natural Science Foundation of Jilin Province [YDZJ202101ZYTS186]. The corresponding author (J. S. Kim) was supported by the Brain Korea 21 FOUR through the National Research Foundation of Korea funded by the Ministry of Education of Korea.

The authors are grateful to the referees whose comments greatly improved the paper.

Conflict of interest

Prof. Junseok Kim is the Guest Editor of special issue "Accurate and Efficient Numerical Methods for the Multi-Component Partial Differential Equations including the Phase-Field Models" for AIMS Mathematics. Prof. Junseok Kim was not involved in the editorial review and the decision to publish this article.

The authors declare there is no conflicts of interest.

Appendix

The following MATLAB source code and data are also available from the corresponding author's webpage: https://mathematicians.korea.ac.kr/cfdkim/open-source-codes/.