An overview of image restoration based on variational regularization

Abstract

Image restoration is a complicated process in which the original information can be recovered from the degraded image model caused by lots of factors. Mathematically, image restoration problems are ill-posed inverse problems. In this paper image restoration models and algorithms based on variational regularization are surveyed. First, we review and analyze the typical models for denoising, deblurring and inpainting. Second, we construct a unified restoration model based on variational regularization and summarize the typical numerical methods for the model. At last, we point out eight diffcult problems which remain open in this field.

Keywords

Regularization / image restoration / inverse problem / total variation / wavelet

1 Introduction

Images are the main medium for humans to acquire and exchange information. In the rapid development of computer technology and the increasingly popularity of digital products today, images are more inseparable with people's production and life. It has become the main tool for computer science, psychology, physiology, and many other disciplines to study visual perception. However, in the process of image acquisition, storage, transmission, and display, due to the influence of various factors, the image quality is often degraded. Usually, the image is contaminated by noise, blurred, reduced resolution, and even some parts of the image are missing. But we know that many areas of production and living need images with high quality. Therefore, eliminating noise in the image, removing blur, and repairing the image are processing very basic and important significance. These image processing techniques are the research content of image restoration.

Image restoration is an image processing technology that uses the prior information of image degradation to establish an image degradation model and relatively accurately reverse the original image information [8]. It is an indispensable complex mathematical process for people to obtain image information and understand images, and also the basis for further processing of images. Therefore, image restoration has always been one of the most basic and important research topics in the field of applied mathematics and image processing, which has important theoretical value and practical significance.

The application of image restoration technology in engineering first appeared in the aerospace field [8]. As early as in the 1950s, humans began to explore space, mainly by launching spacecraft into space and observing them with telescopes. These observations are usually some images. For example, in 1964, the United States spent nearly 10 million dollars to launch the Mariner 4 spacecraft to explore Mars and finally produce 22 images. When acquiring images in space, the motion speed of the spacecraft is much faster than the exposure speed of the camera. In addition, any imaging system is inevitably disturbed by noise. These factors cause the decline in the quality of the image, reduce the scientific value of the image, and cause huge economic losses. It is hoped that there is a technology can recover useful information from reduced images and make up for the loss of imperfect image acquisition systems. Luckily, image restoration technology arises with the need of application. So far, aerospace image is still one of the main application fields of image restoration technology. In the field of medicine, image restoration technology is widely used in X-ray, CT imaging system, and it is used to suppress the noise of medical images and improve the resolution of medical images [109]. In the field of remote sensing [136], image restoration technology can overcome the optical remote sensing of atmospheric disturbance image quality, the imaging system relative ground movement caused by fuzzy, suppress noise and thin cloud, and coherent noise in remote imaging, ghost and tail dragging. In the field of military security, image restoration technology is inseparable from terrain recognition of cruise missiles, face recognition in surveillance video, seal identification, and expired text recognition. Image restoration techniques also have important applications in visual analysis and the film industry [7]. Image restoration techniques also have important applications in visual analysis and the film industry. For example, digital restoration of valuable old paintings and calligraphy, restoration of old photographs and negatives with scratches and defects, removal of occlusion from computer images, removal of text fonts from images for special purposes, etc.

In conclusion, in today's digital age, image restoration has penetrated into many areas of human production and life. Therefore, it is important to study the models and algorithms of image restoration to make it better and faster to serve people's production and life.

Image recovery belongs to the scope of inverse problems, and its mathematical foundation covers generalized functional analysis, variational methods, partial differential equations, function space theory, wavelet analysis, probability statistics, differential geometry, optimization theory, machine learning theory, numerical analysis, etc. The research on image restoration models and algorithms has promoted the development, intersection, and integration of these disciplines those have led to new methods, such as the level set method, the Bregman splitting algorithm, and the graph tangent method. These new methods play an important role in many disciplines such as theory and computation of inverse problems [64, 122].

So far, variational regularized image restoration models and algorithms have become an important research direction in applied mathematics. In order to better the reader's work in this area and to follow the frontiers as soon as possible, this paper reviews the important models and algorithms in the existing literature. Section 2 presents some important mathematical models and theoretical results in image restoration and gives a proper review. Section 3 establishes a unified variational regularized image restoration model and summarizes various typical numerical solution methods. The last section presents some prospects for image restoration research.

2 Models and analysis

Historically, the first problems studied in image restoration were denoising and deblurring [8]. In the last two decades, a large body of literature has studied these two types of problems [33], and the resulting methods such as full-variance regularization have guided techniques such as image restoration and image reconstruction (e.g., compressed perception, photoacoustic laminar imaging [108, 113]). Some typical denoising (decomposition), deblurring, and restoration models are discussed below.

2.1 Decomposition Model

Grayscale noise and sampling noise are two common types of noise that cause image degradation. They are random and are usually assumed to obey some probability distribution [109]. A great a posteriori estimation of the numerical characteristics of these distributions as a priori information can lead to a denoising model based on random methods [58, 122].

• Grayscale Noise Denoising Model

According to the mechanism of grayscale noise formation, it can be divided into additive noise, multiplicative noise, and image-dependent noise. Additive grayscale noise can describe the noise generated by amplifiers, etc. Suppose u^{\delta }=u+\delta, where u^{\delta } is the image contaminated by Gauss white noise \delta with zero expectation and variance of {\sigma }^2. By solving the function

(2.1)\underset{u\in {\mathbb{R}}^{N\times N}}{min}\sum _{1⩽i,j⩽N}((u_{ij}-u_{ij}^{\delta }{)}^2+\lambda |\mathrm{\nabla }{u}_{ij}{|}^2),

the original image u is obtained, where \lambda is a constant proportional to the noise standard deviation \sigma.

Regarding multiplicative grayscale noise, such as optical quantum noise, film grain noise, it is usually processed into additive noise by taking logarithm [93, 107]. Image-dependent noise, such as noise obeying Poisson distribution or pretzel noise generated by CCD sensor, is usually treated by median filtering and other methods [9, 28, 111].

• Sampling Noise Denoising Model

Suppose u^{\delta }=u+\delta \cdot |\mathrm{\nabla }u|, where u^{\delta } is the image with noise, and \delta is Gauss white noise with zero expectation and variance {\sigma }^2. The gradient of the original image u obeys the Laplace distribution, and can be solved by solving

(2.2)\underset{u\in {\mathbb{R}}^{N\times N}}{min}\sum _{1⩽i,j⩽N}(\frac{1}{2}\frac{(u_{ij}-u_{ij}^{\delta }{)}^2}{|\mathrm{\nabla }{u}_{ij}{|}^2}+\frac{\lambda }{2}|\mathrm{\nabla }{u}_{ij}|).

We obtain u, where \lambda is a constant proportional to the noise standard deviation \sigma.

The denoising models based on the probability distribution of noise as the a priori information and the great a posteriori estimation have solved the denoising problem to a certain extent, but the numerical experiments show that these two models also have great shortcomings. The sharp edges of the images are polished. The edge of the image is the most important part of the human visual system to perceive the image [68, 72, 85]. In order to reveal the advantages and disadvantages of (2.1), (2.2) mathematically, they can be expressed in the corresponding continuous form as

(2.3)\underset{u\in X}{min}\left(\frac{1}{p}{\int }_{\mathrm{\Omega }}\frac{(u-u^{\delta }{)}^p}{|\mathrm{\nabla }u{|}^q}+\frac{\lambda }{r}|\mathrm{\nabla }u{|}^r\right).

Here (p,q,r)=(2,0,2) or (2,1,1), corresponding to models (2.1) and (2.2), respectively, where \lambda >0,\mathrm{\Omega }\subset {\mathbb{R}}^2 is a bounded connected open set with Lipschitz boundary \mathrm{\partial }\mathrm{\Omega }. Model (2.3) is called the variational model, and the solution of this problem is often transformed into solving the corresponding Euler-Lagrange equation, i.e., the corresponding partial differential equation (PDE) model, by variationalization, which are equivalent. For more discussions on statistical methods, variational methods, PDE methods, and the relationship between wavelet methods are discussed in [9, 17, 94].

It can be seen that model (2.3) uses the Sobolev space [139] W^{1,2}(\mathrm{\Omega }) with semi-parametric \Vert \mathrm{\nabla }u{\Vert }_{L^2(\mathrm{\Omega })} as the regularizing term. The function in this space has some smoothness, which can filter out the noise in the form of oscillations, but cannot characterize the jumps of the function. So it cannot protect the image. Therefore, it cannot protect the edges of the image.

• ROF Model

