• We want to denoise by replacing a pixel with an average of similar pixels.
• This method works because of the variance law in probability theory.
  • If 9 similar pixels are averaged, then the noise standard deviation of the average will be divided by 3.

How these similar pixels can be founded?

Motivation

For a pixel \(p = (x,y)\) of the image \(\bold{I}\), we compute the value of pixel \(p\) as follows:

$$\bold{I}[p] = \frac{1}{C(p)} \sum_{q \in B(p,r)} \bold{I}[q] \cdot w(p,q), \qquad C(p) = \sum_{q \in B(p,r)} w(p,q),$$

where \(B(p,r)\) is a search window centered at \(p\) with size \((2r+1)\times(2r+1)\) pixels. The value of search size, \(r\), is determined by the standard deviation of noise \(\sigma\).

Pixelwise Implementation

The weight \(w(p,q)\) is computed by using an exponential kernel

$$w(p,q) = e^{-\frac{\max(d^2(p,q)-2\sigma^2,0)}{h^2}},$$

where \(d^2(p,q)\) is the squared Euclidean distance between two patches of size \((2f+1) \times (2f+1)\) centered at \(p\) and \(q\)

$$d^2(p,q) = \frac{1}{(2f+1)^2} \sum_{j \in B(0,f)} (\bold{I}[p+j]-\bold{I}[q+j])^2$$

and \(h\) is the filtering parameter defined as

$$h = k \sigma.$$

The value of \(k\) decreases as the search size, \(r\), increases.

Pixelwise Implementation

For a patch \(B=B(p,f)\) centered at \(p=(x,y)\) with size \((2f+1) \times (2f+1)\) of the image \(\bold{I}\), we compute the value of pixel \(p\) as follows:

$$\bold{I}[p]=\frac{1}{(2f+1)^2} \sum_{Q(q,f) | q \in B(p,f)} \hat{Q}[p],$$

$$\hat{Q} = \frac{1}{C} \sum_{Q(q,f) \in B(p,r)} Q \cdot w(B,Q), \qquad C = \sum_{Q(q,f) \in B(p,r)} w(B,Q),$$

where \(B(p,r)\) is a search window centered at \(p\) with size \((2r+1)\times(2r+1)\) pixels. The value of search size, \(r\), is determined by the standard deviation of noise \(\sigma\).

Patchwise Implementation

The weight \(w(p,q)\) is computed by using an exponential kernel

$$w(p,q) = e^{-\frac{\max(d^2-2\sigma^2,0)}{h^2}},$$

where \(d^2\) is defined to be the squared Euclidean distance 

$$d^2(p,q) = \frac{1}{(2f+1)^2} \sum_{j \in B(0,f)} (\bold{I}[p+j]-\bold{I}[q+j])^2$$

and \(h\) is the filtering parameter defined as

$$h = k \sigma.$$

The value of \(k\) decreases as the search size, \(r\), increases.

Patchwise Implementation

For pixel \(\bold{x} = (x_1, x_2)\), and \(\bold{y} = (y_1, y_2)\) of the image \(\bold{I}\), we consider their corresponding square neighborhoods:

$$S(\bold{x}) = \{(x_1 + s, x_2 + t) | -k \le s, t\le k\}$$

$$S(\bold{y}) = \{(y_1 + s, y_2 + t) | -k \le s, t\le k\}$$

Define \(I_\bold{x} = I|_{S(\bold{x})}\), and \(I_\bold{y} = I|_{S(\bold{y})}\) as two \((2k+1) \times (2k+1)\) patches. We first apply a smoothing Gaussian filter to \(I_\bold{x} \) and \(I_\bold{y} \):

$$\tilde{I }_\bold{x} = g*I_\bold{x} $$

$$\tilde{I }_\bold{y} = g*I_\bold{y} $$

Alternative Formulation - pixelwise

Then, we compute the norm distance between the two image patches \(\tilde{I}_{\bold{x}}\) and \(\tilde{I}_{\bold{y}}\):

$$d(\bold{x}, \bold{y}) = ||\tilde{I}_\bold{x} - \tilde{I}_\bold{y}||^2 $$

For each two pixels, we define the weight:

$$w (\bold{x}, \bold{y}) = \exp^{-d(\bold{x}, \bold{y})/h^2}$$

where h is a parameter (typically you choose it large if the noise is large)

Finally, we define the output at pixel \(\bold{x}\) as follows:

$$J(\bold{x}) = \frac{\sum_\bold{y} w(\bold{x}, \bold{y})I(\bold{y})}{\sum_{\bold{y}} w(\bold{x}, \bold{y})}$$

For pixel \(\bold{x} = (x_1, x_2)\), and \(\bold{y} = (y_1, y_2)\) of the image \(\bold{I}\), we consider their corresponding square neighborhoods:

$$S(\bold{x}) = \{(x_1 + s, x_2 + t) | -k \le s, t\le k\}$$

$$S(\bold{y}) = \{(y_1 + s, y_2 + t) | -k \le s, t\le k\}$$

Define \(I_\bold{x} = I|_{S(\bold{x})}\), and \(I_\bold{y} = I|_{S(\bold{y})}\) as two \((2k+1) \times (2k+1)\) patches. We first apply a smoothing Gaussian filter to \(I_\bold{x} \) and \(I_\bold{y} \):

$$\tilde{I }_\bold{x} = g*I_\bold{x} $$

$$\tilde{I }_\bold{y} = g*I_\bold{y} $$

Alternative Formulation - pixelwise rotation

