雙邊濾波算法原理及實現

雙邊濾波算法原理

雙邊濾波是一種非線性濾波器，它可以達到保持邊緣、降噪平滑的效果。和其他濾波原理一樣，雙邊濾波也是采用加權平均的方法，用周邊像素亮度值的加權平均代表某個像素的強度，所用的加權平均基於高斯分布[1]。最重要的是，雙邊濾波的權重不僅考慮了像素的歐氏距離（如普通的高斯低通濾波，只考慮了位置對中心像素的影響），還考慮了像素范圍域中的輻射差異（例如卷積核中像素與中心像素之間相似程度、顏色強度，深度距離等），在計算中心像素的時候同時考慮這兩個權重。

雙邊濾波的核函數是空間域核與像素范圍域核的綜合結果：在圖像的平坦區域，像素值變化很小，對應的像素范圍域權重接近於1，此時空間域權重起主要作用，相當於進行高斯模糊；在圖像的邊緣區域，像素值變化很大，像素范圍域權重變大，從而保持了邊緣的信息。

雙邊濾波器代碼實現

void cv::bilateralFilter( InputArray _src, OutputArray _dst, int d,
                      double sigmaColor, double sigmaSpace,
                      int borderType )
{
    Mat src = _src.getMat();
    _dst.create( src.size(), src.type() );
    Mat dst = _dst.getMat();

    if( src.depth() == CV_8U )
        bilateralFilter_8u( src, dst, d, sigmaColor, sigmaSpace, borderType );
    else if( src.depth() == CV_32F )
        bilateralFilter_32f( src, dst, d, sigmaColor, sigmaSpace, borderType );
    else
        CV_Error( CV_StsUnsupportedFormat,
        "Bilateral filtering is only implemented for 8u and 32f images" );
}


static void
bilateralFilter_8u( const Mat& src, Mat& dst, int d,
    double sigma_color, double sigma_space,
    int borderType )
{

    int cn = src.channels();
    int i, j, k, maxk, radius;
    Size size = src.size();

    CV_Assert( (src.type() == CV_8UC1 || src.type() == CV_8UC3) &&
              src.type() == dst.type() && src.size() == dst.size() &&
              src.data != dst.data );

    if( sigma_color <= 0 )
        sigma_color = 1;
    if( sigma_space <= 0 )
        sigma_space = 1;

    // 計算顏色域和空間域的權重的高斯核系數, 均值 μ = 0;  exp(-1/(2*sigma^2))  
    double gauss_color_coeff = -0.5/(sigma_color*sigma_color);
    double gauss_space_coeff = -0.5/(sigma_space*sigma_space);

    // radius 為空間域的大小: 其值是 windosw_size 的一半    
    if( d <= 0 )
        radius = cvRound(sigma_space*1.5);
    else
        radius = d/2;
    radius = MAX(radius, 1);
    d = radius*2 + 1;

    Mat temp;
    copyMakeBorder( src, temp, radius, radius, radius, radius, borderType );

    vector<float> _color_weight(cn*256);
    vector<float> _space_weight(d*d);
    vector<int> _space_ofs(d*d);
    float* color_weight = &_color_weight[0];
    float* space_weight = &_space_weight[0];
    int* space_ofs = &_space_ofs[0];

    // 初始化顏色相關的濾波器系數: exp(-1*x^2/(2*sigma^2))  
    for( i = 0; i < 256*cn; i++ )
        color_weight[i] = (float)std::exp(i*i*gauss_color_coeff);

    // 初始化空間相關的濾波器系數和 offset:  
    for( i = -radius, maxk = 0; i <= radius; i++ )
    {
        j = -radius;

        for( ;j <= radius; j++ )
        {
            double r = std::sqrt((double)i*i + (double)j*j);
            if( r > radius )
                continue;
            space_weight[maxk] = (float)std::exp(r*r*gauss_space_coeff);
            space_ofs[maxk++] = (int)(i*temp.step + j*cn);
        }
    }

    // 開始計算濾波后的像素值  
    for( i = 0; i < 0, size.height; i++ )
    {
        const uchar* sptr = temp->ptr(i+radius) + radius*cn;  // 目標像素點 
        uchar* dptr = dest->ptr(i);

        if( cn == 1 )
        {
            // 按行開始遍歷    
            for( j = 0; j < size.width; j++ )
            {
                float sum = 0, wsum = 0;
                int val0 = sptr[j];
                
                // 遍歷當前中心點所在的空間鄰域  
                for( k = 0; k < maxk; k++ )
                {
                    int val = sptr[j + space_ofs[k]];
                    float w = space_weight[k]*color_weight[std::abs(val - val0)];
                    sum += val*w;
                    wsum += w;
                }
                
                // 這里不可能溢出, 因此不必使用 CV_CAST_8U. 
                dptr[j] = (uchar)cvRound(sum/wsum);
            }
        }
        else
        {
            assert( cn == 3 );
            for( j = 0; j < size.width*3; j += 3 )
            {
                float sum_b = 0, sum_g = 0, sum_r = 0, wsum = 0;
                int b0 = sptr[j], g0 = sptr[j+1], r0 = sptr[j+2];
                k = 0;
                
                for( ; k < maxk; k++ )
                {
                    const uchar* sptr_k = sptr + j + space_ofs[k];
                    int b = sptr_k[0], g = sptr_k[1], r = sptr_k[2];
                    float w = space_weight[k]*color_weight[std::abs(b - b0) +
                                                            std::abs(g - g0) + std::abs(r - r0)];
                    sum_b += b*w; sum_g += g*w; sum_r += r*w;
                    wsum += w;
                }
                wsum = 1.f/wsum;
                b0 = cvRound(sum_b*wsum);
                g0 = cvRound(sum_g*wsum);
                r0 = cvRound(sum_r*wsum);
                dptr[j] = (uchar)b0; 
                dptr[j+1] = (uchar)g0; 
                dptr[j+2] = (uchar)r0;
            }
        }
    }
}

參考

[1]: Bilateral Filters（雙邊濾波算法）原理及實現
[2]: 雙邊濾波算法介紹與實現