In 1992, Rudin et al. [106] introduced the full variational semi-parametrization as a regularization term into the image decomposition model from the perspective of energy minimization, and obtained the well-known ROF Model. This model overcomes the defects of the previous denoising models which cannot protect the edges, and has been very successful. Since then, the full-variance regularization method has been widely used in many problems of image processing, and the research related to it is still in progress [108].

Definition 2.1 [52]　Set u\in L^1(\mathrm{\Omega }). u is a bounded variational function defined on \mathrm{\Omega }, if its weak derivative can be expressed as a finite Radon measure on \mathrm{\Omega }, i.e., if for a certain constant value of Radon on \mathrm{\Omega } measure Du=(D_{1}u,D_{2}u), there is

(2.4){\int }_{\mathrm{\Omega }}u\mathrm{\nabla }\varphi =-{\int }_{\mathrm{\Omega }}\varphi Du,\mathrm{\forall }\varphi \in C_c^1(\mathrm{\Omega }).

Here, the full variance of u is noted as {\int }_{\mathrm{\Omega }}|Du|.

By the Riesz representation theorem, the continuous linear generalized function on C_c(\mathrm{\Omega }) can be expressed as a Radon measure on \mathrm{\Omega } [130]. Therefore, the full variance can also be equivalently defined as

(2.5){\int }_{\mathrm{\Omega }}|Du|=sup\left\{{\int }_{\mathrm{\Omega }}u\mathrm{\nabla }\cdot g:\dot{g}=(g_1,g_2)\in (C_c^1(\mathrm{\Omega }){)}^2,|g|=\sqrt{g_1^2+g_2^2}⩽1\right\}.

The full variance finite function of all functions on \mathrm{\Omega } by parametric \Vert u{\Vert }_{BV(\mathrm{\Omega })}=\Vert u{\Vert }_{L^1(\mathrm{\Omega })}+{\int }_{\mathrm{\Omega }}|Du| forms the Banach space BV(\mathrm{\Omega }). It is also customary to denote the full variance of u over \mathrm{\Omega }{\int }_{\mathrm{\Omega }}|Du| as |u{|}_{BV(\mathrm{\Omega })}.

ROF Model is

(2.6)\underset{u\in BV(\mathrm{\Omega })}{inf}\left(\lambda \Vert u^{\delta }-u{\Vert }_{L^2(\mathrm{\Omega })}+|u{|}_{BV(\mathrm{\Omega })}\right),

when u has smoothness of W^{1,1}(\mathrm{\Omega }), the full variational semi-parametrization |u{|}_{BV(\mathrm{\Omega })} of u is the semi-parametrization {\int }_{\mathrm{\Omega }}|\mathrm{\nabla }u| of W^{1,1}(\mathrm{\Omega }). But |u{|}_{BV(\mathrm{\Omega })} can allow jumps of the function [2], which better characterizes the edges of the image. Therefore, using |u{|}_{BV(\mathrm{\Omega })} as the regularization term is the fundamental reason for the great success of the ROF Model.

• TV-G Model

In 2001, Meyer [88] pointed out the ROF Model (2.6) uses \Vert \cdot {\Vert }_{L^2(\mathrm{\Omega })} as the fidelity term. The detailed oscillatory component texture in the image is lost while denoising. To this end, he proposes a Banach space G({\mathbb{R}}^2) that can characterize the oscillatory components of the image, and pointed out that the G-parameter of the oscillatory components such as texture is very small. Therefore, when the regularization parameter \lambda in the ROF Model (2.6) obtained is too small, there will be excessive loss of texture by denoising. In view of this, he proposed the well-known TV-G Model, decomposes the image into parts with geometric structure and oscillatory components (texture or noise).

Definition 2.2　Let G=\mathrm{\nabla }\cdot (L^{\mathrm{\infty }}(\mathrm{\Omega },{\mathbb{R}}^2)),v\in G. Then the G-parametrization of v is

(2.7)\begin{array}{l}\Vert v{\Vert }_{G(R^2)}=inf\left\{|g{|}_{L^{\mathrm{\infty }}({\mathbb{R}}^2)}:v=\mathrm{\nabla }\cdot g,g=(g_1,g_2)\right\}.\end{array}

Here, (G,|\cdot {\Vert }_G) forms the Banach space.

Theorem 2.1 [88]　 Let v\in G\cap L^2(\mathrm{\Omega }),u\in BV(\mathrm{\Omega }), Tr(u)=0. Then

(2.8)⟨u,v⟩={\int }_{\mathrm{\Omega }}uv⩽\Vert v{\Vert }_G|u{|}_{BV(\mathrm{\Omega })}.

It follows from the theorem that (G,\Vert \cdot {\Vert }_G) can be regarded as the dual of BV(\mathrm{\Omega }) in some sense [52]. In fact, (G,\Vert \cdot {\Vert }_G) is isometric isomorphic to the dual W^{-1,\mathrm{\infty }}(\mathrm{\Omega }) of W_0^{1,1}(\mathrm{\Omega }), and therefore larger than the dual space of BV(\mathrm{\Omega }). From this theorem, it also follows that when \Vert u^{\delta }{\Vert }_G⩽\frac{1}{2\lambda }. The solution of the ROF Model (2.6) is u=0, which mathematically explains the loss of texture.

Meyer replaces the L^2 parametrization in the ROF Model with the G parametrization as the approximation measure to obtain the TV-G Model

(2.9)\underset{(u,v)\in BV\times G:u^{\delta }=u+v}{inf}\left(|u{|}_{BV({\mathbb{R}}^2)}+\lambda |v{\Vert }_G\right).

Since TV-G Model (2.9) is defined over the whole {\mathbb{R}}^2, it is difficult to solve numerically, and many scholars have proposed to approximate or replace it with the corresponding model. The more successful ones are Osher et al. [99, 121] and Aujol et al. [3, 4].

Aujol followed Meyer's method to extend the G({\mathbb{R}}^2) to G(\mathrm{\Omega }), giving a deep and practical inscription of G(\mathrm{\Omega }) as

(2.10)G(\mathrm{\Omega })=\left\{v\in L^2(\mathrm{\Omega }):{\int }_{\mathrm{\Omega }}v=0\right\},

and the model is proposed as follows:

(2.11)\underset{(u,v)\in BV(\mathrm{\Omega })\times G(\mathrm{\Omega }),\Vert v{\Vert }_{G(\mathrm{\Omega })}⩽\mu }{inf}\left(\Vert u^{\delta }-u-v{\Vert }_{L^2(\mathrm{\Omega })}+\lambda |u{|}_{BV(\mathrm{\Omega })}\right).

Here, it is understood that the image u^{\delta } is decomposed into the geometric structure u, the texture v, and the noise u^{\delta }-u-v. When discussing the denoising problem, taking u+v as the desired result, it serves to protect the texture. Aujol et al. also discussed that the solution of the discrete form of the model converges to the solution of the discrete form corresponding to Meyer's TV-G Model (2.9) under certain conditions.

• O-S-V Model

Osher et al. [99] then approximated model (2.9) from the perspective of \mathrm{\nabla }\cdot L^p\to \mathrm{\nabla }\cdot L^{\mathrm{\infty }}=G(p\to \mathrm{\infty }), obtaining the model

(2.12)\underset{(u,g)\in BV(\mathrm{\Omega })\times (L^{\mathrm{\infty }}(\mathrm{\Omega }){)}^2}{inf}(|u{|}_{BV(\mathrm{\Omega })}+\Vert u^{\delta }-u-\mathrm{\nabla }\cdot g{\Vert }_{L^2(\mathrm{\Omega })}^2+\lambda \Vert g|{\Vert }_{L^p(\mathrm{\Omega })}).

Here, u^{\delta } is decomposed into the sum of three components u+v+w where v=\mathrm{\nabla }\cdot g,w=u_0-u-v. The first term in the model is the regularization term, which ensures that u\in BV(\mathrm{\Omega }). The second term is the approximation term, which ensures that u^{\delta }\approx u+v. The third term is the term for v under the penalty of approximating G parametrization. When \mu \to \mathrm{\infty },p\to \mathrm{\infty }, model (2.12) converges to Meyer's model (2.9). The paper [99] also shows that the numerical results are insensitive to the value of 1⩽p⩽10. In particular, taking p=2 yields the model

(2.13)\underset{(u,v)\in BV(\mathrm{\Omega })\times H^{-1}(\mathrm{\Omega })}{inf}(|u{|}_{BV(\mathrm{\Omega })}+\mu \Vert u^{\delta }-u-v{\Vert }_{L^2(\mathrm{\Omega })}^2+\lambda \Vert v{\Vert }_{H^{-1}(\mathrm{\Omega })}).

This model simplifies model (2.12), but the model solution needs to be transformed into solving a fourth-order nonlinear partial differential equation, and the numerical implementation is slow. For this reason, Daubechies et al. [45] proposed to replace the space BV(\mathrm{\Omega }) in model (2.13) with the Besov space B_1^1(L^1(\mathrm{\Omega })) from the viewpoint of improving computational efficiency, obtaining

