# Optimizing Shrinkage Curves and Application in Image Denoising

Authors: Hongyao Deng; Qingxin Zhu; Jinsong Tao; Xiuli Song
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

A shrinkage curve optimization is proposed for weighted nuclear norm minimization and is adapted to image denoising. The proposed optimization method employs a penalty function utilizing the difference between a latent matrix and its observation and uses odd polynomials to shrink the singular values of the observation matrix. As a result, the coefficients of polynomial characterize the shrinkage operator fully. Furthermore, the Frobenius norm of the penalty function is converted into the corresponding spectral norm, and thus the parameter optimization problem can be easily solved by using off-and-shelf plain least-squares. In the practical application, the proposed denoising method does not work on the whole image at once, but rather a series of matrix termed Rank-Ordered Similar Matrix (ROSM). Simulation results on 256 noisy images demonstrate the effectiveness of the proposed algorithms.


---

## Body

## 1. Introduction

Low rank matrix approximation has been attracting significant research interest in recent years. This approach aims to reconstruct the latent data from its degraded observation matrix and is frequently applied in many fields, such as machine learning [1], computer version [2], recommendation system [3], and image processing [4]. As a branch of this research, a regularized nuclear norm minimization problem is widely considered over matrices,(1)minX∈Rn×m⁡fX+μX∗,where X denotes a matrix, scalar μ>0 is a parameter, and fX and X∗ are the data fidelity term and the data regularization term, respectively. In this formula, X∗=∑i=1min⁡m,nσiX is called the nuclear norm of X, where σiX denotes the ith largest singular value of X.IffX is convex, then fX+μX∗ is a convex function because nuclear norm X∗ is also convex. Thus, problem (1) is a convex optimization problem and can be treated by using various classic iterative optimization algorithms including steepest-descent, conjugate-gradient, and interior-point algorithms. When the data fidelity term fX=1/2X-YF2, where Y is an observation matrix and ·F denotes the Frobenius norm operator, this is the well-known nuclear norm minimization (NNM) [5]. The NNM problem was proved that it can be solved by applying a soft threshold operation on the singular values of Y, and the solution can be achieved using a Singular Value Thresholding (SVT) algorithm [6].Despite the success of NNM, it is not flexible enough to handle more complex issues. To pursue the convex property, NNM treats each singular value equally. As a result, the soft-thresholding operator shrinks each singular value by the same amount [6]. In principal component analysis however, different principal directions quantify different information. For example, the large singular value delivers the major feature information such as edges and texture. This implies that, in image denoising, the larger the singular value is, the lesser the amount shrinks. Obviously, the NNM model and the corresponding solver cannot handle this issue.To overcome this limitation, a regularized nuclear norm minimization with weights was put forward. These weights may enhance the representation capability of the original nuclear norm. Its form is as follows:(2)minX∈Rn×m⁡fX+∑i=1min⁡m,ngiσiX,where gi is the weight designed to σiX.Problem (2) is a nonconvex nonsmooth low rank minimization problem. Of course, if gi is replaced with μ, problem (2) reverts to problem (1). Solving problem (2) is challenging, or even NP-hard. To solve this problem, researchers presented some assumption to handle it. Gu et al. [7] assumed that g is nondescending on [0,∞), 0≤g1≤⋯≤gmin⁡(m,n), and thus problem (2) becomes convex and can be solved by a soft-thresholding operation. Moreover, the authors devised a solver in which gi is inversely proportional to σiX. Gasso et al. [8] argued that if both f and -g are convex, problem (2) can be solved by DC (Difference of Convex functions) programming. Lu et al. [9] assumed that g is concave increasing monotonically on [0,∞) and f satisfies Lipschitz continuous condition; the weights are achieved at the super gradient point of the concave function g. Based on this assumption, the authors proposed Iteratively Reweighted Nuclear Norm (IRNN) method. In addition, Hu et al. [10] reported their Truncated Nuclear Norm Regularization (TNNR) method, based on the same assumptions as in [9].By applying these low rank matrix approximation theories, different image denoising methods have been reported. For example, a method of coupling sparse denoising and unmixing with low rank constraint is proposed for hyperspectral image in [11]; a scheme of incorporating iterative support detection into TNNR is presented to reduce white Gaussian noise in [12]; the eigenvectors of the Laplacian are considered to suppress Gaussian noise in [13]; a weighted nuclear norm minimization model is presented [14] and is used in three applications, that is, image denoising, background subtraction, and image completion. These methods achieve high quality results. A main reason is that all of them employ a powerful patch-based technique.Inspired by weighted nuclear norm minimization and patch-based technique, a parameter optimization method is proposed in this paper. The proposed method utilizes the difference between a latent matrix and its observation to design a penalty function and employs odd polynomials to shrink the singular values of the observation matrix. As a result, the coefficients of polynomial fully characterize the shrinkage operator. Furthermore, for the penalty function, its Frobenius norm is converted into a spectral norm. Thereby, the parameter optimization can be easily solved by using plain least-squares.To validate the effectiveness, the optimization theory is applied in image denoising. Since the proposed method is to optimize shrinkage curves, it is called OSC method. In the practical application, the OSC method does not work on the whole image at once, but rather a series of matrix termed Rank-Ordered Similar Matrix (ROSM, see Definition2). Thirty-two images were tested. Experimental results show that the OSC method achieves better results than the Bilateral Filter; when the noisy standard deviation is less than 20, the results achieved by OSC are better than those by BM3D, and when the noisy standard deviation varies from 20 to 40, the results by OSC are weaker than by BM3D.The contribution of this paper is twofold. Firstly, in the penalty function devised for weighted nuclear norm minimization, the weight representation is replaced by odd polynomials. So, the coefficients of polynomial characterize the role of the weights fully. Furthermore, the Frobenius norm of the penalty function is converted into a spectral norm. Secondly, the proposed optimization method is adapted to image denoising. Experimental results show that the proposed OSC method outperforms the Bilateral Filter and also is superior to the BM3D method on the case of low noise.The rest of paper is organized as follows. In Section2, the shrinkage curve optimization is formulated. Section 3 describes the image denoising algorithm, and the corresponding analysis is followed in Section 4. Section 5 reports the experimental results, and conclusions are drawn in Section 6.


## 2. Optimizing Shrinkage Curves

In this section, the problem to be discussed is formulated, and the method of optimizing shrinkage curves is followed.


### 2.1. Problem Formulation

LetX be a unknown square matrix in Rn×n, and let Y be its observation. The observed matrix is corrupted by white Gaussian noise N with deviation σ2. This is expressed as(3)Y=X+N.To reconstruct the original square matrixX from its noisy version, the following weighted nuclear norm minimization with a constraint is considered:(4)minX∗,w,subject toY-XF2<ε,where ε is a threshold, X∗,w=∑iwiσiX1, and Y-XF2=∑i∑jyij-xij2. In this formula, σiX denotes the ith largest singular value of X, and X∗,w and Y-XF are the weighted nuclear norm of X and the Frobenius norm of Y-X, respectively. Our aim is to use polynomial coefficients to characterize the weights and obtain the solution of the problem.


### 2.2. Optimization Method

As in [7], it is proven that the weighted nuclear norm minimization (4) can be solved by imposing a soft threshold operation on the singular values of observation matrix. The form is as follows:(5)X=USΣVT,where Y=UΣVT is the Singular Value Decomposition (SVD) of Y and S is the soft-thresholding operator. This operator with weight vector w shrinks the singular values; SΣii=max⁡(Σii-wi,0). To obtain the thresholds, the following penalty function is employed:(6)fS=X-USΣVTF2.Assuming that shrinkage operation is applied differently to every singular value in matrixΣ, and thus it can be broken into(7)SΣ=s1Σ10⋯00s2Σ2⋯0⋮⋮⋱⋮00⋯snΣn,where Σ[i] is the ith largest singular value of Y and Si denotes the ith shrinkage operator. In our work, odd polynomials are taken to represent these shrinkage operators, where the coefficients of polynomial characterize the shrinkage operator fully. Thus, the shrinkage operator is expressed as(8)siΣ=∑j=1JcijΣ2×j-1,fori=1,2,…,n.Substituting (8) into (7), the shrinkage operator can be rewritten as(9)ScΣ=∑j=1Jc1jΣ12×j-10⋯00∑j=1Jc2jΣ22×j-1⋯0⋮⋮⋱⋮00⋯∑j=1JcnjΣn2×j-1.Substituting (9) into (6) and considering V is a unitary matrix, the penalty function can be rewritten as(10)fc1,c2,…,cn=XV-UScΣF2.The focus now shifts to optimization of the penalty function. To obtain the optimal solution by using plain least-squares, the Frobenius norm in (10) is converted into a spectral norm. For ease of formulation, a vector c∈RnJ is defined that contains all the coefficients series: {cij:j=1,…,J}, fori=1,2,…,n; that is,(11)cT=c11,c12,…,c1J,c21,…,cn1,cn2,…,cnJ.A matrixZ∈Rn2×nJ is also defined that is a block-diagonal matrix with n blocks,(12)Z=z1z2⋱zn,where each of the blocks z[i] is size of n×J, with the content(13)zi=11⋮1∗Σi1Σi3⋯Σi2×J-1,fori=1,2,…,n.In addition, two operators,vec and diag, are introduced. The vec operator returns a vector by concatenating the columns in a matrix; the diag operator returns a diagonal matrix by putting the elements of a vector on the main diagonal. For example, if M=[m11,m12;m21,m22] is a matrix, vecM returns the vector m=m11,m21,m12,m22T; if m=m1,m2,m3,m4T is a vector, diagm returns the diagonal matrix where the main diagonal is m1,m2,m3,m4.Using these notations and operators, the penalty function (10) can be rewritten as(14)fc=vecXV-diagvecUZc22=b-AZc22,where b=vec(XV) and A=diagvecU.The optimal set of parameters that define the shrinkage curves is(15)∂fc∂c=0=AZTAZc-b,which leads to(16)copt=ZTATAZ-1ZTATb.


