双边滤波器解析与代码

81 阅读 0 评论 54 点赞

我是靠谱客的博主火星上斑马，最近开发中收集的这篇文章主要介绍双边滤波器解析与代码，觉得挺不错的，现在分享给大家，希望可以做个参考。

概述

大家熟知的滤波器，如高斯滤波器、中值滤波器都是为了平滑图像以抑制图形噪声，可是这带来一个问题：边缘也被平滑掉了。那么要想在一定程度上避免这个问题，我们就可以利用双边滤波器。

双边滤波器（Bilateral Filtering）是一种可以保边去噪的滤波器。

双边滤波这个概念最初由Tomasi和Manduchi在文献[1]提出，在处理相邻各像素值的灰度值或彩色信息时，不仅考虑到几何上的邻近关系，也考虑到了亮度上的相似性，通过对二者的非线性组合，自适应滤波后得到平滑图像。这样处理过的图像在滤除噪声的同时还能够很好地保持图像的边缘信息。

简单地讲：双边滤波器类似于高斯滤波器，它也是给每一个邻域像素分配一个加权系数。不过，这些加权系数包含两个部分, 第一部分加权方式与高斯滤波一样，第二部分的权重则取决于该邻域像素与当前像素的灰度差值。

有些同学可能忘了高斯滤波器，猛击这里。高斯滤波是将输入数组的每一个像素点与高斯内核卷积，将卷积和当作输出像素。

上面那个图就是一维高斯函数的模样，可以发现中间像素的加权系数是最大的，周边像素的加权系数随着它们远离中间像素的距离增大而逐渐减小。

不过，图像是二维的，所以应用的是二维高斯函数，公式如下所示：

现在回归正题：

双边滤波器不像普通的高斯滤波器/卷积低通滤波器，只考虑了位置对中心像素的影响，它还考虑了卷积核中像素与中心像素之间相似程度的影响。

看如下公式：

上式中f是经过双边滤波器的输出图像，g是原始输入图像，w是权重系数。到这里为止，其实跟高斯滤波器的思路一样，就是通过加权平均求值。

双边滤波器的主要区别在于它的权重系数，它的权重系数由两部分组成。一个是空间邻近度因子ws、一个是亮度相似度因子wr。它是ws和wr的乘积。

即w的公式如下：

在这个权重系数中同时考虑了空间域和值域的差别。

ws随着像素点与中心点之间欧几里德距离的增大而减小，wr随着两像素亮度值之差的增大而减小。

在图像平滑的区域，邻域内像素亮度值相差不大，双边滤波器转化为高斯低通滤波器；

在图像变化剧烈的区域，滤波器利用边缘点附近亮度值相近的像素点的亮度值平均代替原亮度值。

一个简单的示例如下：

上图中(a)表示一个边缘，(b)表示权重系数，(c)是滤波后的结果。

【PS:这个双边滤波器的OpenCV代码看得很痛苦啊，源代码文件放在：opencvsourcesmodulesgpusrc里面，自然就是应用了GPU计算，可惜对GPU没有研究。后面我再贴出了PCL库的双边滤波器代码】

还有一种方法同样可以保边滤波：BEEPS，具体查看文献[5].

OpenCV的代码如下：

