BOOSTING MR IMAGE IMPULSE NOISE REMOVAL WITH OVERLAPPING GROUP SPARSE FRACTIONAL-ORDER TOTAL VARIATION AND MINIMAX CONCAVE PENALTY

The acquisition and transmission of magnetic resonance images are susceptible to noise, particularly impulse noise. Although the method based on the ℓ 0 -norm and overlapping group sparse total variation ( ℓ 0 -OGSTV) is effective for impulse noise image restoration, it can only mitigate the staircase artifacts to a certain extent. To boost the impulse noise removal performance of ℓ 0 -OGSTV, we propose a new restoration model that consists of two terms. Specifically, in the first term, we keep using the ℓ 0 -norm as the data fidelity term to eliminate impulse noise. In the second term, we first introduce an overlapping group sparsity fractional-order total variation regularizer to eliminate staircase artifacts while preserving structural information. Then, we adopt the minimax-concave penalty to further accurately estimate the image edges. Finally, we employ an alternate direction method of multipliers to solve the proposed optimization model. Clinical experiments demonstrate its effectiveness in denoising medical images.


INTRODUCTION
Due to the reasons of acquisition technology and system, noise and artifacts are introduced in magnetic resonance (MR) images, and impulse noise is one of the main noises.The impulse noise may be misunderstood as an anomaly of the human health system, and denoising can enrich the visual quality of images.Denoising plays a crucial role in medical imaging and scanning, where even a small amount of noise can be misinterpreted as an anomaly in the human health system.However, removing noise is a complex process because it can sometimes compromise the visual effect and details of the image.Therefore, the restoration of MR images destroyed by impulse noise has become a prominent area of research in the field of image processing.
Many methods have been proposed for medical image denoising, one of the most popular image restoration methods is total variation (TV) model (Rudin et al., 1992), which could preserve edges and remove image noise in homogeneous regions.However, it tended to produce staircase artifacts in smooth regions.To solve this problem, various solutions have been proposed.For example, high-order TV methods (Adam et al., 2021;Ge et al., 2023), the overlapping group sparsity total variation (OGSTV) methods (Liu et al., 2014;Shi et al., 2016), fractional TV methods (Lian and Liu, 2023;Rahman Chowdhury et al., 2020;Zhu et al., 2022) and the minimaxconcave (MC) penalty methods (Du and Liu, 2018;Chen and Zhao, 2023).Among these methods, the OGSTV method excels in image restoration due to its structural sparsity.In (Ge et al., 2023), a method that combined the hyper-Laplacian prior regularization with OGSTV and nonconvex second-order TV, performed well in removing noise from MR images.Ji and Zhao (2023) proposed a novel model with a non-convex penalty that combined the OGS regularizer and the MC penalty (OGS-MCTV).The non-convex MC penalty could preserve edges.Therefore, this method demonstrated a superior image denoising effect compared to other models.Bhutto et al. (2023) proposed an image denoising algorithm that cleverly integrated the advantages of the fractional-order variation domain with an OGS measure, which acted as its regularization component.The fidelity term of these models used the ℓ 2 -norm, which is commonly used to restore images degraded by additive Gaussian noise.In addition, the ℓ 2norm is sensitive to outliers and can easily result in unsatisfactory image restoration.
To effectively remove impulse noise, as described in (Gao et al., 2018), Bayesian statistical rules suggested that the ℓ 1 -norm fidelity was more suitable for restoring images corrupted by impulse noise than the ℓ 2 -norm.In (Yang et al., 2009), the fidelity term of the ℓ 1 -norm was utilized to restore the fuzzy multichannel image damaged by impulse noise.Chan et al. (2010) proposed a two-phase image restoration method based on TV regularization combined with the ℓ 1 -norm data fidelity term for impulse noise removal.Numerical results proved that the proposed method makes good progress in restoration capability.Although the ℓ 1 -norm has demonstrated significant advantages in sparse signal processing and image restoration, it may overly penalize the obtained solution in impulse noise removal (Kuang et al., 2018).
To address the aforementioned issues, a method for removing impulse noise using ℓ 0 total variation (ℓ 0 -TV) was proposed in (Yuan and Ghanem, 2017).It can be depicted as where u ∈ R n×m is the desired original clean image, f ∈ R n×m is the degraded image, λ > 0 is the regularization parameter, o ∈ {0, 1} n is specified by the user, ⊙ denotes an elementwise product, K ∈ R n×n is a linear operator.In this paper, we are concerned with K = I, the identity operator, which constitutes a denoising problem.
Recently, the ℓ 0 -norm data fidelity term has been used to remove impulse noise in (Yin et al., 2022).The OGSTV serves as a regularizer to effectively eliminate staircase artifacts, making this model highly proficient in image restoration tasks even under high impulse noise levels.Sun and Liu (2023) integrated both the ℓ 0norm data fidelity term and the nonconvex generalized regularizer, demonstrating not its remarkable ability to suppress impulse noise and its superior capability in preserving sharp contours while reducing staircase artifacts.These models collectively demonstrate the suitability of ℓ 0 -norm for restoring images corrupted by impulse noise.
According to the literature survey above, the existing methods primarily focus on restoring natural images.However, when it comes to medical image restoration, more emphasis is placed on preserving intricate textures to ensure clinical diagnostic value.Therefore, this paper primarily investigates boosting MR image impulse noise removal.The model comprises an ℓ 0 -norm data fidelity term, a regularizer of overlapping group sparse fractionalorder total variation (OGS-FOTV) and the MC penalty.Our proposed method combines the advantages of the fractional-order variation domain with the OGS measure, which can effectively measure complex texture details and reduce staircase artifacts.Additionally, the MC penalty can improve the sparsity of images in the gradient domain and improve the estimation of high-frequency components, that is, preserve edges.To address the computational challenges stemming from the model's complexity, we employ the alternate direction multiplier algorithm to solve the subproblems.Finally, we conduct numerical experiments to analyze the effectiveness of our proposed model.The rest of this article is organized as follows.Section 2 presents some elementary concepts and preliminaries related to the proposed algorithm.In Section 3, we propose a new model for impulse noise removal and derive an efficient algorithm to solve the corresponding minimization problem.In Section 4, the superiority of the proposed method is proved by numerical experiments.Finally, a conclusion is made in Section 5.

