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}$$
And the quadratic deviation as:
$$\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 radius = 2;
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;
}
The full source code for this article is located at at my github.
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.
No comments:
Post a Comment