天天看點

聯合雙邊濾波器(joint bilateral filter) 代碼及詳細注釋【OpenCV】

原理部分可以參看前一篇 部落格

void jointBilateralFilter(const Mat &src, Mat &dst, int d, double sigma_color,
                          double sigma_space, Mat &joint = Mat(), int borderType =
                              BORDER_REPLICATE)
{
    Size size = src.size();
    if (dst.empty())
        dst = Mat::zeros(src.size(), src.type());
 
    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;
 
    double gauss_color_coeff = -0.5 / (sigma_color * sigma_color);
    double gauss_space_coeff = -0.5 / (sigma_space * sigma_space);
 
    if (joint.empty())
        src.copyTo(joint);
 
    const int cn = src.channels();
    const int cnj = joint.channels();
 
    int radius;
    if (d <= 0)
        radius = cvRound(sigma_space * 1.5);    // 根據 sigma_space 計算 radius
    else
        radius = d / 2;
    radius = MAX(radius, 1);
    d = radius * 2 + 1; // 重新計算 像素“矩形”鄰域的直徑d,確定是奇數
 
    // 擴充 src 和 joint 長寬各2*radius
    Mat jim;
    Mat sim;
    copyMakeBorder(joint, jim, radius, radius, radius, radius, borderType);
    copyMakeBorder(src, sim, radius, radius, radius, radius, borderType);
 
    // cnj: joint的通道數
    vector<float> _color_weight(cnj * 256);
    vector<float> _space_weight(d * d);  // (2*radius + 1)^2
    vector<int> _space_ofs_jnt(d * d);
    vector<int> _space_ofs_src(d * d);
    float *color_weight = &_color_weight[0];
    float *space_weight = &_space_weight[0];
    int *space_ofs_jnt = &_space_ofs_jnt[0];
    int *space_ofs_src = &_space_ofs_src[0];
 
    // initialize color-related bilateral filter coefficients
    // 色差的高斯權重
    for (int i = 0; i < 256 * cnj; i++)
        color_weight[i] = (float) std::exp(i * i * gauss_color_coeff);
 
    int maxk = 0;   // 0 - (2*radius + 1)^2
    // initialize space-related bilateral filter coefficients
    for (int i = -radius; i <= radius; i++)
    {
        for (int j = -radius; 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_jnt[maxk] = (int) (i * jim.step + j * cnj);           // joint 鄰域内的相對坐标 (i, j)【偏移量】, 左上角為(-radius, -radius),右下角為(radius, radius)
            space_ofs_src[maxk++] = (int) (i * sim.step + j * cn);      // src 鄰域内的相對坐标 (i, j)
        }
    }
    #pragma omp parallel for
    for (int i = 0; i < size.height; i++)
    {
        const uchar *jptr = jim.data + (i + radius) * jim.step + radius * cnj;  // &jim.ptr(i+radius)[radius]
        const uchar *sptr = sim.data + (i + radius) * sim.step + radius * cn; // &sim.ptr(i+radius)[radius]
        uchar *dptr = dst.data + i * dst.step;                                              // dst.ptr(i)
 
        // src 和 joint 通道數不同的四種情況
        if (cn == 1 && cnj == 1)
        {
            for (int j = 0; j < size.width; j++)
            {
                float sum = 0, wsum = 0;
                int val0 = jptr[j]; // jim.ptr(i + radius)[j + radius]
 
                for (int k = 0; k < maxk; k++)
                {
                    int val = jptr[j + space_ofs_src[k]];    // jim.ptr(i + radius + offset_x)[j + radius + offset_y]
                    int val2 = sptr[j + space_ofs_src[k]];  // sim.ptr(i + radius + offset_x)[j + radius + offset_y]
 
                    // 根據joint目前像素和鄰域像素的 距離權重 和 色差權重,計算綜合的權重
                    float w = space_weight[k]
                              * color_weight[std::abs(val - val0)];
                    sum += val2 * w;    // 統計 src 鄰域内的像素帶權和
                    wsum += w;          // 統計權重和
                }
                // overflow is not possible here => there is no need to use CV_CAST_8U
                // 歸一化 src 鄰域内的像素帶權和,并賦給 dst對應的像素
                dptr[j] = (uchar) cvRound(sum / wsum);
            }
        }
        else if (cn == 3 && cnj == 3)
        {
            for (int j = 0; j < size.width * 3; j += 3)
            {
                float sum_b = 0, sum_g = 0, sum_r = 0, wsum = 0;
                int b0 = jptr[j], g0 = jptr[j + 1], r0 = jptr[j + 2];   // jim.ptr(i + radius)[j + radius][0...2]
                for (int k = 0; k < maxk; k++)
                {
                    const uchar *sptr_k = jptr + j + space_ofs_src[k];  
                    const uchar *sptr_k2 = sptr + j + space_ofs_src[k];
 
                    int b = sptr_k[0], g = sptr_k[1], r = sptr_k[2];     // jim.ptr(i + radius + offset_x)[j + radius + offset_y][0...2]
                    float w = space_weight[k]
                              * color_weight[std::abs(b - b0) + std::abs(g - g0)
                                             + std::abs(r - r0)];
                    sum_b += sptr_k2[0] * w;    // sim.ptr(i + radius + offset_x)[j + radius + offset_y][0...2]
                    sum_g += sptr_k2[1] * w;
                    sum_r += sptr_k2[2] * 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;
            }
        }
        else if (cn == 1 && cnj == 3)
        {
            for (int j = 0, l = 0; j < size.width * 3; j += 3, l++)
            {
                float sum_b = 0, wsum = 0;
                
                int b0 = jptr[j], g0 = jptr[j + 1], r0 = jptr[j + 2];   // jim.ptr(i + radius)[j + radius][0...2]
                for (int k = 0; k < maxk; k++)
                {
                    int val = *(sptr + l + space_ofs_src[k]);   // sim.ptr(i + radius + offset_x)[l + radius + offset_y]
 
                    const uchar *sptr_k = jptr + j + space_ofs_jnt[k];
                    int b = sptr_k[0], g = sptr_k[1], r = sptr_k[2];        // jim.ptr(i + radius + offset_x)[j + radius + offset_y][0...2]
 
                    float w = space_weight[k]
                              * color_weight[std::abs(b - b0) + std::abs(g - g0)
                                             + std::abs(r - r0)];
                    sum_b += val * w;
                    wsum += w;
                }
                wsum = 1.f / wsum;
                b0 = cvRound(sum_b * wsum);
                dptr[l] = (uchar) b0;
            }
        }
        else if (cn == 3 && cnj == 1)
        {
            for (int j = 0, l = 0; j < size.width * 3; j += 3, l++)
            {
                float sum_b = 0, sum_g = 0, sum_r = 0, wsum = 0;
                int val0 = jptr[l]; // jim.ptr(i + radius)[l + radius]
                for (int k = 0; k < maxk; k++)
                {
                    int val = jptr[l + space_ofs_jnt[k]];   // jim.ptr(i + radius + offset_x)[l + radius + offset_y]
 
                    const uchar *sptr_k = sptr + j + space_ofs_src[k];      // sim.ptr(i + radius + offset_x)[j + radius + offset_y] 
 
                    float w = space_weight[k]
                              * color_weight[std::abs(val - val0)];
 
                    sum_b += sptr_k[0] * w; // sim.ptr(i + radius + offset_x)[j + radius + offset_y] [0...2]
                    sum_g += sptr_k[1] * w;
                    sum_r += sptr_k[2] * w;
                    wsum += w;
                }
 
                // overflow is not possible here => there is no need to use CV_CAST_8U
                wsum = 1.f / wsum;
                dptr[j] = (uchar) cvRound(sum_b * wsum);
                dptr[j + 1] = (uchar) cvRound(sum_g * wsum);
                dptr[j + 2] = (uchar) cvRound(sum_r * wsum);
            }
        }
    }
}      

效果如圖:

聯合雙邊濾波器(joint bilateral filter) 代碼及詳細注釋【OpenCV】

繼續閱讀