Jul 21, 2015

Denoise filters: Mean

Continuing with article series about the denoise filters, we had seen that the Gauss filter had a very high sensitivity to noise due to it's static nature. Now is time to analyze a dynamic denoise filter based in the arithmetic mean of the input pixels, this is called the Mean Denoise Filter. This filter adapts it's convolution kernel according to the values of the input pixels.

The weights for the convolution kernel are calculated in a very similar way as we did with the Gauss filter, but this time, the exponential function depends of the values of the pixels, and the mean ($\mu$) and standard deviation ($\sigma$) between them.

$$w_{ij}=\frac{1}{\sigma\sqrt{2\pi}}\exp^{-\frac{(I_{x+i-r,y+j-r}-\mu)^2}{2\sigma^2}}$$

The mean is calculated as:

$$\mu_{xy}=\frac{1}{(2r+1)^2}\sum_{j=-r}^{r}\sum_{i=-r}^{r}I_{x+i,y+j}$$

$$\sigma_{xy}^2=\frac{1}{(2r+1)^2}\sum_{j=-r}^{r}\sum_{i=-r}^{r}(I_{x+i,y+j}-\mu_{xy})^2$$

The problem with this equations is that we need to loop over all pixels in the scan windows to calculate the mean, then loop over all pixels again to calculate the standard deviation, and finally, loop over again to calculate the convolution, making it very slow for real time image processing. So, I'll make some arrangements to the equations to make the algorithm slightly faster.

$$\sigma_{xy}^2=\frac{1}{(2r+1)^2}\sum_{j=-r}^{r}\sum_{i=-r}^{r}\left(I_{x+i,y+j}^2-2\mu I_{x+i,y+j}+\mu_{xy}^2\right)$$

$$\sigma_{xy}^2=\frac{1}{(2r+1)^2} \left(\sum_{j=-r}^{r}\sum_{i=-r}^{r}I_{x+i,y+j}^2-2\mu_{xy}\sum_{j=-r}^{r}\sum_{i=-r}^{r}I_{x+i,y+j}+\mu_{xy}^2\sum_{j=-r}^{r}\sum_{i=-r}^{r}1\right)$$

And because

$$\sum_{j=-r}^{r}\sum_{i=-r}^{r}1=(2r+1)^2$$

then

$$\sigma_{xy}^2=\frac{1}{(2r+1)^2} \left(\sum_{j=-r}^{r}\sum_{i=-r}^{r}I_{x+i,y+j}^2-2\mu_{xy}^2(2r+1)^2+\mu_{xy}^2(2r+1)^2\right)$$

An finally

$$\sigma_{xy}^2=\left(\frac{1}{(2r+1)^2} \sum_{j=-r}^{r}\sum_{i=-r}^{r}I_{x+i,y+j}^2 \right )-\mu_{xy}^2$$

The mean ($\mu$) and the term

$$\sum_{j=-r}^{r}\sum_{i=-r}^{r}I_{x+i,y+j}^2$$

can be calculated faster with a trick we learned before, Integral Image. The code for calculating the integral image and the quadratic integral image is:
void integralImage(const QImage &image,
QVector<PixelU8> &planes,
QVector<PixelU32> &integral,
QVector<PixelU64> &integral2)
{
int oWidth = image.width() + 1;
int oHeight = image.height() + 1;
planes.resize(image.width() * image.height());
integral.resize(oWidth * oHeight);
integral2.resize(oWidth * oHeight);

for (int y = 1; y < oHeight; y++) {
const QRgb *line = (const QRgb *) image.constScanLine(y - 1);
PixelU8 *planesLine = planes.data()
+ (y - 1) * image.width();

// Reset current line summation.
PixelU32 sum;
PixelU64 sum2;

for (int x = 1; x < oWidth; x++) {
QRgb pixel = line[x - 1];

// Accumulate pixels in current line.
sum += pixel;
sum2 += pow2(pixel);

// Offset to the current line.
int offset = x + y * oWidth;

// Offset to the previous line.
// equivalent to x + (y - 1) * oWidth;
int offsetPrevious = offset - oWidth;

planesLine[x - 1] = pixel;

// Accumulate current line and previous line.
integral[offset] = sum + integral[offsetPrevious];
integral2[offset] = sum2 + integral2[offsetPrevious];
}
}
}


And the code for the Mean Denoise is:
int main(int argc, char *argv[])
{
QCoreApplication a(argc, argv);
Q_UNUSED(a)

QImage inImage("lena.png");
inImage = inImage.convertToFormat(QImage::Format_RGB32);
QImage outImage(inImage.size(), inImage.format());

// Here we configure the denoise parameters.
int mu = 0;
qreal sigma = 1;

// Add noise to the image
qsrand(QTime::currentTime().msec());

for (int i = 0; i < 25000; i++) {
inImage.setPixel(qrand() % inImage.width(),
qrand() % inImage.height(),
qRgb(qrand() % 256,
qrand() % 256,
qrand() % 256));
}

int oWidth = inImage.width() + 1;
QVector<PixelU8> planes;
QVector<PixelU32> integral;
QVector<PixelU64> integral2;
integralImage(inImage, planes, integral, integral2);

QElapsedTimer timer;
timer.start();

for (int y = 0, pos = 0; y < inImage.height(); y++) {
const QRgb *iLine = (const QRgb *) inImage.constScanLine(y);
QRgb *oLine = (QRgb *) outImage.scanLine(y);
int yp = qMax(y - radius, 0);
int kh = qMin(y + radius, inImage.height() - 1) - yp + 1;

for (int x = 0; x < inImage.width(); x++, pos++) {
int xp = qMax(x - radius, 0);
int kw = qMin(x + radius, inImage.width() - 1) - xp + 1;

// Calculate summation and cuadratic summation of the pixels.
PixelU32 sum = integralSum(integral.constData(), oWidth,
xp, yp, kw, kh);
PixelU64 sum2 = integralSum(integral2.constData(), oWidth,
xp, yp, kw, kh);
qreal ks = kw * kh;

// Calculate mean and standard deviation.
PixelReal mean = sum / ks;
PixelReal dev = sqrt(ks * sum2 - pow2(sum)) / ks;

mean = bound(0., mean + mu, 255.);
dev = bound(0., mult(sigma, dev), 127.);

PixelReal sumP;
PixelReal sumW;

for (int j = 0; j < kh; j++) {
const PixelU8 *line = planes.constData()
+ (yp + j) * inImage.width();

for (int i = 0; i < kw; i++) {
// Calculate weighted avverage.
PixelU8 pixel = line[xp + i];
PixelReal d = mean - pixel;
PixelReal h = mult(-2., dev * dev);
PixelReal weight = exp(d * d / h);
sumP += weight * pixel;
sumW += weight;
}
}

// Normalize result.
sumP /= sumW;

oLine[x] = qRgba(sumP.r, sumP.g, sumP.b, qAlpha(iLine[x]));
}
}

qDebug() << timer.elapsed();
outImage.save("mean.png");

return EXIT_SUCCESS;
}


For an input image with a a noise of 25000

The output for a radius = 2 and sigma = 1 is

This time you can also try changing the mu value, this will make the image brighter or darker.