Abstract
This paper proposes a new anisotropic diffusion model in image restoration that is understood from a variational optimization of an energy functional. Initially, a family of new diffusion functions based on cubic Hermite spline is provided for optimal image denoising. After that, the existence and uniqueness of weak solutions for the corresponding Euler–Lagrange equation are proven in an appropriate functional space, and a consistent and stable numerical model is also shown. We complement this work by illustrating some experiments on different actual brain Magnetic Resonance Imaging (MRI) scans, showing the proposed model’s impressive results.
1. Introduction
During the last few decades, the digital image has been significantly utilized as a noninvasive medical diagnosis or treatment tool. For instance, a different method has been developed in the electroencephalogram (EEG) (see [1]). Since a digital image contains relevant information that should be carefully extracted from it, many engineers and mathematicians have established numerous techniques and theories [2–4] to collect quantitative or qualitative data. At that point, by analyzing these reliable data and understanding its contents, one can make the right decisions at the diagnosis time.
Along with the meaningful details, digital images carry false information as noise and other undesirable artifacts due to various reasons. This erroneous information can occur during image formation, transmission, or recording processes. Therefore, one needs good models to identify and eliminate these noises while preserving relevant information and structures.
Image denoising is a necessary preprocessing step for many applications that rely on image quality, such as image segmentation and pattern recognition. Hence, how to restore a degraded image to its original form? Many approaches have been developed to process images and showed great interest in using the variational approach [5–15], which attempts to establish mathematical models strictly related to the physical world by defining a diffusion equation from a particular optimization problem, which generally has the following form:where is an open bounded domain, and are the observed and the reconstructed images defined as functions of that associate each pixel to the gray level or , is a nonnegative increasing function with , and is a weighting parameter that enables to adjust the influence of the data term in the regularizing term.
One common approach to solve the problem (1) is by seeking the steady-state solution of the following PDE, which corresponds to the Euler–Lagrange equation of the energy functional :
Perona and Malik (Perona–Malik) [16] were the first to propose such a model in image processing, which corresponds to nonconvex energy functional. Unfortunately, this property leads to multiple solutions with staircase effects observed in practice [6, 17]. In [8], Charbonnier et al. suggested a strictly convex energy functional to circumvent the Perona–Malik model’s ill-posedness. Besides, many variational approaches have also been developed during the last thirty years. The most famous one was the total variation denoising suggested by Rudin et al. in their seminal paper [5]. Several improved nonlinear diffusion models for image denoising derived from Perona–Malik or TV models have been proposed in the last twenty years; for more detailed information, we refer to [2, 10–15, 18] and the references therein.
Nevertheless, because of the above model’s isotropy, the diffusion equation (2) is controlled by the function and produces the same amount of blurring in all its directions. This means that the process cannot successfully eliminate noises at the edges. Thus, it would be important to consider smoothing along edges by adopting anisotropic diffusion [19].
Furthermore, the type (1) model uses the gradient magnitude as a local descriptor operator for the image edge detector. However, digital images present some difficulties in their discrete structure; they are discrete in space and discrete in intensity value. Consequently, one may need to adapt to the digital image structure by considering differential operators that respond to vertical, horizontal, and diagonal edges and using a consistent approximation. Therefore, motivated by the above reasoning and inspired by Weickert’s anisotropic idea [19], we conceived new convex optimization problems that lead to anisotropic diffusion equations with a novel matrix diffusion tensor.
Therefore, this paper provides a new variational PDE model based on directional edge detectors in Section 2. The construction of a new diffusion function using cubic Hermite spline and the existence and uniqueness of weak solutions are illustrated in Section 3. Next, Section 4 implements a consistent explicit symmetric finite difference approximation that involves the neighborhoods at each point and shows the experiments on different real images. We terminate this work with a conclusion in Section 5.
2. Novel Anisotropic Diffusion Model for Image Restoration
Let us now consider the rectangular domain with being its boundary. This approach’s general idea is to search for a solution among the minima of some functional energy. The proposed functional energy has a regularizing term that depends on the combination of a function and local structures of the image expressed by the four directional derivatives:where is the canonical basis of . Hence, our problem consists of solving the following energy minimization problem:such thatunder the following assumptions:
The Euler–Lagrange equation of (5) is obtained fromwhich leads to
This is equivalent to
Let be a positive function such that . We may reformulate equation (9) and present it as follows.
Proposition 1. The Euler–Lagrange equation (9) is equivalent towhere is a real symmetric positive definite matrix of defined as follows:
Proof. First, it is clear thatOn the other hand, we haveBesides, the matrix has two eigenvalues as follows:Since andthen we deduce that the matrix is symmetric positive definite, which completes the proof.
By applying the steepest descent method, the Euler–Lagrange equation (10) can be solved considering the steady state of the following evolution equation:In [20], the existence and uniqueness of weak solutions for problem (16) were proven if the weighting parameter .
3. Mathematical Investigation of the Proposed Model
3.1. Description of the Proposed Model
Our model’s main objective is to allow strong directional smoothing within the areas where , , , or is small and prevent blurring boundaries, contours, or corners that separate neighboring areas, where one or a combination of these differential operators has substantial value. Some sufficient assumptions have to be made to ensure that our nonlinear diffusion meets acceptable behavior along these lines. Therefore, we can assume, as exhibited by Perona–Malik [16], some minimal hypotheses on :
Besides, due to the decreasing positivity of the function , it is evident that our model’s behavior (16) encourages smoothing along edges in the , , , or directions and ceases across them.
Moreover, we may provide another investigation analysis of our model by examining the eigenvalues-eigenvectors of the matrix . The matrix has two eigenvalues with being the corresponding eigenvectors:whereprovided that . We can then expand the first equation of (16) into
Accordingly, the diffusion caused by (20) is measured by the values and oriented toward . Specifically, it is clear from the expression of that , which means that the diffusion toward is privileged over . Furthermore, we can deduce the following results:(i)In flat areas, we have which means that . Then, the diffusion is isotropic and linear, and the smoothing effect is the same in all directions.(ii)At straight edges: .
For instance, we can assume that and . Then, we obtain and , which implies that
Consequently, the diffusion process is anisotropic and oriented along the direction.
(iii) Corners: .
According to the characteristics of the function , diffusion is anisotropic and oriented along and directions.
In fact, the difference gives insights into our diffusion model’s anisotropic property. In other words, it indicates the isotropic diffusion for the zero value and the anisotropic diffusion for larger values.
3.2. New Adaptive Diffusion Function Using Hermite Spline
For effective image denoising and better control of the diffusion via the process (16), we will approximate numerically the unknown diffusion function that meets the requirements by using Hermite spline [21]. Therefore, we can use cubic Hermite spline to interpolate numeric data specified at 0, , and to obtain a function that meets the requirements . To this end, we propose a function as follows:where and are the coefficients used to define the position and the velocity vector at a specific point, and are two threshold parameters, and is a family of the basis functions composed of polynomials of degree 3 used on the interval such that
And, we may consider
Hence, we can reformulate the expression of as follows:
By considering the continuity of at and , and using , it follows thenwhere
The values of the coefficients and are determined experimentally under the conditions of on the intervals and . Besides, we may introduce some sufficient conditions on the coefficients to ensure that the functions and satisfy on :
Now, we establish a useful growth condition for the function .
Proposition 2. For all , we can prove thatwhere and .
Proof. We know that . Then, : with . : with . : we have . Then, we may find an for an such that , for all .
3.3. Existence and Uniqueness
In this section, we investigate the existence and uniqueness of weak solutions of the Euler–Lagrange equation associated with the energy functional defined in (5). For and , we consider
First, we introduce a new functional space :
Next, we define a weak solution for problem (31).
Definition 1. A function with for is called a weak solution for problem (31) if, for any we haveAnd when is a constant function, we obtainBefore stating our main theorem, we need first to introduce a useful lemma.
Lemma 1 (uniform integrability and weak convergence [22]). Assume that is bounded and let be a sequence of functions in satisfyingSuppose alsoThen, there exist a subsequence and such thatNow, we state our main theorem.
Theorem 1. There exists a unique weak solution for problem (31).
Proof. We consider the variational problemwhereand the functional as defined in (5). It is obvious that knowing that .
Sincethen we can construct a minimizing sequence in such that andBesides,It follows thenOn the other hand, from Proposition 2, we may find and such thatThen, we getMoreover, we haveThen, for there exist and positive constants and (independent of ) such thatOn the other hand, we know thatThen, given , let and choose such that for all , we have . Hence, for , we obtainand this is true for all and arbitrary . It follows then that, for , we haveFrom (43), (45)–(50) and by Lemma 1 and the weak compactness of , we can find a subsequence of and a function such thatTherefore, we haveAdditionally, knowing that the function for is increasing and convex, the function is also convex for all . Therefore, we obtain for Integrating the above inequality over with with , we haveSince and by letting , we getThen, by letting , we deduceIt follows then . Besides, we havetAnd by following the same reasoning as set out above, we know that the function for is increasing and convex. Then, we can easily deduce thatTherefore, by combining (57) with (58), we concludewhich signifies that is a minimizer of the energy functional , i.e.,Furthermore, for all and for all , we have with . Then, , whereHence, we have , which meansBecause of (52), we haveWe conclude then that is a weak solution for problem (31).
Now, assuming that there is another minimizer of and using the fact that is strictly convex (5) and (6), we havea contradiction. Thus, there is only one minimizer, which completes the proof.
4. Numerical Implementation and Experimental Results
4.1. Consistency and Stability of Finite Difference Approximation
For a consistent and stable discretization scheme, one can use the following accurate finite difference scheme approximation: at time , , and mesh points , , we denote by the finite difference approximation of and exhibit the time-space derivatives discretization as follows:
By assuming and denoting
then we may approximate the solution for problem (16) by the above scheme to obtain the following discrete diffusion filter:where is the degraded image, , , and . Furthermore, we use the discrete Neumann boundary condition:
Then, using the filter (67) on every initial image yields a unique sequence [19]. Besides, due to the continuity of the function , for every finite , depends continuously on . Thus, under a specific condition, equation (67) satisfies the following maximum-minimum principle, which describes a stability property for the discrete scheme.
Theorem 2. (discrete extremum principle). For an iteration step satisfyingscheme (67) satisfiesfor all , , and .
Proof. For all , , and , we haveSince , , and , the following inequality holds:It follows thenSimilarly, we haveConsequently, the result (70) is true for becauseTherefore, by induction, it is easy to prove the result (70), which completes the proof.
4.2. Experimental Procedures and Results
4.2.1. Experimental Procedures
This section is devoted to comparing our model (16), as an image denoising algorithm, with the ones proposed by Wang and Zhou (Wang–Zhou) [10] and Maiseli [15]. All the experiments are conducted under Windows 10 and GNU Octave version 6.2.0, running on a laptop with an Intel® CoreTM i7-10510U (8 MB cache, 4 Core, up to 4.9 GHz), 16 GB memory (LPDDR3, Dual-channel, 2133 MHz), and 512 GB storage (PCIe, SSD, ). The experiments were done on three actual MRI scans (Figure 1) affected by different -values of zero-mean white Gaussian noise and restored using our filter (67) with the boundary-initial conditions (68), provided that the diffusion function (23) satisfies the assumptions .
(a)
(b)
(c)
We also used the discrete diffusion filter as described in [15] for Wang-Zhou and Maiseli models with the following diffusion functions: Wang–Zhou [10] diffusion function: Maiseli [15] diffusion function is a combination between the Perona–Malik diffusion function [16] and the Charbonnier diffusion function [8]:where and are positive constants.
Besides, in order to evaluate the quality of the restored images from different image denoising methods, we used two image quality metrics:(i)Peak signal-to-noise ratio [26]: where is the original or the uncorrupted image and is the distorted or the restored image. is one of the oldest image quality metrics evaluating an image’s signal strength relative to noise, and it is always positive. We evaluate the metric by using the Octave built-in function “psnr.” However, due to its limitations and its failure in some circumstances as an adequate measure of visual quality [27], we used another metric that significantly correlates with the human visual perception.(ii)Structural similarity index [28]: where and are tuning parameters. , , and stand for the mean, variance, and covariance, respectively. It is a method for measuring the similarity between a degraded image and a perfect one, and it is bounded between zero and one. For a good similarity between the original and the restored images, we need higher values of -index. In all experiments, we estimate the -index using “ssim.m” with automatic downsampling, downloaded from the website: https://ece.uwaterloo.ca/z70wang/research/ssim/.
In all experiments and to prioritize -index over , we used a combined metricto quantify the quality of the restored images. The iterations are stopped when the value reaches the maximum.
The image denoising process using the proposed model (67) can be implemented as shown in Algorithm 1. First, we consider as one of the three original brain MRI scans (Figure 1) and generate in it a Gaussian white noise with zero-mean and -value, using the Octave built-in function “imnoise.” Next, the initial value is set to be the noisy image in the while-loop. After that, a while-loop is executed, and in each iteration, our filter is used to build a new image with a new value. Finally, when the stopping criterion is confirmed, we obtain the restored image with the corresponding iteration.
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
4.2.2. Experimental Results
Several tests have been performed on the images in Figure 1 to obtain the best possible parameters: The time step: The weighting parameter: The threshold parameters for Maiseli diffusion function (77): The different parameters for our diffusion function (23):
Tables 1 and 2 show quantitative results on real images, corrupted with various white Gaussian noises, filtered by different filters. From Table 1, it is clear that the value of the restored images via Maiseli method has, in most cases, higher values with slight differences over the other methods. However, as shown in Table 2, when it comes to the -index, the proposed method reveals impressive results against the Wang-Zhou and Maiseli methods. Additionally, our method converges more rapidly to the solution than the others with the less iteration number.
From a visual comparison, Figure 2 shows that the denoising process using our diffusion function removes noises more efficiently and preserves essential image features. Besides, it can be seen from different images that all the restored images via Wang–Zhou and Maiseli methods are more affected by blocky artifacts than the proposed one, despite the higher values in . This can be attributed to the fact that the Wang–Zhou and Maiseli methods have lower -indices.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
(m)
(n)
(o)
5. Conclusion
This paper presented a new anisotropic diffusion model for image restoration, eliminating the corruption caused by white Gaussian noise. The existence and uniqueness of weak solutions of functions in Orlicz–Sobolev space that minimize the energy functional have been proven, provided that the functions and satisfy some specific conditions . Besides, we have established a consistent and stable numerical model for the denoising process and used the cubic Hermite spline to approximate the best possible diffusion function that revealed its efficiency regarding the optimal image denoising via (67). We have also proved that our method provides better results compared to the Wang–Zhou and Maiseli methods.
The next stage of the research attempts to prove the existence and uniqueness of weak solutions for the evolution problem (16) and use different numerical methods, such as finite element and mixed finite element methods.
Data Availability
No data were used to support this study.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.