双边滤波(bilateral filter)以及联合双边滤波（joint bilateral filter）

双边滤波

自用备忘，若侵则删。

理论公式

利用二维高斯函数生成空间域核，一维高斯函数生成颜色域核。
- 空间域核：

其中，$(k,l)$为核中心坐标，$(i,j)$为核内邻域坐标。${\sigma }_d$为高斯函数的标准差，我个人感性上认为就是$\sigma$小的时候，高斯函数瘦瘦高高的，反之，矮矮胖胖的。

- 色彩域核（下式仅为灰度值的表示）：

$f(i,j)$代表图像在$(i,j)$处的灰度值，其他标识与空间域一致。

- 双边滤波器模板： 空间域核 色彩域核

- 思考一下，为什么可以保边？
虽然说在边缘处，空间域的权重很大，但是色彩变化大，这就导致了色彩域的权重很小，空间域的大权重优势被色彩域的小权重拉了下来，因此边界邻域的像素作用不大，保住了边缘自身。而对于平坦区域来说，色彩域的作用不大，几乎纯靠空间域核来撑，也就近似于普通的高斯滤波了。

- 示意图：

• 简易表达：
  
  滤波后的y值即为在小模板中加权（融合空间权和色彩权）平均后再除以权重和的值。
  在代码中也是进行了归一化，与公式统一，具体见四重循环内的代码。

代码（C++）

/*
双边滤波器思想：
	1. 在高斯滤波器基础上增加了图像像素之间的相似度的考虑，进而达到保边的效果。

双边滤波器的实现：
	1. 分别求空间域以及颜色域的权重系数
	2. 根据双边公式，移动窗口，处理整张影像
*/

#include "opencv2/opencv.hpp"
#include"opencv2/highgui/highgui.hpp"
#include <math.h>
using namespace cv;
using namespace std;

const int kernel_size = 7;
const int sigma_distance = 3;
const int sigma_color = 20;
#define DEBUG 0

// 空间域权重确定
// @ distance_kernel 空间域核,1D
void GetDistanceWeight(double* distance_kernel)
{
	int k = kernel_size / 2;
	double distance_kernel_2D[kernel_size][kernel_size];
	double delta_square = 2* sigma_distance * sigma_distance; //分母
	for (int i = -k; i <= k; i++) {
		for (int j = -k; j <= k; j++) {
			double distance_numerator = i * i + j * j;
			distance_kernel_2D[i + k][j + k] = exp(-1.0 * distance_numerator / delta_square);
#ifdef DEBUG
			{
				cout << i+k << " " <<j+k <<" "<< distance_kernel_2D[i + k][j + k] << endl;
			}
#endif
		}
	}

	// 将2D kernel 转换为 1D kernel
	for (int i = 0; i < kernel_size; i++) {
		for (int j = 0; j < kernel_size; j++) {
			distance_kernel[kernel_size * i + j] = distance_kernel_2D[i][j];
#ifdef DEBUG
			{
				cout << kernel_size * i + j << " " << distance_kernel[kernel_size * i + j] << endl;
			}
#endif
			
		}
	}
}

// 颜色域权重确定
// @ color_kernel 颜色域核，1D，长度为256
void GetColorWeight(double* color_kernel) {
	for (int i = 0; i < 256; i++) {
		color_kernel[i] = exp(-1.0 * (i * i) / (2 * sigma_color * sigma_color));
#ifdef DEBUG
		cout<< i <<" "<<  color_kernel[i] << endl;
#endif // DEBUG

	}
}

// 双边滤波
// @ src 待滤波的影像
// @ distance_kernel 空间域核
// @ color_kernel 颜色域核
// @ dst 输出的影像

