双边滤波

Posted 2020-09-21 cynchanpin

Qt 平台，双边滤波原理代码例如以下：

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

using namespace cv;
using namespace std;

double Gaussian( double x, double sigma )//高斯核參数计算函数
{
    return exp(-x/(2*sigma*sigma));
}

Mat GaussiCore( int d, double sigma )//高斯核计算
{
    int dHalf = (d-1)/2;
    Mat core( d, d, CV_64FC1, Scalar(0));
    for ( int i = -dHalf; i < dHalf; i ++ )
    {
        double *corePtr = core.ptr<double>(i+dHalf);
        for ( int j = -dHalf; j < dHalf; j ++ )
        {
            corePtr[j+dHalf] = Gaussian((i*i+j*j), sigma);
        }
    }
    return core;
}

int main()
{
    double duration;
    int d = 11;//滤波核直径
    int dHalf = (d-1)/2;//滤波核半径
    double sigma_d = 22;//双边滤波核高斯系数标准差
    double sigma_r = 11;//双边滤波核空间域系数标准差
    double temp1 = 0.0;
    double temp2 = 0.0;
    Mat src = imread("lena.jpg", 0);
    Mat srcBorder( src.rows+d-1, src.cols+d-1, CV_8UC1, Scalar(128));
    int srcRow = src.rows;
    int srcCol = src.cols;
    Mat dst( srcRow, srcCol, CV_8UC1, Scalar(0));
    Mat gaussiCore;
duration = static_cast<double>(getTickCount());
    for ( int i = 0; i < srcRow; i ++ )
    {
        uchar *srcPtr = src.ptr<uchar>(i);
        uchar *srcBorderPtr = srcBorder.ptr<uchar>(i+dHalf);
        for ( int j = 0; j < srcCol; j ++ )
        {
            srcBorderPtr[j+dHalf] = srcPtr[j];
        }
    }

    gaussiCore = GaussiCore( d, sigma_d );

    for ( int i = 0; i < srcRow; i ++ )
    {
        uchar *dstPtr = dst.ptr<uchar>(i);
        uchar *srcBorderPtr1 = srcBorder.ptr<uchar>(i+dHalf);
        for ( int j = 0; j < srcCol; j ++ )
        {
            temp1 = 0.0;
            temp2 = 0.0;
            for ( int n = -dHalf+i, rr = 0; n < dHalf+i; n ++, rr ++ )
            {
                uchar *srcBorderPtr2 = srcBorder.ptr<uchar>(n);
                double *gaussiCorePtr = gaussiCore.ptr<double>(rr);
                for ( int l = -dHalf+j, rl = 0; l < dHalf+j; l ++, rl ++ )
                {
                    temp1 += (double)srcBorderPtr2[l] * gaussiCorePtr[rl]
                            * Gaussian((srcBorderPtr2[l]-srcBorderPtr1[j+dHalf])*(srcBorderPtr2[l]-srcBorderPtr1[j+dHalf]), sigma_r);
                    temp2 +=  gaussiCorePtr[rl]
                            * Gaussian((srcBorderPtr2[l]-srcBorderPtr1[j+dHalf])*(srcBorderPtr2[l]-srcBorderPtr1[j+dHalf]), sigma_r);
                }
            }
            dstPtr[j] = temp1/temp2;
        }
    }

duration = static_cast<double>(getTickCount()) - duration;
duration /= getTickFrequency();
cout << duration << endl;
    namedWindow("src", 0);
    imshow("src", src);
    namedWindow("dst", 0);
    imshow("dst", dst);
    waitKey(0);
    return 0;
}

由于代码单纯介绍双边滤波原理。所以与 OpenCV 自带双边滤波速度还是有非常大差距，有时间会继续研究关于双边滤波的加速算法。下面是双边滤波核计算公式，对于双边滤波为什么能保边降噪。这里简单说明一下：首先降噪是肯定的。由于滤波核有高斯滤波系数。总所周知高斯函数有滤波降噪的作用。而双边滤波系数除了有高斯系数之外。还增加了像素领域的值域差系数，即像素领域内的像素值与像素本身的像素值相差越大，依据最后一个总体公式能够看出，总体系数会变小，从而该像素受到影响也就越小，而仅仅有在图像的边缘区域才会有这样的情况出生，所以在图像的边缘附近差点儿不受滤波系数的影响，这样就起到了保边的效果。

1、离散卷积公式：

2、高斯系数：

3、空间域系数：

4、总体双边滤波系数：