#include "precomp.hpp"
using namespace cv;
using namespace cv::gpu;
using namespace std;
#if !defined (HAVE_CUDA) || defined (CUDA_DISABLER)
cv::gpu::DisparityBilateralFilter::DisparityBilateralFilter(int, int, int) { throw_nogpu(); }
cv::gpu::DisparityBilateralFilter::DisparityBilateralFilter(int, int, int, float, float, float) { throw_nogpu(); }
void cv::gpu::DisparityBilateralFilter::operator()(const GpuMat&, const GpuMat&, GpuMat&, Stream&) { throw_nogpu(); }
#else /* !defined (HAVE_CUDA) */
namespace cv { namespace gpu { namespace device
{
namespace disp_bilateral_filter
{
void disp_load_constants(float* table_color, PtrStepSzf table_space, int ndisp, int radius, short edge_disc, short max_disc);
template<typename T>
void disp_bilateral_filter(PtrStepSz<T> disp, PtrStepSzb img, int channels, int iters, cudaStream_t stream);
}
}}}
using namespace ::cv::gpu::device::disp_bilateral_filter;
namespace
{
const float DEFAULT_EDGE_THRESHOLD = 0.1f;
const float DEFAULT_MAX_DISC_THRESHOLD = 0.2f;
const float DEFAULT_SIGMA_RANGE = 10.0f;
inline void calc_color_weighted_table(GpuMat& table_color, float sigma_range, int len)
{
Mat cpu_table_color(1, len, CV_32F);
float* line = cpu_table_color.ptr<float>();
for(int i = 0; i < len; i++)
line[i] = static_cast<float>(std::exp(-double(i * i) / (2 * sigma_range * sigma_range)));
table_color.upload(cpu_table_color);
}
inline void calc_space_weighted_filter(GpuMat& table_space, int win_size, float dist_space)
{
int half = (win_size >> 1);
Mat cpu_table_space(half + 1, half + 1, CV_32F);
for (int y = 0; y <= half; ++y)
{
float* row = cpu_table_space.ptr<float>(y);
for (int x = 0; x <= half; ++x)
row[x] = exp(-sqrt(float(y * y) + float(x * x)) / dist_space);
}
table_space.upload(cpu_table_space);
}
template <typename T>
void disp_bilateral_filter_operator(int ndisp, int radius, int iters, float edge_threshold,float max_disc_threshold,
GpuMat& table_color, GpuMat& table_space,
const GpuMat& disp, const GpuMat& img, GpuMat& dst, Stream& stream)
{
short edge_disc = max<short>(short(1), short(ndisp * edge_threshold + 0.5));
short max_disc = short(ndisp * max_disc_threshold + 0.5);
disp_load_constants(table_color.ptr<float>(), table_space, ndisp, radius, edge_disc, max_disc);
if (&dst != &disp)
{
if (stream)
stream.enqueueCopy(disp, dst);
else
disp.copyTo(dst);
}
disp_bilateral_filter<T>(dst, img, img.channels(), iters, StreamAccessor::getStream(stream));
}
typedef void (*bilateral_filter_operator_t)(int ndisp, int radius, int iters, float edge_threshold, float max_disc_threshold,
GpuMat& table_color, GpuMat& table_space,
const GpuMat& disp, const GpuMat& img, GpuMat& dst, Stream& stream);
const bilateral_filter_operator_t operators[] =
{disp_bilateral_filter_operator<unsigned char>, 0, 0, disp_bilateral_filter_operator<short>, 0, 0, 0, 0};
}
cv::gpu::DisparityBilateralFilter::DisparityBilateralFilter(int ndisp_, int radius_, int iters_)
: ndisp(ndisp_), radius(radius_), iters(iters_), edge_threshold(DEFAULT_EDGE_THRESHOLD), max_disc_threshold(DEFAULT_MAX_DISC_THRESHOLD),
sigma_range(DEFAULT_SIGMA_RANGE)
{
calc_color_weighted_table(table_color, sigma_range, 255);
calc_space_weighted_filter(table_space, radius * 2 + 1, radius + 1.0f);
}
cv::gpu::DisparityBilateralFilter::DisparityBilateralFilter(int ndisp_, int radius_, int iters_, float edge_threshold_,
float max_disc_threshold_, float sigma_range_)
: ndisp(ndisp_), radius(radius_), iters(iters_), edge_threshold(edge_threshold_), max_disc_threshold(max_disc_threshold_),
sigma_range(sigma_range_)
{
calc_color_weighted_table(table_color, sigma_range, 255);
calc_space_weighted_filter(table_space, radius * 2 + 1, radius + 1.0f);
}
void cv::gpu::DisparityBilateralFilter::operator()(const GpuMat& disp, const GpuMat& img, GpuMat& dst, Stream& stream)
{
CV_DbgAssert(0 < ndisp && 0 < radius && 0 < iters);
CV_Assert(disp.rows == img.rows && disp.cols == img.cols && (disp.type() == CV_8U || disp.type() == CV_16S) && (img.type() == CV_8UC1 || img.type() == CV_8UC3));
operators[disp.type()](ndisp, radius, iters, edge_threshold, max_disc_threshold, table_color, table_space, disp, img, dst, stream);
}
#endif /* !defined (HAVE_CUDA) */