PRELIMINARIES
In order to better describe the model, in this section, we introduce the definition of the ℓ 0 fidelity term, the Moreau envelope, the MC penalty, overlapping group sparsity, discrete fractional-order difference and the ADMM framework.
THE ℓ 0 FIDELITY TERM First, we give some basic definitions and properties related to the ℓ 0 -norm fidelity term.The o ∈ {0, 1} n is specified by the user.More specifically, when o i = 0, the pixel at position i is an outlier, and when o i = 1, the pixel at position i is a potential outlier.For this paper, we set o i = 0, f i = u min or u max 1, otherwise for the salt-andpepper impulse noise.
The following lemma, as delineated in (Yuan and Ghanem, 2017), presents the variational characterization of the ℓ 0 -norm.
Lemma 1 For any given w ∈ R n , it holds that and z * = 1 − sign(|w|) is the unique optimal solution to problem (2).Here, the standard signum function sign is employed in component form, and sign(0) = 0.

THE MINIMAX-CONCAVE PENALTY
The Moreau envelope is convex, continuous, differentiable and real valued (Zhou and Zhao, 2021).
Definition 2 Let b ≥ 0, the minimax-concave penalty of ∥x∥ 2 with parameter b is the function φb (x) : R N → R which is given by (Du and Liu, 2018) where the function φb (x) is non-convex (Ji and Zhao, 2023;Shen et al., 2021).
OVERLAPPING GROUP SPARSITY Liu et al. (2015;2014) expanded the OGSTV function as a new regularizer from the onedimensional signal denoising problem (Selesnick and Chen, 2013) to the general two-dimensional case.This approach effectively reduces staircase artifacts.They also defined a where and ⌊x⌋ denotes the largest integer less than or equal to x.The center of gi, j,K is (i, j).Let g i, j,K be a vector which is obtained by stacking the K columns of the matrix gi, j,K , i.e., g i, j,K = gi, j,K (:).Then the OGS regularizer can be defined by

DEFINITION OF DISCRETE FRACTIONAL-ORDER DIFFERENCE
There exist multiple interpretations of fractional differences, and the Grünwald-Letnikov (G-L) fractional-order derivative (Zhang et al., 2012) is among the prevalent methods.Defined by G-L, let the size of an image u be N × M. Thus, the discrete form of the fractional-order gradient ∇ α u can be evaluated by where α is the fractional order and we set 1 ≤ α < 2 in this paper.The discrete gradients D α x u, D α y u ∈ R N×M along the x-axis and the y-axis are given by Here K is the number of adjacent pixels that are used to calculate the fractionalorder derivative at each pixel.The coefficients with the Gamma function Γ(x).Furthermore, the conjugate operator of the fractional-order gradient operator is (∇ α ) * = (−1) α div α .In the discrete case, the vector p = (p (1) , p (2) ) ∈ R N×M × R N×M discrete fractionalorder divergence is defined as (Rahman Chowdhury et al., 2020;Bhutto et al., 2023) THE ALTERNATING DIRECTION

