图像去噪之非局部均值算法-原理和实现

最新推荐文章于 2024-04-14 22:46:26 发布

不停歇的蚂蚁

最新推荐文章于 2024-04-14 22:46:26 发布

阅读量1.5w

点赞数 4

分类专栏： opencv

本文链接：https://blog.csdn.net/feitian944/article/details/80369385

版权

opencv 专栏收录该内容

2 篇文章 1 订阅

订阅专栏

1、基本原理

是一种空间域滤波，和高斯模糊、均值滤波类似，当前像素点的值时通过周围像素的加权平均得到。不同之处在于权值计算策略不同，也就是下面公式中w计算方式不同；
这里写图片描述

2、权重策略

这里写图片描述

以上图为例，对其去噪；对于帽子边缘的一个像素点，相比于高斯模糊或者均值模糊，用附近帽子边缘的像素点求均值更能好。那么计算机如何知道哪些像素也同样的是帽子边缘的像素？帽子边缘具有相似的结构，如果我们在每个像素周围画一个方框，方框内像素的欧式距离越小，说明结构越相似，中间的像素越是处于相同的结构位置，对应权重要更大，这个框的大小称kernel size。对多大范围内(B(p,r))的像素求加权平均，在这个算法中成为search size;

距离计算公式：
这里写图片描述

距离越小，越相似，权重越大
这里写图片描述

3、实现代码

#include <iostream>
#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/photo.hpp>

using namespace std;
using namespace cv;




int distance2(Mat m1, Mat m2)
{
    Mat dst;
    absdiff(m1, m2, dst);
    int sum = 0;
    for (int j = 0; j < dst.rows; j++) {
        for (int i = 0; i < dst.cols; i++) {
            sum += dst.at<uchar>(j,i) *  dst.at<uchar>(j, i);
        }
    }
    return sum;
}

void nlmean(Mat src, Mat &dst, double h, double sigma, int halfKernelSize, int halfSearchSize)
{
    Mat boardSrc;
    dst.create(src.rows, src.cols, CV_8UC1);
    int boardSize = halfKernelSize + halfSearchSize;
    copyMakeBorder(src, boardSrc, boardSize, boardSize, boardSize, boardSize, BORDER_REFLECT);
    double h2 = h*h;
    double sigma2 = 2 * sigma*sigma;

    int rows = src.rows;
    int cols = src.cols;

    int lookupTableSize = 256 * 256 *(2* halfKernelSize +1)*(2* halfKernelSize+1);
    double *lookupTable = new double[lookupTableSize];

    for (int i = 0; i < lookupTableSize; i++) {
        if(i<sigma2)
            lookupTable[i] =1;
        else
            lookupTable[i] = exp(-i/ h2/(2* halfKernelSize+1) /(2* halfKernelSize+1));
    }

    for (int j = boardSize; j < boardSize+ rows; j++) {
        for (int i = boardSize; i < boardSize+cols; i++) {
            Mat patchA=boardSrc( Range(j - halfKernelSize, j + halfKernelSize), Range (i - halfKernelSize, i + halfKernelSize));
            double w = 0;
            double p = 0;
            double sumw = 0;

            for (int sr = -halfSearchSize; sr <= halfSearchSize; sr++) {
                for (int sc = -halfSearchSize; sc <= halfSearchSize; sc++) {

                    Mat patchB = boardSrc(Range(j+sr - halfKernelSize, j+sr + halfKernelSize), Range(i +sc- halfKernelSize, i+sc + halfKernelSize));
                    int d2 = distance2(patchA, patchB);

                    w = lookupTable[d2];
                    p += boardSrc.at<uchar>(j+sr, i+sc)*w;

                    sumw += w;
                }
            }

            dst.at<uchar>(j - boardSize, i - boardSize) = saturate_cast<uchar>(p / sumw);

        }
    }
    delete[]lookupTable;

}

int main()
{
    Mat src = imread("lena.png", IMREAD_GRAYSCALE);
    Mat dst;
    Mat cvdst;

    fastNlMeansDenoising(src, cvdst, 3, 7, 10);


    namedWindow("src");
    imshow("src", src);
    imshow("cvdst", cvdst);
    nlmean(src, dst, 3, 5, 3, 10);
    imshow("mydst", dst);

    Mat diff;
    absdiff(src, dst, diff);
    imshow("diff", diff);
    waitKey(0);
}