做了一点加速。比原来的快几十毫秒，代码例如以下：

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

using namespace cv;
using namespace std;

double Gaussian( double x, double sigma )//高斯核參数计算函数
{
    return exp(-x/(2*sigma*sigma));
}

Mat GaussiCore( int d, double sigma )//高斯核计算
{
    int dHalf = (d-1)/2;
    int row = d;
    int col = d;
    Mat core( row, col, CV_64FC1, Scalar(0));
    if ( core.isContinuous() )
    {
        double *corePtr = core.ptr<double>(0);
        for ( int i = 0; i < row; i ++ )
        {
            for ( int j = 0; j < col; j ++ )
            {
                corePtr[i*col+j] = Gaussian(((i-dHalf)*(i-dHalf)+(j-dHalf)*(j-dHalf)), sigma);
            }
        }
    }
    else
    {
        for ( int i = -dHalf; i < dHalf; i ++ )
        {
            double *corePtr = core.ptr<double>(i+dHalf);
            for ( int j = -dHalf; j < dHalf; j ++ )
            {
                corePtr[j+dHalf] = Gaussian((i*i+j*j), sigma);
            }
        }
    }
    return core;
}

int main()
{
    double duration;
    int d = 11;//滤波核直径
    int dHalf = (d-1)/2;//滤波核半径
    double sigma_d = 22;//双边滤波核高斯系数标准差
    double sigma_r = 11;//双边滤波核空间域系数标准差
    double temp1 = 0.0;
    double temp2 = 0.0;
    Mat src = imread("lena.jpg", 0);
    Mat srcBorder( src.rows+d-1, src.cols+d-1, CV_8UC1, Scalar(128));
    int srcRow = src.rows;
    int srcCol = src.cols;
    Mat dst( srcRow, srcCol, CV_8UC1, Scalar(0));
    Mat gaussiCore;
duration = static_cast<double>(getTickCount());

    if ( src.isContinuous() && srcBorder.isContinuous() )
    {
        uchar *srcPtr = src.ptr<uchar>(0);
        uchar *srcBorderPtr = srcBorder.ptr<uchar>(0);
        for ( int i = 0; i < srcRow; i ++ )
        {
            for ( int j = 0; j < srcCol; j ++ )
            {
                srcBorderPtr[(i+dHalf)*(srcCol+d-1)+j+dHalf] = srcPtr[i*srcCol+j];
            }
        }
    }
    else
    {
        for ( int i = 0; i < srcRow; i ++ )
        {
            uchar *srcPtr = src.ptr<uchar>(i);
            uchar *srcBorderPtr = srcBorder.ptr<uchar>(i+dHalf);
            for ( int j = 0; j < srcCol; j ++ )
            {
                srcBorderPtr[j+dHalf] = srcPtr[j];
            }
        }
    }

    gaussiCore = GaussiCore( d, sigma_d );

    for ( int i = 0; i < srcRow; i ++ )
    {
        uchar *dstPtr = dst.ptr<uchar>(i);
        uchar *srcBorderPtr1 = srcBorder.ptr<uchar>(i+dHalf);
        for ( int j = 0; j < srcCol; j ++ )
        {
            temp1 = 0.0;
            temp2 = 0.0;
            for ( int n = -dHalf+i, rr = 0; n < dHalf+i; n ++, rr ++ )
            {
                uchar *srcBorderPtr2 = srcBorder.ptr<uchar>(n);
                double *gaussiCorePtr = gaussiCore.ptr<double>(rr);
                for ( int l = -dHalf+j, rl = 0; l < dHalf+j; l ++, rl ++ )
                {
                    temp1 += (double)srcBorderPtr2[l] * gaussiCorePtr[rl]
                            * Gaussian((srcBorderPtr2[l]-srcBorderPtr1[j+dHalf])*(srcBorderPtr2[l]-srcBorderPtr1[j+dHalf]), sigma_r);
                    temp2 +=  gaussiCorePtr[rl]
                            * Gaussian((srcBorderPtr2[l]-srcBorderPtr1[j+dHalf])*(srcBorderPtr2[l]-srcBorderPtr1[j+dHalf]), sigma_r);
                }
            }
            dstPtr[j] = temp1/temp2;
        }
    }

duration = static_cast<double>(getTickCount()) - duration;
duration /= getTickFrequency();
cout << duration << endl;
    namedWindow("src", 0);
    imshow("src", src);
    namedWindow("dst", 0);
    imshow("dst", dst);
    waitKey(0);
    return 0;
}