A recommender system to restore images with impulse noise

TL;DR

This paper introduces a collaborative filtering recommender system for restoring images with impulse noise, utilizing a new color image representation and parameter tuning to improve image quality.

Contribution

The paper proposes a novel collaborative filtering approach for impulse noise removal using a new color image representation and parameter analysis.

Findings

Effective noise restoration demonstrated on a standard image database.

Parameter tuning impacts bias and variance, influencing performance.

Potential extension to image inpainting and super-resolution.

Abstract

We build a collaborative filtering recommender system to restore images with impulse noise for which the noisy pixels have been previously identified. We define this recommender system in terms of a new color image representation using three matrices that depend on the noise-free pixels of the image to restore, and two parameters: $k$, the number of features; and $\lambda$, the regularization factor. We perform experiments on a well known image database to test our algorithm and we provide image quality statistics for the results obtained. We discuss the roles of bias and variance in the performance of our algorithm as determined by the values of $k$ and $\lambda$, and provide guidance on how to choose the values of these parameters. Finally, we discuss the possibility of using our collaborative filtering recommender system to perform image inpainting and super-resolution.

Tables

Table 1. Table 1: Parameter optimization for image I23_06_4

Method	$k$	$\lambda$	PSNR (dB)	$1-\text{MSSIM}$
Empirical	352	11	45.69
	380	13		0.001928
L-curve	352	9	45.67
	380	10		0.002028
$\Sigma$-curve	352	9	45.67
	380	10		0.002028

Equations

