Numerical Comparisons of Different Imaging Algorithms

Soulef Bougueroua^1,2, Nourreddine Daili^1,2

¹Department of Mathematics, Faculty of Sciences, University Ferhat Abbas Setif 1, Setif, Algeria.
²Fundamental and Numerical Mathematics Laboratory, University Ferhat Abbas Setif 1, Setif, Algeria.

DOI: 10.4236/jamp.2023.119175                

Abstract

Image processing is the set of operations performed to extract “information” from the image. An interesting problem in digital image processing is the restoration of degraded images. It often happens that the resulting image is different from the expected image. Our problem will therefore be to recover an image close to the original image from a poor quality image (that has been skewed by Gaussian and additive noise). There are several algorithms on how we can improve the broken image in better quality. We present in this paper our numerical results obtained with the models of Tikhonov regularization, ROF, Vese Osher, anisotropic and isotropic TV denoising algorithms.

Keywords

Comparative, Denoising Algorithms, Image Restoration

1. Introduction

A digital image is composed of basic units (referred to as pixels) that each represent a specific area of the image. The width and height of an image are determined by the practically infinite number of pixels that make up each of those dimensions, as well as the range of grayscale or colors that each pixel can contain (we speak of image dynamics). There are three categories of digital images:

• Binary images: in the simplest images, a pixel can only take the values black or white. When a piece of text just has one color, this form of image is typically used to scan it;

• The grayscale images: images with gray levels typically display 256 shades of gray. Simply put, each of the 256 colors in a 256-color image is defined by the range of gray. According to tradition, 0 represents black (null luminous intensity) and 255 represents white (maximum luminous intensity);

• The color images: in order to represent the colors red, green, and blue, a color image is actually made up of three images. Each of these three images is referred to as a canal. This representation in red, green, and blue mimics how the human visual system works.

The fundamentally ill-posed character of some practical problems is recognized and is manifested in a very large class of problems, called “inverse problems”. There are several types of ill-posed inverse problems, and their applications can be found in many fields such as image processing.

Digital image processing is one of the most crucial components of machine learning or computer vision. A fascinating area of digital image processing is the restoration of images. During the acquisition of image (especially through photography), it is typical for the final image to diverge from the expected image.

Image processing refers to all techniques and methods used to modify, remove Gaussian noise and improve or analyze digital images. It aims to extract useful information, improve visual quality, detect specific patterns or objects, and identify various imaging tasks. Image processing plays a crucial role in many areas of our daily lives, such as digital photography, television, medicine, computer vision, pattern recognition, augmented reality, robotics, surveillance, security and surveillance, geography. It plays an essential role in understanding and exploiting the visual information contained in the images.

There are many examples and fields of applications of image processing. The two main areas that have enabled image processing to develop are:

1) The military domain: missiles of all kinds (self-directed (short range), cruise (long range), …); intelligence (remote sensing from satellite images, the accuracy of which can go up to a few centimeters today, photo-interpretation); real simulators (aircraft, tank, …).

2) The medical field: medical imaging (ultrasound, MRI, tomography; angiography; x-ray; ultrasound; scanner; MRI, …).

Image processing is a discipline that concerns the manipulation, analysis and improvement of images using algorithms and computer techniques. Different techniques can be combined to get better results. There are many algorithms used in imaging to process and analyze images such as resizing, rotation, filtering, as well as more advanced operations such as segmentation, object detection, pattern recognition, and image restoration.

Image restoration algorithms are used to reduce Gaussian noise, remove artifacts, improve image resolution or overall image quality. Common methods include deconvolution techniques, model-based restoration, use of adaptive filters.

Denoising is the opposite problem from removing noise from an image; the outcome would be subpar if noise was left in the image. Noise is parasitic information that is added to the scene. Since noise has a wide range of origins and characteristics, it can be replicated in many different ways. There are many different kinds of noise, but the case study in this article is Gaussian additive noise with grayscale images, or $f=u+\eta$ . The original image is represented by u, the observed noisy image is represented by f, and the Gaussian random fluctuation to zero mean is represented by $\eta$ . Gaussian noise is often referred to as normal noise in a predefined density function. It is a common technique for including noise in images. According to the following definition, this noise can be produced randomly and separately inside the image with

$p\left(z\right)=\frac{1}{\sqrt{2\pi \sigma }}{\text{e}}^{-\frac{1}{2{\sigma }^2}{\left(z-\bar{z}\right)}^2}\mathrm{,}$ (1)

where z stands for intensity, $\bar{z}$ represents the mean value of z, and $\sigma$ stands for standard deviation. To return the image to a better level of visual quality, descalation techniques are necessary. The investigation of various image restoration models with Gaussian noises will be covered in this paper.

2. Some Models of Image Restoration

This work as an introduction to image recovery, which is an interesting and ill-positioned problem and is of crucial importance to the idea of image processing. The process of recovering an image nearly identical to the original from an observation, usually a fuzzy or blurred image of an real image, is known as Restore name image. Several recovery models have been applied to many problems that are poorly posed in the mathematical literature. Among them are the following models.

2.1. Tikhonov Regularization

The oldest regular method still in used to address inverse problems is the Tikhonov regularization method. In other words, we replace the ill-posed original problem with a well-posed alternative approximation problem. It is one of the most well-known methods of regularization in both statistical and digital analysis.

The Tikhonov regularization is a very commonplace yet overly simplistic regularization method for image processing. If we assume that the additive noise v is Gaussian and that f represents the observed image, then we attempt to reconstruct or restore the image u.

Let $V=H_0^1\left(\Omega \right)$ and $H=L^2\left(\Omega \right)$ , we take the original minimization problem (adjustment to the data):

$\left(\mathcal{P}\right){\alpha }^{\ast }\mathrm{<mo>\; :\; </mo>}=\underset{u\in V}{\min }{\Vert u-f\Vert }_H^2\mathrm{,}$ (2)