(2.14)\underset{(u,v)\in B_1^1(L^1(\mathrm{\Omega }))\times H^{-1}(\mathrm{\Omega })}{inf}(\Vert u^{\delta }-K(u+v){\Vert }_{L^2(\mathrm{\Omega })}^2+\mu |u{|}_{B_1^1(L^1(\mathrm{\Omega }))}+\lambda \Vert v{\Vert }_{H^{-1}(\mathrm{\Omega })}^2).

Here, when K is the convolution operator, model (2.14) is the deblurring model; when K is the constant operator, model (2.14) is the denoising model.

The reason this works so well is that [46]

B_1^1(L^1(\mathrm{\Omega }))\subset BV(\mathrm{\Omega })\subset (B_1^1(L^1(\mathrm{\Omega })),wk).

At this point, since the spaces B_1^1(L^1(\mathrm{\Omega })), H^{-1}(\mathrm{\Omega }), L^2(\mathrm{\Omega }) have the corresponding wavelet inscriptions [118], so model (2.14) can be transformed to the wavelet domain for solution, and and the objective function in the wavelet domain is differentiable, which can be transformed into a one-dimensional function of the minimax problem. It is easy to solve and improve the speed of operation. Regarding the inscription of the function space, a series of results are given in [79–81].

Regarding Meyer TV-G model (2.9), further theoretical analysis and applied research are also available in [57, 60, 77]

• TV-L¹ Model

In 2005, Chan et al. [6] proposed the TV-L^1 Model [6] under the influence of [1, 93],

(2.15)\underset{u\in BV(\mathrm{\Omega })}{inf}(\Vert u^{\delta }-u{\Vert }_{L^1(\mathrm{\Omega })}^2+\lambda |u{|}_{BV(\mathrm{\Omega })}).

This model changes the L^2 parametrization in ROF Model (2.6) to L^1 parametrization as the approximation metric, which can remove the pretzel noise [28]. Paper [6] also gives a fast numerical method for solving model (2.15). In [29], the mathematical properties of the model are investigated in depth.

Although these above models can solve the image denoising and image decomposition problems well, due to the oscillatory nature of both texture and noise are oscillatory in nature. It is a challenging problem to distinguish texture from noise, especially for texture-based image denoising. Peyré et al. [86] proposed an adaptive model for anisotropic locally parallel textures, using a weighted l^2 parametrization of the windowed Fourier coefficients to measure texture components. The weights can be generated by solving a variational molecular problem iteratively based on the local orientation and frequency of the texture, which reflects the geometric structure information of the local parallel texture. In [59], the concept of local variance measure is proposed, and pointed out that the local variance of noise is small compared to other components. In view of this, wavelet inscriptions in Hardy space H^1 and binary BMO space and their relationship with the local variance of images are investigated, and a denoising method that can distinguish texture from noise is given in [70, 133, 134]. Both theoretical analysis and numerical experiments show the effectiveness of the described methods.

These Denoising (decomposition) Models have been studied throughout the development of image restoration theory and techniques in the last 20 years. The ideas and methods have contributed to the development of image restoration and other disciplines and fields.

2.2 Deblurring Model

There are many causes of image blurring, which can be usually classified as linear and nonlinear [24]. Motion blur is a common blur, which is caused by the relative motion of the camera and the shooting target in the image smoothing, is a linear process, mathematically can be described as

u^{\delta }=Ku.

Here, u^{\delta } is the blurred image, u is the image to be recovered, K is the operator with translation invariance, usually taken as the convolution operator. Since noise is unavoidable in practice, we consider the image with zero expectation and variance of {\sigma }^2 with additive Gauss white noise \delta interference.

It can be described by u^{\delta }=Ku+\delta. It is well known that the convolution operator is a bounded tight operator and its inverse is unbounded, which leads to the discomfort of solving the model [90]. Therefore, the regularization [117] is the preferred method for deblurring, and in order to protect the edges of the image, {\int }_{\mathrm{\Omega }}|Du| is chosen as the regularization term and then minimized to obtain the TV deblurring model [34].

(2.16)\underset{u\in BV(\mathrm{\Omega })}{inf}(\Vert Ku-u^{\delta }{\Vert }_{L^2(\mathrm{\Omega })}^2+\lambda |u{|}_{BV(\mathrm{\Omega })}),

from which u can be solved.

If K is known, replacing u in the denoising model discussed earlier with Ku, then (2.16) is the corresponding common deblurring model. For example, in [44], we propose to use Besov space B_1^1(L^1(\mathrm{\Omega })) instead of BV(\mathrm{\Omega }) in model (2.16) in BV(\mathrm{\Omega }), to obtain the deblurring model based on the iterative wavelet shrinkage algorithm

(2.17)\underset{u\in B_1^1(L^1(\mathrm{\Omega }))}{inf}(\Vert Ku-u^{\delta }{\Vert }_{L^2(\mathrm{\Omega })}^2+\lambda |u{|}_{B_1^1(L^1(\mathrm{\Omega }))}).

References [18, 19, 35, 53, 102, 136] discussed the use of spatial domains in the K unknown case using the spatial domain iterative method, the cubic correlation method, the simulated annealing method, least squares method, and the fact that the convolution operator K can be sparsely represented by a wavelet framework and thus solved in alternating directions, etc. The numerical simulations show the effectiveness of these methods. However, these methods do not completely solve the nonlinearities of the problem due to the unknown K. References [24, 103] discussed nonlinear defuzzification problems based on statistical models and pointed out that building more reasonable nonlinear models is an effective way to study these problems.

2.3 Restoration Model

The term digital image restoration (inpainting) was first introduced in 2000 at an academic conference [11], which is an important element in image restoration research. It is a process of filling in the information-deficient areas of an image so that the observer does not notice that the image was once deficient or has been restored. Nowadays, scientists and engineers are working on automatic restoration techniques, in which the user selects the restored area and the computer does the rest of the work automatically. Research on image restoration models can be divided into two broad (and major) categories.

• Non-texture image restoration

Let the image u^{\delta } be defined on \mathrm{\Omega }, D\subseteq \mathrm{\Omega } denote the defective part of the image, and \mathrm{\Omega }\setminus D denote the undeficient region of the image. Earlier image restoration techniques were based on PDE. For example, Sapiro et al. in [11] proposed a method to use the edge \mathrm{\partial }D of the region to be repaired, to estimate the direction of iso-illumination from coarse to fine, and uses a diffusion mechanism to propagate the information to the region to be repaired. Oliveira et al. [95] proposed a fast patching technique in \mathrm{\partial }D. They propose to determine the direction of the iso-illumination on \mathrm{\partial }D and connect the corresponding iso-illumination with a straight line. The information adjacent to the region to be repaired is diffused and filled in the area of the iso-illumination line. The disadvantage of these methods is that they are relatively complicated to compute. Nowadays, the most studied model is the variational image restoration model. This method is based on the assumption of some a priori information of the original image u, such as u which function space it belongs to, combined with the known information of u^{\delta }{|}_{\mathrm{\Omega }\setminus D} (which in practice is the information disturbed by noise and ambiguity), restores the image by minimizing the corresponding energy generalization. For example, Chan et al. [32] proposed a TV restoration model. Considering the noise and blur interference, u^{\delta } on the region \mathrm{\Omega }\setminus D can be described as

(2.18)u^{\delta }{|}_{\mathrm{\Omega }\setminus D}=(Ku+\delta ){|}_{\mathrm{\Omega }\setminus D}.

Here, \delta is the additive Gauss white noise with zero mean and {\sigma }^2 variance, and K is the convolution operator that causes motion blur. Considering the edge protection, given the a priori assumption of u\in BV(\mathrm{\Omega }) by minimizing

(2.19)\begin{array}{l}\underset{u\in BV(\mathrm{\Omega })}{inf}(\Vert Ku-u^{\delta }{\Vert }_{L^2(\mathrm{\Omega }\setminus D)}^2+\lambda |u{|}_{BV(\mathrm{\Omega })}),\end{array}

u can be fixed. Here, the regularization parameter \lambda is proportional to the noise variance prior information. For the convenience of solving, model (2.19) is rewritten as

(2.20)\begin{array}{l}\underset{u\in BV(\mathrm{\Omega })}{inf}(\Vert (Ku-u^{\delta }){\chi }_{\mathrm{\Omega }\setminus D}{\Vert }_{L^2(\mathrm{\Omega })}^2+\lambda |u{|}_{BV(\mathrm{\Omega })}).\end{array}

