Denoise filters: Gauss

We have seen before about convolution filters, and at the end of that article we have seen a kernel used for removing noise in an image, that is random spurious pixels in a image, the process of removing that random pixels is called Denoise. In this case I will show you the general form of the denoising kernel based in gaussian distribution.

The gaussian kernel is generated with this equation:

$$w_{ij}=\frac{1}{\sigma\sqrt{2\pi}}\exp^{-\frac{(i - radius)^2 + (j - radius)^2}{2\sigma^2}}$$

By definition, the gaussian kernel tends to create a slightly blurred image, but where the pivot/central pixel has more importance than the surrounding pixels. If the value of $\sigma$ tends to 0 (zero) the filter will output the original image, If the value of $\sigma$ tends to $\infty$ we will have just the blur filter we had seen before.

The output pixel in the coordinates (x, y) will be the weighted average of all pixels in the scan window:

$$C_{xy}=\frac{1}{W_{xy}}\sum_{j=0}^{2r}\sum_{i=0}^{2r}w_{ij}I_{x+i-r,y+j-r}$$

with:

$$W_{xy}=\sum_{j=0}^{2r}\sum_{i=0}^{2r}w_{ij}$$

And we can see the weighted average as a kernel convolution in the form:

$$ kernel=\frac{1}{W_{xy}}\begin{vmatrix} w_{0,0} & w_{1,0} & \cdots & w_{2r,0} \\ w_{0,1} & w_{1,1} & \cdots & w_{2r,1} \\ \cdots & \cdots & \cdots & \cdots \\ w_{0,2r} & w_{1,2r} & \cdots & w_{2r,2r} \end{vmatrix} $$

The source code for the kernel is:

inline QVector<qreal> gaussKernel(int radius, qreal sigma, int *kl)
{
    int kw = 2 * radius + 1;
    QVector<qreal> kernel(kw * kw);
    qreal sum = 0;

    // Create convolution matrix according to the formula:
    for (int j = 0; j < kw; j++)
        for (int i = 0; i < kw; i++) {
            int x = i - radius;
            int y = j - radius;
            qreal sigma2 = -2 * sigma * sigma;
            qreal weight = exp((x * x + y * y) / sigma2);
            kernel[i + j * kw] = weight;
            sum += weight;
        }

    // Normalize weights.
    for (int i = 0; i < kernel.size(); i++)
        kernel[i] /= sum;

    *kl = kw;

    return kernel;
}

In the code, I just get rid of the constant:

$$\frac{1}{\sigma\sqrt{2\pi}}$$

It was not necessary because it was multiplying both numerator and denominator. The source code for the convolution is the same as we seen before, I just added the code for generate some random noise, added time metrics, and fixed the convolution result in the cases in which the kernel is smaller than required, for example when the convolution is being calculated in a corner or border.

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 radius = 3;
    qreal sigma = 1000;

    // Create gaussian denoise kernel.
    int kw;
    QVector<qreal> kernel = gaussKernel(radius, sigma, &kw);

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

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

    QElapsedTimer timer;
    timer.start();

    for (int y = 0; y < inImage.height(); y++) {
        const QRgb *iLine = (const QRgb *) inImage.constScanLine(y);
        QRgb *oLine = (QRgb *) outImage.scanLine(y);

        for (int x = 0; x < inImage.width(); x++) {
            PixelReal sum;
            qreal sumW = 0;

            // Apply kernel.
            for (int j = 0, pos = 0; j < kw; j++) {
                const QRgb *line = (const QRgb *) inImage.constScanLine(y + j - radius);

                if (y + j < radius
                    || y + j >= radius + inImage.height())
                    continue;

                for (int i = 0; i < kw; i++, pos++) {
                    if (x + i < radius
                        || x + i >= radius + inImage.width())
                        continue;

                    PixelU8 pixel(line[x + i - radius]);
                    qreal weight = kernel[pos];
                    sum += weight * pixel;
                    sumW += weight;
                }
            }

            // We normallize the kernel because the size of the kernel is not
            // fixed.
            sum /= sumW;

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

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

    return EXIT_SUCCESS;
}

If for example, the noisy image is:

The output for a noise of 5000, it will output:

As you can see, the output image is slightly blurred and with a lower contrast, this is because the gaussian filter is by definition a light blur, and because this is a static filter it is very sensitive to the noise, and because the probability of a value is 50/50, the noise tends to a gray image. You can go to mi github and download the full source code for this article.