where $f\mathrm{<mo>\; :\; </mo>}\Omega \subset {ℝ}^N\to ℝ$ is the observed image and the following regularized problem:

$\left({\mathcal{P}}_{\alpha }\right)\beta \mathrm{<mo>\; :\; </mo>}=\underset{u\in V}{\min }\left\{{\Vert u-f\Vert }_H^2+\alpha {\Vert \nabla u\Vert }_H^2\right\}\mathrm{,}\alpha >0.$ (3)

The gradient must be “very minimal” in order for us to merely adjust u to the data f (it depends on the parameter). A slight gradient, in an image, is “smoothed”. The restoration will provide a blurry image because the margins are eroded.

Proposition 1. [1] Assume that $\left(\mathcal{P}\right)$ requires at least one answer $\tilde{u}$ . The problem $\left({\mathcal{P}}_{\alpha }\right)$ requires a one-of-a-kind solution $u_{\alpha }$ . When $\alpha \to 0$ , one can extract a subsequence from the family ${\left(u_{\alpha }\right)}_{\alpha }$ that converges (possibly) in V to a solution $u^{\ast }$ of $\left(\mathcal{P}\right)$ .

Proof. The solution exists, unique and therefore converges the well-posed regularized Tichonov problem. The proof of this proposition will be found in detail in [1] .

The restored image u is far less clear (in particular, the edges are eroded), which makes the problem of image restoration incompatible with the common expression for image restoration, $L\left(u\right)={\Vert \nabla u\Vert }_2^2$ (Tikhonov regularization). Think about the overall variation, or consider $L\left(u\right)={\displaystyle \int }\left|Du\right|$ . This strategy is significantly more successful. With respect to the problem of functions of bounded variation spaces, this leads to a functional minimization in a particular Banach space.

2.2. The Continuous Model of Rudin-Osher-Fatemi

Rudin, Osher, and Fatemi (ROF) proposed the first image restoration model from a given noisy image having additive noise using regularization (TV), which is defined by

$TV\left(u\left(x\mathrm{,}y\right)\right)\mathrm{<mo>\; :\; </mo>}=\left\{{\displaystyle {\int }_{\Omega }}\left|\nabla u\left(x\mathrm{,}u\right)\right|\text{d}x\text{d}y\mathrm{,}\text{with}\left|\nabla u\right|=\sqrt{u_x^2+u_y^2}\right\}\mathrm{.}$

The regularization of total variation (TV) approach of image processing is used to reduce noise from digital images. (TV) is a technique that was originally developed by ROF, it has since been applied to a multitude of other imaging problems.

Rudin, Osher, and Fatemi developed the method known as (TV) to address the problem of visual degradation. Now, it has been used to solve numerous additional image problems.

In [2] , a model has been proposed by Rudin, Osher, and Fatemi and in which the image is divided into two parts $f=u+v$ , where u is an unknown image and v is the noise. f is a brilliant measure that is usual at the beginning of a clean image, and is an agreement parameter. We will thus try to solve the problem and simply apply the regularization to the “ noise” portion using the $u+v$ formula with $u\in BV\left(\Omega \right)$ and $v\in L^2\left(\Omega \right)$ . If $f\in L^2$ is correct, the ROF problem is well-posed and the minimizer u exists, unique, and stable in $L^2\left(\Omega \right)$ . ROF proposed the following minimization problem:

$\left({\mathcal{P}}_{ROF}\right){\alpha }_{RoF}\mathrm{<mo>\; :\; </mo>}=\underset{u}{\inf }\left\{J\left(u\right)+\frac{1}{2\lambda }{\Vert v\Vert }_2^2\mathrm{<mo>\; :\; </mo>}u\in BV\left(\Omega \right)\mathrm{,}v\in L^2\left(\Omega \right)\mathrm{,}f=u+v\right\}\mathrm{.}$ (4)

This results in a $\left(BV\left(\Omega \right)\mathrm{,}{L}^2\left(\Omega \right)\right)$ decomposition of the image f.

$J\left(u\right)$ denotes the total variance of u and $\lambda >0$ .

$J\left(u\right)=\sup \left\{{\displaystyle {\int }_{\Omega }}u\left(x\right)div\left(\phi \left(x\right)\right)\text{d}x\mathrm{<mo>\; :\; </mo>}\phi \in C_c^1\left(\Omega \mathrm{,}{ℝ}^2\right)\mathrm{,}{\Vert \phi \Vert }_{\infty }\le 1\right\}\mathrm{.}$

Also known as BV, or functions of bounded variation space, according to

$BV\left(\Omega \right)=\left\{u\in L^1\left(\Omega \right),J\left(u\right)<+\infty \right\}.$

Here $J\left(u\right)$ denotes the TV of u and $\lambda >0$ is a weight parameter.

Theorem 2. [2] [3] The problem $\left({\mathcal{P}}_{ROF}\right)$ requires a single solution, which is provided by

$u=f-\lambda {\Pi }_{\lambda K}\left(f\right)\mathrm{,}$ (5)

where $\Pi$ is the orthogonal projector on $\lambda K$ (dilatation of K by $\lambda$ ), and K is the overall closure in $L^2$ .

$K:=\left\{div\left(\phi \right):\phi \in C_c^1\left(\Omega ,{ℝ}^2\right),{\Vert \phi \Vert }_{\infty }\le 1\right\}.$

2.3. Meyer’s Model

In [4] , Yves Meyer shows that if $\lambda$ is small enough, the ROF model will erase the texture. Yves Meyer suggests the use of a space of functions, which is in some ways the dual of the BV space, to extract both the u component in BV and the v component as an oscillating function (texture or noise) from f. The following definition is given by Meyer.

Definition 1. [4] $G\left({ℝ}^2\right)$ is a Banach space made of v distributions that may be written