Then, we rotate \(I_\bold{y}\) by \(\theta \in (0, \frac{\pi}{2}, \pi, \frac{3\pi}{2})\), compute the norm distance between the two image patches \(\tilde{I}_{\bold{x}}\) and \(\tilde{I}_{\bold{y}}\):

$$d(\bold{x}, \bold{y})_{\theta} = ||\tilde{I}_\bold{x} - \tilde{I}_\bold{y}|\theta||^2 $$

and find the angle that minimize the norm distance:

$$\theta^*=argmin_{\theta}(d(\bold{x}, \bold{y})_{\theta})$$

For each two pixels, we define the weight:

$$w (\bold{x}, \bold{y}) = \exp^{-d(\bold{x}, \bold{y})_{\theta^*}/h^2}$$

Finally, we define the output at pixel \(\bold{x}\) as follows:

$$J(\bold{x}) = \frac{\sum_\bold{y} w(\bold{x}, \bold{y})I(\bold{y})}{\sum_{\bold{y}} w(\bold{x}, \bold{y})}$$

Alternative Formulation - patchwise

For pixel \(\bold{x} = (x_1, x_2)\), and \(\bold{y} = (y_1, y_2)\) of the image \(\bold{I}\), we consider their corresponding square neighborhoods:

$$S(\bold{x}) = \{(x_1 + s, x_2 + t) | -k \le s, t\le k\}$$

$$S(\bold{y}) = \{(y_1 + s, y_2 + t) | -k \le s, t\le k\}$$

Define \(I_\bold{x} = I|_{S(\bold{x})}\), and \(I_\bold{y} = I|_{S(\bold{y})}\) as two \((2k+1) \times (2k+1)\) patches. We first apply a smoothing Gaussian filter to \(I_\bold{x} \) and \(I_\bold{y} \):

$$\tilde{I }_\bold{x} = g*I_\bold{x} $$

$$\tilde{I }_\bold{y} = g*I_\bold{y} $$

Then, we compute the norm distance between the two image patches \(\tilde{I}_{\bold{x}}\) and \(\tilde{I}_{\bold{y}}\):

$$d(\bold{x}, \bold{y}) = ||\tilde{I}_\bold{x} - \tilde{I}_\bold{y}||^2 $$

For each two pixels, we define the weight:

$$w (\bold{x}, \bold{y}) = \exp^{-d(\bold{x}, \bold{y})/h^2}$$

Next, we calculate the weighted result for a patch centered at \(\bold{x}\) as follows:

$$J_\bold{x} = J|S(\bold{x}) = \frac{\sum_\bold{y} w(\bold{x}, \bold{y}) I_\bold{y}}{\sum_{\bold{y}} w(\bold{x}, \bold{y})}$$

The previous step will result in multiple estimates for each pixel, hence, we calculate the final estimate of each pixel as follows:

$$J(\bold{x}) = \frac{1}{(2k + 1)^2} \sum_{i} B_i(\bold{x})$$

where \(B_i(x)\) is a \((2k+1)\times (2k+1)\) patch in which x is included.

Note: we use \((2k+1)^2\) in the equation above for most of the pixels in the central region of an image. This would not be the case for boundary pixels.

Alternative Formulation - patchwise rotation

For pixel \(\bold{x} = (x_1, x_2)\), and \(\bold{y} = (y_1, y_2)\) of the image \(\bold{I}\), we consider their corresponding square neighborhoods:

$$S(\bold{x}) = \{(x_1 + s, x_2 + t) | -k \le s, t\le k\}$$

$$S(\bold{y}) = \{(y_1 + s, y_2 + t) | -k \le s, t\le k\}$$

Define \(I_\bold{x} = I|_{S(\bold{x})}\), and \(I_\bold{y} = I|_{S(\bold{y})}\) as two \((2k+1) \times (2k+1)\) patches. We first apply a smoothing Gaussian filter to \(I_\bold{x} \) and \(I_\bold{y} \):

$$\tilde{I }_\bold{x} = g*I_\bold{x} $$

$$\tilde{I }_\bold{y} = g*I_\bold{y} $$

Then, we rotate \(I_\bold{y}\) by \(\theta \in (0, \frac{\pi}{2}, \pi, \frac{3\pi}{2})\), compute the norm distance between the two image patches \(\tilde{I}_{\bold{x}}\) and \(\tilde{I}_{\bold{y}}\):

$$d(\bold{x}, \bold{y})_{\theta} = ||\tilde{I}_\bold{x} - \tilde{I}_\bold{y}|\theta||^2 $$

and find the angle that minimize the norm distance:

$$\theta^*=argmin_{\theta}(d(\bold{x}, \bold{y})_{\theta})$$

For each two pixels, we define the weight:

$$w (\bold{x}, \bold{y}) = \exp^{-d(\bold{x}, \bold{y})_{\theta^*}/h^2}$$

Next, we calculate the weighted result for a patch centered at \(\bold{x}\) as follows:

$$J_\bold{x} = J|S(\bold{x}) = \frac{\sum_\bold{y} w(\bold{x}, \bold{y}) I_\bold{y}|\theta^*}{\sum_{\bold{y}} w(\bold{x}, \bold{y})}$$

The previous step will result in multiple estimates for each pixel, hence, we calculate the final estimate of each pixel as follows:

$$J(\bold{x}) = \frac{1}{(2k + 1)^2} \sum_{i} B_i(\bold{x})$$

where \(B_i(x)\) is a \((2k+1)\times (2k+1)\) patch in which x is included.

Note: we use \((2k+1)^2\) in the equation above for most of the pixels in the central region of an image. This would not be the case for boundary pixels.

Results

Results

Results

Results

Results

Results

Results