双边滤波原理浅析

前言

双边滤波（Bilateral filtering）是一种结合了图像像素空间位置关系和像素值大小的非线性组合的保边滤波算法，无需迭代计算，操作简单。从原理上讲，图像领域的滤波通常意义是希望去除噪声，使图像更加平滑。而传统的均值滤波、高斯滤波等只考虑空域，认为临近像素应该拥有相近的特性，这种假设在图像中的Edge处则不满足，如果仍用这种方式滤波则必定会导致边缘被模糊掉。为了达到保边效果，双边滤波通过结合空域和值域，在保证滤波效果的同时保持了边缘。对于双边滤波实际的应用场景，可用于noise reduction、 HDR compression等。部分情况下可以用于人脸磨皮美颜。

原论文《Bilateral Filtering for Gray and Color Images》摘要：
Bilateral filtering smooths images while preserving edges, by means of a nonlinear combination of nearby image values. The method is noniterative, local, and simple. It combines gray levels or colors based on both their geometric closeness and their photometric similarity, and prefers near values to distant values in both domain and range. In contrast with filters that operate on the three bands of a color image separately, a bilateral filter can enforce the perceptual metric underlying the CIE-Lab color space, and smooth colors and preserve edges in a way that is tuned to human perception. Also, in contrast with standard filtering, bilateral filtering produces no phantom colors along edges in color images, and reduces phantom colors where they appear in the original image.

原理

原论文在Gray和Color Image情况下均进行了讨论，这里只对Gray Image进行说明。Color Image一般会转为 CIE-Lab 颜色空间处理。

双边滤波实际上是由一个low-pass Domain filter和一个Range filter组合而成，其中Domain filter表示的是像素点空域相关性，即某像素点的临近点对该像素点的影响。定义如下：

其中$c(\xi ,$ x $)$是度量输入图像中某像素点 x的临近点ξ与x之间的空间关系。f (ξ)是输入图像的像素值，h (x)是输出的像素值。

$k_d$ ( x $)$是归一化参数，如果这个Domain filter是一个Shift-invariant system（平移不变系统，可以简单的认为是一个固定数值的kernel），那么$c(\xi ,$ x $)$只与 ξ - x 有关，$k_d$ ( x $)$是一个常数。

Range filter（值域滤波器）拥有与Domain filter相似的定义：

需要注意的是$c(\xi ,$ x $)$变成了$s($ f(ξ), f(x) $)$，也就是与像素值差值有关。

归一化参数：

将 Domain filter和Range filter组合起来形成Combined filtering，实现局部像素的空域和值域相关性。这就是bilateral filtering的含义。

原论文中带边缘情况的权重及最终滤波效果示意图如下：

论文《A Gentle Introduction to Bilateral Filtering and its Applications》中对该部分的描述彩图更加便于理解：

可以看到作用在X及相邻像素的权重在边界有很明显的分界，从而几乎只对X所属的边缘一侧的像素点进行加权。

实现

OpenCV中对双边滤波有较为高效的CUDA实现，代码如下：

__global__ void bilateral_kernel(const PtrStepSz<T> src, PtrStep<T> dst, const B b, const int ksz, const float sigma_spatial2_inv_half, const float sigma_color2_inv_half)
        {
            typedef typename TypeVec<float, VecTraits<T>::cn>::vec_type value_type;

            int x = threadIdx.x + blockIdx.x * blockDim.x;
            int y = threadIdx.y + blockIdx.y * blockDim.y;

            if (x >= src.cols || y >= src.rows)
                return;

            value_type center = saturate_cast<value_type>(src(y, x));

            value_type sum1 = VecTraits<value_type>::all(0);
            float sum2 = 0;

            int r = ksz / 2;
            float r2 = (float)(r * r);

            int tx = x - r + ksz;
            int ty = y - r + ksz;

            if (x - ksz/2 >=0 && y - ksz/2 >=0 && tx < src.cols && ty < src.rows)
            {
                for (int cy = y - r; cy < ty; ++cy)
                    for (int cx = x - r; cx < tx; ++cx)
                    {
                        float space2 = (x - cx) * (x - cx) + (y - cy) * (y - cy);
                        if (space2 > r2)
                            continue;

                        value_type value = saturate_cast<value_type>(src(cy, cx));

                        float weight = ::exp(space2 * sigma_spatial2_inv_half + sqr(norm_l1(value - center)) * sigma_color2_inv_half);
                        sum1 = sum1 + weight * value;
                        sum2 = sum2 + weight;
                    }
            }
            else
            {
                for (int cy = y - r; cy < ty; ++cy)
                    for (int cx = x - r; cx < tx; ++cx)
                    {
                        float space2 = (x - cx) * (x - cx) + (y - cy) * (y - cy);
                        if (space2 > r2)
                            continue;

                        value_type value = saturate_cast<value_type>(b.at(cy, cx, src.data, src.step));

                        float weight = ::exp(space2 * sigma_spatial2_inv_half + sqr(norm_l1(value - center)) * sigma_color2_inv_half);

                        sum1 = sum1 + weight * value;
                        sum2 = sum2 + weight;
                    }
            }
            dst(y, x) = saturate_cast<T>(sum1 / sum2);
        }

主要有两个参数sigma_spatial和sigma_color需要设置，对应在CUDA函数输入中则是sigma_spatial2_inv_half 和sigma_color2_inv_half。

sigma_spatial表示坐标空间中滤波器的sigma值，如果该参数的值较大，则说明颜色相近的较远的像素将相互影响，从而使更大的区域中足够相似的颜色获取相同的颜色。

sigma_color表示值域空间中滤波器的sigma值，决定了多少差值之内的像素会被计算进来，如果该参数的值越大，表明该像素邻域内有越宽广的颜色会被混合到一起，产生较大的半相等颜色区域。

核心权重计算代码：

float weight = ::exp(space2 * sigma_spatial2_inv_half + sqr(norm_l1(value - center)) * sigma_color2_inv_half);

使用

OpenCV的双边滤波CUDA实现的Python调用方式如下：

# 参数设置
kernel_size=0
sigma_color=10
sigma_spatial=3

cv2.cuda.bilateralFilter(
            src=src_gmat,
            dst=dst_gmat,
            kernel_size=kernel_size,
            sigma_color=sigma_color,
            sigma_spatial=sigma_spatial)

扩展讨论

Bilateral Filter可能存在梯度翻转现象（ 图片引自参考资料[9] ）：
关于该问题在论文《Guided Image Filtering》中的解释：

参考资料

[1] Bilateral Filtering for Gray and Color Images
[2] 双边滤波原理（Bilateral Filtering）
[3] GitHub - opencv_contrib/modules/cudaimgproc/src/bilateral_filter.cpp
[4] GitHub - opencv_contrib/modules/cudaimgproc/src/cuda/bilateral_filter.cu
[5] GitHub - aashikgowda / Bilateral-Filter-CUDA
[6] opencv 中的双边滤波用法总结(10)
[7] opencv学习(二十二)之双边滤波bilateralFilter
[8] Bilateral Filtering: Theory and Applications
[9] 导向滤波原理（Guided Filter）
[10] A Gentle Introduction to Bilateral Filtering and its Applications
[11] 论文《Guided Image Filtering》