MR Image Restoration by Utilizing Fractional-Order TV and Recursive Filtering

Authors: nana  wei
wei  xue
Venue: IEEE_org_ICIST_2024_Conference
Published: 2024-08-14 16:54:20

Abstract: Total variation based methods are effective models in magnetic resonance image restoration. For eliminating impulse noise, an effective way is to use the l0-norm total variation model. However, the TV image restoration consistently produces staircase artifacts, particularly at noise levels with high density. In this paper, we propose a novel MR image restoration model incorporating fractional-order regularization and filtering methods. Specifically, the first term employs the l0-norm as the data fidelity term to effectively eliminate impulse noise. The second term introduces a fractional-order total variation regularizer, which preserves structural information while mitigating staircase artifacts during deblurring. Given its suboptimal performance in texture detail recovery, we utilize recursive filtering for high-quality edge-preserving filtering. Finally, we solve the corresponding optimization model by using the alternating direction method of multipliers. Experimental results demonstrate the effectiveness of our method in restoring medical images.

Main paper: 
 i. introduction: 
 Magnetic Resonance (MR) Image is a powerful tool used in the field of medical imaging for the diagnosis and monitoring of various diseases. However, the quality of MR images can often be compromised by the presence of impulse noise, also known as impulse noise. This type of noise appears as sparsely occurring white and black pixels, and can significantly degrade the image quality, making it difficult for medical professionals to interpret the images accurately. Impulse noise in MR images can be caused by a variety of factors, including hardware faults, transmission errors, and other interferences. The challenge lies in developing effective methods to remove this noise while preserving the important details and edges in the images.To remove the impulse noise, a spectrum of approaches has been explored. Among them, regularization methods and filtering methods are widely recognized. A cornerstone in the realm of medical image restoration is the Total vriation (TV) model, which has demonstrated remarkable efficacy in denoising, deblurring, segmentation, and superresolution tasks. However, the TV model encounters limitations when This work was supported in part by the Natural Science Foundation of the Anhui Higher Education Institutions of China (Grant Nos. 2023AH040149 and 2022AH050310) and the Anhui Provincial Natural Science Foundation (Grant No. 2208085MF168). confronted with noise in piecewise affine regions, often manifesting as undesirable staircase artifacts. To overcome this issue, a series of regularization functions has been proposed. Such as, the High-Order TV method [3] , Non-Convex High-Order TV [4] , Fractional-Order TV (FOTV). [5] and [6] are among the notable advancements, their numerical experiments showed that their methods alleviated the staircase effect effectively. Besides, [2] proposed the Overlapping Group Sparse Fractional-Order TV (OGSTV), which has garnered attention for its refined approach. Besides, [10] proposed an enhanced ℓ 1 -NCFOTV method, which promises to elevate the restoration precision of images afflicted by impulse noise and to diminish the staircase effect.In the filtering domain, the Median Filter (MF) and Adaptive Median Filter (AMF) have been mainstays since the early days of denoising endeavors. Despite their popularity, these filters falter under the weight of high noise intensity, prompting the evolution of more sophisticated filtering techniques such as Recursive Filtering (RF), Recursive Mean Filtering, A nonlocal adaptive mean filter(NAMF) [8] and improved median filter [9] . Additionally, recursive filtering demonstrates effective performance in restoring edges.To effectively remove impulse noise, the ℓ 2 -norm [1] is sensitive to outliers and can easily result in unsatisfactory image restoration, which is commonly used to restore images degraded by additive Gaussian noise. Furthermore, there has been a growing interest in using an ℓ 1 -fidelity term instead of the ℓ 2 -fidelity term for image restoration in many literature such as [11] , [12] . The numerical evidence clearly indicates that the proposed method has made substantial advancements in its capacity for restoration. While the ℓ 1 -norm has been recognized for its considerable benefits in the realms of sparse signal processing and image restoration, it has been noted that it might excessively penalize the resulting solution in the context of impulse noise elimination, as highlighted by [7] .To circumvent this issue, a method for removing impulse noise using ℓ 0 total variation (ℓ 0 -TV) was proposed in [13] . This method employs the ℓ 0 -norm for the data fidelity and uses the proximal alternating direction method of multipliers for image restoration. The empirical evidence show better performance than ℓ 1 -norm based methods. Consequently, this indicates that ℓ 0 -norm is suitable for the restoration of images distorted by impulse noise. The image restoration of the optimization problem can be depicted as:EQUATIONwhere u ∈ R n×m is the desired original clean image, b ∈ 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, H ∈ R n×n is a linear operator. Another circumstance, when H = I, the identity operator, it translates into a denoising issue.In recent work, the ℓ 0 -norm has been employed as a data fidelity criterion for the elimination of impulse noise. [2] proposed the ℓ 0 -norm with an OGSTV and MC penalty model, which can significantly enhance the model's proficiency in image restoration, particularly under conditions characterized by high levels of impulse noise. [14] proposed ℓ 0 -OGSTV model effectively diminishes staircasing artifacts and precisely restores images marred by impulse noise. [15] has successfully integrated the ℓ 0 -norm data fidelity term with a nonconvex generalized regularization approach. This integration enhances the model's ability to preserve sharp image contours while simultaneously reducing staircase artifacts. Collectively, these models showcase the effectiveness of the ℓ 0 -norm in restoring images that have been corrupted by impulse noise.To the best of our knowledge, the combination of ℓ 0 -norm data fidelity, fractional-order total variation and recursive filtering has not been addressed by any of the existing optimization models. Therefore, this has motivated us to combine them for the restoration of MR images in medical imaging. In this article, we focus on addressing the challenge of restoring MR image in medical imaging that have been degraded by blurred boundaries and residual noise. The model comprises an ℓ 0norm data fidelity term to restore images from impulse noise, a regularizer of Fractional-Order total variation (FOTV) to eliminate the staircase effects and the recursive filtering (RF) to improve edge quality significantly. Our method capitalizes on the strengths of the Forward-Backward Total Variation to provide superior denoising capabilities, effectively mitigating the staircasing artifacts. Additionally, the recursive filtering technique employed not only refines the texture details within medical images, but also reduces the computational process. 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. In Section 2 introduces fundamental ideas and preliminary information pertinent to the proposed algorithm. In Section 3, We propose a novel framework aimed at eliminating blur and impulse noise 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.In this section, we briefly review a few key components of the ℓ 0 -norm fidelity term, the discrete fractional-order difference, the recursive filtering and the Alternating Direction Method of Multipliers framework.A. The ℓ 0 Fidelity Term Initially, we provide an overview of the essential definitions and characteristics associated with ℓ 0 -norm fidelity term. In the expression ∥o ⊙ (Hu -f )∥ 0 , the o ∈ {0, 1} n is specified by the user. Specifically, when o i = 0, it implies the pixel in position i is an outlier, while when o i = 1, it implies the pixel in position i is a potential outlier. For this paper, we set o i = 0, t i = u min or u max 1, otherwise for the salt-and-pepper impulse noise.The following lemma, as meticulously outlined in the work of [13] , elucidates the variational formulation that underpins the ℓ 0 -norm.Lemma 1: For any given w ∈ R n , it holds thatEQUATIONand 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 result of Lemma 1 implies that the ℓ 0 -norm minimization problem in (1) is equivalent toEQUATIONif u * is an overall optimal solution of (1), then (u * , 1 -sign (|Hu * -b|)) is overall optimal solution of (3). In a similar manner, if u * is a global Ideal solution of (3), then (u * , 1 -sign (|Hu * -b|)) is an overall optimal solution of (1). The ability of ℓ 0 -norm has proven to be remarkably effective in restoring images, offering superior performance.Given an image domain Ω ⊂ R 2 , we discretize it as a rectangular grid {(x i , y j ) :1 ≤ i ≤ M, 1 ≤ j ≤ N },Consequently, the image discretized on the Euclidean plane R M ×N , denoted as u i,j = u (x i , y j ). Let C α 0 (Ω) with α > 0 be the space of α-order continuously differentiable functions defined on Ω with compact support. Based on the GL fractional-order derivatives, Thus, the discrete form of the fractional-order gradient ∇ a u can be evaluated byEQUATIONwhere α is the fractional order and we set 1 ≤ α < 2 in this paper. where the fractional-order derivatives (D α x u), (D α y u) ∈ R M ×N along the x-axis and the y-axis are approximated by theEQUATIONIn this context, K denotes the quantity of adjacent pixels employed to calculate the fractional-order derivative for each pixel. The coefficients {C α k } K-1 k=0 are given byC α k = Γ(α+1) Γ(k+1)Γ(α+1-k)with the Gamma function Γ(x). Hence, the discrete fractional order TV of u is expressed in terms ofEQUATIONwhere (∇ α ) * = (-1) α div α is the conjugate operator of the fractional order gradient operator. In the discrete case, the vector div α p = (p (1) , p (2) ) ∈ R N ×M × R N ×M discrete fractional-order divergence is defined as [6] , [17] EQUATIONi+k,j + pi,j+k ). (7) Observe that the divergence ( 7) is the adjoint of the gradient (4).The first-order recursive filtering (RF) was initially introduced by [16] . Denoting I [u] and J [u] are the noisy image and the denoising image, respectively. RF computes J [u] recursively.EQUATIONwhere a ∈ [0, 1] is a feedback coefficient, and its implemen-EQUATIONwhere σ s and σ r represent the spatial and range parameters, respectively. The rapid iteration of the recursion in (8) is attributed to the pre-computed values from (9) and the independent computation of J [u] for each row. For gray images, one can sequentially perform (8) forward and backward in two directions. With color images, it's necessary to iterate through this operational sequence for each color channel. As d increases, a d goes to zero, stopping the propagation chain, thus preserving edges. Furthermore, the expanded recursion of Eq. (8) also elucidates this issueJ[u] = n ℓ=0 ℓ k=0 a d[u-k+1] (1 -a d[u-ℓ] )I[u -ℓ]. (10)The alternative direction method of multipliers (ADMM) is a computational framework for solving optimization problems , which is to solve the following constrained separable optimization problems:EQUATIONwhere f (•), g(•) : χ i → R are closed convex functions, A, B ∈ R l×n are linear transforms, χ i → R are nonempty closed convex sets, and d ∈ R l is a given vector. For problem (11) , we establish the augmented Lagrangian functionEQUATION)where µ ∈ R l is the Lagrange multiplier and λ > 0 is a penalty parameter which controls the linear constraint. The objective is to find the saddle point of L A by alternatively minimizing L A with respect to u, w and µ.The problem (11) is addressed by presenting the ADMM algorithm as Algorithm 1.Algorithm 1 ADMM for minimizing the problem (11) . Input: penalty parameter λ > 0, number of iterations. Initialize: Initial image u 0 = b, counter k = 0, Lagrange multipliers µ. Output: Restored image u.1: For k = 0, compute u k+1 , w k+1 , µ k+1 2: u k+1 = arg min u f (u) + λ 2 ||Au + Bw k -d + µ k λ || 2 2 , 3: w k+1 = arg min w g(w) + λ 2 ||Au k+1 + Bw -d + µ k λ || 2 2 , 4: µ k+1 = µ k + λ(Au k+1 + Bw k+1 -d), 5: k = k + 1 6:until a stopping criterion is satisfied.In this part, we first introduce the proposed MR image restoration model and the corresponding solution methods. Finally, an ADMM solution framework is provided.The proposed MR image restoration model is as followsEQUATIONwhere ∇ α u denotes the fractional-order TV. ϕ denotes the recursive filtering. λ 1 > 0 and λ 2 > 0 represent regularization parameters.Using variable splitting, the problem is rephrased as a constrained optimization problem that followsEQUATION)The corresponding augmented Lagrangian functional is given by LA(u, v, x, y, z, µv, µx, µy, µz) =< 1, 1 -v > +λ1||x|| + λ2ϕ(z)EQUATIONwhere variables µ v , µ x , µ y , and µ z are the Lagrange multipliers associated with the constraints of ( 14). β v , β x , β y , and β z > 0 are the corresponding penalty parameters. We utilize the alternating direction method of multipliers [18] to solve the proposed model (15) . According to the ADMM scheme, we can alternately solve for the following problems.1) u-subproblem: The u-subproblemEQUATION)Based on the first-order optimality conditions, we are tasked with resolving a system of linear equationEQUATIONwhere µ = µz βz , considering u with periodic boundary constraints. Due to the circulant and circulant blocks (BCCB) structure, matrices (∇ α ) T ∇ α and H T H can be diagonalized by 2D discrete fast Fourier transforms (F F T ). Therefore, solving for u can be efficiently solved using 2D F F T and 2D F F T inverse operations. The process for acquiring the optimal u is outlined as followsEQUATIONwhere κ = H T (β y (b+y)-µ y )+(∇ α ) T (β x x+µ x )+β z (z-µ), F and F -1 represent the Fourier transform and its inverse.2) v-subproblem: The v subproblem can be written asEQUATION)The v-subproblem in (19) is equivalent toEQUATIONtherefore, projection method is engaged to find the solutionEQUATION)This subproblem mentioned is a projection onto a convex set, ensuring that the pixel values of the restored image remain within the range of 0 to 1.3) y-subproblem: Solving the y-subproblem involves employing a soft thresholding technique coupled with a shrink operator. The formula is presented as followsEQUATIONEq. ( 22) holds equivalence to Eq. ( 23)EQUATIONBy expanding(23) and discarding the constant terms, we can rewrite (23) asEQUATIONThe minimizer is determined through the application of a fourdimensional shrinkage operator Simplify Eq. ( 24), we getEQUATIONwhere shrink(s, γ) = sgn(s)⊙max {∥s∥ 1 -γ, 0} , and sgn(•) denotes the signum function.x-subproblem: Variable x in (26) isupdated by solving the following problem:EQUATION)With the aid of the previously mentioned shrinkage operator, the x-subproblem can be directly tackled.The z-subproblem is a recursive filtering issue, with the specific form as followsEQUATIONwe define ϑ = λ2 βz and z (k) = u + µ, substitute it into Eq. ( 27). where z (k) donates the 'blurry' image, (28) minimizes the residue between z (k) and the 'clean' image z by employing the prior ϕ(z).EQUATIONExpanding on this insight, [19] proposed that the ADMM algorithm can be executed without the prior speciffication of ϕ. Consequently, we may utilize recursive filtering for the resolution of (28). denoted by ϕEQUATIONFinally, the Lagrange multipliers are updated by the followingEQUATIONAlgorithm 2 solving the minimization problem (13) .Input: Regularization parameter λ 2 > 0, α, penalty parameters β x , β y , β z , β v > 0, number of iterations. Initialize: Initial image u 0 = b, counter k = 0, Lagrange multipliers µ x , µ y , µ, µ v . Output: Restored image u.1: for k = 0 to number of iterations do 2:// Solve the subproblems: 3:Compute u k+1 according to Eq. (Compute v k+1 according to Eq. (Compute y k+1 according to Eq. (Compute x k+1 according to Eq. (26) 7:Compute z k+1 according to Eq. (// Update the Lagrange multipliers:9:Lagrange multipliers according to Eq. (30) 10:// Check stopping criterion:11: if ∥u k+1 -u k ∥ 2 ∥u k ∥ 2 ≤ 1 × 10 -4 then 12: break 13:end if 14: end for Our method is systematically presented in Algorithm 2. Additionally, there are two remarks that are worth noting regarding this algorithm.Remark 1: When α = 1, The fractional-order TV degrades to the conventional TV, then the Eq. ( 26) reduces to the conventional total variation regularization issue. In the experimental section, we provide an in-depth analysis of how the value of α impacts the noise reduction capabilities of our proposed model.Remark 2: When ϕ (u) = ∥u∥ T V (the total variation norm), then the Eq. ( 28) is a canonical total variation-based image denoising task [20] .In this section, we exhibit experimental results that verify the effectiveness of our proposed method for image restoration. The real MR test images are shown in Fig 1 . The experiment is under Windows 10 and MATLAB R2019b operating system, and the CPU is i5-8250U. The restored image quality was assessed through the peak signal-to-noise ratio (PSNR) and the structural similarity index (SSIM) [19] . Higher PSNR and SSIM values indicate better image quality.Our experiments employed a relative error-based stopping criterion for the algorithm.EQUATIONwhere u k+1 and u k are the restored image at the current iterate and previous iterate respectively.The performance of the model is influenced by a number of primary parameters, including the fractional order α, the parameters λ 1 , λ 2 , and penalty parameters of β v , β x , β y , and β z , additionally, these parameters must be meticulously adjusted to achieve higher precision in the outcomes.Here, we primarily discuss the FOTV term parameter α and the RF parameter λ 2 . Firstly, the value selected for parameter α is crucial importance and the fractional order α is 1 ≤ α < 2. In Fig. 2 the 'head' image was processed with a 5 × 5 Gaussian blur kernel (σ=5, impulse=70%) while the 'abdomen' image was subjected to a 7×7 Gaussian blur kernel (σ=10, impulse=50%). The figures show the PSNR and SSIM values increase with the α value. Therefore, for Gaussian blur, a value of α = 1.9 can obtain the best PSNR and SSIM results.In Fig. 3 , the 'ankle' image was processed with a 5 × 5 average blur kernel ( impulse=70%) while the 'head' image was subjected to a 7 × 7 average blur kernel ( impulse=50%). Our analysis of the images show that the optimal values for PSNR and SSIM are achieved at α = 1.3 value, thus our preference for this value in the context of average blur kernel. In this paper, the parameters λ 2 control the weight and also manage the magnitude of the value ϑ within the RF. Therefore, the selection of the λ 2 parameter is crucial. In Fig. 4 , we conducted experiments with the Gaussian and average blur kernel. The results show that the optimal effect is achieved with λ 2 = 0.001 when adding a Gaussian blur kernel, and with λ 2 = 0.005 when adding an average blur kernel.In this experimental section, This process involves adjusting one parameter at a time while maintaining the remaining parameters at their default values.In order to verify the effectiveness of our method in removing staircasing artifacts through fractional-order TV and further noise reduction and texture enhancement through RF, the subsequent sections will validate these capabilities. Our proposed method is called FOTVF ( 13), and the model that removes the third term ϕ(u)) from FOTVF is named as FOTV. The model is also solved using the ADMM algorithm. Fig. 5 and Fig. 6 show the PSNR and SSIM values of FOTVF and FOTV at 5 × 5 Gaussian blur kernel (std=5 and impulse=40%), and Fig. 6 is the restoration results of FOTVF and FOTV for 'head' at 5 × 5 average blur kernel (impulse=60%). It can be seen that the PSNR and SSIM values of FOTV, without the third term decrease a lot. Furthermore, it can be observed from these figures that the restoration effect of the FOTV is not ideal. In contrast, the FOTVF recovery image not only smooths out the staircase artifacts but also achieves a higher SSIM value.EQUATIONThe experimental results of our proposed model are compared with three related methods: ℓ 0 -OGSTV [14] , HNHOTV-OGS [3] , and Median Filter (built-in noise reduction technology in MATLAB).In this experiment, we assumed that the blurring kernel is known. To simulate a noisy blurred image, we blur the original images with a 5×5 Gaussian blur kernel with standard deviation σ = 5 and a 7×7 Gaussian blur kernel with standard deviation σ = 10, or two different sizes of average blur kernel 5 × 5 and 7 × 7 to the original images. Then, we introduce impulse noise with different densities to the blurred images.For parameters, we fixed λ 1 = 1. The other parameters were manually selected to obtain the most satisfactory restoration quality. The parameters of HNHOTV-OGS and ℓ 0 -OGSTV are consistent with the original text. Three different impulse noise levels of 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 I to Table IV .In each table, we can observe when using various intensities of noise and different blur kernels, our method almost always achieves higher PSNR and SSIM values compared to the other methods. Only in a few specific cases does ℓ 1 -OGSTV perform slightly better. From Table .I to Table .IV, it can be observed that our method performs well when relatively low levels of noise, Despite a modest reduction in denoising efficacy compared to ℓ 0 -OGSTV under conditions of high-density noise, it continues to deliver superior PSNR and SSIM values.For the Fig. 7 , we present the outcomes of three distinct denoising models applied to MR images. The image 'sacroiliac' is treated with average blur and 30%, 50%, and 70% impulse noise. In the results, the HNHOTV-OGS image denoising effect is not satisfactory; there are slight blocky artifacts, and the denoising performance is poor at higher noise densities.The main competition for our technique is the ℓ 0 -OGSTV model, which effectively eliminates impulse noise and mitigates staircase effects through the application of the ℓ 0 -norm. Despite the addition of varying degrees of Gaussian kernels and average kernels, our method consistently outperformed others, maintaining its advantage in producing better results.We proposed a MR image restoration model that combines regularization and filtering methods, aiming at effectively removing impulse noise and staircasing artifacts present in MR images. We demonstrated the effectiveness of using the ℓ 0 -norm as a data fidelity term to eliminate impulse noise, while incorporating fractional-order total variation and recursive filtering as penalty terms to mitigate staircasing artifacts and preserve important edges. We solved the proposed model using the alternating direction method of multipliers.In experiments involving Gaussian and average blurring, our method outperformed three other methods in terms of PSNR and SSIM across various levels of blur and noise. 
 
 Appendix:  