Numerical simulations show that model (2.20) does not guarantee the connectivity principle in vision theory and produces a step effect. Chan et al. [31] proposed a CDD repair model to overcome this drawback. The basic idea is to control the diffusion strength \frac{1}{|\mathrm{\nabla }u|} of the Euler-Lagrange equation corresponding to (2.20), so that the diffusion strength depends on the geometric information of the iso-illuminated lines, such as curvature. So \frac{1}{|\mathrm{\nabla }u|} is replaced by \frac{g|\kappa |}{|\mathrm{\nabla }u|}, and the CDD repair model as

(2.21)\{\begin{array}{l}-\mathrm{\nabla }\cdot \left({\displaystyle \frac{\mathrm{\nabla }ug(|\kappa |)}{|\mathrm{\nabla }u|}}\right)+\lambda K^{\ast }(Ku-u^{\delta }){\chi }_{\mathrm{\Omega }\setminus D}=0,\\ {{\displaystyle \frac{\mathrm{\partial }u}{\mathrm{\partial }\overrightarrow{n}}}|}_{\mathrm{\partial }\mathrm{\Omega }}=0.\end{array}

Here, \kappa =\mathrm{\nabla }\cdot (\frac{\mathrm{\nabla }u}{|\mathrm{\nabla }u|}) is the iso-illuminated line curvature, K is the constant operator in the original, \overrightarrow{n} is the outer normal vector of \mathrm{\partial }\mathrm{\Omega }, and the function g is defined as

This choice can make the diffusion stronger at the corresponding large curvature and weaker at the small curvature, which can compensate the deficiency of model (2.20). However, solving the restoration model (2.21) requires solving a third-order nonlinear PDE, and the numerical implementation is slow. Costanzino [41] builds the reconciliation model by making the a priori assumption that u\in W^{1,2}(\mathrm{\Omega }):

(2.22)\begin{array}{l}\underset{u\in BV(\mathrm{\Omega })}{inf}\left(\Vert (Ku-u^{\delta }){\chi }_{\mathrm{\Omega }\setminus D}{\Vert }_{L^2(\mathrm{\Omega })}^2+\lambda {\int }_{\mathrm{\Omega }}|\mathrm{\nabla }u{|}^2\right).\end{array}

This model can avoid the defect of step effect in model (2.20) and is faster than the CDD restoration model (2.21), but model (2.22) uses \Vert \cdot {\Vert }_{W^{1,2}(\mathrm{\Omega })} parametrization as the regularization term, which leads to blurred image edges, which is a fatal flaw. Based on this, Zhang et al. in [50] proposed the p-tuning restoration model:

(2.23)\underset{u\in BV(\mathrm{\Omega })}{inf}\left(\Vert (Ku-u^{\delta }){\chi }_{\mathrm{\Omega }\setminus D}{\Vert }_{L^2(\mathrm{\Omega })}^2+\frac{\lambda }{p}{\int }_{\mathrm{\Omega }}|\mathrm{\nabla }u{|}^p\right).

It overcomes the above shortcomings. Here, , K is the constant operator in the original text.

Other common restoration models are the Mumford-Shah restoration model [32, 119]. This model has low computational complexity, but it is not as effective as in the traditional applications of image segmentation and denoising. Esedoglua and Shen [50] obtained the Mumford-Shah-Euler model by introducing Euler's elastica, which improved the Mumford-Shah restoration model. Another class is the TV-Stokes model that regularizes the tangential vector field of the iso-illumination line by taking the incompressibility of the normal direction of the iso-illumination line as a constraint [105, 115].

Another important model for non-texture restoration is the model combining wavelet (frame) and variational in wavelet domain restoration [16, 30, 37, 131], which are here unified into a single model without considering the interference of fuzziness

(2.24)\underset{{\mathrm{\Psi }}^u}{inf}\left(\frac{1}{p}{\int }_{\mathrm{\Omega }}|\mathrm{\nabla }u{|}^p+\frac{\lambda }{2}\Vert {\chi }_I({\mathrm{\Psi }}^u-{\mathrm{\Psi }}^{u^{\delta }}){\Vert }_{l^2(\mathrm{\Gamma })}^2\right),

where \mathrm{\Gamma }=Z\times Z^2,u={\sum }_{j,k}⟨u,{\mathrm{\Psi }}_{j,k}⟩{\mathrm{\Psi }}_{j,k},(j,k)\in \mathrm{\Gamma },{\mathrm{\Psi }}^u=\{⟨u,{\mathrm{\Psi }}_{j,k}⟩{\}}_{(j,k)\in \mathrm{\Gamma }}, and \mathrm{\Psi } is the wavelet (frame) of the parent function [84, 87], I=\{(j,k)\in \mathrm{\Gamma }, the wavelet coefficients {\mathrm{\Psi }}_{j,k}^{u^{\delta }} are corrupted by the absence of D is corrupted}, 1⩽p⩽2.

• Texture image repair

For the restoration of richly textured images, a common approach is to perform structural texture decomposition according to the denoising (decomposition) model, repair the structural part according to the above non-texture restoration model, and use the existing texture synthesis algorithm for the textured part [91, 110, 127] to synthesize the texture; another common method is to synthesize the texture directly using a multi-scale geometric analysis tool [100, 114, 125]. Recently, Shen et al. proposed the following structure-texture repair model using frame coefficients for structure and overlapping local cosine transform coefficients for texture [16, 22]:

(2.25)\underset{u,v}{inf}({\lambda }_1\Vert Fu{\Vert }_{L^1}+{\lambda }_2\Vert Lv{\Vert }_{L^1}),\mathrm{s}.\mathrm{t}.P_{\mathrm{\Omega }\setminus D}(u+v)=P_{\mathrm{\Omega }\setminus D}{u}^{\delta }.

Here, F,L,P denote the frame operator, the overlapping local cosine transform and the projection operator, respectively. The model is simple, easy to implement numerically, and works well.

3 Numerical algorithm

First, we analyze and compare the features of the models described in the previous section and propose a unified model that can be used to handle various image restoration as well as image sparse approximation problems (continuous form) based on l^1 regularization as follows:

(3.1)\underset{u\in X}{min}\left(F_{p,q,r,\lambda }^{u^{\delta }}(K,T,u)=\frac{1}{p}{\int }_{{\mathrm{\Omega }}^{\prime }}\frac{(Ku-u^{\delta }{)}^p}{|Tu{|}^q}+\frac{\lambda }{r}{\int }_{\mathrm{\Omega }}|Tu{|}^r\right).

Here, p⩾1,r⩾1,q⩾0,\lambda ⩾0,{\mathrm{\Omega }}^{\prime }=\mathrm{\Omega } or \mathrm{\Omega }\setminus D.

When p,q,r,K,T take different values, different recovery models corresponding to those in the previous section and some models in image reconstruction can be obtained.

For example:

1) Take p=2, q=0, r=1, {\mathrm{\Omega }}^{\prime }=\mathrm{\Omega }, K is the constant operator, T is the gradient operator (frame operator), then it is the ROF Model (frame-let-based denoising model).

2) Take p=2, q=0, , {\mathrm{\Omega }}^{\prime }=\mathrm{\Omega }, K as the constant operator and T as the gradient operator, then it is a generalized ROF Model [112].

3) Take p=2, q=0, r=1, {\mathrm{\Omega }}^{\prime }=\mathrm{\Omega }\setminus D, K is the projection operator, T is the gradient operator (frame operator), then it is a TV (framelet-based) repair model.

4) Take p=2, q=0, r=1, {\mathrm{\Omega }}^{\prime }=\mathrm{\Omega }, K is the convolution operator, T is the gradient operator (frame operator), then it is a TV(framelet-based) defuzzification model.

5) If p=1, q=0, r=1, {\mathrm{\Omega }}^{\prime }=\mathrm{\Omega }, K is the constant operator, and T is the gradient operator, then it is a TV-L^1 model.

Moreover, if p=2,q=0,r=1,T is a constant operator and K is an operator with high compression ratio, then the corresponding discrete model is applied to the compression perception [23, 39, 47] problem in the Basis Pursuit model. This is a class of l^1-minimization models, and the corresponding operator K can be chosen according to practical needs, which is now widely used in compressive imaging [56, 116], MRI, CT imaging [82, 128], Biosensing [42], image decomposition and machine vision [38, 129], and multi-sensing networks [104, 123], and other fields.

For model (3.1), there are two main numerical solution methods. One is to optimize first and then discretize, i.e., to obtain (3.1) by the variational method first. This infinite-dimensional optimization problem satisfies the Euler-Lagrange equation [49], usually called the gradient flow equation (a linear or nonlinear PDE satisfying the chi-square Neumann boundary condition), is then solved by the finite difference [78], finite element [71], level set [97], and so on. Another approach is to first discretize (3.1), and then solving it numerically by a finite dimensional optimization technique. Since the condition numbers of the coefficient matrices of the discretized linear equations obtained by the former method are usually large, which leads to slow convergence, more and more scholars have been working on this method in recent years. Therefore, in recent years, more and more scholars have adopted the latter method. Here, we only describe this approach. Without loss of generality, the following model is solved as:

(3.2)min(F_{2,0,1,\lambda }^{u^{\delta }}(K,\mathrm{\nabla },u)),

where K takes the convolution operator (constant operator), i.e., the TV defuzzification model (ROF Model).

3.1 Chan-Golub-Mulet’s original-pairing method

The method was proposed by Chan et al. [26, 36]. The Euler-Lagrange equation corresponding to model (3.2) is (here and later understood as its corresponding discrete form)

(3.3)\{\begin{array}{l}\mathrm{\nabla }\cdot \left({\displaystyle \frac{\mathrm{\nabla }u}{|\mathrm{\nabla }u|}}\right)+\lambda K^{\ast }(Ku-u^{\delta })=0,\\ {{\displaystyle \frac{\mathrm{\partial }u}{\mathrm{\partial }\overrightarrow{n}}}|}_{\mathrm{\partial }\mathrm{\Omega }}=0.\end{array}

To overcome the singularity of the above equation, introduce the auxiliary variable

(3.4)w=\frac{\mathrm{\nabla }u}{|\mathrm{\nabla }u|}.

We get

(3.5)\{\begin{array}{l}|\mathrm{\nabla }u|w-\mathrm{\nabla }u=0,\\ -\mathrm{\nabla }\cdot w+\lambda K^{\ast }(Ku-f)=0,\\ {{\displaystyle \frac{\mathrm{\partial }u}{\mathrm{\partial }\overrightarrow{n}}}|}_{\mathrm{\partial }\mathrm{\Omega }}=0.\end{array}

Then using the Newton's method of inexact linear search with Armijo steps [12] to solve (3.5) numerically. This method is fast and has global convergence.

3.2 Chambolle’s pairwise method

Solving the second equation from (3.5)

(3.6)u=(\lambda K^{\ast }K{)}^{-1}(\lambda K^{\ast }{u}^{\delta }+\mathrm{\nabla }\cdot w).

Substituting into the first equation of (3.5), we get

(3.7)\mu w=H(w).

Here,

\mu =|H(w)|,H(w)=\mathrm{\nabla }((\lambda K^{\ast }K{)}^{-1}(\lambda K^{\ast }{u}^{\delta }+\mathrm{\nabla }\cdot w)).

From (3.5), an iterative format for immobile points can be devised [27]

(3.8)w^0=0,w^{n+1}=\frac{w^n+\tau H(w^n)}{w^n+\tau |H(w^n)|},

where \tau is the parameter that controls the step size, the iterations converge if and only if \tau ⩽\frac{1}{8}\Vert (K^{\ast }K{)}^{-1}\Vert. The algorithm is easy to implement. It is theoretically sound and has been widely used [4, 5, 13, 14].

3.3 Semi-smooth Newton method

The KKT condition for the Chambolle pairwise method is obtained from (3.4), (3.5):

(3.9)\{\begin{array}{l}\mu w=H(w),\\ \mu (|w{|}^2-1)=0,\\ \mu ⩾0,\\ |w|⩽1.\end{array}

Using the Fisher-Burmeister function

\phi (a,b)=\sqrt{a^2+b^2}-a-b,

combine the last three equations of (3.9) into one equation, we get

(3.10)\{\begin{array}{l}\mu w=H(w),\\ \phi (\mu ,1-|w{|}^2)=0.\end{array}

Ng et al. [92] note that the semismoothness of (3.10), using a semismooth Newton method with superlinear convergence rate [69] for the solution.

3.4 Original-pairwise tight constraint set method

The method proposed by Kunisch et al. [75] is an algorithm that deals with (3.9) in a more general framework algorithm for linearly tight and non-tight constraints. Subsequently, Krishnan et al. [73, 74] also extended it to the case of nonlinear tight and non-tight constraints.

3.5 Bregman iteration

Bregman iterations include primitive and split iterations were proposed by Osher et al. [63, 96] in 2005 and 2009. It can efficiently handle many constrained and unconstrained optimization problems, and is easy to implement numerically. The Bregman iterative algorithm and its applications are still a hot research topic, see [15, 20–22, 98, 135] and its cited literature.

Let J(u) be a convex generic function. The Bregman distance on J(u) is defined as

(3.11)B_J^P(u,v)=J(u)-J(v)-⟨u-v,p⟩,p\in \mathrm{\partial }J(v).

Taking J(u)=\lambda {\int }_{\mathrm{\Omega }}|Du|, the Bregman primal iteration algorithm of (3.2) is

(3.12)\{\begin{array}{l}{u}_0^{\delta }=u^{\delta },\\ u_{j+1}=u\left({\displaystyle \frac{1}{2}}||u-u^{\delta }|{|}^2+B_J^{p_j}(u,u_j)\right),\\ u_{j+1}^{\delta }=u_j^{\delta }+(u^{\delta }-u_{j+1}),\\ p_{j+1}=(u_{j+1}^{\delta }-u^{\delta }).\end{array}

The third equation in (3.12) shows that the method makes up for the lost information (u^{\delta }-u_{j+1}) when the next iteration is performed. This reduces the loss of image details in the ROF Model. The algorithm does not need to iterate to convergence in order to obtain the desired result in denoising. If the algorithm is iterated to convergence, then u_{j+1}\to u^{\delta } is obtained. Obviously, this is not the desired result, but it is based on this convergence that Bregman primitive iteration can be effectively applied to problems such as compressed perception [20, 98].

The Bregman split iteration algorithm solves the ROF Model by introducing intermediate variables v=\mathrm{\nabla }u to transform the unconstrained problem into a constrained problem, and then solving it by the penalty method

\underset{u,v}{min}\left(\frac{1}{2}\Vert u-u^{\delta }{\Vert }_{L^2(\mathrm{\Omega })}^2+\lambda \Vert |v\Vert |{|}_{L^1(\mathrm{\Omega })}\right).

However, the difference from the common penalty method is that for a fixed penalty factor \lambda, the Bregman primitive iterative method to solve the above equation and iterate until convergence, we can obtain the ROF Model and the constraint v=\mathrm{\nabla }u holds automatically. This avoids the numerical instability caused by the penalty factor \lambda \to \mathrm{\infty }, which is usually chosen to ensure that the constraint v=\mathrm{\nabla }u holds. And it is possible to change the linearity in the numerical solution by appropriately selecting the penalty factor \lambda to improve the efficiency of the solution.

3.6 Lagrange method of enlargement

Tai et al. [126] pointed out that the Bregman method for solving (3.2) is effective because it is essentially an incremental Lagrange method [61]. Consider the general case, and let J(u) be a convex generic function. For the equationally constrained optimization problem

\underset{u}{min}J(u),\mathrm{s}.\mathrm{t}.G(u)=0,

whose increasing Lagrange function is

(3.13)L(u,\mu )=J(u)+\mu G(u)+\frac{\lambda }{2}\Vert G(u){\Vert }^2.

Here, \mu is the Lagrange multiplier, and \lambda is the penalty parameter. Solve the saddle point problem corresponding to (3.13)

\{\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }L(u,\mu )}{\mathrm{\partial }u}}=0,\\ {\displaystyle \frac{\mathrm{\partial }L(u,\mu )}{\mathrm{\partial }\mu }}=0.\end{array}

Then u is obtained.

3.7 Original-parity hybrid gradient method

Solving the original problem and the dual problem of (3.2) have their own advantages and disadvantages. Chen et al. [137] proposed to mix them together in the original-parity hybrid gradient method, i.e., given the initial guess value u_0 to solve

(3.14)\{\begin{array}{l}{w}_{k+1}=P_{|w|⩽1}(w_k-{\tau }_k\mathrm{\nabla }{u}_k),\\ u_{k+1}=u_k-{\theta }_k(K^{\ast }(K{u}_k-u^{\delta })+\lambda \mathrm{\nabla }\cdot w_{k+1}).\end{array}

Here, P is the projection operator, \lambda is the parameter controlling the step size, [138] gives a number of selection criteria. Numerical simulations show that this algorithm converges faster than the Bregman split iteration method and the Chambolle pairwise algorithm. Chan et al. also studied the relationship between the original-dual hybrid gradient method and some commonly used algorithms [51], for example, the forward-backward splitting approximation method [40, 89], alternating direction minimization algorithm [120], alternating direction multiplier method [55], Douglas-Rachford splitting method [48], the non-strict splitting Uzawa method [83], the mean mean gradient method [124], etc.

Other common algorithms include the graph cut method [43], the quadratic programming method [54], the second-order cone programming method [62], Extremely small method [76], splitting algorithm [67], etc.

4 Summary and outlook

This paper reviews the important models and various numerical solutions in variational regularized image restoration techniques. However, there are still many challenging problems in the field of image restoration, which are summarized as follows.

1) For texture-based images, both noise and texture are oscillatory components, and many models cannot effectively distinguish between them. Although Zhang and Fan in [133] achieved good results by using the local variance information of both as constraints to build variational optimization models, they did not solve the essential difficulty of modeling texture-based images. The adaptive texture model proposed by Peyré et al. in [86] for modeling anisotropic local parallel textures can well reflect the geometric structure information of local parallel textures, but the model is a nonconvex variational model, the convergence of the algorithm is not guaranteed [101]. Therefore, further research is needed to find a model that can differentiate texture and noise adaptively for denoising.

2) For some images, such as medical images, the noise generated by the imaging system is image-dependent, and such noise is neither additive nor multiplicative. Therefore, modeling and algorithms for such noise need to be investigated.

3) In practical problems, the convolution operator corresponding to linear fuzzy is often unknown. Therefore, the study of blind deconvolution is more challenging. The modeling of nonlinear blurred images is usually simplified to approximate linear blur. Although this is computationally convenient, but the image recovery effect is often poor, so it is worthwhile to explore the appropriate nonlinear model.

4) A large number of numerical experiments have shown that the original-dual hybrid gradient method and Bregman split iteration method are faster than other common algorithms. The convergence speed of these two algorithms has yet to be analyzed theoretically.