void BilateralFilter(Mat& src, double* distance_kernel, double* color_kernel, Mat& dst) {
	dst = src.clone();
	int n_rows = dst.rows;
	int n_cols = dst.cols;
	int n_channels = dst.channels();
	int n_cols_with_channels = n_cols * n_channels;
	int half_kernel_size = kernel_size / 2;

	int index;
	double pixel_sum;
	double weight_sum = 0;
	double temp_bilateral_weight = 0;
	// 边界不做处理
	for (int i = half_kernel_size; i < (n_rows - half_kernel_size); i++) {
		uchar* pt_dst = dst.ptr<uchar>(i);
		uchar* pt_src = src.ptr<uchar>(i);
		for (int j = n_channels * half_kernel_size; j < (n_cols_with_channels - n_channels * half_kernel_size); j++) {
			index = 0;
			pixel_sum = weight_sum = 0;

			// 内层kx,ky循环,空间域内滤波
			for (int kx = i - half_kernel_size; kx <= i + half_kernel_size; kx++) {
				uchar* pt_k_src = src.ptr<uchar>(kx);
				for (int ky = j - n_channels * half_kernel_size; ky <= (j + n_channels * half_kernel_size); ky += n_channels) {
					temp_bilateral_weight = distance_kernel[index++] * color_kernel[(int)abs(pt_src[j] - pt_k_src[ky])];
					weight_sum += temp_bilateral_weight;
					pixel_sum += (pt_k_src[ky] * temp_bilateral_weight); // 邻域某像素与中心点的双边权重乘积
				}
			}

			pixel_sum /= weight_sum; // 归一化
			pt_dst[j] = saturate_cast<uchar>(pixel_sum); //加权赋值
		}
	}
}

int main() {
	Mat src, dst, dst_cv;
	src = imread("test.png");
	if (!src.data) {
		cout << "loading failed" << endl;
		system("pause");
		return -1;
	}

	double distance_kernel[kernel_size * kernel_size];
	double color_kernel[256];
	GetDistanceWeight(distance_kernel);
	GetColorWeight(color_kernel);
	// 测试自己实现的双边滤波函数
	BilateralFilter(src, distance_kernel, color_kernel, dst);
	namedWindow("src");
	imshow("src", src);
	namedWindow("my bf");
	imshow("my bf", dst);
	// 测试opencv中的bilateralFilter， 观察与自实现的函数的异同
	bilateralFilter(src, dst_cv, 7, 50, 3);
	namedWindow("cv bf");
	imshow("cv bf", dst_cv);
	waitKey(0);
	return 0;
}

数学辅助理解

• 一维高斯分布：

高斯函数的$\sigma$设置越小则越窄。

• 二维高斯分布的公式为：
  
  在二维空间生成的曲面可视化图为：

$\sigma$太小的时候平滑效果不明显。

联合双边滤波（joint bilateral filter）

联合双边滤波相对于双边滤波来说，最大的特点就是引入了一幅引导影像。
具体可使用以下数学表达：
$$J_p=\frac{1}{k_p}\sum _{q\in \Omega }{I}_{q}f(\Vert p-q\Vert )g\left(\Vert {\tilde{I}}_p-{\tilde{I}}_q\Vert \right)$$
其中${\tilde{I}}_p$以及${\tilde{I}}_q$就是引导影像上的像素灰度值。
在opencv的contrib模块中，也提供了联合双边滤波的API：jointBilateralFilter(引导图，待滤波的图，滤波后的图，像素邻域直径，灰度域sigma，空间域sigma)。
简单的代码实现思路可以为将双边滤波中灰度域权重的计算过程进行修改，取引导影像上的灰度而非待处理影像的灰度，在此不赘述。至于opencv中API调用的代码则为：

#include <opencv2/opencv.hpp>
#include <ximgproc.hpp>
 
int main()
{
	cv::Mat src = cv::imread("data/dp.png", 1); // 原始带噪声的深度图
	cv::Mat joint = cv::imread("data/teddy.png", 0);
	cv::Mat dst;
	cv::ximgproc::jointBilateralFilter(joint, src, dst, -1, 3, 9);
  
	imshow("src", src);
	imshow("joint", joint);
	imshow("jointBilateralFilter", dst);
	cv::waitKey(0);
 
    return 0;
}

参考链接

https://blog.csdn.net/panda1234lee/article/details/52839205

写在最后

我们最近创建了一个“三维重建技术动向与商业落地”的知识星球，这个星球汇聚了来自985和国际顶级学府的专家和学者，他们分享了最新的三维重建技术和商业应用的前沿知识和经验。如果你对三维重建领域感兴趣，那么这个知识星球是你不可错过的。通过加入这个知识星球，你可以学习到最新的三维重建技术和商业应用，提高自己的技能和能力。同时，如果你是一个三维重建领域的专家，你也可以在这个知识星球上分享自己的知识和经验，让更多的人受益。我们会追踪最新的AIGC与3D的技术，并试图从投资人、技术人、产品人以及用户的视角提出一些看法。加入知识星球，让我们一起探索三维重建领域的商业落地想法和前沿知识！如果你想加入这个知识星球，可以添加我的微信号（zhuanshanqwer），我可以免费为你提供名额。