## 2.1. Problem Formulation

LetX be a unknown square matrix in Rn×n, and let Y be its observation. The observed matrix is corrupted by white Gaussian noise N with deviation σ2. This is expressed as(3)Y=X+N.To reconstruct the original square matrixX from its noisy version, the following weighted nuclear norm minimization with a constraint is considered:(4)minX∗,w,subject toY-XF2<ε,where ε is a threshold, X∗,w=∑iwiσiX1, and Y-XF2=∑i∑jyij-xij2. In this formula, σiX denotes the ith largest singular value of X, and X∗,w and Y-XF are the weighted nuclear norm of X and the Frobenius norm of Y-X, respectively. Our aim is to use polynomial coefficients to characterize the weights and obtain the solution of the problem.


## 2.2. Optimization Method

As in [7], it is proven that the weighted nuclear norm minimization (4) can be solved by imposing a soft threshold operation on the singular values of observation matrix. The form is as follows:(5)X=USΣVT,where Y=UΣVT is the Singular Value Decomposition (SVD) of Y and S is the soft-thresholding operator. This operator with weight vector w shrinks the singular values; SΣii=max⁡(Σii-wi,0). To obtain the thresholds, the following penalty function is employed:(6)fS=X-USΣVTF2.Assuming that shrinkage operation is applied differently to every singular value in matrixΣ, and thus it can be broken into(7)SΣ=s1Σ10⋯00s2Σ2⋯0⋮⋮⋱⋮00⋯snΣn,where Σ[i] is the ith largest singular value of Y and Si denotes the ith shrinkage operator. In our work, odd polynomials are taken to represent these shrinkage operators, where the coefficients of polynomial characterize the shrinkage operator fully. Thus, the shrinkage operator is expressed as(8)siΣ=∑j=1JcijΣ2×j-1,fori=1,2,…,n.Substituting (8) into (7), the shrinkage operator can be rewritten as(9)ScΣ=∑j=1Jc1jΣ12×j-10⋯00∑j=1Jc2jΣ22×j-1⋯0⋮⋮⋱⋮00⋯∑j=1JcnjΣn2×j-1.Substituting (9) into (6) and considering V is a unitary matrix, the penalty function can be rewritten as(10)fc1,c2,…,cn=XV-UScΣF2.The focus now shifts to optimization of the penalty function. To obtain the optimal solution by using plain least-squares, the Frobenius norm in (10) is converted into a spectral norm. For ease of formulation, a vector c∈RnJ is defined that contains all the coefficients series: {cij:j=1,…,J}, fori=1,2,…,n; that is,(11)cT=c11,c12,…,c1J,c21,…,cn1,cn2,…,cnJ.A matrixZ∈Rn2×nJ is also defined that is a block-diagonal matrix with n blocks,(12)Z=z1z2⋱zn,where each of the blocks z[i] is size of n×J, with the content(13)zi=11⋮1∗Σi1Σi3⋯Σi2×J-1,fori=1,2,…,n.In addition, two operators,vec and diag, are introduced. The vec operator returns a vector by concatenating the columns in a matrix; the diag operator returns a diagonal matrix by putting the elements of a vector on the main diagonal. For example, if M=[m11,m12;m21,m22] is a matrix, vecM returns the vector m=m11,m21,m12,m22T; if m=m1,m2,m3,m4T is a vector, diagm returns the diagonal matrix where the main diagonal is m1,m2,m3,m4.Using these notations and operators, the penalty function (10) can be rewritten as(14)fc=vecXV-diagvecUZc22=b-AZc22,where b=vec(XV) and A=diagvecU.The optimal set of parameters that define the shrinkage curves is(15)∂fc∂c=0=AZTAZc-b,which leads to(16)copt=ZTATAZ-1ZTATb.


## 3. Application in Image Denoising

In this section, the OSC method is introduced for image denoising, containing denoising modeling data and the denoising algorithm.


### 3.1. Denoising Modeling Data

The proposed OSC method does not work on the whole image at once, but rather a matrix-set in which each matrix contains a fixed number of similar patches extracted from the original noisy images. A patch is first defined as follows.Definition 1.
Y denotes an image, sized N1×N2 pixels. Let yr=Yi,j be a reference pixel. A block of size n×n is extracted from Y, where yr resides at the top-left corner. By applying vec to the block, the n×n  block is identified with a vector in Rn. The corresponding patch is defined by(17)pr=vec⁡Yi,j⋯Yi,j+n-1⋮⋱⋮Yi+n-1,j⋯Yi+n-1,j+n-1=yr⋮yr+n-1,where i∈{1,2,…,N1-n+1} and j∈{1,2,…,N2-n+1}. When all pixels are complete, a patch-set in Rn can be built, denoted by P=pr:r=1,2,…,N1-n+1×N2-n+1.Definition 2.
Letpr∈P be a reference patch and · denote the Euclidean norm; the distance between pr and pj can be calculated, using the metric(18)drj=pr-pj2,forpj∈P.
These scalar distances are then sorted and the patches inP are correspondingly ordered,(19)dr1≤dr2≤⋯≤drM⟶p1≺p2≺⋯≺pM,where dr(j) denotes the jth smallest distance value, ≺ denotes the order relation between patches, and M=P=N1-n+1×N2-n+1 denotes the number of patches in P. The denoising modeling data, termed Rank-Ordered Similar Matrix (ROSM), are defined as(20)Yr=ROSMr=p1,p2,…,pn.
Obviously, ROSM is a square matrix of sizen×n and p(1)=pr. When all patches in P are complete, a matrix-set in Rn×n can be built, denoted by M~=Yr:r=1,2,…,M.


### 3.2. Denoising Algorithm