5) The localization of the damaged region in image restoration is currently considered as a separate problem, and almost all of the literature is based on a template to mark the location of the damaged region before the restoration process. This is extremely inconvenient in practical applications. Therefore, it is important to investigate how to automatically detect and locate the damaged area, or to perform blind repair without locating the damaged area.

6) The research on video image restoration techniques is still in the preliminary stage. The existing approaches usually use static image restoration techniques for video image restoration based on the spatial correlation of each frame, but rarely consider the principle of temporal correlation of each frame. Therefore, the video image restoration model and its corresponding algorithm combining the principles of temporal and spatial correlation should be studied.

7) The image restoration results are sensitive to the regularization parameters, and the selection of the regularization parameters often depends on empirical debugging, and there are few quantitative analysis results available. Therefore, theoretically, the quantitative study of the regularization parameters is a difficult problem in the inverse problem. In addition, from the development of image restoration models, new models often originate from some improvement of the regularization term or approximation term of the original model. The model proposed in this paper (3.1) unifies the study of denoising, deblurring, restoration, compressed perception, etc., and can be used to select the Omega model that can reasonably describe the mathematical mechanisms of image acquisition, storage and display systems in specific problems. We can select the \mathrm{\Omega },{\mathrm{\Omega }}^{\prime },D,p,q,r,K,T that can reasonably describe these processes and unify the study of their suitability, the selection of regularization parameters and the algorithms from the perspective of solving the inverse problems by regularization methods. The problem of the selection of the regularization parameters and the design of the algorithm are also studied.

8) The existing image restoration evaluation system is not perfect, and people's judgment of the processing effect is based on visual observation to a large extent. The existing objective evaluation indexes, such as PSNR indexes, do not take into account the visual characteristics of human eyes, and there is no original image to refer to in practical applications. Therefore, it remains to be studied how to establish quantitative objective evaluation criteria to guide modeling and computation in accordance with people's visual cognitive system.

References

[1] Alliney S. A property of the minimum vectors of a regularizing functional defined by means of the absolute norm. IEEE Trans Signal Processing 1997; 45(4): 913–917

[2] AmbrosioLFuscoNPallaraD. Functions of Bounded Variation and Free Discontinuity Problems. New York: Clarendon Press, 2000

[3] Aubert G, Aujol J F. Modeling very oscillating signals: application to image processing. Appl Math Optim 2005; 51(2): 163–182

[4] Aubert G, Aujol J F. A variational approach to remove multiplicative noise. SIAM J Appl. Math. 2008; 68(4): 925–946

[5] Aujol J F. Dual norms and image decomposition models. Int. J. Comput. Vision 2005; 63(1): 85–104

[6] Aujol J F, Gilboa G, Chan T F, Osher S J. Structure-texture image decomposition modeling, algorithms and parameter selection. Int J Comput Vis 2005; 67(1): 111–136

[7] Ballester C, Bertalmio M, Caselles V, Sapiro G, Verdera J. Filling-in by joint interpolation of vector fields and grey levels. IEEE Trans Image Processing 2001; 10(8): 1200–1211

[8] Banham M R, Katsaggelos A K. Digital image restoration. IEEE Signal Processing Mag 1997; 14(2): 24–41

[9] Bardsley J M, Luttman A. Total variation-penalized Poisson likelihood estimation for ill-posed problems. Advances in Computational Mathematics, Special Volume on Mathematical Methods for Image Processing 2009; 31(13): 35–59

[10] Bardsley J M. An efficient computational method for total variation-penalized Poisson likelihood estimation. Inverse Problems and Imaging 2008; 2(2): 167–185

[11] BertalmioMSapiromGCasellesVBallesterC. Image inpainting. In: Proceedings of the 27th Annual Conference on Computer Graphics and Interactive Techniques, New Orleans, USA, July 2000, 417‒424

[12] BoydSVandenbergheL. Convex Optimization. Cambridge University Press, 2004

[13] Bresson X, Chan T F. Fast dual minimization of the vectorial total variation norm and applications to color image processing. Inverse Problems and Imaging 2008; 2(4): 455–484

[14] BressonXChanT F. Non-local unsupervised variational image segmentation models. UCLA CAM Report, 2008

[15] Cai J F, Osher S J, Shen Z W. Convergence of the linearized Bregman iteration for l₁-norm minimization. Math Comput 2009; 78(268): 2127–2136

[16] Cai J F, Chan R H, Shen Z W. A framelet-based image inpainting algorithm. Appl Comput Harmon Anal 2008; 24(2): 131–149

[17] CaiJ FDongBOsherS JShenZ W. Image restoration: total variation; wavelet frames; and beyond. UCLA CAM Report, 2011

[18] CaiJ FJiHLiuC QShenZ W. Blind motion deblurring from a single image using sparse approximation. In: IEEE Conf Computer Vision and Pattern Recognition, CVPR 2009

[19] Cai J F, Ji H, Liu C Q, Shen Z W. Blind motion deblurring using multiple images. J Computational Physics 2009; 228(14): 5057–5071

[20] Cai J F, Osher S J, Shen Z W. Linearized Bregman iterations for compressed sensing. Math Comput 2009; 78(267): 1515–1536

[21] Cai J F, Osher S J, Shen Z W. Linearized Bregman Iterations for frame based image deblurring. SIAM J Imaging Sciences 2009; 2(1): 226–252

[22] Cai J F, Osher S J, Shen Z W. Split Bregman methods and frame based image restoration. Multiscale Model Simul 2009; 8(2): 337–369

[23] Candès E, Romberg J, Tao T. Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information. IEEE Trans Information Theory 2006; 52(2): 489–509

[24] Carasso A S. Linear and nonlinear image deblurring: a documented study. SIAM J Numer Anal 1999; 36(6): 1659–1689

[25] Carasso A S. Singular integrals, image smoothness, and the recovery of texture in image deblurring. SIAM J Appl Math 2004; 64(5): 1749–1774

[26] CarterJ L. Dual Methods for Total Variation Based Image Restoration, Ph D thesis. Los Angeles, CA: UCLA, 2001

[27] Chambolle A. An algorithm for total variation minimization and applications. Journal of Mathematical Imaging and Vision 2004; 20(12): 89–97

[28] Chan R H, Ho C W, Nikolova M. Salt-and-pepper noise removal by median type noise detectors and detail preserving regularization. IEEE Trans Image Processing 2005; 14(10): 1479–1485

[29] Chan T F, Esedoglu S. Aspects of total variation regularized function approximation. SIAM J Appl Math 2005; 65(5): 1817–1837

[30] Chan T F, Shen J. Variational image inpainting. Communications on Pure and Applied Mathematics 2005; 58(5): 579–619

[31] Chan T F, Shen J H. Nontexture inpainting by curvature-driven diffusions. J Visual Commun Image Rep 2001; 12(4): 436–449

[32] Chan T F, Shen J H. Mathematical models for local nontexture inpaintings. SIAM J Appl Math 2002; 62(3): 1019–1043

[33] ChanT FShenJ H. Image Processing and Analysis: Variational, PDE, Wavelet, and Stochastic Methods. Philadelphia: SIAM, 2005