$G\left({ℝ}^2\right)\mathrm{<mo>\; :\; </mo>}=\left\{v\left(x\mathrm{,}y\right)={\partial }_x{g}_1\left(x\mathrm{,}y\right)+{\partial }_y{g}_2\left(x\mathrm{,}y\right)\mathrm{<mo>\; /\; </mo>}{g}_1\mathrm{,}{g}_2\in L^{\infty }\left({ℝ}^2\right)\right\}\mathrm{.}$

We will see that the space G allows for oscillating functions v, as justified by Meyer, and that the oscillations are well measured by the norm

$\begin{array}{l}{\Vert v\Vert }_G\mathrm{<mo>\; :\; </mo>}=\inf \{{\Vert g\Vert }_{L^{\infty }\left({ℝ}^2\right)}=\text{ess}\underset{x\in {ℝ}^2}{\sup }\left|g\left(x\right)\right|\mathrm{<mo>\; /\; </mo>}v=divg\mathrm{,}\\ \stackrel{}{}g=\left(g_1\mathrm{,}{g}_2\right)\in L^{\infty }\left({ℝ}^2\right)\times L^{\infty }\left({ℝ}^2\right)\mathrm{,}\left|g\right|=\sqrt{g_1^2+g_2^2}\}\mathrm{.}\end{array}$

Meyer suggests the following new image restoration model:

$\left({\mathcal{P}}_{Meyer}\right){\alpha }_{Meyer}\mathrm{<mo>\; :\; </mo>}=\underset{u}{\inf }\left\{J\left(u\right)+\alpha {\Vert v\Vert }_G\mathrm{<mo>\; /\; </mo>}u\in BV\left(\Omega \right)\mathrm{,}v\in G\left(\Omega \right)\mathrm{<mo>\; ;\; </mo>}f=u+v\right\}\mathrm{.}$ (6)

$J\left(u\right)={\displaystyle \int }\left|\nabla u\right|$ denotes the total variation of u and $\alpha >0$ , while $G\left({ℝ}^2\right)$ denotes the space of oscillating functions.

Description of the Model

The interest in this space of oscillating functions stems from the fact that a strongly oscillating image with a small average norm in $G\left({ℝ}^2\right)$ , can have large oscillations but a small average norm, and that the $L^2\left(\Omega \right)$ norm is not the best choice for capturing the oscillating portion of an image. That is why he created a new space that was better suited from the start the G oscillating functions space. We have ${\Vert \mathrm{.}\Vert }_G$ low for oscillating functions with a null average and high for geometric functions.

We have arrived at the following conclusion based on a close approximation of the $L^{\infty }$ standard:

${\Vert \sqrt{g_1^2+g_2^2}\Vert }_{L^{\infty }}=\underset{p\to \infty }{lim}{\Vert \sqrt{g_1^2+g_2^2}\Vert }_{L^p}\mathrm{,}{g}_1\mathrm{,}{g}_2\in L^{\infty }\left({ℝ}^2\right)\mathrm{.}$

Then, if $\lambda ,\mu >0$ are tuning parameters, $\lambda \to \infty$ and $p\to \infty$ the approximation of Meyer model is given by

$\underset{u\mathrm{,}{g}_1\mathrm{,}{g}_2}{\inf }\left\{G_p\left(u\mathrm{,}{g}_1\mathrm{,}{g}_2\right)={\displaystyle \int }\left|\nabla u\right|+\lambda {\displaystyle \int }{\left|f-u-{\partial }_x{g}_1-{\partial }_y{g}_2\right|}^2\text{d}x\text{d}y+\mu {\left[{\left(\sqrt{g_1^2+g_2^2}\right)}^P\right]}^{\frac{1}{P}}\right\}$ (7)

where

${\displaystyle \int }\left|\nabla u\right|$ insures that $u\in BV\left({ℝ}^2\right)$ ,

${\displaystyle \int }{\left|f-u-{\partial }_x{g}_1-{\partial }_y{g}_2\right|}^2\text{d}x\text{d}y$ insures that $f\approx u+div\left(g\right)$ ,

$\mu {\left[{\left(\sqrt{g_1^2+g_1^2}\right)}^p\right]}^{\frac{1}{p}}$ is a penalty on the norm $v=div\left(g\right)$ in G.

As a result, the form of the Euler-Lagrange equation is given here.

$u=f-{\partial }_x{g}_1-{\partial }_y{g}_2+\frac{1}{2\lambda }div\left(\frac{\nabla u}{\left|\nabla u\right|}\right)\mathrm{.}$ (8)

$\mu {\left({\Vert \sqrt{g_1^2+g_2^2}\Vert }_p\right)}^{1-p}{\left(\sqrt{g_1^2+g_2^2}\right)}^{p-2}{g}_1=2\lambda \left[\frac{\partial }{\partial x}\left(u-f\right)+{\partial }_{xx}^2{g}_1+{\partial }_{xy}^2{g}_2\right]\mathrm{.}$ (9)

$\mu {\left({\Vert \sqrt{g_1^2+g_2^2}\Vert }_p\right)}^{1-p}{\left(\sqrt{g_1^2+g_2^2}\right)}^{p-2}{g}_2=2\lambda \left[\frac{\partial }{\partial y}\left(u-f\right)+{\partial }_{xy}^2{g}_1+{\partial }_{yy}^2{g}_2\right]\mathrm{.}$ (10)

2.4. Vese-Osher Model

Vese and Osher, who were the first to propose an approach to solve Meyer’s problem numerically; that is to say to realize the program, they used the approximation of Meyer model as follows:

$\begin{array}{l}\left({\mathcal{P}}_{Vese\mathrm{<mo>\; \_\; </mo>}Osher}\right){\alpha }_{Vese\mathrm{<mo>\; \_\; </mo>}Osher}\mathrm{<mo>\; :\; </mo>}=\underset{\left(u\mathrm{,}v\right)\in BV\times G\left(\Omega \right)}{\inf }\{J\left(u\right)+\lambda {\Vert f-u-v\Vert }_2^2+\mu {\Vert v\Vert }_G\mathrm{<mo>\; /\; </mo>}\\ {{\displaystyle }}^{}u\in BV\left(\Omega \right)\mathrm{,}v\in G\left(\Omega \right)\mathrm{<mo>\; ;\; </mo>}u+v=f\}\mathrm{.}\end{array}$ (11)

In our numerical calculations, the steps to calculate the solution of this problem are:

1) Replace the term ${\Vert v\Vert }_G$ by ${\Vert \sqrt{g_1^2+g_2^2}\Vert }_p$ with $v=div\left(g_1\mathrm{,}{g}_2\right)$ ;

2) $p=1$ is used for digital resonances because it allows for faster calculations each iteration;

3) Gives the equation of Euler-Lagrange;

4) We apply a fixed point iterative technique with a finite difference semi-implicit scheme.

2.4.1. The Numerical Discretization of Meyer’s Model

The numerical discretization of Equations (8), (9) and (10) is performed using the semi-implicit method of difference and the iterative algorithm based on the fixed point. We used the following initial values for the iterative algorithm:

$h=\mathrm{1,}p=1\text{and}n=100.$