The proposed denoising algorithm consists of Algorithms1 and 2. The former trains the polynomial coefficients copt, and the latter reduces noise.Algorithm 1:Training the coefficients copt.
Input:N paired images, Xi,Yi∣i=1,2,…,N
(1)   Initialize the table W=0;
(2)   for  t=1:2
(3)    Initialize a parameter coptt=[0]; // of size nJ×1.
(4)    for  i=1:N
(5)      Initialize S1=[0];
(6)      Initialize s2=[0];
(7)      Get a paired images, X,Y←Xi,Yi;
(8)      Based on Definition 1, extract all patches from image Y and build the patch-set P=pr
(9)      for each patch pr in P  do
(10)       Based on Definition 2, obtain a ROSM Yr;
(11)       Obtain the matrix Xr corresponding to Yr;
(12)       Singular value decomposition, UrΣrVrT=SVDYr;
(13)       Map XrVr to a vector, br=vecXrVr;
(14)       Map Ur to a diagonal matrix, Ar=diagvecUr;
(15)       Map Σr to a diagonal block matrix Zr, according to Eq. (12) and (13);
(16)       Accumulation, S1=S1+ZrTArTArZr;
(17)       Accumulation, s2=s2+ZrTArTbr;
(18)     end for
(19)     Obtain a optimized parameter, coptt=S1-1s2;
(20)     if  t=1  do
(21)        for each patch  pr in Y  do
(22)         Obtain the estimation Y^r=UrScΣrVrT;
(23)         Plug Y^r into the image canvas of the noisy image Y;
(24)        end for
(25)        Obtain a new Yi←Average the pixels for fixed position in the image canvas;
(26)     end if
(27)     Save the coptt to Table W;
(28)    end for
(29)    end for
Output: TableW that containing the parameters copt.Algorithm 2:Optimizing shrinkage curves based denoising algorithm.
Input: Noisy imageY and the shrinkage parameters ct
(1)    Initialize the intermediate image X^0=Y-127/128 and adjusted image Y^0=Y-127/128;
(2)      MIN=mix(Y-127/128);
(3)      MAN=max(Y-127/128);
(4)    for  t=1:3
(5)    Set a counter to zeros, D=0;
(6)    Obtain an adjusted estimate, Y^t=X^t-1+γY-X^(t-1);
(7)    According to Definition 1, build the patch-set P=pr, associated with image Y;
(8)    for each patch pr in P  do
(9)      According to Definition 2, obtain a ROSM Yr;
(10)      Singular value decomposition, UrΣrVrT=SVDYr;
(11)      Find the parameter ct corresponding to Yr from Table W;
(12)      Obtain the estimation Y^r=UrScΣrVrT;
(13)      if  (Y^r<MIN), then Y^r=Yr;
(14)      if  (Y^r>MAN), then Y^r=Yr;
(15)      Put Y^r back into the image canvas;
(16)      In counter D, the entries associated with pixels in Y^r is added by 1;
(17)    end for
(18)    Obtain an estimated image, X^t← obtained canvas image is entry-wisely divided by D;
(19)  end for
Output: The final estimated imageY^3×128+127.SupposeY is a noisy image and X is its noise-free version. Based on Definition 2, the matrix-set M~ of Y is built, and the corresponding noise-free matrix-set is also obtained, denoted by M=Xr:r=1,…,M. Thereby, there exist M paired samples Xr,Yr:r=1,…,M to train copt. The corresponding penalty function is as follows:(21)fc=∑r=1MXrVr-UrScΣrF2,where UrΣrVrT=SVDYr. The optimum parameters are(22)copt=∑r=1MZrTArTArZr-1∑r=1MZrTArTbr,where br=vec(XrVr), Ar=diag(vec(Ur)), and Zr can be obtained by using (12) and (13).Next, the two-iteration denoising algorithm is introduced. SupposeY is a noisy image to be processed, of size N1×N2. Employing Definition 2, the matrix-set M~ of Y is built. For any ROSM Yr in M~, the estimate of its noise-free version can be obtained by(23)X^r=UrScΣrVrT,where Yr=UrΣrVrT is the SVD of Yr and ScΣk is expressed as (10). When all ROSMs in M~ are complete, an estimated-set can be built, denoted by M^=X^r:r=1,2,…,M. Every estimate X^r is put back into the original image canvas, then intensities of pixels falling in the same position in the canvas are averaged. And thus, the first-estimate of the noise-free image is obtained, denoted by X^(t=1). Repeating the above operations again, the second-estimate X^(t=2) is also yielded.In addition, a scaled version of the residual error, the difference between the noisy and estimated image, is considered. LetRr denote the operation putting the rth estimate X^r back into the original image canvas; the denoising method can be expressed as the following tuples formula:(24)X^t=∑r=1MRrUrSctΣrVrT⊘DY^t=X^t+γY-X^t,fort=1,2,where X^t is the tth image estimate, Y^t is the tth adjustment, γ is the scaled factor of residual error, ⊘ denotes entry-wise division, and D denotes the matrix with the same size as the noisy image, in which each entry records the number that the pixel on the same location is processed.


## 3.1. Denoising Modeling Data

The proposed OSC method does not work on the whole image at once, but rather a matrix-set in which each matrix contains a fixed number of similar patches extracted from the original noisy images. A patch is first defined as follows.Definition 1.
Y denotes an image, sized N1×N2 pixels. Let yr=Yi,j be a reference pixel. A block of size n×n is extracted from Y, where yr resides at the top-left corner. By applying vec to the block, the n×n  block is identified with a vector in Rn. The corresponding patch is defined by(17)pr=vec⁡Yi,j⋯Yi,j+n-1⋮⋱⋮Yi+n-1,j⋯Yi+n-1,j+n-1=yr⋮yr+n-1,where i∈{1,2,…,N1-n+1} and j∈{1,2,…,N2-n+1}. When all pixels are complete, a patch-set in Rn can be built, denoted by P=pr:r=1,2,…,N1-n+1×N2-n+1.Definition 2.
Letpr∈P be a reference patch and · denote the Euclidean norm; the distance between pr and pj can be calculated, using the metric(18)drj=pr-pj2,forpj∈P.
These scalar distances are then sorted and the patches inP are correspondingly ordered,(19)dr1≤dr2≤⋯≤drM⟶p1≺p2≺⋯≺pM,where dr(j) denotes the jth smallest distance value, ≺ denotes the order relation between patches, and M=P=N1-n+1×N2-n+1 denotes the number of patches in P. The denoising modeling data, termed Rank-Ordered Similar Matrix (ROSM), are defined as(20)Yr=ROSMr=p1,p2,…,pn.
Obviously, ROSM is a square matrix of sizen×n and p(1)=pr. When all patches in P are complete, a matrix-set in Rn×n can be built, denoted by M~=Yr:r=1,2,…,M.


## 3.2. Denoising Algorithm

The proposed denoising algorithm consists of Algorithms1 and 2. The former trains the polynomial coefficients copt, and the latter reduces noise.Algorithm 1:Training the coefficients copt.
Input:N paired images, Xi,Yi∣i=1,2,…,N
(1)   Initialize the table W=0;
(2)   for  t=1:2
(3)    Initialize a parameter coptt=[0]; // of size nJ×1.
(4)    for  i=1:N
(5)      Initialize S1=[0];
(6)      Initialize s2=[0];
(7)      Get a paired images, X,Y←Xi,Yi;
(8)      Based on Definition 1, extract all patches from image Y and build the patch-set P=pr
(9)      for each patch pr in P  do
(10)       Based on Definition 2, obtain a ROSM Yr;
(11)       Obtain the matrix Xr corresponding to Yr;
(12)       Singular value decomposition, UrΣrVrT=SVDYr;
(13)       Map XrVr to a vector, br=vecXrVr;
(14)       Map Ur to a diagonal matrix, Ar=diagvecUr;
(15)       Map Σr to a diagonal block matrix Zr, according to Eq. (12) and (13);
(16)       Accumulation, S1=S1+ZrTArTArZr;
(17)       Accumulation, s2=s2+ZrTArTbr;
(18)     end for
(19)     Obtain a optimized parameter, coptt=S1-1s2;
(20)     if  t=1  do
(21)        for each patch  pr in Y  do
(22)         Obtain the estimation Y^r=UrScΣrVrT;
(23)         Plug Y^r into the image canvas of the noisy image Y;
(24)        end for
(25)        Obtain a new Yi←Average the pixels for fixed position in the image canvas;
(26)     end if
(27)     Save the coptt to Table W;
(28)    end for
(29)    end for
Output: TableW that containing the parameters copt.Algorithm 2:Optimizing shrinkage curves based denoising algorithm.
Input: Noisy imageY and the shrinkage parameters ct
(1)    Initialize the intermediate image X^0=Y-127/128 and adjusted image Y^0=Y-127/128;
(2)      MIN=mix(Y-127/128);
(3)      MAN=max(Y-127/128);
(4)    for  t=1:3
(5)    Set a counter to zeros, D=0;
(6)    Obtain an adjusted estimate, Y^t=X^t-1+γY-X^(t-1);
(7)    According to Definition 1, build the patch-set P=pr, associated with image Y;
(8)    for each patch pr in P  do
(9)      According to Definition 2, obtain a ROSM Yr;
(10)      Singular value decomposition, UrΣrVrT=SVDYr;
(11)      Find the parameter ct corresponding to Yr from Table W;
(12)      Obtain the estimation Y^r=UrScΣrVrT;
(13)      if  (Y^r<MIN), then Y^r=Yr;
(14)      if  (Y^r>MAN), then Y^r=Yr;
(15)      Put Y^r back into the image canvas;
(16)      In counter D, the entries associated with pixels in Y^r is added by 1;
(17)    end for
(18)    Obtain an estimated image, X^t← obtained canvas image is entry-wisely divided by D;
(19)  end for
Output: The final estimated imageY^3×128+127.SupposeY is a noisy image and X is its noise-free version. Based on Definition 2, the matrix-set M~ of Y is built, and the corresponding noise-free matrix-set is also obtained, denoted by M=Xr:r=1,…,M. Thereby, there exist M paired samples Xr,Yr:r=1,…,M to train copt. The corresponding penalty function is as follows:(21)fc=∑r=1MXrVr-UrScΣrF2,where UrΣrVrT=SVDYr. The optimum parameters are(22)copt=∑r=1MZrTArTArZr-1∑r=1MZrTArTbr,where br=vec(XrVr), Ar=diag(vec(Ur)), and Zr can be obtained by using (12) and (13).Next, the two-iteration denoising algorithm is introduced. SupposeY is a noisy image to be processed, of size N1×N2. Employing Definition 2, the matrix-set M~ of Y is built. For any ROSM Yr in M~, the estimate of its noise-free version can be obtained by(23)X^r=UrScΣrVrT,where Yr=UrΣrVrT is the SVD of Yr and ScΣk is expressed as (10). When all ROSMs in M~ are complete, an estimated-set can be built, denoted by M^=X^r:r=1,2,…,M. Every estimate X^r is put back into the original image canvas, then intensities of pixels falling in the same position in the canvas are averaged. And thus, the first-estimate of the noise-free image is obtained, denoted by X^(t=1). Repeating the above operations again, the second-estimate X^(t=2) is also yielded.In addition, a scaled version of the residual error, the difference between the noisy and estimated image, is considered. LetRr denote the operation putting the rth estimate X^r back into the original image canvas; the denoising method can be expressed as the following tuples formula:(24)X^t=∑r=1MRrUrSctΣrVrT⊘DY^t=X^t+γY-X^t,fort=1,2,where X^t is the tth image estimate, Y^t is the tth adjustment, γ is the scaled factor of residual error, ⊘ denotes entry-wise division, and D denotes the matrix with the same size as the noisy image, in which each entry records the number that the pixel on the same location is processed.