[34] ChanT FShenJ H. Theory and computation of variational image deblurring. In: IMS Lecture Notes, 2006, 93‒131

[35] Chan T F, Wong C K. Total variation blind deconvolution. IEEE Trans Image Processing 1998; 7(3): 370–375

[36] Chan T F, Golub G H, Mulet P. A nonlinear primal-dual method for total variation based image restoration. SIAM J Sci Comp 1999; 20(6): 1964–1977

[37] Chan T F, Shen J H, Zhou H M. Total variation wavelet inpainting. Journal of Mathematical Imaging and Vision 2006; 25(1): 107–125

[38] ChenTHuangT SYinW TZhouX S. A new coarse-to-fine framework for 3D brain MR image registration. In: Computer Vision for Biomedical Image, Lecture Notes in Computer Science, 2005, 3765: 114‒124

[39] Cohen A, Dahmen W, DeVore R. Compressed sensing and best k-term approximation. J American Mathematical Society 2009; 22(1): 211–231

[40] Combettes P L, Wajs V R. Signal recovery by proximal forward backward splitting. SIAM J Multiscale Modeling and Simulation 2005; 4(4): 1168–1200

[41] CostanzinoN. EN161 Project Presentation Ⅲ: Structure Inpainting via Variational Methods, 2004

[42] Dai W, Sheikh M, Milenkovic O, Baraniuk R G. Compressed sensing DNA microarrays. EURASIP J Bioinformatics and Systems Biology 2009; 162824

[43] Darbon J, Sigelle M. Image restoration with discrete constrained total variation part I: fast and exact optimization. J Math Imaging Vis 2006; 26(3): 261–276

[44] Daubechies I, Defrise M, De Mol C. An iterative thresholding algorithm for linear inverse problems with a sparsity constraint. Comm Pure and Applied Math 2004; 57(11): 1413–1457

[45] Daubechies I, Teschke G. Variational image restoration by means of wavelets: simultaneous decomposition, deblurring and denoising. Applied and Computational Harmonic Analysis 2005; 19(1): 1–16

[46] DeVore R A. Nonlinear approximation. Acta Numerica 1998; 7: 51–150

[47] Donoho D. Compressed sensing. IEEE Trans Information Theory 2006; 52(4): 1289–1306

[48] Eckstein J, Bertsekas D P. On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators. Mathematical Programming 1992; 55(13): 293–318

[49] EkelandITemamR. Convex Analysis and Variational Problems. Philadelphia: SIAM, 1999

[50] Esedogl S, Shen J H. Digital inpainting based on the Mumford-Shah-Euler image model. European J Appl Math 2002; 13(4): 353–370

[51] EsserEZhangX QChanT F. A general framework for a class of first order primal-dual algorithms for TV minimization. UCLA CAM Report (09–67), 2009

[52] EvansL CGariepR F. Measure Theory and Fine Properties of Functions. Boca Raton, FL: CRC Press, 1992

[53] Fish D A, Brinicombe A M, Pike E R, Walker J G. Blind deconvolution by means of the Richardson-Lucy algorithm. J Opt Soc Am A 1996; 12(1): 58–65

[54] Fu H Y, Ng M K, Nikolova M, Barlow J L. Efficient minimization methods of mixed l₂-l₁ and l₁-l₁ norms for image restoration. SIAM J Scientific Computing 2006; 27(6): 1881–1902

[55] Gabay D, Mercier B. A dual algorithm for the solution of non-linear variational problems via finite element approximations. Comp Math Appl 1976; 2(1): 17–40

[56] GanL. Block compressed sensing of natural images. In: Proc Int Conf Digital Signal Processing, Cardiff, UK, 2007

[57] Garnett J B, Le T M, Meyer Y, Vese L A. Image decompositions using bounded variation and generalized homogeneous Besov spaces. Appl Comput Harmon Anal 2007; 23(1): 25–56

[58] Geman S, Geman D. Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images. IEEE Trans Pattern Anal Mach Intell 1984; 6(6): 721–741

[59] Gilboa G, Sochen N, Zeevi Y. Texture preserving variational denoising using an adaptive fidelity term. IEEE Trans Image Processing 2006; 15(8): 2281–2289

[60] Gilles J, Meyer Y. Properties of BV-G structures + textures decomposition models: application to road detection in satellite images. IEEE Trans Image Processing 2010; 19(11): 2793–2800

[61] GlowinskiRLe TallecP. Augmented Lagrangians and Operator-Splitting Methods in Nonlinear Mechanics. Philadelphia: SIAM, 1989

[62] Goldfarb D, Yin W T. Second-order cone programming methods for total variation based image restoration. SIAM J Scientific Computing 2006; 27(2): 622–645

[63] Goldstein T, Osher S J. The split Bregman method for L₁-regularized problems. SIAM J Imaging Sci 2009; 2(2): 323–343

[64] Grasmair M. Generalized Bregman distances and convergence rates for non-convex regularization methods. Inverse Problems 2010; 26(11): 115014

[65] HamzaA BKrimH. A variational approach to maximum a posteriori estimation for image denoising. In: Energy Minimization Methods in Computer Vision and Pattern Recognition. Lecture Notes in Computer Science, 2001, 19–34

[66] Hamza A B. Krim H, Unal B. Unifying probabilistic and variational estimation. IEEE Signal Process Mag 2002; 19(5): 37–47

[67] Huang Y, Ng M K, Wen Y W. A fast total variation minimization method for image restoration. SIAM J Multiscale Modeling and Simulation 2008; 7(2): 774–795

[68] Hubel D H, Wiesel T N. Receptive fields, binocular intersection and functional architecture in the cat’s visual cortex. J Physiology 1962; 160(1): 106–154

[69] Ito K, Kunisch K. Semi-smooth Newton methods for variational inequalities of the first kind. Math Modelling and Numerical Analysis 2003; 37(1): 41–62

[70] JiaoY LFanQ B. An Efficient Model for Distinguishing Between Texture and Noise. Tech Report, School of Math, Stat, Wuhan University, Nov 2010

[71] JohnsonC. Numerical Solution of Partial Differental Equations by the Finite Element Method. UK: Cambridge University Press, 1987

[72] KanizsaG. Organization in Vision. New York: Praeger, 1979

[73] Krishnan D, Lin P, Yip A M. A primal-dual active-set method for non-negativity constrained total variation deblurring problems. IEEE Trans Image Processing 2007; 16(11): 2766–2777

[74] Krishnan D, Pham Q V, Yip A M. A primal-dual active-Set algorithm for bilaterally constrained total variation deblurring and piecewise constant Mumford-Shah segmentation problems. Advances in Comput Math 2009; 31(13): 237–266

[75] Kunisch K, Hintermüller M. Total bounded variation regularization as a bilaterally constrained optimization problem. SIAM J Appl Math 2004; 64(4): 1311–1333

[76] LangeK L. Optimization. New York: Springer-Verlag, 2004

[77] Le T M, Vese L A. Image decomposition using total variation and div(BMO), Multiscale Modeling and Simulation. SIAM Interdis Journal 2005; 4(2): 390–423

[78] LeVequeR J. Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-state and Time-dependent Problems (Classics in Applied Mathematics). SIAM, 2007

[79] Lu D Y, Fan Q B. Characterizations of L_p(R) using tight wavelet frames. Wuhan University J of Natural Sciences 2010; 15(6): 461–466

[80] Lu D Y, Fan Q B. Gabor frames generated by multiple generators. Scientia Sinica Mathematica 2010; 40(7): 693–708

[81] Lu D Y, Fan Q B. A class of tight framelet packets. Czechoslovak Mathematical Journal 2011; 61(3): 623–639

[82] Lustig M, Donoho D, Pauly J M. Sparse MRI: the application of compressed sensing for rapid MR imaging. Magnetic Resonance in Medicine 2007; 58(6): 1182–1195

[83] Malgouyres F, Zeng T. A predual proximal point algorithm solving anon negative basis pursuit denoising model. Int Journal of Computer Vision 2009; 83(3): 294–311

[84] MallatS. A Wavelet Tour of Signal Processing, 3rd Ed. Academic Press, 2008

[85] Marr D, Hildreth E. Theory of edge detection. Proc Roy Soc London, Ser B 1980; 207(1167): 187–217

[86] MaurelPAujolJ FPeyréG. Locally parallel textures modeling with adapted Hilbert spaces. In: Lecture Notes in Computer Science, 2009, 5681: 429‒442

[87] MeyerY. Wavelets and Operators. UK: Cambridge University Press, 1992

[88] MeyerY. Oscillating Patterns in Image Processing and Nonlinear Evolution Equations. University Lecture Series, Vol 22, Providence, RI: American Mathematical Society, 2001

