原理部分可以参看前一篇博客
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); } } } }
效果如图：