## 4. Algorithm Analysis

In this section, the influence of different levels of noise for the singular vectors is first discussed; then the shrunken scales of singular values are analyzed.


### 4.1. Influence of Different Levels of Noise for Singular Vectors

Two matrices,Yα and Yβ, are chosen arbitrarily from the noisy matrix-set M~, and yα=vecYα and yβ=vecYβ. Thus, yα=xα+nα and yβ=xβ+nβ, where x denotes the noise-free vector and n is white noise. The squared Euclidean distance between the two vectors yα and yβ is(25)yα-yβ22=xα-xβ22+2xα-xβ,nα-nβ+nα-nβ22,where the symbol · is the inner-product operator.The vector difference,nα-nβ, is the Gaussian noise with covariance matrix 2σ2n2I, where symbol I denotes identity matrix. Based on the argument in [15], if the dimension of y is large, the norm squared, nα-nβ22, is concentrated around its mean 2σ2n2, with high probability. Similarly, the projection xα-xβ,nα-nβ of the noise on the deterministic vector xα-xβ is concentrated around 0 with high probability. Therefore, (17) can be interpreted as a translation-invariant procedure.On the other hand, as the noise level increases, the singular vectors ofYr∈M~ are no longer aligned with the original singular vectors of noise-free matrix Xr. As a result, the number of reliable singular vectors of Yr becomes smaller. For example, when Xr is corrupted by Gaussian noise at very small magnitude, all the singular vectors of Yr track the noise; when corrupted at a large magnitude, the number of the singular vectors tracking the noise becomes smaller. Therefore, at increasing noise level, the singular vectors of Yr become more and more unreliable and result in worse restored feature.


### 4.2. Shrunken Scales of Different Singular Values

LetYr,Xr be a paired samples and Σi be the ith largest singular value of Yr∈Rn×n. The coefficients vector copt, of size nJ×1, is trained on the paired samples and then is divided into n parts; ciT=ci1,ci2,…,ciJfori=1,2,…,n, where ci corresponds to the ith largest singular value of Yr. And thus, the estimate of the ith largest singular value of Yr is as follows:(26)Σ^i=∑j=1JcijΣi2∗j-1,fori=1,2,…,n.To observe the shrunken scales of the differentΣi visually, the shrinkage curves are shown in Figure 1. Test results show that (i) different singular values Σi have different shrunk scales; (ii) the larger the singular value, the smaller the shrunk scales; (iii) the shrunken scales for fixed index singular value are almost the same although Yr is corrupted by different levels of noise.Figure 1
Illustration of the shrinkage scale of theith largest singular value. From left to right and top to bottom, the ith inset is associated with the ith largest singular value. In each inset, the blue diagonal line is the ground singular value, and red line is its estimate. The parameters: patch-size n×n=9×9 and J=1.


## 4.1. Influence of Different Levels of Noise for Singular Vectors

Two matrices,Yα and Yβ, are chosen arbitrarily from the noisy matrix-set M~, and yα=vecYα and yβ=vecYβ. Thus, yα=xα+nα and yβ=xβ+nβ, where x denotes the noise-free vector and n is white noise. The squared Euclidean distance between the two vectors yα and yβ is(25)yα-yβ22=xα-xβ22+2xα-xβ,nα-nβ+nα-nβ22,where the symbol · is the inner-product operator.The vector difference,nα-nβ, is the Gaussian noise with covariance matrix 2σ2n2I, where symbol I denotes identity matrix. Based on the argument in [15], if the dimension of y is large, the norm squared, nα-nβ22, is concentrated around its mean 2σ2n2, with high probability. Similarly, the projection xα-xβ,nα-nβ of the noise on the deterministic vector xα-xβ is concentrated around 0 with high probability. Therefore, (17) can be interpreted as a translation-invariant procedure.On the other hand, as the noise level increases, the singular vectors ofYr∈M~ are no longer aligned with the original singular vectors of noise-free matrix Xr. As a result, the number of reliable singular vectors of Yr becomes smaller. For example, when Xr is corrupted by Gaussian noise at very small magnitude, all the singular vectors of Yr track the noise; when corrupted at a large magnitude, the number of the singular vectors tracking the noise becomes smaller. Therefore, at increasing noise level, the singular vectors of Yr become more and more unreliable and result in worse restored feature.


## 4.2. Shrunken Scales of Different Singular Values

LetYr,Xr be a paired samples and Σi be the ith largest singular value of Yr∈Rn×n. The coefficients vector copt, of size nJ×1, is trained on the paired samples and then is divided into n parts; ciT=ci1,ci2,…,ciJfori=1,2,…,n, where ci corresponds to the ith largest singular value of Yr. And thus, the estimate of the ith largest singular value of Yr is as follows:(26)Σ^i=∑j=1JcijΣi2∗j-1,fori=1,2,…,n.To observe the shrunken scales of the differentΣi visually, the shrinkage curves are shown in Figure 1. Test results show that (i) different singular values Σi have different shrunk scales; (ii) the larger the singular value, the smaller the shrunk scales; (iii) the shrunken scales for fixed index singular value are almost the same although Yr is corrupted by different levels of noise.Figure 1
Illustration of the shrinkage scale of theith largest singular value. From left to right and top to bottom, the ith inset is associated with the ith largest singular value. In each inset, the blue diagonal line is the ground singular value, and red line is its estimate. The parameters: patch-size n×n=9×9 and J=1.


## 5. Experiments

In this section, the OSC denoising performance is evaluated. In Section5.1, setting parameters are formulated. Section 5.2 introduces two metrics, the Peak Signal-to-Noise-Ratio (PSNR) and the Mean Structural Similarity Index Measure (MSSIM). Section 5.3 reports the experimental results, and the comparisons with other methods are discussed.


### 5.1. Setting Parameters

A total of three parameters are set in the proposed method, a scale factorγ of residual error, an integer J regulating odd polynomial order, and a patch-size n×n. The large γ may lead to oscillatory effects. For all noise levels, it is desirable for γ to be set to 0.1. To study the effect of parameters (J,n) on denoising, the coefficients vector copt was trained on the Lena image, and the trained copt is used to reduce noise for the Barbara image. The influences of J and n for PSNR results are shown in the left and right insets in Figure 2, respectively. By experience, these parameter values used in our experiments are shown in Table 1. During implementation process, because the behavior of coefficients copt on large entries may be completely distorted, thereby an input value outside the training range may be violently magnified. For this case, the algorithm checks whether the input is within the learning data range, and if not, it is not modified at all.Table 1
Parameters used in training and denoising procedures.
σ
5
10
15
20
25
30
35
40
t
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2
γ
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
J
4
4
4
4
4
4
4
4
4
4
4
4
4
4
4
4
n
6
5
7
6
8
6
9
6
10
7
11
7
12
8
12
9Figure 2
Illustration of influence of different parameters for denoising performance. The left inset illustrates the dependence of the parameterJ on PSNR results when the size of patches is fixed to 7∗7. The right inset illustrates the dependence of the patch sizes on PSNR results when the parameter J is fixed to 4.


### 5.2. Two Metrics

The PSNR measurement is based on pixel intensity errors between the noise-free and the restored images, given in decibels (dB). Higher PSNR means better denoising capability. LetX and X^ be the noise-free and the restored image, respectively, X denote the cardinality of X, and ·F denote the Frobenius norm; the calculation of PSNR is as follows:(27)PSNR=20log10⁡255X-1∑q=13Xq-X^qF2.The MSSIM measurement is the mean SSIM that yields mean value of structural similarities. MSSIM values are bounded in the range0,1, and the closer the value is to 1, the better the denoising scheme is implied. The calculation of SSIM involves paired blocks extracted from the noise-free image X and the restored X^, respectively. Let y1 and y2 denote paired blocks, μy1 and μy2 be the mean values of y1 and y2, respectively, σy1 and σy2 be respective variances, and σy1y2 be their covariance. Assuming c1 and c2 are two stabilization variables, the SSIM value can be calculated by(28)SSIMy1,y2=2μy1μy2+c12σy1y2+c2μy12+μy22+c1σy12+σy22+c2.