[89] Micchelli C A, Shen L.X, Xu Y S. Proximity algorithms for image models: denoising. Inverse Problems 2011; 27(4): 045009

[90] MorozovV A. Methods for Solving Incorrectly Posed Problems. New York: Springer-Verlag, 1984

[91] NealenAAlexaM. Hybrid texture synthesis. In: Proceedings of the Eurographics Symposium on Rendering, 2003, 97‒105

[92] Ng M K, Qi L Q, Tang Y F, Huang Y M. On semismooth Newton’s methods for total variation minimization. Journal of Mathematical Imaging and Vision 2007; 27(3): 265–276

[93] Nikolova M. A variational approach to remove outliers and impulse noise. Journal of Mathematical Imaging and Vision 2004; 20(1-2): 99–120

[94] Nikolova M. Model distortions in Bayesian map reconstruction. Inverse Problems and Imaging 2007; 1(2): 399–422

[95] OliveiraM MBowenBMcKennaRChangY S. Fast digital image inpainting. In: Proceedings of the International Conference on Visualization, Imaging and Image Processing (VⅡP 2001), Marbella, Spain, Sep 3‒5, 2001

[96] Osher S, Burger M, Goldfarb D, Xu J J, Yin W T. An iterative regularization method for total variation based image restoration. Multiscale Model Simul 2005; 4(2): 460–489

[97] OsherS JFedkiwR P. Level Set Methods and Dynamic Implicit Surfaces. New York: Springer-Verlag, 2003

[98] Osher S J, Mao Y, Dong B, Yin W T. Fast linearized Bregman iteration for compressed sensing and sparse denoising. Communications in Mathematical Sciences 2010; 8(1): 93–111

[99] Osher S J, Solé A, Vese L A. Image decomposition and restoration using total variation minimization and the H₋₁ norm. Multiscale Modeling and Simulation 2003; 1(3): 349–370

[100] Peyré G. Texture synthesis with grouplets. IEEE Trans on Pattern Analysis and Machine Inte 2010; 32(4): 733–746

[101] Peyré G. A review of adaptive image representations. IEEE Journal of Selected Topics in Signal Processing 2011; 5(5): 896–911

[102] Pruessner A, O’Leary D P. Blind deconvolution using a regularized structured total least norm algorithm. SIAM J Matrix Anal Appl 2003; 24(4): 1018–1037

[103] PuetterR C. Pixons and Bayesian image reconstruction. In: Proceedings of SPIE, 1994, 23(2): 112‒131

[104] RabbatMHauptJSinghANowakRD. Decentralized compression and predistribution via randomized gossiping. In: Proc Int Conf on Information Processing in Sensor Networks, Nashville, Tennessee, 2006

[105] RahmanTTaiX COsherS J. A TV-Stokes denoising algorithm. In: Lecture Notes in Computer Science, Berlin: Springer-Verlag, 2007, 4485: 473‒483

[106] Rudin L I, Osher S J, Fatemi E. Nonlinear total variation based noise removal algorithms. Physica D 1992; 60: 259–268

[107] Scherzer O. Scale space methods for denoising and inverse problem. Advances in Imaging and Electron Physics 2003; 128: 445–530

[108] ScherzerO. Handbook of Mathematical Methods in Imaging. Berlin: Springer-Verlag, 2010

[109] ScherzerOGrasmairMGrossauerHHaltmeierMLenzenF. Variational Methods in Imaging. New York: Springer-Verlag, 2009

[110] ShikhminA. Synthesizing natural textures. In: Proceedings of 2001 ACM Symposium on Interactive 3D Graphics, 2001, 217‒226

[111] Snyder D L, Hammoud A M, White R L. Image recovery from data acquired with a charged-coupled-device camera. J Opt Soc Am A 1993; 10(5): 1014–1023

[112] SongB. Topics in Variational PDE Image Segmentation, Inpainting and Denoising. Ph D thesis, Los Angeles, CA: University of California, 2003

[113] StarckJ LMurtaghFFadiliJ M. Compressed sensing. In: Sparse Image and Signal Processing: Wavelets, Curvelets, Morphological Diversity, Cambridge University Press, 2010

[114] StoecklerJWellandG V. Beyond Wavelet. New York: Academic Press, 2003

[115] TaiX COsherS JHolmR. Image inpainting using a TV-Stokes equation. In: Image Processing Based on Partial Differential Equations, Heidelberg: Springer-Verlag, 2007, 3‒22

[116] TakharDLaskaJWakinM BDuarteM FBaronDSarvothamSKellyK FBaraniukR G. A new compressive imaging camera architecture using optical-domain compression. In: Proc Computational Imaging IV at SPIE Electronic Image, San Jose, California, 2006

[117] TikhonovA NArseninV Y. Solutions of Ill-Posed Problems. Washington D C: John Wiley and Sons, 1977

[118] TriebelH. Theory of Function Spaces Ⅱ. Monographs in Mathematics, Vol 84, Boston: Birkhäuser, 1992

[119] Tsai A, Yezzi A R, Willsky A S. Curve evolution implementation of the Mumford-Shah functional for image segmentation, denoising, interpolation, and magnification. IEEE Trans Image Processing 2001; 10(8): 1169–1186

[120] Tseng P. Applications of a splitting algorithm to decomposition in convex programming and variational inequalities. SIAM J Control and Optimization 1991; 29(1): 119–138

[121] Vese L A, Osher S J. Modeling textures with total variation minimization and oscillating patterns in image processing. J Sci Comput 2003; 19: 553–572

[122] VogelC R. Computational Methods for Inverse Problems. Frontiers Appl Math, Vol 23, Philadelphia: SIAM, 2002

[123] WangWGarofalakisM NRamchandranK. Distributed sparse random projections for refinable approximation. In: Proc Int Conf Information Processing in Sensor Networks, Cambridge, Massachusetts, 2007

[124] Weiss P, Blanc-Féraud L, Aubert G. Efficient schemes for total variation minimization under constraints in image processing. SIAM J Scientific Computing 2009; 31(3): 2047–2080

[125] Woiselle A, Starck J L, Fadili J. 3-D data denoising and inpaint- ing with the low redundancy fast curvelet transform. Journal of Mathematical Imaging and Vision 2011; 39(2): 121–139

[126] Wu C L, Tai X C. Augmented Lagrangian method, dual methods, and split Bregman iteration for ROF, vectorial TV, and high order models. SIAM Journal on Imaging Sciences 2010; 3(3): 300–313

[127] YamauchiHHaberJSeidelH P. Image restoration using multiresolution texture synthesis and image inpainting. In: Proceedings of Computer Graphics International Conference, Japan, 2003, 120‒125

[128] Ye J C. Compressed sensing shape estimation of star-shaped objects in Fourier imaging. IEEE Signal Processing Letters 2007; 14(10): 750–753

[129] Yin W T, Goldfarb D, Osher S J. A comparison of three total variation based texture extraction models. J Visual Communication and Image Representation 2007; 18(3): 240–252

[130] YosidaK. Functional Analysis: Reprint of the 1980 Edition. Berlin: Springer-Verlag, 1995

[131] Zhang H Y, Peng Q C, Wu Y D. Wavelet inpainting based on p-Laplace operator. Acta Automatica Sinica 2007; 33(5): 546–549

[132] Zhang H Y, Wu B, Peng Q C, Wu Y D. Digital image inpainting based on p-harmonic energy minimization. Chinese Journal of Electronics 2007; 3: 525–530

[133] Zhang T, Fan Q B. Wavelet characterization of dyadic BMO norm and its application in image decomposition for distinguishing between texture and noise. Int Journal of Wavelet, Multiresolution and Information Processing 2011; 9(3): 445–457

[134] Zhang T, Fan Q B, Gao Q L. Wavelet characterization of Hardy space and its application in variational image decomposition. Int Journal of Wavelet, Multiresolution and Information Processing 2010; 8(1): 71–87

[135] Zhang X Q, Burger M, Osher S J. A unified primal-dual algorithm gramework based on Bregman iteration. J Sci Comput 2011; 46(1): 20–46

[136] ZhouM Y. Deconvolution and Signal Recovery. Beijing: National Defence Industry Press, 2001 (in Chinese)

[137] ZhuM QChanT F. An efficient primal-dual hybrid gradient algorithm for total variation image restoration. UCLA CAM Report (08–34), 2008

[138] Zhu M Q, Wright J, Chan T F. Duality based algorithms for total variation regularized image restoration. Comput Optim Appl 2008; 47(3): 377–400

[139] ZiemerW P. Weakly Differentiable Functions: Sobolev Spaces and Functions of Bounded Variation. Graduate Texts in Mathematics, Vol 120, New York: Springer-Verlag, 1989

RIGHTS & PERMISSIONS

Higher Education Press 2024