METHOD OF MULTIPLIERS
The alternating direction method of multipliers (ADMM) is to solve the following constrained separable optimization problems: where ξ i (•) : χ i → R are closed convex functions, A, B ∈ R l×n and d ∈ R l is a given vector.The augmented Lagrangian function (Hestenes, 1969) for the problem (10) is where µ ∈ R l is the Lagrange multiplier and δ > 0 is a penalty parameter.The objective is to find the saddle point of L A by alternatively minimizing L A with respect to x, y and µ.The ADMM algorithm to solve problem ( 10) is presented as Algorithm 1.

THE PROPOSED ALGORITHM
In this section, we first introduce the proposed MR Image denoising model, and then solve it in the ADMM framework.

MODEL
The proposed MR Image denoising model is as follows min where Φ K b (u) is a new non-convex penalty based on the MC penalty (Ji and Zhao, 2023) and the overlapping group sparse fractional-order total variation (OGS-FOTV) (Bhutto et al., 2023), which is defined as follows where φ (•) is the OGS regularizer as defined by Eq. ( 6).And φ (∇ α u) represents the OGS-FOTV regularizer, ∇ α u is the fractional-order gradient as defined by Eq. ( 7).
Next, according to Eq. ( 13), we proposed model ( 12) can be reformulated as

u-subproblem
The u-subproblem Subproblem ( 17) is a least squares problem, we can solve the equivalent normal equation, For the periodic boundary condition of u, (∇ α ) T ∇ α is the block circulant with circulant blocks, which can be diagonalized by 2D discrete Fourier transform (Wu and Tai, 2010).Therefore, through a fast Fourier transform (FFT) operation and a FFT inverse operation, we can obtain the optimal u as follows x-subproblem The x-subproblem is the overlapping group sparse problem, we have The problem ( 12) can be solved iteratively by the majorization-minimization (MM) algorithm, and the process of solving the related problem is discussed in detail in (Yin et al., 2022).Here, we express it in lemma 2, as follows Lemma 2 we consider a minimization problem of the form , where ρ is a positive parameter and the functional Φ(v) = ∑ n i, j=1 g (i, j),K 2 .In order to minimize P(v), the MM algorithm is continuously iteratively solved to obtain

v-subproblem
The v subproblem can be written as The algorithm for v-subproblem is the same as the algorithm for x-subproblem.

z-subproblem
The z-subproblem is of the following form The solution z k+1 can be computed using projection as: y-subproblem The y-subproblem can be solved by a soft thresholding with the shrink operatorction.Besides, this subproblem shares the following form which can be simplified as Subproblem ( 26) can be calculated by where shrink(x, y) = sgn(x) • max {∥x∥ 1 − y, 0} , and sgn(•) denotes the signum function.

Updating Lagrangian multipliers
Finally, the Lagrange multipliers are updated by the following The proposed method is presented as Algorithm 2. And for this algorithm, we have two remarks.
Remark 1 When α = 1, the fractional-order TV is reduced to the standard TV, so the Eq. ( 13) is degraded into the OGS-MCTV penalty in (Ji and Zhao, 2023).In the experimental section, we discuss in detail the effect of the value of α on the denoising performance of the proposed model.
Remark 2 When K = 1, the Eq. ( 13) is degraded into the MC penalty of the fractional-order TV term in (Chen and Zhao, 2023).Different K values lead to different denoising results, specific experiments are shown in the following section.

NUMERICAL EXPERIMENTS
In this section, we present several experimental results to verify the effectiveness of the proposed method for image denoising.The test image is real MR images, as shown in Fig. 1.All the test images are from the Department of Radiology, Maanshan People's Hospital, Maanshan, China.The experiment is under Windows 10 and MATLAB R2021a operating system, and the CPU is AMD R7 5800H 3.20GHz and 16GB RAM.MR images used for analysis were corrupted by salt-and-peppe noise levels of 30%, 50%, 70%.All denoising methods measure the quality of recovered images by Peak Signal-to-Noise Ratio (PSNR) and Structural Similarity (SSIM) indices, which reflect human subjective sensory and visual perception quality respectively (Irum et al., 2015;Kuang et al., 2018).
The stopping criterion for all tested algorithms is set to