### 5.3. Experimental Results and Comparisons

A test set was built to evaluate the OSC method, which was a combination of two groups, denoted byΓ=Γ1,Γ2. Each group Γk contained noisy images, with eight different levels of white noise. The Γ1 group is associated with the sixteen noise-free images widely used, shown in Figure 3; the Γ2 group is associated with the sixteen noise-free biomedical images, shown in Figure 5. These noise-free images were taken from the CVG-Granada database.Figure 3
Sixteen widely used noise-free images associated with the groupΓ1.The corresponding training set was also built, which was a subset of the testΓ, denoted by Ψ=Ψ1,Ψ2. The group Ψ1 contained the noisy versions of the five noise-free images, the Einstein, Elaine, Barbara, Fingerprint, and the Man image; the group Ψ2 included the noisy versions of the six noise-free images, the C05c, Celulas, Cromo2, Mr2, Mr032, and the Mr034 image. The OSC method was applied to the test set Γ. Using the parameter values shown in Algorithm 1, the coefficients copt were trained by Algorithm 1. Algorithm 2 used the same parameter values and the trained copt. The PSNR and MSSIM results associated with both groups are reported in Tables 2 and 3, respectively. Moreover, the visual results and zoom-in for three widely used images are shown in Figure 4 and are shown in Figure 6 for three biomedical images.Table 2
PSNR/SSIM comparisons with bilateral filter and BM3D.
Methods
Noise levelsσ
5
10
15
20
25
30
35
40
(1) Baboon 256×256
B. Filter
28.39/0.845
27.71/0.842
26.43/0.811
24.67/0.746
22.77/0.665
20.96/0.585
19.32/0.512
17.87/0.449
BM3D
34.74/0.964
29.62/0.899
26.91/0.830
25.13/0.761
23.84/0.693
22.89/0.633
22.16/0.580
21.56/0.530
OSC
37.84/0.942
33.95/0.846
32.49/0.762
31.26/0.695
30.41/0.650
29.72/0.627
29.09/0.621
28.59/0.639
(2) Cameraman 256×256
B. Filter
32.16/0.895
31.63/0.887
29.87/0.800
27.15/0.636
24.37/0.488
21.96/0.385
19.95/0.315
18.27/0.266
BM3D
38.29/0.962
34.18/0.932
31.91/0.901
30.48/0.875
29.45/0.854
28.64/0.837
27.93/0.822
27.18/0.806
OSC
37.02/0.963
31.94/0.892
29.72/0.831
27.88/0.785
26.62/0.748
25.59/0.720
24.81/0.690
24.16/0.654
(3) Einstein 256×256
B. Filter
31.48/0.802
31.19/0.804
29.73/0.739
27.17/0.598
24.45/0.452
22.05/0.340
20.04/0.261
18.36/0.206
BM3D
37.49/0.934
34.19/0.871
32.55/0.832
31.43/0.803
30.56/0.778
29.84/0.757
29.10/0.736
28.43/0.717
OSC
37.71/0.981
33.51/0.935
31.78/0.890
30.35/0.846
29.39/0.814
28.62/0.776
28.02/0.747
27.57/0.710
(4) Elaine 256×256
B. Filter
31.16/0.825
30.63/0.822
29.05/0.767
26.63/0.651
24.10/0.522
21.83/0.412
19.90/0.330
18.26/0.268
BM3D
36.32/0.940
32.60/0.858
31.03/0.810
30.05/0.782
29.30/0.761
28.69/0.743
28.10/0.726
27.51/0.710
OSC
37.32/0.991
32.77/0.952
30.70/0.910
29.05/0.867
27.83/0.820
26.96/0.781
26.32/0.750
25.71/0.715
(5) F16 256×256
B. Filter
32.23/0.906
31.39/0.894
29.42/0.810
26.71/0.656
24.05/0.513
21.74/0.411
19.81/0.339
18.18/0.287
BM3D
37.53/0.957
33.43/0.926
31.19/0.899
29.70/0.876
28.57/0.855
27.70/0.836
26.96/0.818
26.22/0.799
OSC
37.47/0.994
32.90/0.962
30.79/0.924
29.20/0.892
27.98/0.858
27.07/0.823
26.37/0.793
25.72/0.754
(6) Lena 256×256
B. Filter
31.98/0.893
31.03/0.878
29.06/0.800
26.48/0.661
23.93/0.526
21.68/0.421
19.78/0.345
18.17/0.288
BM3D
38.16/0.961
34.49/0.932
32.40/0.906
30.95/0.883
29.86/0.863
28.99/0.844
28.23/0.825
27.43/0.805
OSC
37.06/1.000
32.14/0.973
29.83/0.939
28.05/0.906
26.84/0.878
25.89/0.861
25.13/0.840
24.49/0.798
(7) Monarch 256×256
B. Filter
31.53/0.929
30.69/0.911
28.85/0.837
26.33/0.710
23.79/0.584
21.55/0.486
19.66/0.412
18.06/0.355
BM3D
38.21/0.975
34.12/0.956
31.86/0.936
30.35/0.918
29.25/0.900
28.36/0.882
27.58/0.865
26.72/0.845
OSC
37.06/0.998
31.97/0.997
29.60/0.980
27.69/0.961
26.29/0.939
25.21/0.916
24.33/0.890
23.65/0.859
(8) House 256×256
B. Filter
33.76/0.878
32.88/0.870
30.59/0.777
27.51/0.602
24.58/0.442
22.10/0.331
20.06/0.259
18.37/0.210
BM3D
39.83/0.957
36.71/0.922
34.94/0.891
33.77/0.873
32.86/0.859
32.09/0.848
31.38/0.837
30.65/0.826
OSC
39.04/0.950
35.25/0.868
33.44/0.790
31.87/0.713
30.78/0.660
29.95/0.619
29.27/0.585
28.67/0.555
(9) Balloon 256×256
B. Filter
35.11/0.918
34.12/0.897
31.46/0.780
28.02/0.579
24.90/0.402
22.32/0.283
20.21/0.209
18.48/0.160
BM3D
40.82/0.960
37.84/0.937
36.03/0.922
34.60/0.906
33.48/0.892
32.57/0.877
31.77/0.862
31.07/0.849
OSC
40.17/0.993
36.31/0.975
34.28/0.953
32.47/0.930
31.04/0.911
29.91/0.885
29.04/0.870
28.35/0.874
(10) Peppers 256×256
B. Filter
32.68/0.907
31.67/0.889
29.51/0.803
26.72/0.654
24.04/0.510
21.74/0.402
19.80/0.325
18.18/0.269
BM3D
38.12/0.956
34.68/0.928
32.70/0.907
31.29/0.887
30.16/0.868
29.28/0.850
28.52/0.834
27.70/0.816
OSC
37.43/0.970
32.81/0.904
30.74/0.838
29.05/0.786
27.84/0.734
26.91/0.700
26.19/0.664
25.57/0.625
(11) Barbara 512×512
B. Filter
31.28/0.886
30.29/0.872
28.47/0.799
26.12/0.667
23.75/0.537
21.62/0.436
19.78/0.361
18.20/0.304
BM3D
38.31/0.965
34.98/0.942
33.11/0.923
31.78/0.905
30.72/0.887
29.81/0.869
28.98/0.848
27.99/0.823
OSC
39.82/0.954
35.87/0.893
34.00/0.847
32.43/0.814
31.34/0.792
30.53/0.772
29.90/0.743
29.42/0.727
(12) Boat 512×512
B. Filter
31.06/0.819
30.61/0.820
29.11/0.759
26.72/0.629
24.18/0.494
21.90/0.388
19.97/0.311
18.32/0.255
BM3D
37.28/0.939
33.92/0.888
32.14/0.854
30.88/0.826
29.91/0.801
29.12/0.780
28.43/0.759
27.74/0.739
OSC
37.74/0.996
33.19/0.967
31.08/0.939
29.36/0.904
28.11/0.876
27.17/0.857
26.40/0.840
25.79/0.826
(13) Couple 512×512
B. Filter
30.67/0.825
30.16/0.823
28.70/0.763
26.44/0.641
24.02/0.512
21.81/0.408
19.91/0.330
18.29/0.272
BM3D
37.52/0.951
34.04/0.909
32.11/0.877
30.76/0.848
29.72/0.820
28.87/0.795
28.15/0.771
27.48/0.747
OSC
38.07/0.991
33.71/0.952
31.81/0.912
30.23/0.877
29.10/0.848
28.28/0.813
27.62/0.770
27.12/0.739
(14) Fingerprint 512×512
B. Filter
29.07/0.926
27.75/0.908
25.95/0.872
24.01/0.816
22.13/0.748
20.43/0.675
18.91/0.604
17.56/0.539
BM3D
36.51/0.988
32.46/0.969
30.28/0.949
28.81/0.930
27.70/0.911
26.83/0.894
26.09/0.877
25.30/0.858
OSC
38.05/0.998
33.27/0.998
31.00/0.989
29.16/0.974
27.83/0.961
26.79/0.943
25.95/0.928
25.26/0.912
(15) Hill 512×512
B. Filter
31.03/0.804
30.66/0.804
29.27/0.745
26.91/0.618
24.35/0.481
22.05/0.371
20.08/0.290
18.41/0.231
BM3D
37.13/0.943
33.62/0.883
31.86/0.839
30.72/0.804
29.85/0.775
29.16/0.750
28.56/0.728
27.99/0.707
OSC
37.79/0.964
33.66/0.888
31.97/0.817
30.60/0.758
29.68/0.712
29.00/0.689
28.46/0.661
27.98/0.653
(16) Man 512×512
B. Filter
31.40/0.843
30.93/0.836
29.38/0.765
26.91/0.626
24.31/0.485
21.99/0.376
20.02/0.298
18.37/0.241
BM3D
37.82/0.954
33.98/0.908
31.93/0.867
30.59/0.833
29.62/0.805
28.86/0.780
28.22/0.758
27.65/0.737
OSC
37.69/0.990
33.24/0.957
31.35/0.907
29.80/0.897
28.70/0.896
27.88/0.868
27.21/0.870
26.66/0.834Table 3
PSNR/SSIM results for sixteen biomedical images.
Images
σ
5
10
15
20
25
30
35
40
Algae
39.37/0.898
36.66/0.818
35.61/0.776
34.47/0.741
33.50/0.703
32.68/0.660
32.01/0.620
31.49/0.586
c05c
39.31/0.944
36.11/0.876
34.71/0.822
33.43/0.790
32.43/0.760
31.60/0.734
30.93/0.714
30.38/0.697
celulas
38.19/0.978
34.40/0.926
32.88/0.880
31.62/0.835
30.72/0.802
30.01/0.771
29.41/0.755
28.91/0.726
crm04280
40.44/0.981
36.81/0.942
35.31/0.905
34.00/0.867
32.98/0.841
32.17/0.818
31.50/0.796
30.95/0.769
crm05210
39.49/0.984
35.36/0.940
33.48/0.904
31.98/0.871
30.85/0.845
30.01/0.822
29.30/0.800
28.74/0.776
cromo2
39.78/0.973
35.37/0.916
33.27/0.869
31.52/0.827
30.21/0.793
29.21/0.748
28.44/0.713
27.82/0.667
Fluocel
41.00/0.953
36.85/0.882
34.90/0.816
33.26/0.758
32.05/0.720
31.08/0.698
30.32/0.680
29.67/0.651
Heart
39.56/0.975
35.87/0.937
34.30/0.903
32.85/0.869
31.83/0.841
31.05/0.814
30.46/0.780
29.97/0.748
mr2
38.86/0.973
35.08/0.923
33.42/0.876
31.89/0.831
30.69/0.796
29.72/0.765
28.90/0.734
28.18/0.708
mr030
37.74/0.998
32.89/0.969
30.63/0.936
28.84/0.912
27.56/0.885
26.60/0.860
25.84/0.829
25.22/0.801
mr031
37.89/0.987
33.13/0.939
30.97/0.893
29.23/0.845
27.98/0.801
27.04/0.768
26.32/0.738
25.73/0.720
mr032
37.98/0.984
33.37/0.931
31.30/0.881
29.63/0.843
28.41/0.802
27.49/0.769
26.75/0.736
26.13/0.706
mr033
37.99/0.987
33.42/0.942
31.40/0.897
29.75/0.851
28.56/0.808
27.65/0.772
26.93/0.745
26.35/0.725
mr034
38.02/0.988
33.48/0.941
31.50/0.894
29.84/0.847
28.62/0.804
27.68/0.771
26.96/0.735
26.38/0.703
Muscle
38.02/0.966
33.55/0.905
31.48/0.851
29.71/0.807
28.40/0.768
27.41/0.738
26.60/0.723
25.91/0.698
X-ray
38.26/0.955
34.15/0.870
32.39/0.793
30.90/0.729
29.80/0.676
29.00/0.635
28.36/0.600
27.86/0.565Figure 4
Visual results comparison. Rows are associated with denoised images and columns correspond to methods, where Columns 1–4 are the visual results and Columns 5–8 are their respective zoom-in.Figure 5
Sixteen noise-free biomedical images associated with the groupΓ2.Figure 6
Visual results for four biomedical images. Row 1 corresponds to noisy images corrupted by noise (σ=35). Row 2 corresponds to denoised images. Row 3 corresponds to zoomed-in views. Columns correspond to different images.To augment the performance evaluations, the OSC method was compared with Bilateral Filter [16] and BM3D [17]. For the Bilateral Filter there are three parameters needed to be set, containing a sliding window of size d×d, a geometric spread standard deviation σd, and a photometric spread standard deviation σr. The three parameters, d, σd, and σr, were set to 5, 5, and 0.1, respectively. The source codes of BM3D were taken from the original authors, and the parameter values used in the experiments were those recommended by the authors. The PSNR and MSSIM results from the two methods are reported in Table 2, and the visual results for three images are shown in Figure 4.In addition, the mean PSNR and mean MSSIM results for fixed noise are calculated for noisy images, the Bilateral Filter, BM3D, and the OSC method, respectively. The corresponding calculations are as follows:(29)PSNR-σ=1Γσ∑k∈ΓσPSNRσk,MSSIM-σ=1Γσ∑k∈ΓσMSSIMσk.In these two formulae,Γσ denotes all the images with the same level noise σ in the test set Γ, and thus Γσ=32; and PSNRσk and MSSIMσk denote the PSNR and MSSIM values of the kth image with the noise σ, respectively. For example, if the OSC method employs formula (29) to calculate, PSNR-5 denotes the arithmetical mean PSNR for all images with the noise σ=5, associated with the OSC method. The mean PSNR and MSSIM values for fixed noise for the four methods are plotted in Figure 7.Figure 7
The mean PSNR and SSIM results for fixed noise from four methods.Some observations and conclusions can be drawn from the quantitative measurements and visual results. Firstly, the OSC method achieved the desirable results. The OSC method achieved the same results as the BM3D method on average and significantly outperformed the noisy method and the Bilateral Filter by 8.49 dB and 5.26 dB on average, respectively, in the PSNR results. The OSC method achieved the same results as the BM3D method on average and significantly outperformed the noisy method and the Bilateral Filter by 0.417 and 0.275 on average, respectively, in the MSSIM results. Secondly, the OSC method has a strong capability to preserve detail. The OSC method reconstructed more image details from noisy images than the Bilateral Filter and was almost the same as the BM3D. As seen in Figure4, the fine-details of the moustache for the Baboon image were preserved, while the Bilateral Filter introduced excessive smoothing and resulted in visual feature blurring. Thirdly, in contrast to BM3D, the OSC method achieved decreasing results at increasing levels noise, as in Figure 7.