$\{\begin{array}{l}{u}_0=f;\\ g_1^0=-\frac{1}{2\lambda }\frac{f_x}{\left|\nabla f\right|};\\ g_2^0=-\frac{1}{2\lambda }\frac{f_y}{\left|\nabla f\right|}.\end{array}$

The following concepts are used: $u_{i\mathrm{,}j}=u\left(ih\mathrm{,}jh\right)$ , $f_{i\mathrm{,}j}=f\left(ih\mathrm{,}jh\right)$ , $g_{\mathrm{1,}i\mathrm{,}j}=g_1\left(ih\mathrm{,}jh\right)$ , with the step $h>0$ and the point $\left(ih\mathrm{,}jh\right)$ for all $0\le i\mathrm{,}j\le M$ , and the variable change is taken.

$H\left(g_1\mathrm{,}{g}_2\right)={\left({\Vert \sqrt{g_1^2+g_2^2}\Vert }_p\right)}^{1-p}{\left(\sqrt{g_1^2+g_2^2}\right)}^{p-2}\mathrm{.}$

So the discretization of Equations (8), (9) and (10) is given by

$\begin{array}{l}{u}_{i\mathrm{,}j}^{n+1}=\frac{1}{1+\frac{1}{2\lambda h^2}\left(c_1+c_2+c_3+c_4\right)}(f_{i\mathrm{,}j}-\frac{g_{\mathrm{1,}i+\mathrm{1,}j}^n-g_{\mathrm{1,}i-\mathrm{1,}j}^n}{2h}-\frac{g_{\mathrm{2,}i\mathrm{,}j+1}^n-g_{\mathrm{2,}i\mathrm{,}j-1}^n}{2h}\\ +\frac{1}{2\lambda h^2}\left(c_1{u}_{i+\mathrm{1,}j}^n+c_2{u}_{i-\mathrm{1,}j}^n+c_3{u}_{i\mathrm{,}j+1}^n+c_4{u}_{i\mathrm{,}j-1}^n\right))\mathrm{<mo>\; ;\; </mo>}\end{array}$ (12)

$\begin{array}{l}{g}_{\mathrm{1,}i\mathrm{,}j}^{n+1}=\frac{2\lambda }{\mu H\left(g_{\mathrm{1,}i\mathrm{,}j}\mathrm{,}{g}_{\mathrm{2,}i\mathrm{,}j}\right)g_{\mathrm{1,}i\mathrm{,}j}}(\frac{u_{i+\mathrm{1,}j}^n-u_{i-\mathrm{1,}j}^n}{2h}-\frac{f_{i+\mathrm{1,}j}-f_{i-\mathrm{1,}j}}{2h}\\ +\frac{g_{\mathrm{1,}i+\mathrm{1,}j}^n-2g_{\mathrm{1,}i\mathrm{,}j}^{n+1}+g_{\mathrm{1,}i-\mathrm{1,}j}^n}{h^2}\\ +\frac{1}{4h^2}\left(g_{\mathrm{2,}i+\mathrm{1,}j+1}^n+g_{\mathrm{2,}i-\mathrm{1,}j-1}^n-g_{\mathrm{2,}i+\mathrm{1,}j-1}^n-g_{\mathrm{2,}i-\mathrm{1,}j+1}^n\right))\mathrm{<mo>\; ;\; </mo>}\end{array}$ (13)

$\begin{array}{l}{g}_{\mathrm{2,}i\mathrm{,}j}^{n+1}=\frac{2\lambda }{\mu H\left(g_{\mathrm{1,}i\mathrm{,}j}\mathrm{,}{g}_{\mathrm{2,}i\mathrm{,}j}\right)g_{\mathrm{2,}i\mathrm{,}j}}(\frac{u_{i\mathrm{,}j+1}^n-u_{i\mathrm{,}j-1}^n}{2h}-\frac{f_{i\mathrm{,}j+1}-f_{i\mathrm{,}j-1}}{2h}\\ +\frac{g_{\mathrm{2,}i\mathrm{,}j+1}^n-2g_{\mathrm{2,}i\mathrm{,}j}^{n+1}+g_{\mathrm{2,}i\mathrm{,}j-1}^n}{h^2}\\ +\frac{1}{4h^2}\left(g_{\mathrm{1,}i+\mathrm{1,}j+1}+g_{\mathrm{1,}i-\mathrm{1,}j-1}^n-g_{\mathrm{1,}i+\mathrm{1,}j-1}^n-g_{\mathrm{1,}i-\mathrm{1,}j+1}^n\right))\mathrm{.}\end{array}$ (14)

The following notations are used:

$\begin{array}{l}{c}_1=\frac{1}{\sqrt{{\left(\frac{u_{i+1,j}^n-u_{i,j}^n}{h}\right)}^2+{\left(\frac{u_{i,j+1}^n-u_{i,j-1}^n}{2h}\right)}^2}};\\ c_2=\frac{1}{\sqrt{{\left(\frac{u_{i,j}^n-u_{i-1,j}^n}{h}\right)}^2+{\left(\frac{u_{i-1,j+1}^n-u_{i-1,j-1}^n}{2h}\right)}^2}};\\ c_3=\frac{1}{\sqrt{{\left(\frac{u_{i+1,j}^n-u_{i-1,j}^n}{2h}\right)}^2+{\left(\frac{u_{i,j+1}^n-u_{i,j}^n}{h}\right)}^2}};\\ c_4=\frac{1}{\sqrt{{\left(\frac{u_{i+1,j-1}^n-u_{i-1,j-1}^n}{2h}\right)}^2+{\left(\frac{u_{i,j}^n-u_{i,j-1}^n}{h}\right)}^2}}.\end{array}$ (15)

2.4.2. Solution of Vese-Osher Problem

In order to solve the Vese-Osher problem, we will study this final problem in the discriminating case, when the image is a vector with two dimensions of size $N\times N$ , the Eulidian space $X={ℝ}^{N\times N}$ and $Y=X\times X$ .

If $u\in X$ then $\nabla u\in Y$ is defined by ${\left(\nabla u\right)}_{i\mathrm{,}j}=\left({\left(\nabla u\right)}_{i\mathrm{,}j}^1\mathrm{,}{\left(\nabla u\right)}_{i\mathrm{,}j}^2\right)$ , where

${\left(\nabla u\right)}_{i\mathrm{,}j}^1=\{\begin{array}{l}{u}_{i+\mathrm{1,}j}-u_{i\mathrm{,}j}\text{if}i<N\mathrm{,}\\ 0\text{if}i=N\mathrm{,}\end{array}\text{and}{\left(\nabla u\right)}_{i\mathrm{,}j}^2\mathrm{<mo>\; =\; </mo>}\{\begin{array}{l}{u}_{i\mathrm{,}j+1}-u_{i\mathrm{,}j}\text{if}j<N\mathrm{,}\\ 0\text{if}j=N\mathrm{.}\end{array}$ (16)

In the discriminating case, the total variance (TV) of u is defined as

$J_d\left(u\right)={\displaystyle {\sum }_{1\le i\mathrm{,}j\le N}}\left|{\left(\nabla u\right)}_{i\mathrm{,}j}\right|\mathrm{.}$

The divergence operator is $div=-{\nabla }^{\ast }$ (the adjoint of $\nabla$ ) where

${\left(div\left(p\right)\right)}_{i\mathrm{,}j}=\{\begin{array}{l}{p}_{i\mathrm{,}j}^1-p_{i-\mathrm{1,}j}^1\text{if}1<i<N\\ p_{i\mathrm{,}j}^1\text{if}i=1\\ -p_{i-\mathrm{1,}j}^1\text{if}i=N\end{array}+\{\begin{array}{l}{p}_{i\mathrm{,}j}^2-p_{i\mathrm{,}j-1}^2\text{if}1<j<N\\ p_{i\mathrm{,}j}^2\text{if}j=1\\ -p_{i\mathrm{,}j-1}^2\text{if}j=N\mathrm{,}\end{array}$ (17)

and the space $G^d$ is defined by

$G^d\mathrm{<mo>\; :\; </mo>}=\left\{v\in X\mathrm{<mo>\; /\; </mo>}\exists g\in Y\text{such}\text{that}v=div\left(g\right)\right\}\mathrm{.}$

Note that

$G_{\mu }\left(\Omega \right)\mathrm{<mo>\; :\; </mo>}=\left\{v\in G\left(\Omega \right)\text{such}\text{that}{\Vert v\Vert }_G\le \mu \right\}\mathrm{,}$

$G_{\mu }^d\left(\Omega \right)\mathrm{<mo>\; :\; </mo>}=\left\{v\in G^d\left(\Omega \right)\text{such}\text{that}{\Vert v\Vert }_{G^d}\le \mu \right\}\mathrm{.}$

As $J_d^{\ast }$ is the indicator function of $G_1^d\left(\Omega \right)$ defined by

$J_d^{\ast }\left(v\right)={\chi }_{G_1^d}\left(v\right)=\{\begin{array}{l}0\text{if}v\in G_1^d\\ +\infty \text{else}\text{.}\end{array}$

So to solve the Vese-Osher problem, we propose the following algorithm.

Algorithm 1. The algorithm for solving Vese Osher problem

${\inf }_{\left(u\mathrm{,}v\right)\in BV\left(\Omega \right)\times G_{\mu }\left(\Omega \right)}\left\{F_{\lambda \mathrm{,}\mu }\left(u\mathrm{,}v\right)=\left\{\begin{array}{c}J\left(u\right)+\frac{1}{2\lambda }{\Vert f-u-v\Vert }_2^2\text{if}v\in G_{\mu }\left(\Omega \right)\\ +\infty \text{if}v\in G\left(\Omega \right)\mathrm{<mo>\; /\; </mo>}{G}_{\mu }\left(\Omega \right)\end{array}\right\}\right\}\mathrm{.}$

By description;

${\inf }_{\left(u\mathrm{,}v\right)\in X\times X}\left\{F_{\lambda \mathrm{,}\mu }\left(u\mathrm{,}v\right)=\left\{\begin{array}{c}{J}_d\left(u\right)+\frac{1}{2\lambda }{\Vert f-u-v\Vert }_X^2\text{if}v\in G_{\mu }^d\left(\Omega \right)\\ +\infty \text{if}v\in X\mathrm{<mo>\; /\; </mo>}{G}_{\mu }^d\left(\Omega \right)\end{array}\right\}\right\}\mathrm{.}$

${\inf }_{\left(u\mathrm{,}v\right)\in X\times X}F\left(u\mathrm{,}v\right)={\inf }_{\left(u\mathrm{,}v\right)\in X\times X}\left\{J_d\left(u\right)+\frac{1}{2\lambda }{\Vert f-u-v\Vert }_X^2+J_d^{\ast }\left(\frac{v}{\mu }\right)\right\}\mathrm{.}$

We divide the problem into two sub-problems:

Pbm1 u solution, v fixed: ${\inf }_{v\in X}\left\{J_d\left(u\right)+\frac{1}{2\lambda }{\Vert f-u-v\Vert }_X^2\right\}\mathrm{.}$

${\curvearrowright }_{\begin{array}{c}\text{accordingtoROF}\\ \text{changeofvariable}\\ f=f-v\end{array}}^{\text{thesolutionis}}\hat{u}=f-v-P_{G_{\mu }^d}\left(f-v\right)\mathrm{.}$

Pbm2 v solution, u fixed: ${\inf }_{v\in G_{\mu }^d\left(\Omega \right)}\left\{{\Vert f-u-v\Vert }_X^2\right\}\mathrm{.}$

${\curvearrowright }_{\begin{array}{c}\text{changeofvariable}\\ f=f-u\end{array}}^{\begin{array}{c}\text{thesolutionis}\\ \text{accordingtoROF}\end{array}}\hat{v}=P_{G_{\mu }^d}\left(f-u\right)\mathrm{.}$

Lemma 3. There is an unique solution $\left(\hat{u}\mathrm{,}\hat{v}\right)\in X\times G_{\lambda }^d$ that minimizes $F_{\lambda \mathrm{,}\mu }\left(u\mathrm{,}v\right)$ in $X\times G_{\lambda }^d$ .

2.5. The Split Bregman Algorithm

Goldstein and Osher first proposed the split Bregman algorithm in [5] to handle more general form optimization problems:

$\varpi \mathrm{<mo>\; :\; </mo>}=\underset{u\in X}{\min }\left\{H\left(u\right)+{\Vert \Phi \left(u\right)\Vert }_1\right\}\mathrm{,}$ (18)

where X is a closed convex set and $\Phi \mathrm{<mo>\; :\; </mo>}X\to ℝ$ , $H\mathrm{<mo>\; :\; </mo>}X\to ℝ$ are convex functions. This problem is the same as the stress minimization problem as below:

$\bar{\varpi }\mathrm{<mo>\; :\; </mo>}=\underset{u\in X\mathrm{,}d\in ℝ}{\min }\left\{H\left(u\right)+{\Vert d\Vert }_1\right\}\text{such}\text{that}d=\Phi \left(u\right)\mathrm{.}$ (19)

Goldstein and Osher introduced the split Bregman algorithm, which was written as follows:

Algorithm 2. The split Bregman algorithm

Initialization: $k=0,u^0=0,b^0=0$ .

While $\Vert u^k-u^{k-1}\Vert >tol$ do,

$u^{k+1}={\min }_{u}H\left(u\right)+\frac{\lambda }{2}{\Vert d^k-\Phi \left(u\right)-b^k\Vert }_2^2$ ,

$d^{k+1}={\min }_d\left|d\right|+\frac{\lambda }{2}{\Vert d-\Phi \left(u^{k+1}\right)-b^k\Vert }_2^2$ ,

$b^{k+1}=b^k+\left(\Phi \left(u^{k+1}\right)-d^{k+1}\right)$ ,

$k=k+1$ ,

End while.

The split Bregman algorithm is used to solve some of the most common form optimization problems:

$\hat{\varpi }\mathrm{<mo>\; :\; </mo>}=\underset{u\in X}{\min }\left\{z\left(u\right)+\frac{1}{2}{\Vert u-f\Vert }_2^2\right\}\mathrm{.}$ (20)

Anisotropic and isotropic TV denoising problems are solved using the split Bregman method.

2.5.1. Anisotropic TV Denoising Problem

The problem of anisotropic TV denoising is considered in [6] .

$\left({\mathcal{P}}_1\right){\tau }_1\mathrm{<mo>\; :\; </mo>}=\underset{u}{\min }\left\{{\Vert \frac{\partial u}{\partial x}\Vert }_1+{\Vert \frac{\partial u}{\partial y}\Vert }_1+\frac{\mu }{2}{\Vert u-f\Vert }_2^2\right\}\mathrm{,}$ (21)

where f is the noisy image, $\frac{\partial u}{\partial x}$ and $\frac{\partial u}{\partial y}$ will be noted by $u_x$ and $u_y$ respectively. The problem is solved using a constraint equivalent to a problem $\left({\mathcal{P}}_1\right)$ . We answer the problem $\left({\mathcal{P}}_2\right)$ as follows:

$\left({\mathcal{P}}_2\right)\{\begin{array}{l}{\tau }_2\mathrm{<mo>\; :\; </mo>}=\underset{u}{\min }{\Vert d_x\Vert }_1+{\Vert d_y\Vert }_1+\frac{\mu }{2}{\Vert u-f\Vert }_2^2\\ \text{subject}\text{to}{d}_x=u_x\mathrm{,}{d}_y=u_y\mathrm{.}\end{array}$

The split Bregman algorithm can be used to tackle this last problem:

$\left({\mathcal{P}}_3\right){\tau }_3\mathrm{<mo>\; :\; </mo>}=\underset{u\mathrm{,}{d}_x\mathrm{,}{d}_y}{\min }\left\{{\Vert d_x\Vert }_1+{\Vert d_y\Vert }_1+\frac{\mu }{2}{\Vert u-f\Vert }_2^2+\frac{\lambda }{2}{\Vert d_x-u_x\Vert }_2^2+\frac{\lambda }{2}{\Vert d_y-u_y\Vert }_2^2\right\}\mathrm{.}$ (22)

We use

$shrink\left(x\mathrm{,}a\right)=\{\begin{array}{l}x-a\text{if}x>a\mathrm{,}\\ x+a\text{if}x<-a\mathrm{,}\\ 0\text{else}\text{.}\end{array}$ (23)

The Gauss-Seidel function is also useful.

$\begin{array}{l}{G}_{i\mathrm{,}j}^k=\frac{\lambda }{\mu +4\lambda }(u_{i+\mathrm{1,}j}^k+u_{i-\mathrm{1,}j}^k+u_{i\mathrm{,}j+1}^k+u_{i\mathrm{,}j-1}^k+d_{x\mathrm{,}i-\mathrm{1,}j}^k+d_{x\mathrm{,}i\mathrm{,}j}^k\\ +d_{y\mathrm{,}i\mathrm{,}j-1}^k+d_{y\mathrm{,}i\mathrm{,}j}^k+b_{x\mathrm{,}i-\mathrm{1,}j}^k+b_{x\mathrm{,}i\mathrm{,}j}^k+b_{y\mathrm{,}i\mathrm{,}j-1}^k+b_{y\mathrm{,}i\mathrm{,}j}^k)+\frac{\mu }{\mu +4\lambda }{f}_{i\mathrm{,}j}\mathrm{.}\end{array}$ (24)

Algorithm 3. The split Bregman algorithm of anisotropic TV denoising

Initialization: $k=0,u^0=0,b^0=0$ .

While $\Vert u^k-u^{k-1}\Vert >tol$ do,

$u^{k+1}=G^k$ , where G is the Gauss-Seidel function,

$d_x^{k+1}=shrink\left({\nabla }_x{u}^{k+1}+b_x^k,\frac{1}{\lambda }\right)$ ,

$d_y^{k+1}=shrink\left({\nabla }_y{u}^{k+1}+b_y^k,\frac{1}{\lambda }\right)$ ,

$b_x^{k+1}=b_x^k\left({\nabla }_x{u}^{k+1}-d_x^{k+1}\right)$ ,

$b_y^{k+1}=b_y^k\left({\nabla }_y{u}^{k+1}-d_y^{k+1}\right)$ ,

$k=k+1$ ,

End while.

2.5.2. Isotropic TV Denoising Problem

The problem of isotropic TV denoising is considered in [6] ,

$\left({{\mathcal{P}}^{\prime }}_1\right)I{s}_1\mathrm{<mo>\; :\; </mo>}=\underset{u}{\min }\left\{{\Vert \nabla u\Vert }_2+\frac{\mu }{2}{\Vert u-f\Vert }_2^2\right\}\mathrm{.}$ (25)

The problem $\left({{\mathcal{P}}^{\prime }}_1\right)$ is solved using a constraint equivalent problem $\left({{\mathcal{P}}^{\prime }}_2\right)$ :

$\left({{\mathcal{P}}^{\prime }}_2\right)\{\begin{array}{l}\widetilde{I{s}_1}\mathrm{<mo>\; :\; </mo>}=\underset{u}{\min }{\Vert \left(d_x\mathrm{,}{d}_y\right)\Vert }_2+\frac{\mu }{2}{\Vert u-f\Vert }_2^2\\ \text{subjectto}{d}_x=u_x,d_y=u_y\mathrm{.}\end{array}$ (26)

To solve the problem $\left({{\mathcal{P}}^{\prime }}_2\right)$ , we solve the following problem without constraint:

$\left({{\mathcal{P}}^{\prime }}_3\right)\widetilde{I{s}_3}\mathrm{<mo>\; :\; </mo>}=\underset{u\mathrm{,}{d}_x\mathrm{,}{d}_y}{\min }\left\{{\Vert \left(d_x\mathrm{,}{d}_y\right)\Vert }_2+\frac{\mu }{2}{\Vert u-f\Vert }_2^2+\frac{\lambda }{2}{\Vert d_x-u_x\Vert }_2^2+\frac{\lambda }{2}{\Vert d_y-u_y\Vert }_2^2\right\}\mathrm{.}$ (27)

The split Bregman algorithm can be used to tackle this last difficulty.

We give the following definitions:

$s^k\mathrm{<mo>\; :\; </mo>}=\sqrt{{\left|u_x^k-b_x^k\right|}^2+{\left|u_y^k-b_y^k\right|}^2}\mathrm{.}$ (28)

Algorithm 4. The split Bregman algorithm of isotropic TV denoising

Initialization: $k=0,u^0=0,b^0=0$ .

While $\Vert u^k-u^{k+1}\Vert >tol$ do,

$u^{k+1}=G^k$ , where G is the Gauss-Seidel function,

$d_x^{k+1}=\frac{s^k\lambda \left(u_x^k+b_x^k\right)}{s^k\lambda +1}$ ,

$d_y^{k+1}=\frac{s^k\lambda \left(u_y^k+b_y^k\right)}{s^k\lambda +1}$ ,

$b_x^{k+1}=b_x^k+\left(u_x^{k+1}-d_x^{k+1}\right)$ ,

$b_y^{k+1}=b_y^k+\left(u_y^{k+1}-b_y^{k+1}\right)$ ,

$k=k+1$ ,

End while.

3. Numerical Experimental Results

We present in this section our numerical results obtained with the following models of: Tikhonov regularization, ROF, anisotropic and isotropic TV denoising. Let X be the matrices that depict an image of size $m\times n$ . We then used Matlab $f=imnoise\left(X\mathrm{,}\text{'}gaussian\text{'}\mathrm{,}sigma\right)$ command to define our noise image f, where sigma is a version of the Gaussian noise level. We used the values $\mu =0.1$ , $\lambda =0.2$ and the tolerance $Tol={10}^{-5}$ in our studies. In our experience, we have tried to implement several models of rehabilitation. Each model aims to produce a better solution to remove noise from the image. However, we are going to implement iast in the script. By calculating Performance metrics as well different sigma values, we try to present the best result. The results for Tikhonov regularization, ROF, anisotropic and isotropic TV denoising algorithms are in Tables 1-8.

In addition, in Tables 9-16, we evaluate quality of images restored by the image restoration models of Tikhonov regularization, ROF, anisotropic and isotropic TV denoising algorithms, we use square error (MSE), signal noise rate (SNR), peak signal to noise ratio (PSNR), image quality index (IQI), normalized cross-correlation (NK), average difference (AD), structural content (SC), maximum difference (MD), and normalized absolute error (NAE), So that the definitions of quality measures are in the Figure 1.

In Figure 2, we did an experiment by taking the original image of Barbara (image without noise), then we added white Gaussian noise (sigma 0.08).

A comparative numerical was carried out between the Tikhonov regularization restoration model and the ROF; TV anistropic and isotropic denoising algorithm for the same parameter sigma 0.08 is shown in Figure 3 and Figure 4.

Table 1. Results for the ROF algorithm, sigma = 0.01.

Table 2. Results for the ROF algorithm, sigma = 0.08.

Table 3. Results for the ROF algorithm, sigma = 0.2.

Table 4. Results for the Tichonov Regularization algorithm, sigma = 0.01.

Table 5. Results for the Tichonov Regularization algorithm, sigma = 0.08.

Table 6. Results for the Tichonov Regularization algorithm, sigma = 0.2.

Table 7. Results for the anisotropic TV denoising algorithm, sigma = 0.08.

Table 8. Results for the isotropic TV denoising algorithm, sigma = 0.08.

Table 9. Performance metrics for ROF algorithm, sigma = 0.01.

Table 10. Performance metrics for ROF algorithm, sigma = 0.08.

Table 11. Performance metrics for ROF algorithm, sigma = 0.2.

Table 12. Performance metrics for Tichonov algorithm, sigma = 0.01.

Table 13. Performance metrics for Tichonov algorithm, sigma = 0.08.

Table 14. Performance metrics for Tichonov algorithm, sigma = 0.2.

Table 15. Performance metrics for the anisotropic TV denoising algorithm, sigma = 0.08.

Table 16. Performance metrics for the isotropic TV denoising algorithm, sigma = 0.08.

Figure 1. Quality measures.

In Figure 5, we did an experiment by taking the original image of girl (image without noise), then we added white Gaussian noise (sigma 0.08).

Finally, in Figure 6 and Figure 7, we show comparisons and numerical results for the Tikhonov regularization restoration model and the ROF model; TV anistropic and isotropic denoising algorithm with a noisy image of girl for the same parameter sigma 0.08.

Figure 2. The original and noisy image Barbara for sigma = 0.08.

Figure 3. Denoised image barbar by Tikhonov and ROF for sigma = 0.08.

Figure 4. Denoised image Barbara by TV anistropic and isotropic for sigma = 0.08.

Figure 5. The original and noisy image girl for sigma = 0.08.

Figure 6. Denoised image girl by Tikhonov and ROF for sigma = 0.08.

Figure 7. Denoised image girl by TV anistropic and isotropic for sigma = 0.08.

Remark

To quantify the restoration quality for a noisy image, we use sometimes measures. On note $x_{j\mathrm{,}k}$ is original image and ${x^{\prime }}_{j\mathrm{,}k}$ is restored image with [M N] this is the size of the images.

Program for calculating image quality measurements in MATLAB

4. Conclusion

In this paper, we have presented and compared theoretical and numerical of different imaging algorithms for solving optimization problems. We are looking for an image that is near to the original as possible among images that have been skewed by Gaussian and additive noise. Image deconstruction is a technique for restoring a noisy image after it has been captured. According to our experimentation, and by calculating performance metrics as well different sigma values, we can conclude that the ROF model is better image quality compared to the Tichonov regularization, because the space BV ensures continuity and allows for a stairway effect in restoring smooth images in applications where edges are not the main feature. We can conclude that the anisotropic TV and isotropic TV denoising algorithms work in a direct correlation relationship. In other words, regardless of how little the sigma value is that we get better and more old image quality results. Finally, it should be mentioned that all the methods that are common for removing parasitic information are added from an image.

Acknowledgements

We would like to thank the referees for some corrections which greatly improved the presentation of this paper.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

References

[1]	Bergounioux, M. (2008-2009) Quelques méthodes mathématiques pour le traite- ment d’image. https://cursus.edu/fr/18910/quelques-methodes-mathematiques-pour-le-traitement-dimage
[2]	Rudin, L., Osher, S. and Fatemi, E. (1992) Total Variation Based Noise Removal Algorithms. Physica D: Nonlinear Phenomena, 60, 259-268. https://doi.org/10.1016/0167-2789(92)90242-F
[3]	Chambolle, A. (2004) An Algorithm for Total Variation Minimization and Applications. Journal of Mathematical Imaging and Vision, 20, 89-97. https://doi.org/10.1023/B:JMIV.0000011321.19549.88
[4]	Meyer, Y. (March 2001) Oscillating Patterns in Image Processing and in Some Nonlinear Evolution Equations: The Fifteenth Dean Jacquelines B. Lewis Memorial Lectures. American Mathematical Society, Boston, USA.
[5]	Goldstein, T. and Osher, S. (2009) The Split Bregman Method for -Regularized Problems. SIAM Journal on Imaging Sciences, 2, 323-343. https://doi.org/10.1137/080725891
[6]	Bush, J. (June 10, 2011) Bregman Algorithms. Senior Thesis, University of California, Santa Barbara, Santa Barbara.