$$y_{i,j}=\{\begin{array}{ll}{a}_{i,j},& \text{with probability }1-p\\ n_{i,j},& \text{with probability }p\end{array}$$

$$\mathbf{X}=\left[\begin{array}{cccc}{\mathbf{x}}^{(1)}& {\mathbf{x}}^{(2)}& \cdots & {\mathbf{x}}^{(m)}\end{array}\right]\text{and}\Theta =\left[\begin{array}{cccc}{\theta }^{(1)}& {\theta }^{(2)}& \cdots & {\theta }^{(n)}\end{array}\right],$$

$$\mathbf{Z}={\mathbf{X}}^{\text{T}}\Theta .$$

$$J(\mathbf{X},\Theta )=J\left({\mathbf{x}}^{(1)},\dots ,{\mathbf{x}}^{(m)},{\theta }^{(1)},\dots ,{\theta }^{(n)}\right)=\frac{1}{2}\sum _{(i,j):r_{i,j}=1}{\left({\mathbf{x}}^{(i)\text{T}}{\theta }^{(j)}-y_{i,j}\right)}^2,$$

$$\underset{\mathbf{X},\Theta }{\min }J(\mathbf{X},\Theta ).$$

$$y=d+x_k{\theta }_k\iff x_k{\theta }_k=c,\text{where }c=y-d.$$

$${\theta }_k=\frac{c}{x_k}\text{and}{\theta }_k^2+x_k^2$$

$$=\frac{x_k^4+c^2}{x_k^2}.$$

$$J_{\lambda }(\mathbf{X},\Theta )=J_{\lambda }\left({\mathbf{x}}^{(1)},\dots ,{\mathbf{x}}^{(m)},{\theta }^{(1)},\dots ,{\theta }^{(n)}\right)=$$

$$\frac{\lambda }{2}\left(\sum _{i=1}^m\Vert {\mathbf{x}}^{(i)}{\Vert }_2^2+\sum _{j=1}^n\Vert {\theta }^{(j)}{\Vert }_2^2\right),$$

$$\underset{\mathbf{X},\Theta }{\min }{J}_{\lambda }(\mathbf{X},\Theta ).$$

$$\mathbf{Y}=\mathbf{A}\circ \mathbf{R}+\mathbf{N}\circ \neg \mathbf{R},$$

$$\mathbf{C}=\mathbf{Y}\circ \mathbf{R}+\left({{\mathbf{X}}_{\mathbf{B},\mathbf{R},\lambda }}^{\text{T}}{\Theta }_{\mathbf{B},\mathbf{R},\lambda }+\mu \cdot 1^{\text{T}}\right)\circ \neg \mathbf{R},$$

$$\displaystyle E\left[\big{(}y-\tilde{f}\big{)}^{2}\right]$$

$$=E[y^2]+E[{\tilde{f}}^2]-2E[y\tilde{f}],$$

$$=\text{Var}(y)+E[y{]}^2+\text{Var}(\tilde{f})+E[\tilde{f}{]}^2-2E[y\tilde{f}],$$

$$=\text{Var}(f)+{\sigma }^2+\text{Var}(\tilde{f})+E[y{]}^2+E[\tilde{f}{]}^2-2E[y\tilde{f}],$$

$$\displaystyle=\text{Var}(f)+\sigma^{2}+\text{Var}(\tilde{f})+\big{(}E[\tilde{f}]-E[f]\big{)}^{2}+2E[f]E[\tilde{f}]-2E[y\tilde{f}],$$

$$\displaystyle=\text{Var}(f)+\sigma^{2}+\text{Var}(\tilde{f})+E[\tilde{f}-f]^{2}-2\big{(}E[y\tilde{f}]-E[y]E[\tilde{f}]\big{)},$$

$$=\text{Var}(f)+{\sigma }^2+\text{Var}(\tilde{f})+\text{Bias}(\tilde{f}{)}^2-2\text{Cov}(y,\tilde{f}),$$

$$=\text{Var}(f)+{\sigma }^2+\text{Var}(\tilde{f})+\text{Bias}(\tilde{f}{)}^2-2\text{Cov}(f,\tilde{f}),$$

$$E[f(x_0)]=\frac{1}{m}\sum _{j=1}^m{f}^{(j)}(x_0)=f(x_0),$$

$$\text{Var}(f(x_0))=E\left[{\left(f(x_0)-E[f(x_0)]\right)}^2\right]=\frac{1}{m}\sum _{j=1}^m{\left(f^{(j)}(x_0)-f(x_0)\right)}^2=0.$$

$$E\left[\big{(}y(x_{0})-\tilde{f}(x_{0})\big{)}^{2}\right]=\text{Var}\big{(}\tilde{f}(x_{0})\big{)}+\text{Bias}\big{(}\tilde{f}(x_{0})\big{)}^{2}+\sigma^{2},$$

$$\text{Var}(\tilde{f}(x_0))=\underset{m\to \infty }{\lim }\frac{1}{m}\sum _{j=1}^m{\left({\tilde{f}}^{(j)}(x_0)-\frac{1}{m}\sum _{i=1}^m{\tilde{f}}^{(i)}(x_0)\right)}^2,$$

$$\text{Bias}(\tilde{f}(x_0))=\underset{m\to \infty }{\lim }\frac{1}{m}\sum _{j=1}^m\left({\tilde{f}}^{(j)}(x_0)-f(x_0)\right).$$

$$xm+xn$$

$$\iff x$$

$$\sqrt{\Vert {\mathbf{X}}_{\mathbf{B},\mathbf{R},\lambda }{\Vert }_F^2+\Vert {\Theta }_{\mathbf{B},\mathbf{R},\lambda }{\Vert }_F^2},$$

$${\Vert \left(\mathbf{Y}-\left({{\mathbf{X}}_{\mathbf{B},\mathbf{R},\lambda }}^{\text{T}}{\Theta }_{\mathbf{B},\mathbf{R},\lambda }+\mu \cdot 1^{\text{T}}\right)\right)\circ \mathbf{R}\Vert }_F,$$

$$\mathbf{D}={\mathbf{X}}^{\text{T}}\Theta +\mathbf{M}.$$

Full text

1 Introduction

Impulse noise in a grayscale image is characterized by the appearance of random single pixel errors on the image. This type of error occurs when there are malfunctioning pixel elements in the camera sensors, faulty memory locations, or timing errors in analog-to-digital conversion [7, 13].

We shall model impulse noise in the following way. Let $a_{i,j}$ and $y_{i,j}$ be the pixel values at position $(i,j)$ in the original and noisy images, respectively. Let $[s_{\min},s_{\max}]$ be the dynamic range of the image, i.e., the set of possible pixel values. Let $p$ be a number, called the noise ratio, such that $p\in[0,1]$, then

[TABLE]

where $n_{i,j}\neq a_{i,j}$ and $n_{i,j}\in[s_{\min},s_{\max}]$. If $n_{i,j}=s_{\min}\text{ or }s_{\max}$, we say that the image has salt-and-pepper noise [2]; and if $n_{i,j}$ is a (feasible) value chosen at random from $\mathcal{U}(s_{\min},s_{\max})$, we say that the image has random-valued impulse noise [5].

In color images, impulse noise is modeled by considering its presence as defined above in any of its color components taken individually as grayscale images.

There are results to detect and correct impulse noise in images, see for example the various techniques described in [5, 13]. Current approaches divide the task in two steps, firstly by identifying the likely pixels that have been corrupted by impulse noise, and secondly by applying a local filter to restore their values, see for example [4, 6]. Here we will concern ourselves with the restoration of color images for which corrupted pixels have already been “detected” by an impulse noise detection algorithm, say like ROLD [5], a detection statistic for random-valued impulse noise. In particular, for our study, we will use the TID2008 image database [14], which provides reference images and corresponding random-valued impulse noise counterpart examples. Unlike in the works referenced above, our goal will be to restore the noisy images with a recommender system which will substitute the pixels labeled as noise with a suggested value for each of them, as opposed to applying a local filter in each case.

This work is divided in the following way. In section 2 we introduce recommender systems in general and the collaborative filtering system that we propose in particular. We then lay out the experimental methodology that we will use to test our collaborative filtering system in section 3 and report on the results obtained in section 4. We follow with a discussion in section 5, where we talk about the bias vs variance trade-off, see section 5.1; how to choose the two main parameters of our recommender system, the number of features $k$ and the regularization factor $\lambda$ that control that trade-off, see section 5.2; and a brief note on the new family of image decompositions that our recommender system is based upon and the implications for the system performance, see section 5.3. Finally, we present our conclusions and future work directions in section 6.

2 A collaborative filtering recommender system

Recommender systems are algorithms and techniques that provide suggestions for their users on a given topic or class of objects. For example, a recommender system could suggest what items to buy, what music to hear, what movie to watch, or what news to read [15]. In our case, we are going to implement a recommender system that will suggest what value to give to a missing or corrupted pixel in an image, a novel approach.

The particular flavor of recommender system that we are going to use is a collaborative filtering algorithm [15]. In this modality, “systems require recommendation seekers to express preferences by rating a dozen or two items, thus merging the roles of recommendation seeker and preference provider. These systems focus on algorithms for matching people based on their preferences and weighting the interests of people with similar tastes to produce a recommendation for the information seeker” [3]. For example, consider a recommender system where the items are movies and the users rank them in a scale from 1 to 5. The recommender system would then compute a ranking for a movie that you have not yet ranked by taking into consideration the rankings given to that movie by other users and how similar your tastes are to the tastes of those users. Then, if the ranking is above certain threshold, the system would recommend that movie to you [12].

In our case, we are going to think of the rows in an image as the movies in the example above, and its columns as the users. Then our recommender collaborative filtering algorithm will provide a ranking for a missing pixel, which in this case will translate into its suggested value.

But why should a system like this work? Imagine that two people with very similar tastes in wine decide to rank 50 wines on a predetermined scale. Now suppose that the first person has scored all 50 wines and that the second one has done so for the first 49 wines in the list. If their rankings are very similar for the first 49 wines, we would expect that the second person would give a very similar ranking to the ranking that the first person gave to the 50th bottle of wine. That is intuitively why we expect the collaborative filtering algorithm to work.

2.1 Collaborative filtering cost function

Consider a grayscale image $\mathbf{Y}$ of size $m\times n$ with integer entries $y_{i,j}\in[0,2^{b}-1]$, where $b\in\mathbb{N}$ is the image’s bit depth. Typically $b=8$. In this case $[0,2^{b}-1]$ is what we called the dynamic range of the image in section 1. Consider a matrix $\mathbf{R}$ also of size $m\times n$ but with entries $r_{i,j}\in\{0,1\}$, where $r_{i,j}=1$ if we know the pixel value $y_{i,j}$ at position $(i,j)$ to be correct and 0 otherwise, i.e., the pixel value is missing or has been labeled as noise by our impulse noise detector. We will call such matrix $\mathbf{R}$ a mask for $\mathbf{Y}$.

Now imagine that for each row $i$ of $\mathbf{Y}$ we associate a vector $\mathbf{x}^{(i)}\in\mathbb{R}^{k}$, where $k\in\mathbb{N}$ is a fixed number of features that we are going to associate with that row. Similarly we associate a vector $\theta^{(j)}\in\mathbb{R}^{k}$ to column $j$ of $\mathbf{Y}$. Hence we will have $m$ vectors $\{\mathbf{x}^{(1)},\mathbf{x}^{(2)},\ldots,\mathbf{x}^{(m)}\}\subset\mathbb{R}^{k}$ and $n$ vectors $\{\theta^{(1)},\theta^{(2)},\ldots,\theta^{(n)}\}\subset\mathbb{R}^{k}$ associated with our image $\mathbf{Y}$. The collaborative filtering algorithm then considers these two sets of vectors and will model the prediction of the pixel value at position $(i,j)$ by computing $z_{i,j}=\mathbf{x}^{(i)\text{T}}\theta^{(j)}$, i.e., the inner product of $\mathbf{x}^{(i)}$ and $\theta^{(j)}$. In other words, if we form the matrices

[TABLE]

then our prediction (recommendation) matrix is

[TABLE]

Define the collaborative filtering cost function for $\mathbf{Y}$ as

[TABLE]

where implicitly we have provided a mask $\mathbf{R}$ for $\mathbf{Y}$ as well.

To link matrices $\mathbf{X}$ and $\Theta$ meaningfully to $\mathbf{Y}$, we require that they solve the minimization problem

[TABLE]

To see why this requirement is meaningful, let $\mathbf{R}$ be a mask for $\mathbf{Y}$ whose entries are all ones and assume that $J(\mathbf{X},\Theta)=0$, then it is easy to see that $y_{i,j}=\mathbf{x}^{(i)\text{T}}\theta^{(j)}$ for all $(i,j)$. That is, $\mathbf{X}$ and $\Theta$ offer a perfect reconstruction (factorization) of image $\mathbf{Y}$. In other words $\mathbf{Y}=\mathbf{X}^{\text{T}}\Theta$.

2.2 Regularization and normalization

In section 2.1 we defined the collaborative cost function $J$ for an image $\mathbf{Y}$ that measures how good two matrices $\mathbf{X}$ and $\Theta$ approximate the known good (non noisy) values of the pixels of $\mathbf{Y}$ when the approximation is given by $\mathbf{X}^{\text{T}}\Theta$.

Two remarks are in order. Assume that $\mathbf{x},\theta\in\mathbb{R}^{k}$ exactly approximate the pixel value $y$ at position $(i,j)$ in image $\mathbf{Y}$. In that case, $y=\mathbf{x}^{\text{T}}\theta$. Also assume that $k>1$ and write the inner product of these two vectors as $\mathbf{x}^{\text{T}}\theta=d+x_{k}\theta_{k}$, where we focus our attention on the product of their last corresponding $k$-th vector components. Then,

[TABLE]

Without loss of generality, assume that $c\neq 0$. Then

[TABLE]

Observe that the real function of real variable $f(x)=\frac{x^{4}+c^{2}}{x^{2}}$ has the following properties:

$\lim_{x\rightarrow\pm\infty}f(x)=+\infty,$ and $\lim_{x\rightarrow 0}f(x)=+\infty$. 2. 2.

$f$ has two minima. One at $x_{1}=\sqrt{|c|}$ and another one at $x_{2}=-\sqrt{|c|}$.

The first remark then is that without any further restrictions, we can find an infinite number of vector pairs $(\mathbf{x},\theta)$ for which their combined $\ell_{2}$-norm is arbitrarily large and for which $y=\mathbf{x}^{\text{T}}\theta$; and that there are at most two vector pairs $(\mathbf{x},\theta)$ with minimal combined $\ell_{2}$-norm satisfying $y=\mathbf{x}^{\text{T}}\theta$. This prompts us to regularize the cost function $J$ defined in equation (4) and propose instead the regularized collaborative filtering cost function for $\mathbf{Y}$ as

[TABLE]

where $\lambda\geq 0$ is the regularization factor. Consequently, the minimization problem (5) is equivalently transformed into the regularized minimization problem for $\mathbf{Y}$

[TABLE]

The introduction of regularization reduces the search space for a solution of problem (9), and also addresses the issue of model overfitting, of which we will talk in section 5.

The second remark comes from the following observation. Assume that image $\mathbf{Y}\in\mathbb{R}^{m\times n}$ is not trivially zero and that a whole column $j$ of pixels is tagged as noise. If we set $\lambda>0$, a solution to minimization problem (9) will necessarily return a matrix $\Theta$ for which $\theta^{(j)}=\mathbf{0}$, the zero vector. Then for any row $i$, the recommender system will return a value of $a=0$ as the suggested value for the pixel at position $(i,j)$, given that $\mathbf{x}^{(i)\text{T}}\theta^{(j)}=0$ in this case. This would be a very poor recommendation since at least there is a row $i$ for which the average value of the noise free pixels in that row is not zero, and assigning that average to the missing pixel would be a more meaningful suggestion than zero. To have an intuition as to why that is the case, think back to the scenario where users rank movies. For a new user that hasn’t ranked any movies yet and wants to know what other users think of a particular movie, the recommender system would do a better job if it offered the average recommendation that other users have given to that movie than offering a ranking value of zero.

One way to address this problem would be to normalize the original image $\mathbf{Y}$ by subtracting from each row the average value of the noise-free pixel values in that row, and then solve problem (9) for this new matrix. For a noisy pixel at position $(i,j)$ we would then recommend the value obtained by computing the corresponding inner product and adding to it the average value aforementioned. We describe this procedure in detail in the following section.

3 Experimental setup

As mentioned in section 1, we use the photographs in the TID2008 image database [14] to conduct our experiments. This database consists of 25 reference images, 17 types of distortions for each reference image, 4 different levels of each type of distortion. From these, we work with the images affected with random-valued impulse noise. The generative model for this distortion is described by equation (1). The 4 levels of distortion correspond to a noise ratio of about $p,2p,4p$ and $8p$ for levels 1 through 4, respectively, with $p\approx 0.0085$.

Each color image in our database is a 24-bit depth color image of size $384\times 512$ pixels which we manipulate into a single matrix of size $384\times 1536$ where each of the 8-bit red, green, and blue color channels are concatenated column wise, see figure 1 for examples. From this point onward, when we mention a color image we will think of it this way, unless we say otherwise, i.e., as a single matrix where all three color channels have been combined column wise.

With this treatment for color images in place, given a noisy image $\mathbf{Y}$ in our dataset, like the one depicted in figure 1(b), the experiment consists in following the steps described in the last paragraph of section 2. In particular, we have to first obtain a mask $\mathbf{R}$ for $\mathbf{Y}$ that identifies the pixels that have been corrupted with random-valued impulse noise. Here is where we can use, for example, the ROLD statistic [5] to mark pixels that are to be considered impulse noise.

The mask $\mathbf{R}=(r_{i,j})\in\{0,1\}^{m\times 3n}$, with $m=384$ and $n=512$ in this case, is a matrix such that $r_{i,j}=1$ if the pixel at position $(i,j)$ is considered noise-free, and zero otherwise, see figure 2.

In our case we count with the noise-free reference image $\mathbf{A}$, see figure 1(a), and this means that we can obtain a perfect mask $\mathbf{R}$ for $\mathbf{Y}$ in the sense that it will identify correctly all pixels that are noise, and only those. This is not necessarily the case when we don’t have a reference image that represents the grown truth. The relationships between the noise-free image $\mathbf{A}$, its related noisy image $\mathbf{Y}$, the mask $\mathbf{R}$ for $\mathbf{Y}$ (with respect to $\mathbf{A}$), and the noise matrix $\mathbf{N}$ is given by the following equation

[TABLE]

where $\circ$ is the Hadamard product of matrices (entry-wise multiplication), and $\neg$ is the negation operator that changes a 1 for a 0 and vice versa. Note that the noise matrix $\mathbf{N}$ is assumed to satisfy the conditions of equation (1), and it can be made unique if we impose the extra condition that it have the least number of nonzero entries.

Next we compute $\mathbf{B}$, the normalized version of the noisy image $\mathbf{Y}$, as described in the last paragraph of section 2, namely, let $\mu\in\mathbb{R}^{m}$ be the vector with entries $\mu_{i}=\frac{1}{\#\{j\,:\,r_{i,j}=1\}}\sum_{j\,:\,r_{i,j}=1}y_{i,j}$, then define $\mathbf{B}\in\mathbb{R}^{m\times 3n}$ to be the matrix with entries $b_{i,j}=y_{i,j}-\mu_{i}$ if $r_{i,j}=1$, and 0 otherwise. We now proceed to obtain a solution to problem (9) for the normalized color image $\mathbf{B}$. For this, we have to select a regularization factor $\lambda$ as well. We chose $\lambda=0$ and $\lambda=10$ to compare non-regularized and regularized solutions.

Since this is a minimization problem, a method such as gradient descent or any flavor of conjugate gradient can be used. Either of these two approaches requires choosing an initial value for matrices $\mathbf{X}\in\mathbb{R}^{k\times m}$ and $\Theta\in\mathbb{R}^{k\times 3n}$, see equation (2), from which the minimization algorithm then iterates to pursue a minima of the regularized collaborative filtering cost function $J_{\lambda}(\mathbf{X},\Theta)$, see equation (2.2). We chose both matrices to have entries taken from the normal distribution with zero mean and variance one, $N(0,1)$, as initial values. The choice of the number of features $k$ will be discussed in section 5. For our experiments we selected $k=1,8,308,352,380,384$.

With this experimental setup, let $(\mathbf{X_{\mathbf{B},\mathbf{R},\lambda},\Theta_{\mathbf{B},\mathbf{R},\lambda}})={\arg\min}_{\mathbf{X},\Theta}J_{\lambda}(\mathbf{X},\Theta)$ be a solution pair for the regularized minimization problem for $\mathbf{B}$. We then propose image

[TABLE]

as the reconstruction from the noisy image $\mathbf{Y}$, from mask $\mathbf{R}$ for $\mathbf{Y}$, and regularization factor $\lambda$, where $\mathbf{1}\in\mathbb{R}^{3n}$ is the vector with all ones, and $n=512$ in this case. Recall that $\mu\in\mathbb{R}^{m}$ and hence $\mathbf{M}=\mu\cdot\mathbf{1}^{\text{T}}\in\mathbb{R}^{m\times 3n}$. Basically, if $\mathbf{R}$ is a perfect mask, as is the case in our experiments, we have that $\mathbf{Y}\circ\mathbf{R}=\mathbf{A}\circ\mathbf{R}$, and comparing equations (10) and (11) gives that we are substituting the noise matrix $\mathbf{N}$ with the sum of the prediction matrix $\mathbf{Z}_{\mathbf{B},\mathbf{R},\lambda}={\mathbf{X}_{\mathbf{B},\mathbf{R},\lambda}}^{\text{T}}\Theta_{\mathbf{B},\mathbf{R},\lambda}$, see equation (3), with the normalization term $\mathbf{M}$, to restore image $\mathbf{A}$ from its noisy version $\mathbf{Y}$. Observe that, by design, $\|(\mathbf{A}-(\mathbf{Z}_{\mathbf{B},\mathbf{R},\lambda}+\mathbf{M}))\circ\mathbf{R}\|_{2}$ is small, which justifies why we would make that substitution of $\mathbf{N}$. We can therefore call the matrix $\mathbf{D}_{\mathbf{B},\mathbf{R},\lambda}=\mathbf{Z}_{\mathbf{B},\mathbf{R},\lambda}+\mathbf{M}={\mathbf{X}_{\mathbf{B},\mathbf{R},\lambda}}^{\text{T}}\Theta_{\mathbf{B},\mathbf{R},\lambda}+\mathbf{M}$ the denoising matrix.

4 Results

In this section we present the results that we obtain when we follow the experimental protocol, described in section 3, for all of the reference photographs in the TID2008 image database [14] that are affected with each of the four impulse noise levels available therein. The noise levels are controlled by the noise ratio parameter $p$, see equation (1), which has in this case, on average, the values of 0.0085, 0.0169, 0.0339, and 0.0678.

To measure the quality of the reconstructions we obtain, we utilize the peak signal-to-noise ratio (PSNR) [2, 17], measured in decibels (dB), and the mean structural similarity index (MSSIM), a number from 0 to 1, where 1 corresponds to two identical images. This image quality statistic is designed to correlate with the way humans perceive image quality [18].

The first set of results we present can be seen in figure 3. For each image in the TID2008 database affected with impulse noise we plot its PSNR and MSSIM image quality metrics. Then we process each image applying to it the experimental protocol described in section 3, and plot on the same figure both image quality metrics for the reconstruction. We can observe that for all images and all four levels of noise ratio $p$ the improvements in quality are very satisfactory, especially as measured by MSSIM. The choice of the number of features $k=352$ and regularization factor $\lambda=11$ in our experiments is justified in section 5.

The next set of results that we present pertain to regularization, see sections 2.2 and 5 for definitions and justification of its use. Figures 4 and 5 show the effects of regularization on PSNR and MSSIM for $k=352$ number of features and a regularization factor of $\lambda=11$. The choices of these two numbers is discussed in section 5.

The results confirm that regularizing the solutions is a good idea that yields better reconstruction results, if marginally in some instances, but dramatic in others. The figures show the results for all 25 images contained in the TID2008 database.

Finally, we show in figure 6 a selection of images affected by random-valued impulse noise for a noise ratio of $p=0.0678$, their reconstruction with our methodology, and the originals from which the noisy versions were generated for comparison purposes.

5 Discussion

5.1 Bias vs variance trade-off

When designing or choosing a model to represent a data set, we are interested in quantifying to what degree the model appropriately captures its essence. To help us quantify the goodness of a model we will show how the concepts of bias and variance play a role.

Suppose that we conduct the following experiment. Given a set $\{x_{1},x_{2},\ldots,x_{n}\}\subseteq X$ we associate to each $x_{i}$ a sample value $y_{i}\in\mathbb{R}$, forming a set of observations $\{(x_{1},y_{1}),(x_{2},y_{2}),\ldots,(x_{n},y_{n})\}$, which we call the training set. Assume that $y_{i}=f(x_{i})+\epsilon$, where $\epsilon$ is a random variable with zero mean and variance $\sigma^{2}>0$, which we call noise, and $f:X\rightarrow\mathbb{R}$ is an unknown deterministic function which we wish to determine. We typically try to do this by training a function $\tilde{f}:X\rightarrow\mathbb{R}$ over the training set that minimizes the mean square error of the approximation by $\tilde{f}$ to the sample values $\{y_{i}\}$ in it. Ideally, we would like to come up with a function $\tilde{f}$ that would also minimize the mean square error for sample values of other possible points $x_{0}$ for which we haven’t gotten any data yet, i.e., $x_{0}\in X\setminus\{x_{i}\}$, so that we could make good predictions $\tilde{f}(x_{0})$ with our model $\tilde{f}$. Essentially, we want to minimize the mean square error over $X$, which we decompose as follows

[TABLE]

where $\text{Bias}(\tilde{f})=E[\tilde{f}-f]$. Now suppose, as mentioned above, that we want to know what the expected test mean square error $E[(y(x_{0})-\tilde{f}(x_{0}))^{2}]$ will be for a point $x_{0}\in X\setminus\{x_{1},x_{2},\ldots,x_{n}\}$. This expectation is the average of the square error for all possible functions $\tilde{f}$ derived from all possible training sets when we perform the experiment described above repeatedly. This identifies $\tilde{f}(x_{0})$ to a random variable.

In the finite case, let $\{(x_{1}^{(j)},y_{1}^{(j)}),(x_{2}^{(j)},y_{2}^{(j)}),\ldots,(x_{n}^{(j)},y_{n}^{(j)})\}$ be the $j$-th training set that is obtained after performing the experiment $j$ times. Assume that we have $m$ training sets, then

[TABLE]

since $f^{(j)}(x_{0})=f(x_{0})$ because $f$ is deterministic, i.e., a regular function. From this,

[TABLE]

Similarly, we can prove that $\text{Cov}(f(x_{0}),\tilde{f}(x_{0}))=0$, and together with equations (13) and (14), equation (5.1) becomes

[TABLE]

with

[TABLE]

and

[TABLE]

The overall expected test mean square error can be computed by averaging E\big{[}\big{(}y(x_{0})-\tilde{f}(x_{0})\big{)}^{2}\big{]} over all possible values of $x_{0}$ in the test set [10].

Equation (15) says that in order to minimize the expected test mean squared error, we must find a function $\tilde{f}$ that simultaneously yields low variance and low bias. Also note that since $\text{Var}(\tilde{f}(x_{0}))\geq 0$ and $\text{Bias}(\tilde{f}(x_{0}))^{2}\geq 0$, this error will always be at least as large as the variance $\sigma^{2}$ of the noise component $\epsilon$ in our model.

The variance speaks to how much $\tilde{f}$ would change if we changed our training data set, and bias refers to the error introduced by approximating a real-life problem by a much simpler model. Often these two quantities are mutually exclusive and therefore we speak of the bias-variance trade-off [10].

In our recommender system model, the bias-variance trade-off is controlled by the number of features $k$, and the regularization factor $\lambda$. In figure 7 we show the error when reconstructing image I23 from noisy image I23_06_4 for values of $k$ equal to 1, 8, and 352 using a variety of values of $\lambda$ between 0 and 80. How to determine a good value of $k$? And how about $\lambda$?

5.2 How to choose $k$ and $\lambda$ without a ground truth reference

For an $m\times n$ grayscale image we need $mn$ values to describe each pixel when no compression is used. For our recommender system model, we need to find two matrices $\mathbf{X}$ and $\Theta$ to propose a recommendation. If $\mathbf{X}\in\mathbb{R}^{k\times m}$ and $\Theta\in\mathbb{R}^{k\times n}$, see equation (2), then we need $km+kn$ numbers to represent those matrices. It is therefore reasonable to assume that we would need a value of $k$ such that $km+kn$ is at least as big as $mn$ to represent the image with fidelity. We consequently try to find the smallest $k\in\mathbb{N}$ such that $km+kn\geq mn$. Let $x\in\mathbb{R}$ and $m,n\in\mathbb{N}$, then

[TABLE]

hence if $k_{min}\in\mathbb{N}$ is the smallest natural number such that $k_{min}\geq\frac{mn}{m+n}$, we would have at least as many numbers necessary to define our recommender system as there are numbers describing the image that we are trying to reconstruct. In our case, since all the images in our database are $384\times 512$ color images, we have that $m=384$ and $n=3\times 512=1536$, see section 3. This gives a minimum $k_{min}=308$ for a system that, with the criterion outlined above, will have enough variance to capture any of the pictures in our dataset.

On the other extreme, since $\frac{mn}{m+n}<\min\{m,n\}$ when $m,n\neq 0$, if we set $k_{max}=\min\{m,n\}$ we should also have enough features in our recommender system to capture the images in the database with fidelity. In our case $k_{max}=\min\{384,1536\}=384$. Finally, through exploration, see figure 8, we find that $k=352$ is the number of features that maximizes the PSNR in our experiments with image I23 and noisy image I23_06_4 for a regularization factor of $\lambda=0$, and $k=380$ maximizes MSSIM under the same circumstances.

To explore the impact of the regularization factor $\lambda$ we performed a numerical exploration for a range of values of $\lambda$ between 0 and 80. When reconstructing image I23 from image I23_06_4, we found that, on average, for $k=352$ with $\lambda=11$ we obtain the largest PSNR, and for $k=380$ with $\lambda=13$ we get the smallest error as measured by $1-\text{MSSIM}$. For $k=308,384$ both image quality measures are not as good for their respective optimal values of $\lambda$. The experiments are summarized in figure 9. They reflect the trade-off between bias and variance that we talked about since the parameter $\lambda$ controls, to a degree, for the variance in the model. For the experiments with all other images in our database we chose a regularization factor $\lambda=11$ and we set $k=352$ for the number of features. We want to emphasize that we do not expect theses choices to produce optimal results for all images and all noise ratios $p$, but our experiments show that they produced very satisfactory results, as discussed in section 4.

It is important to notice that for a given number of features $k$ we were able to find by exploration the optimal value for the regularization factor $\lambda$ because we counted with both the original and distorted versions of the particular image we chose to fine tune these parameters. But what if we didn’t have that luxury, i.e., what if we only counted with the distorted version of the image to figure out a reasonable value for $\lambda$, as is the real world scenario? We propose in that case to use the L-curve method [8, 9] to select the value of $\lambda$. In figure 10(a) we plot the norm of the solution as given by,

[TABLE]

versus the norm of the residue computed as,

[TABLE]

both averaged from ten different random seed reconstructions of image I23 from I23_06_4.

We observe that the value of the regularization factor $\lambda$ at the “elbow” in the L-shaped curve is very close to the value of $\lambda$ for which we experimentally obtain the best reconstruction results, see table 1. That is, in the case of $k=352$ number of features, the L-shaped curve has an elbow at $\lambda=9$; whereas for $k=380$, the elbow appears at $\lambda=10$. These values of $\lambda$ are best seen when we plot the sum of the average residue and the average solution norms, and they correspond to the minima of what we call the $\Sigma$-curves in figure 10(b). What is remarkable of this approach is that we can find a very good approximation to the empirically obtained optimal value of $\lambda$ for a given $k$ via the L or $\Sigma$-curves without knowledge of the ground truth.

5.3 Image decomposition and error: many choices

Another topic to address is the image decomposition that we have proposed. In particular, the additive term $\mathbf{M}$ on the right hand side in the representation

[TABLE]

In our case we set $\mathbf{M}=\mu\cdot\mathbf{1}^{\text{T}}\in\mathbb{R}^{m\times 3n}$, with $m=384$ and $n=512$, where $\mu$ is the vector with entries $\mu_{i}=\frac{1}{\#\{j\,:\,r_{i,j}=1\}}\sum_{j\,:\,r_{i,j}=1}y_{i,j}$, and $\mathbf{1}\in\mathbb{R}^{3n}$ is the vector with all ones, see equation (11). Matrix $\mathbf{M}$ is the term that provides normalization in our recommender system, see section 2.2. This term averages the available pixel values in each row and proposes that average for any missing pixel in that row as a starting point, which is then further enhanced by the contribution coming from $\mathbf{X}^{\text{T}}\Theta$ for that particular pixel. This approach has the caveat that, in images, the value of a pixel may be more likely to be correlated with the values of nearby pixels, than with those that are farther away. This can be corrected by proposing different ways to build the normalizing matrix $\mathbf{M}$.

One approach could be to assign a weight to the value of a contributing pixel that takes into account its distance to the reference pixel, for example; or to just average the pixel values of the pixels in the same class that it belongs to, which could be defined in a variety of ways, such as belonging in the neighborhood of a $K$-means partition of the image, as another example, see figure 11.

Finally, we want to mention how to address the slight asymmetry on how the rows and columns of the image are treated. We can transpose the image, process it in the same way the original image was processed, revert to the original orientation, and finally add and average both representations. We do this for our test image I23 and its noisy representation I23_06_4 and obtain that, over five different random initializations, we obtain an average improvement of $\approx 2.68$ dB in PSNR, and $7.2860\times 10^{-4}$ units in MSSIM.

6 Conclusions and future work

We have proposed a collaborative filtering recommender system to restore images with impulse noise. This system assumes that we know the location of the pixels that have been corrupted by impulse noise and it uses this knowledge to build a recommendation for the values that those pixels should have. In the process we have proposed an image decomposition that consists in the product of two matrices plus the sum of a third one.

To test and calibrate this collaborative recommender system, we have used a well known image database. There are two parameters—$k$, the number of features; and $\lambda$, the regularization factor—which we have explored and optimized for our particular image data set to obtain the best image quality. We have measured the performance of our algorithm using two well known image quality metrics, PSNR and MSSIM. Our experiments show that this system performs very well at reconstructing the original images from their noisy counterparts provided we give a complete list of the distorted pixels in each of the originals.

We have provided a rationale on how to choose the parameters aforementioned, but this is an area of further exploration. Also, we have pointed out that the image decomposition proposed is not unique and we have hinted at other possible decompositions that fit our recommender system. Exploring these could lead to new insights, improvements, and applications of our algorithm.

As described above, the bulk of our work focused on reconstructing images with isolated pixels deemed to be noisy, or missing. What happens when we consider missing segments of an image? Figures 12(a), 12(b), and 12(c) show the performance of our system when presented with an image that has been excised of a $32\times 32$ square portion of its pixels. The reconstruction obtained compares remarkably well with the original and goes from a 29.21 dB mutilated image to a reconstruction at 41.81 dB. Figures 12(d) and 12(e) showcase an image whose original is not available that has been defaced with a grid of dots and for which our recommender system has provided a reconstruction that could be the original undisturbed image by visual inspection. These examples fall under the purview of the general techniques known as inpainting [1, 16], and although the reconstructions shown are far from the current state-of-the-art in inpainting, the fact that our algorithm recovers the geometry and texture simultaneously without having been designed implicitly to solve this holy grail of inpainting makes exploring this application of our algorithm worthy of pursuit. One can envision as well applications of our work to superresolution.

I would like to thank the support of the leadership at the Institute for Physical Science and Technology, without which this work would not have been possible. I also want to personally thank John J. Benedetto, Director of the Norbert Wiener Center in the Mathematics Department at the University of Maryland, College Park; and Mike Dellomo, Associate Director and Advisor of the Masters in Telecommunications Program in the Department of Electrical and Computer Engineering also at the University of Maryland, College Park; for their feedback and valuable input in the writing of this manuscript. Finally, I want to acknowledge and thank the support of ARO Grants W911NF-15-1-0112 and W911NF-17-1-0014.