## 5.1. Setting Parameters

A total of three parameters are set in the proposed method, a scale factorγ of residual error, an integer J regulating odd polynomial order, and a patch-size n×n. The large γ may lead to oscillatory effects. For all noise levels, it is desirable for γ to be set to 0.1. To study the effect of parameters (J,n) on denoising, the coefficients vector copt was trained on the Lena image, and the trained copt is used to reduce noise for the Barbara image. The influences of J and n for PSNR results are shown in the left and right insets in Figure 2, respectively. By experience, these parameter values used in our experiments are shown in Table 1. During implementation process, because the behavior of coefficients copt on large entries may be completely distorted, thereby an input value outside the training range may be violently magnified. For this case, the algorithm checks whether the input is within the learning data range, and if not, it is not modified at all.Table 1
Parameters used in training and denoising procedures.
σ
5
10
15
20
25
30
35
40
t
1
2
1
2
1
2
1
2
1
2
1
2
1
2
1
2
γ
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
J
4
4
4
4
4
4
4
4
4
4
4
4
4
4
4
4
n
6
5
7
6
8
6
9
6
10
7
11
7
12
8
12
9Figure 2
Illustration of influence of different parameters for denoising performance. The left inset illustrates the dependence of the parameterJ on PSNR results when the size of patches is fixed to 7∗7. The right inset illustrates the dependence of the patch sizes on PSNR results when the parameter J is fixed to 4.


## 5.2. Two Metrics