下面的是PCL中的代码，比上面opencv直观，不用多说：

代码网址：http://www.pointclouds.org/documentation/tutorials/writing_new_classes.php#example-a-bilateral-filter

 #include <pcl/point_types.h>
#include <pcl/io/pcd_io.h>
#include <pcl/kdtree/kdtree_flann.h>
typedef pcl::PointXYZI PointT;
float G (float x, float sigma)
{
return exp (- (x*x)/(2*sigma*sigma));
}
int main (int argc, char *argv[])
{
std::string incloudfile = argv[1];
std::string outcloudfile = argv[2];
float sigma_s = atof (argv[3]);
float sigma_r = atof (argv[4]);
// Load cloud
pcl::PointCloud<PointT>::Ptr cloud (new pcl::PointCloud<PointT>);
pcl::io::loadPCDFile (incloudfile.c_str (), *cloud);
int pnumber = (int)cloud->size ();
// Output Cloud = Input Cloud
pcl::PointCloud<PointT> outcloud = *cloud;
// Set up KDTree
pcl::KdTreeFLANN<PointT>::Ptr tree (new pcl::KdTreeFLANN<PointT>);
tree->setInputCloud (cloud);
// Neighbors containers
std::vector<int> k_indices;
std::vector<float> k_distances;
// Main Loop
for (int point_id = 0; point_id < pnumber; ++point_id)
{
float BF = 0;
float W = 0;
tree->radiusSearch (point_id, 2 * sigma_s, k_indices, k_distances);
// For each neighbor
for (size_t n_id = 0; n_id < k_indices.size (); ++n_id)
{
float id = k_indices.at (n_id);
float dist = sqrt (k_distances.at (n_id));
float intensity_dist = abs (cloud->points[point_id].intensity - cloud->points[id].intensity);
float w_a = G (dist, sigma_s);
float w_b = G (intensity_dist, sigma_r);
float weight = w_a * w_b;
BF += weight * cloud->points[id].intensity;
W += weight;
}
outcloud.points[point_id].intensity = BF / W;
}
// Save filtered output
pcl::io::savePCDFile (outcloudfile.c_str (), outcloud);
return (0);
}

例程一个：

#include <iostream>
#include <opencv2opencv.hpp>
using namespace std;
using namespace cv;
int main()
{
Mat dst;
Mat input_img = imread("C:\Users\Scien_Lin\Desktop\IMG\cat1.jpeg");
namedWindow("Ori");
imshow("Ori",input_img);
if (!input_img.data)
{
cerr<<"Error: image read problem!"<<endl;
return -1;
}
int KERNEL_SIZE = 31;
for (int i = 1; i < KERNEL_SIZE; i = i + 2)
{
bilateralFilter(input_img,dst,i,i*2,i/2);
}
namedWindow("Bil");
imshow("Bil",dst);
waitKey(0);
return 0;
}

更多图像处理、机器视觉资源请关注博客：LinJM-机器视觉微博：林建民-机器视觉

参考文献：

[1] Tomasi C, Manduchi R. Bilateral filtering for gray and color images[C]//Computer Vision, 1998. Sixth International Conference on. IEEE, 1998: 839-846

[2] 靳明, 宋建中. 一种自适应的图像双边滤波方法[J]. 光电工程, 2004, 31(7): 65-68.

[3] 张志强, 王万玉. 一种改进的双边滤波算法[J]. 中国图象图形学报, 2009 (3): 443-447.

[4] Rachel Zhang 的博客http://blog.csdn.net/abcjennifer/article/details/7616663

[5] Thévenaz P, Sage D, Unser M. Bi-exponential edge-preserving smoother[J]. Image Processing, IEEE Transactions on, 2012, 21(9): 3924-3936.

[6] http://www.pointclouds.org/documentation/tutorials/writing_new_classes.php#example-a-bilateral-filter