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 [24] 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 [515], 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, 1015, 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 havewhich 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 then