The PSNR measurement is based on pixel intensity errors between the noise-free and the restored images, given in decibels (dB). Higher PSNR means better denoising capability. LetX and X^ be the noise-free and the restored image, respectively, X denote the cardinality of X, and ·F denote the Frobenius norm; the calculation of PSNR is as follows:(27)PSNR=20log10⁡255X-1∑q=13Xq-X^qF2.The MSSIM measurement is the mean SSIM that yields mean value of structural similarities. MSSIM values are bounded in the range0,1, and the closer the value is to 1, the better the denoising scheme is implied. The calculation of SSIM involves paired blocks extracted from the noise-free image X and the restored X^, respectively. Let y1 and y2 denote paired blocks, μy1 and μy2 be the mean values of y1 and y2, respectively, σy1 and σy2 be respective variances, and σy1y2 be their covariance. Assuming c1 and c2 are two stabilization variables, the SSIM value can be calculated by(28)SSIMy1,y2=2μy1μy2+c12σy1y2+c2μy12+μy22+c1σy12+σy22+c2.


## 5.3. Experimental Results and Comparisons

A test set was built to evaluate the OSC method, which was a combination of two groups, denoted byΓ=Γ1,Γ2. Each group Γk contained noisy images, with eight different levels of white noise. The Γ1 group is associated with the sixteen noise-free images widely used, shown in Figure 3; the Γ2 group is associated with the sixteen noise-free biomedical images, shown in Figure 5. These noise-free images were taken from the CVG-Granada database.Figure 3
Sixteen widely used noise-free images associated with the groupΓ1.The corresponding training set was also built, which was a subset of the testΓ, denoted by Ψ=Ψ1,Ψ2. The group Ψ1 contained the noisy versions of the five noise-free images, the Einstein, Elaine, Barbara, Fingerprint, and the Man image; the group Ψ2 included the noisy versions of the six noise-free images, the C05c, Celulas, Cromo2, Mr2, Mr032, and the Mr034 image. The OSC method was applied to the test set Γ. Using the parameter values shown in Algorithm 1, the coefficients copt were trained by Algorithm 1. Algorithm 2 used the same parameter values and the trained copt. The PSNR and MSSIM results associated with both groups are reported in Tables 2 and 3, respectively. Moreover, the visual results and zoom-in for three widely used images are shown in Figure 4 and are shown in Figure 6 for three biomedical images.Table 2
PSNR/SSIM comparisons with bilateral filter and BM3D.
Methods
Noise levelsσ
5
10
15
20
25
30
35
40
(1) Baboon 256×256
B. Filter
28.39/0.845
27.71/0.842
26.43/0.811
24.67/0.746
22.77/0.665
20.96/0.585
19.32/0.512
17.87/0.449
BM3D
34.74/0.964
29.62/0.899
26.91/0.830
25.13/0.761
23.84/0.693
22.89/0.633
22.16/0.580
21.56/0.530
OSC
37.84/0.942
33.95/0.846
32.49/0.762
31.26/0.695
30.41/0.650
29.72/0.627
29.09/0.621
28.59/0.639
(2) Cameraman 256×256
B. Filter
32.16/0.895
31.63/0.887
29.87/0.800
27.15/0.636
24.37/0.488
21.96/0.385
19.95/0.315
18.27/0.266
BM3D
38.29/0.962
34.18/0.932
31.91/0.901
30.48/0.875
29.45/0.854
28.64/0.837
27.93/0.822
27.18/0.806
OSC
37.02/0.963
31.94/0.892
29.72/0.831
27.88/0.785
26.62/0.748
25.59/0.720
24.81/0.690
24.16/0.654
(3) Einstein 256×256
B. Filter
31.48/0.802
31.19/0.804
29.73/0.739
27.17/0.598
24.45/0.452
22.05/0.340
20.04/0.261
18.36/0.206
BM3D
37.49/0.934
34.19/0.871
32.55/0.832
31.43/0.803
30.56/0.778
29.84/0.757
29.10/0.736
28.43/0.717
OSC
37.71/0.981
33.51/0.935
31.78/0.890
30.35/0.846
29.39/0.814
28.62/0.776
28.02/0.747
27.57/0.710
(4) Elaine 256×256
B. Filter
31.16/0.825
30.63/0.822
29.05/0.767
26.63/0.651
24.10/0.522
21.83/0.412
19.90/0.330
18.26/0.268
BM3D
36.32/0.940
32.60/0.858
31.03/0.810
30.05/0.782
29.30/0.761
28.69/0.743
28.10/0.726
27.51/0.710
OSC
37.32/0.991
32.77/0.952
30.70/0.910
29.05/0.867
27.83/0.820
26.96/0.781
26.32/0.750
25.71/0.715
(5) F16 256×256
B. Filter
32.23/0.906
31.39/0.894
29.42/0.810
26.71/0.656
24.05/0.513
21.74/0.411
19.81/0.339
18.18/0.287
BM3D
37.53/0.957
33.43/0.926
31.19/0.899
29.70/0.876
28.57/0.855
27.70/0.836
26.96/0.818
26.22/0.799
OSC
37.47/0.994
32.90/0.962
30.79/0.924
29.20/0.892
27.98/0.858
27.07/0.823
26.37/0.793
25.72/0.754
(6) Lena 256×256
B. Filter
31.98/0.893
31.03/0.878
29.06/0.800
26.48/0.661
23.93/0.526
21.68/0.421
19.78/0.345
18.17/0.288
BM3D
38.16/0.961
34.49/0.932
32.40/0.906
30.95/0.883
29.86/0.863
28.99/0.844
28.23/0.825
27.43/0.805
OSC
37.06/1.000
32.14/0.973
29.83/0.939
28.05/0.906
26.84/0.878
25.89/0.861
25.13/0.840
24.49/0.798
(7) Monarch 256×256
B. Filter
31.53/0.929
30.69/0.911
28.85/0.837
26.33/0.710
23.79/0.584
21.55/0.486
19.66/0.412
18.06/0.355
BM3D
38.21/0.975
34.12/0.956
31.86/0.936
30.35/0.918
29.25/0.900
28.36/0.882
27.58/0.865
26.72/0.845
OSC
37.06/0.998
31.97/0.997
29.60/0.980
27.69/0.961
26.29/0.939
25.21/0.916
24.33/0.890
23.65/0.859
(8) House 256×256
B. Filter
33.76/0.878
32.88/0.870
30.59/0.777
27.51/0.602
24.58/0.442
22.10/0.331
20.06/0.259
18.37/0.210
BM3D
39.83/0.957
36.71/0.922
34.94/0.891
33.77/0.873
32.86/0.859
32.09/0.848
31.38/0.837
30.65/0.826
OSC
39.04/0.950
35.25/0.868
33.44/0.790
31.87/0.713
30.78/0.660
29.95/0.619
29.27/0.585
28.67/0.555
(9) Balloon 256×256
B. Filter
35.11/0.918
34.12/0.897
31.46/0.780
28.02/0.579
24.90/0.402
22.32/0.283
20.21/0.209
18.48/0.160
BM3D
40.82/0.960
37.84/0.937
36.03/0.922
34.60/0.906
33.48/0.892
32.57/0.877
31.77/0.862
31.07/0.849
OSC
40.17/0.993
36.31/0.975
34.28/0.953
32.47/0.930
31.04/0.911
29.91/0.885
29.04/0.870
28.35/0.874
(10) Peppers 256×256
B. Filter
32.68/0.907
31.67/0.889
29.51/0.803
26.72/0.654
24.04/0.510
21.74/0.402
19.80/0.325
18.18/0.269
BM3D
38.12/0.956
34.68/0.928
32.70/0.907
31.29/0.887
30.16/0.868
29.28/0.850
28.52/0.834
27.70/0.816
OSC
37.43/0.970
32.81/0.904
30.74/0.838
29.05/0.786
27.84/0.734
26.91/0.700
26.19/0.664
25.57/0.625
(11) Barbara 512×512
B. Filter
31.28/0.886
30.29/0.872
28.47/0.799
26.12/0.667
23.75/0.537
21.62/0.436
19.78/0.361
18.20/0.304
BM3D
38.31/0.965
34.98/0.942
33.11/0.923
31.78/0.905
30.72/0.887
29.81/0.869
28.98/0.848
27.99/0.823
OSC
39.82/0.954
35.87/0.893
34.00/0.847
32.43/0.814
31.34/0.792
30.53/0.772
29.90/0.743
29.42/0.727
(12) Boat 512×512
B. Filter
31.06/0.819
30.61/0.820
29.11/0.759
26.72/0.629
24.18/0.494
21.90/0.388
19.97/0.311
18.32/0.255
BM3D
37.28/0.939
33.92/0.888
32.14/0.854
30.88/0.826
29.91/0.801
29.12/0.780
28.43/0.759
27.74/0.739
OSC
37.74/0.996
33.19/0.967
31.08/0.939
29.36/0.904
28.11/0.876
27.17/0.857
26.40/0.840
25.79/0.826
(13) Couple 512×512
B. Filter
30.67/0.825
30.16/0.823
28.70/0.763
26.44/0.641
24.02/0.512
21.81/0.408
19.91/0.330
18.29/0.272
BM3D
37.52/0.951
34.04/0.909
32.11/0.877
30.76/0.848
29.72/0.820
28.87/0.795
28.15/0.771
27.48/0.747
OSC
38.07/0.991
33.71/0.952
31.81/0.912
30.23/0.877
29.10/0.848
28.28/0.813
27.62/0.770
27.12/0.739
(14) Fingerprint 512×512
B. Filter
29.07/0.926
27.75/0.908
25.95/0.872
24.01/0.816
22.13/0.748
20.43/0.675
18.91/0.604
17.56/0.539
BM3D
36.51/0.988
32.46/0.969
30.28/0.949
28.81/0.930
27.70/0.911
26.83/0.894
26.09/0.877
25.30/0.858
OSC
38.05/0.998
33.27/0.998
31.00/0.989
29.16/0.974
27.83/0.961
26.79/0.943
25.95/0.928
25.26/0.912
(15) Hill 512×512
B. Filter
31.03/0.804
30.66/0.804
29.27/0.745
26.91/0.618
24.35/0.481
22.05/0.371
20.08/0.290
18.41/0.231
BM3D
37.13/0.943
33.62/0.883
31.86/0.839
30.72/0.804
29.85/0.775
29.16/0.750
28.56/0.728
27.99/0.707
OSC
37.79/0.964
33.66/0.888
31.97/0.817
30.60/0.758
29.68/0.712
29.00/0.689
28.46/0.661
27.98/0.653
(16) Man 512×512
B. Filter
31.40/0.843
30.93/0.836
29.38/0.765
26.91/0.626
24.31/0.485
21.99/0.376
20.02/0.298
18.37/0.241
BM3D
37.82/0.954
33.98/0.908
31.93/0.867
30.59/0.833
29.62/0.805
28.86/0.780
28.22/0.758
27.65/0.737
OSC
37.69/0.990
33.24/0.957
31.35/0.907
29.80/0.897
28.70/0.896
27.88/0.868
27.21/0.870
26.66/0.834Table 3
PSNR/SSIM results for sixteen biomedical images.
Images
σ
5
10
15
20
25
30
35
40
Algae
39.37/0.898
36.66/0.818
35.61/0.776
34.47/0.741
33.50/0.703
32.68/0.660
32.01/0.620
31.49/0.586
c05c
39.31/0.944
36.11/0.876
34.71/0.822
33.43/0.790
32.43/0.760
31.60/0.734
30.93/0.714
30.38/0.697
celulas
38.19/0.978
34.40/0.926
32.88/0.880
31.62/0.835
30.72/0.802
30.01/0.771
29.41/0.755
28.91/0.726
crm04280
40.44/0.981
36.81/0.942
35.31/0.905
34.00/0.867
32.98/0.841
32.17/0.818
31.50/0.796
30.95/0.769
crm05210
39.49/0.984
35.36/0.940
33.48/0.904
31.98/0.871
30.85/0.845
30.01/0.822
29.30/0.800
28.74/0.776
cromo2
39.78/0.973
35.37/0.916
33.27/0.869
31.52/0.827
30.21/0.793
29.21/0.748
28.44/0.713
27.82/0.667
Fluocel
41.00/0.953
36.85/0.882
34.90/0.816
33.26/0.758
32.05/0.720
31.08/0.698
30.32/0.680
29.67/0.651
Heart
39.56/0.975
35.87/0.937
34.30/0.903
32.85/0.869
31.83/0.841
31.05/0.814
30.46/0.780
29.97/0.748
mr2
38.86/0.973
35.08/0.923
33.42/0.876
31.89/0.831
30.69/0.796
29.72/0.765
28.90/0.734
28.18/0.708
mr030
37.74/0.998
32.89/0.969
30.63/0.936
28.84/0.912
27.56/0.885
26.60/0.860
25.84/0.829
25.22/0.801
mr031
37.89/0.987
33.13/0.939
30.97/0.893
29.23/0.845
27.98/0.801
27.04/0.768
26.32/0.738
25.73/0.720
mr032
37.98/0.984
33.37/0.931
31.30/0.881
29.63/0.843
28.41/0.802
27.49/0.769
26.75/0.736
26.13/0.706
mr033
37.99/0.987
33.42/0.942
31.40/0.897
29.75/0.851
28.56/0.808
27.65/0.772
26.93/0.745
26.35/0.725
mr034
38.02/0.988
33.48/0.941
31.50/0.894
29.84/0.847
28.62/0.804
27.68/0.771
26.96/0.735
26.38/0.703
Muscle
38.02/0.966
33.55/0.905
31.48/0.851
29.71/0.807
28.40/0.768
27.41/0.738
26.60/0.723
25.91/0.698
X-ray
38.26/0.955
34.15/0.870
32.39/0.793
30.90/0.729
29.80/0.676
29.00/0.635
28.36/0.600
27.86/0.565Figure 4
Visual results comparison. Rows are associated with denoised images and columns correspond to methods, where Columns 1–4 are the visual results and Columns 5–8 are their respective zoom-in.Figure 5
Sixteen noise-free biomedical images associated with the groupΓ2.Figure 6
Visual results for four biomedical images. Row 1 corresponds to noisy images corrupted by noise (σ=35). Row 2 corresponds to denoised images. Row 3 corresponds to zoomed-in views. Columns correspond to different images.To augment the performance evaluations, the OSC method was compared with Bilateral Filter [16] and BM3D [17]. For the Bilateral Filter there are three parameters needed to be set, containing a sliding window of size d×d, a geometric spread standard deviation σd, and a photometric spread standard deviation σr. The three parameters, d, σd, and σr, were set to 5, 5, and 0.1, respectively. The source codes of BM3D were taken from the original authors, and the parameter values used in the experiments were those recommended by the authors. The PSNR and MSSIM results from the two methods are reported in Table 2, and the visual results for three images are shown in Figure 4.In addition, the mean PSNR and mean MSSIM results for fixed noise are calculated for noisy images, the Bilateral Filter, BM3D, and the OSC method, respectively. The corresponding calculations are as follows:(29)PSNR-σ=1Γσ∑k∈ΓσPSNRσk,MSSIM-σ=1Γσ∑k∈ΓσMSSIMσk.In these two formulae,Γσ denotes all the images with the same level noise σ in the test set Γ, and thus Γσ=32; and PSNRσk and MSSIMσk denote the PSNR and MSSIM values of the kth image with the noise σ, respectively. For example, if the OSC method employs formula (29) to calculate, PSNR-5 denotes the arithmetical mean PSNR for all images with the noise σ=5, associated with the OSC method. The mean PSNR and MSSIM values for fixed noise for the four methods are plotted in Figure 7.Figure 7
The mean PSNR and SSIM results for fixed noise from four methods.Some observations and conclusions can be drawn from the quantitative measurements and visual results. Firstly, the OSC method achieved the desirable results. The OSC method achieved the same results as the BM3D method on average and significantly outperformed the noisy method and the Bilateral Filter by 8.49 dB and 5.26 dB on average, respectively, in the PSNR results. The OSC method achieved the same results as the BM3D method on average and significantly outperformed the noisy method and the Bilateral Filter by 0.417 and 0.275 on average, respectively, in the MSSIM results. Secondly, the OSC method has a strong capability to preserve detail. The OSC method reconstructed more image details from noisy images than the Bilateral Filter and was almost the same as the BM3D. As seen in Figure4, the fine-details of the moustache for the Baboon image were preserved, while the Bilateral Filter introduced excessive smoothing and resulted in visual feature blurring. Thirdly, in contrast to BM3D, the OSC method achieved decreasing results at increasing levels noise, as in Figure 7.


## 6. Conclusions

In this paper, a shrinkage curves optimization is proposed for weighted nuclear norm minimization. Based on the theory [7] that the weighted nuclear norm minimization can be solved by imposing a soft threshold operation on the singular values, the proposed optimization method designs a penalty function using the difference between a latent matrix and its observation and employs odd polynomials to shrink the singular values of the observation matrix. As a result, the coefficients of polynomial fully characterize the shrinkage operator. Furthermore, the Frobenius norm of the penalty function is converted into the corresponding spectral norm, so the parameter optimization can be easily solved by off-and-shelf plain least-squares.To validate the effectiveness, the proposed parameter optimization method is adapted to image denoising. In our experiments, a total of 256 noisy images were tested. They contained the noisy version of thirty-two noise-free images, corrupted by eight different levels of noise. Experimental results show that the proposed denoising method is effective in terms of the comparative results. The proposed method achieves better results than the Bilateral Filter; when the noisy standard deviation is less than 20, the proposed method slightly outperforms the BM3D method and is weaker than the BM3D method when the noisy standard deviation varies from 20 to 40.



---

*Source: 4965262-2017-04-30.xml*