Reviews: review: title: title: This is a good paper.
main_review: main_review: This paper proposes a novel MR image restoration model incorporating fractional-order regularization and filtering methods, which eliminate impulse noise and preserves structural information while mitigating staircase artifacts during deblurring. Finally, this paper solves the corresponding optimization model by using the alternating direction method of multipliers. Experimental results demonstrate the effectiveness of the proposed method in restoring medical images. 
1). Please compare with existing papers to highlight the innovative of this article.
2). Since the value chosen for parameter $alpha$  is crucial, provide a guide to finding the optimal parameter $alpha$.
title: title: Manuscript Accept
main_review: main_review: Does the proposed MR image restoration model address the effect of the noise levels with high density?
What is the specific formula of the stopping criterion?  Algorithms 1 and 2 should keep a uniform format. 
What is the difficulty when the authors make the combination of ℓ0-norm data fidelity, fractional-order total variation and recursive filtering?
title: title: Well-founded approach
main_review: main_review: The paper presents a novel and well-founded approach to MR image restoration, combining fractional-order TV regularization with recursive filtering. This topic is interesting, the following comments need to further consider: (1) Simplify the language in the abstract and introduction, and provide brief explanations for technical terms such as "fractional-order TV" and "recursive filtering." This will help readers from different backgrounds understand the significance of the research. (2) Provide a brief summary of how the proposed method compares to existing approaches in terms of PSNR and SSIM. Mentioning the specific methods used for comparison and highlighting key results would provide context and reinforce the paper's contributions. (3) Include a brief discussion on the assumptions made in the model, such as the conditions under which the method is expected to perform best. Additionally, acknowledging any limitations of the approach would present a more balanced view of the research. (4) Discuss the potential practical applications of the method in more detail.