RESULTS AND ANALYSIS
The experimental results of our proposed model are compared with three related methods: ℓ 0 -OGSTV (Yin et al., 2022), ℓ 1 -OGSTV (Liu et al., 2015) and HNHOTV-OGS (Adam et al., 2021).Their regularizers are closely related to the proposed method in this paper.
In the whole experiment, we fixed α = 1.9,K = 3, N = 5, λ 1 = 20, λ 2 = 30, b=200, and other parameters were manually selected to obtain the most satisfactory restoration quality.For better comparison, HNHOTV-OGS and ℓ 0 -OGSTV keep within the scope of the authors suggest the parameters.For ℓ 1 -OGSTV, we manually select the parameters to get the best PSNR or SSIM value.
The effectiveness of this method on salt-and-peppe noise noise denoising is verified by experiments.Three different noise levels 30%, 50% and 70% are added to the test image respectively to generate each observed image.The obtained PSNR and SSIM values are shown in Table 1, Table 2 and Table 3.In each table, we observed that even at different noise levels, compared with the denoising results of the other three methods, the PSNR and SSIM values of the proposed method were almost higher than those of the other three methods.Only about the picture "sacroiliac" performs slightly worse than ℓ 0 -OGSTV.
In Figs.2-4, we show the comparison of the denoised images of the three methods for the MR Images "abdomen" and "pelvic" at 30% noise level, "sacroiliac" and "head" at 50% noise level, and "ankle" and "mrcp" at 70% noise level, respectively.From the results, ℓ 1 -OGSTV image denoising effect is not ideal, there are still slight blocky artifacts, and some important texture edges are blurred.The closest competitor to our method is ℓ 0 -OGSTV, whose recovered images remove impulse noise well and reduce staircase artifacts due to the use of ℓ 0 -norm.However, as the noise level increased, our method still achieved better results.Because our method similarly uses the ℓ 0 -norm and combines the advantages of the overlapping group sparse fractional-order total variation and MC penalty term, it can well remove staircase artifacts and preserve detail edges even at high noise levels.

PARAMETER SENSITIVITY ANALYSIS
The parameters that affect the performance mainly include the fractional order α, the group size K, the number of MM iterations N, MC penalty parameter b and parameters λ 1 , λ 2 , and these parameters need to be carefully tuned to get more accurate results.Therefore, we selected some images to test under different noise conditions to demonstrate the sensitivity of the proposed model to these parameters.
Firstly, in order to test the sensitivity of iteration number N, other parameters are fixed.Table 4 shows the influence of different MM iteration number N on PSNR, SSIM, overall algorithm iteration number and time, and the best result is obtained when N=5.
The variation of the value of the group size K is important for the quality of the recovered image, and it can be obtained from Fig. 5 that when the group size K=3, the performance is the best, and the image quality Image Anal Stereol 2024;43:53-66 Fig. 2: The first and third lines are the recovered results for "abdomen" and "pelvic" with the 30% noise level, respectively, while the second and fourth lines show the fragments corresponding to the zoomed images.( Fig. 3: The first and third lines are the recovered results for "sacroiliac" and "head" with the 50% noise level, respectively, while the second and fourth lines show the fragments corresponding to the zoomed images.does not continue to improve if the K value continues to increase.In our experiment, the range of values about fractional order α is 1 ≤ α < 2. In Fig. 6, we test two images at 30% and 50% noise levels respectively.The figures show the PSNR and SSIM values increase with the α value.Therefore, α = 1.9 can obtain the best PSNR and SSIM results.In this paper, the parameters λ 1 and λ 2 respectively control the weight of the overlapping group sparse fractional-order total variation and MC penalty terms.Because the new regularizer is based on the overlapping group sparse fractional-order total variational minus the MC penalty term, the new regularizer is more successful in generating sparsity.According to the experimental test, when λ 1 and λ 2 are [20,30], the denoising effect of the model is relatively ideal.In Fig. 7, the recovery results of "abdomen" at 50% noise level and "head" at 70% noise level are shown for different values of λ 1 and λ 2 , respectively, where (c1)-(c2) perform the best, and we set λ 1 =20 and Finally, the analysis is about the parameter b in the MC penalty term.As shown in Fig. 8, we test the effect of different b values on PSNR and SSIM values at two noise levels, and it is observed that b=200 is the optimal choice.

CONCLUSION
In this paper, we propose a novel denoising model for MR images aimed at effectively removing impulse noise and staircase artifacts.We demonstrate the effectiveness of using the ℓ 0 -norm as a data fidelity term to eliminate impulse noise, while the combination of the overlapping group sparse fractional-order total variation and MC penalty as regularizers can mitigate staircase artifacts and preserve important edges.To solve the proposed model, we employ the alternating direction method of multipliers.Numerical experiments validate that our proposed method outperforms three alternative methods in terms of PSNR and SSIM across varying levels of noise.In future work, we aim to extend our proposed method to address other types of noise and reduce computational time.

Table 1 :
The PSNR and SSIM values for denoised images by different methods when 30% noise level.

Table 2 :
The PSNR and SSIM values for denoised images by different methods when 50% noise level.

Table 3 :
The PSNR and SSIM values for denoised images by different methods when 70% noise level.

Table 4 :
Denoising results of different MM iterations (N).