图像处理之双边滤波介绍与源码实现
1 双边滤波简介
双边滤波（Bilateral filter）是一种非线性的滤波方法，是结合图像的空间邻近度和像素值相似度的一种折衷处理，同时考虑空域信息和灰度相似性，达到保边去噪的目的。具有简单、非迭代、局部的特点。
双边滤波器的好处是可以做边缘保存（edge preserving），一般过去用的维纳滤波或者高斯滤波去降噪，都会较明显地模糊边缘，对于高频细节的保护效果并不明显。双边滤波器顾名思义比高斯滤波多了一个高斯方差sigma－d，它是基于空间分布的高斯滤波函数，所以在边缘附近，离的较远的像素不会太多影响到边缘上的像素值，这样就保证了边缘附近像素值的保存。但是由于保存了过多的高频信息，对于彩色图像里的高频噪声，双边滤波器不能够干净的滤掉，只能够对于低频信息进行较好的滤波。
2 双边滤波原理
滤波算法中，目标点上的像素值通常是由其所在位置上的周围的一个小局部邻居像素的值所决定。在2D高斯滤波中的具体实现就是对周围的一定范围内的像素值分别赋以不同的高斯权重值，并在加权平均后得到当前点的最终结果。而这里的高斯权重因子是利用两个像素之间的空间距离（在图像中为2D）关系来生成。通过高斯分布的曲线可以发现，离目标像素越近的点对最终结果的贡献越大，反之则越小。其公式化的描述一般如下所述：
其中的c即为基于空间距离的高斯权重，而用来对结果进行单位化。
高斯滤波在低通滤波算法中有不错的表现，但是其却有另外一个问题，那就是只考虑了像素间的空间位置上的关系，因此滤波的结果会丢失边缘的信息。这里的边缘主要是指图像中主要的不同颜色区域（比如蓝色的天空，黑色的头发等），而Bilateral就是在Gaussian blur中加入了另外的一个权重分部来解决这一问题。Bilateral滤波中对于边缘的保持通过下述表达式来实现：
其中的s为基于像素间相似程度的高斯权重，同样用来对结果进行单位化。对两者进行结合即可以得到基于空间距离、相似程度综合考量的Bilateral滤波：
上式中的单位化分部综合了两种高斯权重于一起而得到，其中的c与s计算可以详细描述如下：
且有
且有
上述给出的表达式均是在空间上的无限积分，而在像素化的图像中当然无法这么做，而且也没必要如此做，因而在使用前需要对其进行离散化。而且也不需要对于每个局部像素从整张图像上进行加权操作，距离超过一定程度的像素实际上对当前的目标像素影响很小，可以忽略的。限定局部子区域后的离散化公就可以简化为如下形式：
上述理论公式就构成了Bilateral滤波实现的基础。为了直观地了解高斯滤波与双边滤波的区别，我们可以从下列图示中看出依据。假设目标源图像为下述左右区域分明的带有噪声的图像（由程序自动生成），蓝色框的中心即为目标像素所在的位置，那么当前像素处所对应的高斯权重与双边权重因子3D可视化后的形状如后边两图所示：
左图为原始的噪声图像；中间为高斯采样的权重；右图为Bilateral采样的权重。从图中可以看出Bilateral加入了相似程度分部以后可以将源图像左侧那些跟当前像素差值过大的点给滤去，这样就很好地保持了边缘。为了更加形象地观察两者间的区别，使用Matlab将该图在两种不同方式下的高度图3D绘制出来，如下：
上述三图从左到右依次为：双边滤波，原始图像，高斯滤波。从高度图中可以明显看出Bilateral和Gaussian两种方法的区别，前者较好地保持了边缘处的梯度，而在高斯滤波中，由于其在边缘处的变化是线性的，因而就使用连累的梯度呈现出渐变的状态，而这表现在图像中的话就是边界的丢失（图像的示例可见于后述）。
3 双边滤波实现步骤
有了上述理论以后实现Bilateral Filter就比较简单了，其实它也与普通的Gaussian Blur没有太大的区别。这里主要包括3部分的操作: 基于空间距离的权重因子生成；基于相似度的权重因子的生成；最终filter颜色的计算。
3.1Spatial Weight
这就是通常的Gaussian Blur中使用的计算高斯权重的方法，其主要通过两个pixel之间的距离并使用如下公式计算而来：
3.1Similarity Weight
与基于距离的高斯权重计算类似，只不过此处不再根据两个pixel之间的空间距离，而是根据其相似程度（或者两个pixel的值之间的距离）。
其中的表示两个像素值之间的距离，可以直接使用其灰度值之间的差值或者RGB向量之间的欧氏距离。
3.4 Color Filtering
其中的相似度分部的权重s主要根据两个Pixel之间的颜色差值计算面来。对于灰度图而言，这个差值的范围是可以预知的，即[-255, 255]，因而为了提高计算的效率我们可以将该部分权重因子预计算生成并存表，在使用时快速查询即可。使用上述实现的算法对几张带有噪声的图像进行滤波后的结果如下所示：
3 双边滤波源码实现
1 #include <time.h> 2 #include <stdio.h> 3 #include <stdlib.h> 4 #include <math.h> 5 #include "bf.h" 6 7 8 //-------------------------------------------------- 9 /** 10 双边滤波器图像去噪方法 11 \param pDest 去噪处理后图像buffer指针 12 \param pSrc 原始图像buffer指针 13 \paran nImgWid 图像宽 14 \param nImgHei 图像高 15 \param nSearchR 搜索区域半径 16 \param nSigma 噪声方差 17 \param nCoff DCT变换系数个数 18 19 return void 20 */ 21 //-------------------------------------------------- 22 void BilateralFilterRGB(float *pDest, BYTE *pSrc,int nImgWid, int nImgHei, int nSearchR, int nSigma, int nCoff) 23 { 24 25 BYTE *pSrcR = new BYTE[nImgHei * nImgWid]; // 新建图像缓冲区R 26 BYTE *pSrcG = new BYTE[nImgHei * nImgWid]; // 新建图像缓冲区G 27 BYTE *pSrcB = new BYTE[nImgHei * nImgWid]; // 新建图像缓冲区B 28 29 30 for(int i = 0;i < nImgHei; i++) 31 { 32 for (int j = 0; j < nImgWid; j++) 33 { 34 pSrcR[i * nImgWid + j] = pSrc[i * nImgWid * 3 + j * 3 + 0]; 35 pSrcG[i * nImgWid + j] = pSrc[i * nImgWid * 3 + j * 3 + 1]; 36 pSrcB[i * nImgWid + j] = pSrc[i * nImgWid * 3 + j * 3 + 2]; 37 } 38 } 39 40 41 float* pDestR = new float[nImgHei * nImgWid]; // 去噪后的图像数据R通道 42 float* pDestG = new float[nImgHei * nImgWid]; // 去噪后的图像数据G通道 43 float* pDestB = new float[nImgHei * nImgWid]; // 去噪后的图像数据B通道 44 memset(pDestR,255,nImgHei * nImgWid * sizeof(float)); // 传递变量初始化 45 memset(pDestG,255,nImgHei * nImgWid * sizeof(float)); 46 memset(pDestB,255,nImgHei * nImgWid * sizeof(float)); 47 48 49 float *pII = new float[nImgHei * nImgWid ]; // 积分图 50 float *pWeights = new float[nImgHei * nImgWid ]; // 每一个窗的权重之和，用于归一化 51 52 53 float mDctc[MAX_RANGE]; // DCT变换系数 54 float mGker[MAX_RANGE]; // 高斯核函数的值 55 56 for (int i = 0; i< MAX_RANGE; ++i) 57 { 58 mGker[i] = exp(-0.5 * pow(i / nSigma, 2.0)); // 首先计算0-255的灰度值在高斯变换下的取值，然后存在mGker数组中 59 } 60 61 CvMat G = cvMat(MIN_RANGE, MAX_RANGE, CV_32F, mGker); // 转换成opencv中的向量 62 CvMat D = cvMat(MIN_RANGE, MAX_RANGE, CV_32F, mDctc); // 把mDctc系数数组转换成opencv中的向量形式 63 64 cvDCT(&G,&D,CV_DXT_ROWS); // 然后将高斯核进行离散的dct变换 65 mDctc[0] /= sqrt(2.0); // 离散余弦变换的第一个系数是1/1.414; 66 67 68 bf(pSrcR, pDestR, mDctc, pII, pWeights, nImgHei, nImgWid, nCoff, nSearchR); 69 bf(pSrcG, pDestG, mDctc, pII, pWeights, nImgHei, nImgWid, nCoff, nSearchR); 70 bf(pSrcB, pDestB, mDctc, pII, pWeights, nImgHei, nImgWid, nCoff, nSearchR); 71 72 73 for(int i=0; i < nImgHei;i++) 74 { 75 for (int j=0; j < nImgWid; j++) 76 { 77 pDest[i * nImgWid * 3 + j * 3 + 0] = pDestR[i * nImgWid + j]; 78 pDest[i * nImgWid * 3 + j * 3 + 1] = pDestG[i * nImgWid + j]; 79 pDest[i * nImgWid * 3 + j * 3 + 2] = pDestB[i * nImgWid + j]; 80 } 81 } 82 83 84 free(pSrcR); // 释放临时缓冲区 85 free(pSrcG); 86 free(pSrcB); 87 free(pDestR); 88 free(pDestG); 89 free(pDestB); 90 91 }
1 //-------------------------------------------------- 2 /** 3 双边滤波器图像去噪方法 4 \param pDest 去噪处理后图像buffer指针 5 \param pSrc 原始图像buffer指针 6 \paran nImgWid 图像宽 7 \param nImgHei 图像高 8 \param pII 积分图缓冲区 9 \param pWeights 归一化权重系数 10 \param nCoff DCT变换系数 11 \param nPatchWid 图像块宽度 12 13 return void 14 */ 15 //-------------------------------------------------- 16 void bf (BYTE *pSrc, float *pDest, float *pDctc, float *pII, float *pWeights, int nImgHei, int nImgWid, int nCoff, int nPatchWid) 17 { 18 if (pDest == NULL || pSrc ==NULL) 19 { 20 return; 21 } 22 23 float *cR = (float*) malloc((nCoff - 1) * MAX_RANGE * sizeof(float)); 24 float *sR = (float*) malloc((nCoff - 1) * MAX_RANGE * sizeof(float)); 25 float *dcR = (float*) malloc((nCoff - 1) * MAX_RANGE * sizeof(float)); 26 float *dsR = (float*) malloc((nCoff - 1) * MAX_RANGE * sizeof(float)); 27 28 29 // 余弦变换查找表 30 int fx = 0; 31 int ck = 0; 32 int r2 = pow(2.0*nPatchWid + 1.0, 2.0); // 这个是图像块的像素点的总个数 33 34 float c0 = pDctc[0]; // 这是DCT的第一个系数 35 float c0r2 = pDctc[0] * r2; // 这个用于初始化归一化权重之和 36 37 38 for (ck = 1; ck < nCoff; ck++) // 注意到这里面的系数并不包含C0，所以后面才有把C0额外的加上 39 { // 一个矩阵，用于存放各种系数和数值，便于查找 40 int ckr = (ck - 1) * MAX_RANGE; // 数组是从0开始的，所以这里减1 41 42 for (fx = 0; fx<MAX_RANGE; fx++) // fx就是图像的灰度值 43 { 44 float tmpfx = PI * float(fx) * float(ck) / MAX_RANGE; // ck其实就相当于那个余弦函数里面的t，fx相当于u，PI/MAX相当于前面的那个系数，这都是余弦变换相关的东西 45 46 cR[ckr + fx] = cos(tmpfx); // 存储余弦变换，这个用在空间域上 47 sR[ckr + fx] = sin(tmpfx); // 存储正弦变换 48 dcR[ckr + fx] = pDctc[ck] * cos(tmpfx); // 存储余弦变换和系数的乘积，这个则用在强度范围上 49 dsR[ckr + fx] = pDctc[ck] * sin(tmpfx); // 存储正弦变换和系数的乘积 50 } 51 } 52 53 54 float *pw = pWeights; // 新建一个归一化权重的中间变量进行数据的初始化 55 float *pwe = pWeights + nImgWid * nImgHei; // 限定最大位置 56 57 while (pw < pwe) 58 { 59 *pw++ = c0r2; // 赋初值，让它们都等于第一个DCT系数乘以图像块中总的像素点的个数 60 } 61 62 for (int ck = 1; ck < nCoff; ck++) 63 { 64 int ckr = (ck-1)*MAX_RANGE; // 数组是从0开始的，所以这里减1 65 66 add_lut(pWeights, pSrc, pII,cR + ckr, dcR + ckr, nImgHei, nImgWid, nPatchWid, nPatchWid); // add cos term to pWeights 67 add_lut(pWeights, pSrc, pII,sR + ckr, dsR + ckr, nImgHei, nImgWid, nPatchWid, nPatchWid); // add cos term to pWeights 68 69 } 70 71 ImgRectSumO(pDest, pSrc, pII, c0, nImgHei, nImgWid, nPatchWid, nPatchWid); //加上C0的变换之后，再初始化滤波后的函数 72 73 for (int ck = 1; ck < nCoff; ck++) 74 { 75 int ckr = (ck-1)*MAX_RANGE; 76 77 add_f_lut(pDest, pSrc, pII, cR + ckr, dcR + ckr, nImgHei, nImgWid, nPatchWid, nPatchWid); // add cosine term to dataf 78 add_f_lut(pDest, pSrc, pII, sR + ckr, dsR + ckr, nImgHei, nImgWid, nPatchWid, nPatchWid); // add sine term to dataf 79 80 } 81 82 83 float *pd = pDest + nPatchWid * nImgWid; 84 pw = pWeights + nPatchWid * nImgWid; 85 86 float *pdie = pd + nImgWid - nPatchWid; 87 float *pdend = pDest + (nImgHei - nPatchWid) * nImgWid; 88 89 while (pd < pdend) 90 { 91 pd += nPatchWid; //把边都给去掉了 92 pw += nPatchWid; //这个也是把边去掉了 93 94 while (pd < pdie) 95 { 96 *pd++ /= ( (*pw++)*255); // 之所以要除以255，是为了把它化到[0,1]之间，便于显示 97 } 98 99 pd += nPatchWid; // 把边都给去掉了 100 pw += nPatchWid; // 这个也是把边去掉了 101 pdie += nImgWid; // 过渡到下一行 102 } 103 104 free(cR); 105 free(sR); 106 free(dcR); 107 free(dsR); // 释放缓冲区 108 }
//-------------------------------------------------- /** 余弦函数首系数积分图 \param pDest 去噪处理后图像buffer指针 \param pSrc 原始图像buffer指针 \paran nImgWid 图像宽 \param nImgHei 图像高 \param pII 积分图缓冲区 \param c0 DCT变换第一个系数 \param nPatchWid 图像块宽度 \param nPatchHei 图像块高度 return void */ //-------------------------------------------------- void ImgRectSumO (float* pDest, const BYTE* pSrc, float* pII, float c0,const int nImgHei, const int nImgWid, const int nPatchWid, const int nPatchHei) { if (pDest == NULL || pSrc ==NULL) { return; } int ri1 = nPatchWid + 1; // 中心像素点 int rj1 = nPatchHei + 1; int ri21 = 2 * nPatchWid + 1; // 图像块的宽度 int rj21 = 2 * nPatchHei + 1; const BYTE *pSrcImg = pSrc; // 原始图像数据 const BYTE *pSrcImgEnd = pSrc + nImgHei * nImgWid; // 最大的像素点位置 const BYTE *pSrcImgWid = pSrcImg + nImgWid; // 第一行的最大的像素点位置，下面循环用的的变量 float *pIITmp = pII; // 积分图数据 float *ppII = pIITmp; // 积分图中间变量 *pIITmp++ = c0 * float( *pSrcImg++ ); // 第一个系数乘以图像数据 while (pSrcImg < pSrcImgWid) // 遍历第一行进行求积分图，其实这个是图像数据的积分图 { *pIITmp++ = (*ppII++) + c0 * float(*pSrcImg++); } while (pSrcImg < pSrcImgEnd) // 遍历所有的像素点的变量 { pSrcImgWid += nImgWid; // 进行第二行的循环 ppII = pIITmp; // 积分图中间变量 *pIITmp++ = c0 * float(*pSrcImg++); // 第二行第一个像素点的处理 while (pSrcImg < pSrcImgWid) { (*pIITmp++) = (*ppII++) + c0 * float(*pSrcImg++); // 求第二行的积分图 } float *pIIWid = pIITmp; // 当前位置像素点 pIITmp = pIITmp - nImgWid; // 上一行的起始点位置 float *pIIup = pIITmp - nImgWid; // 上上一行的起始点位置 while (pIITmp < pIIWid) // 遍历整个行的每一个数据 { (*pIITmp++) = (*pIITmp) + (*pIIup++); // 行与行之间的积分图 } } float *pres = pDest + ri1 * nImgWid; // 最小的行位置 float *pend = pDest + (nImgHei - nPatchWid) * nImgWid; // 最大的行位置 float *pii1 = NULL; float *pii2 = NULL; float *pii3 = NULL; float *pii4 = NULL; // 求积分图的四个点 pii1 = pII + ri21 * nImgWid + rj21; // 右下角 pii2 = pII + ri21 * nImgWid; // 左下角 pii3 = pII + rj21; // 右上角 pii4 = pII; // 左上角 while (pres < pend) { float *pe = pres + nImgWid - nPatchHei; pres += rj1; while (pres < pe) // 这个只是求图像数据的积分图 { (*pres++) = (*pii1++) - (*pii2++) - (*pii3++) + (*pii4++); } pres += nPatchHei; pii1 += rj21; pii2 += rj21; pii3 += rj21; pii4 += rj21; } }
//-------------------------------------------------- /** 余弦积分图加函数 \param pNomalWts 像素点的归一化权重矩阵 \param pSrc 原始图像buffer指针 \param pII 积分图 \paran nImgWid 图像宽 \param nImgHei 图像高 \param pII 积分图缓冲区 \param pCosMtx 余弦函数矩阵 \param pCosDctcMtx 余弦函数与DCT变换系数乘积矩阵 \param nPatchWid 图像块宽度 \param nPatchHei 图像块高度 return void */ //-------------------------------------------------- void add_lut (float* pNomalWts, const BYTE* pSrc, float* pII, float* pCosMtx, float* pCosDctcMtx, const int nImgHei, const int nImgWid, const int nPatchWid, const int nPatchHei) { int ri1 = nPatchWid + 1; // kernel相关 int rj1 = nPatchHei + 1; int ri21 = 2 * nPatchWid + 1; // 这是kernel的宽度和长度 int rj21 = 2 * nPatchHei + 1; const BYTE *pi = pSrc; // 这个是图像数据 float *pii = pII; // 这个是中间的缓冲变量，积分图的缓冲区 const BYTE *piw = pi + nImgWid; // 也是赋初值的时候使用的中间变量 float *pii_p = pii; // 直线积分图指针的指针 *pii++ = pCosMtx[*pi++]; //图像数据值所对应的查找表中的余弦或者正弦变换 while (pi < piw) { *pii++ = (*pii_p++) + pCosMtx[*pi++]; // 第一行的积分图 } const BYTE *piend = pSrc + nImgHei * nImgWid; // 限定循环位置 while (pi < piend) // 遍历整个图像数据 { piw += nImgWid; // 第二行 pii_p = pii; // 指向积分图指针的指针 *pii++ = pCosMtx[*pi++]; // 获取第一个图像数据值所对应的查找表中的余弦或者正弦变换 while (pi < piw) { (*pii++) = (*pii_p++) + pCosMtx[*pi++]; // 求这一行的积分图 } float *piiw = pii; pii = pii - nImgWid; float *pii_p1 = pii - nImgWid; while (pii < piiw) { (*pii++) = (*pii) + (*pii_p1++); // 行与行之间求积分图 } } float *pres = pNomalWts + ri1 * nImgWid; // 定位要处理点的像素的那一行 float *pend = pNomalWts + (nImgHei - nPatchWid) * nImgWid; // 限定了边界 float *pii1 = NULL; float *pii2 = NULL; float *pii3 = NULL; float *pii4 = NULL; pii1 = pII + ri21 * nImgWid + rj21; // 获得积分图中那个最大的数据右下角 pii2 = pII + ri21 * nImgWid; // 积分图左下角 pii3 = pII + rj21; // 积分图右上角 pii4 = pII; // 积分图左上角 pi = pSrc + ri1 * nImgWid; // 定位要处理的像素的位置 while (pres < pend) // 设定高度做循环 { float *pe = pres + nImgWid - nPatchHei; // 限定了宽度的范围 pres = pres + rj1; // 定位了要处理的像素点的归一化系数矩阵的位置 pi = pi + rj1; // 定位了要处理的像素点 while (pres < pe) // 设定宽度做循环 { (*pres++) = (*pres) + pCosDctcMtx[*pi++] * ( (*pii1++) - (*pii2++) - (*pii3++) + (*pii4++) ); // 得到的就是归一化系数矩阵之和 } pres += nPatchHei; pi += nPatchHei; pii1 += rj21; pii2 += rj21; pii3 += rj21; pii4 += rj21; // 略过不处理的边界区域 } }
//-------------------------------------------------- /** 积分图 \param pDest 像素点的归一化权重矩阵 \param pSrc 原始图像buffer指针 \param pII 积分图 \paran nImgWid 图像宽 \param nImgHei 图像高 \param pII 积分图缓冲区 \param pCosMtx 余弦函数矩阵 \param pCosDctcMtx 余弦函数与DCT变换系数乘积矩阵 \param nPatchWid 图像块宽度 \param nPatchHei 图像块高度 return void */ //-------------------------------------------------- void add_f_lut (float* pDest,const BYTE* pSrc,float* pII,float* pCosMtx,float* pCosDctcMtx,const int nImgHei, const int nImgWid,const int nPatchWid, const int nPatchHei) { int ri1 = nPatchWid + 1; // kernel相关，目标是指向中心像素点 int rj1 = nPatchHei + 1; // 目标是指向中心像素点 int ri21 = 2 * nPatchWid + 1; // 这是kernel的宽度和长度 int rj21 = 2 * nPatchHei + 1; // 其实就是指向搜索区域的右下角，也就是最大位置的那个像素点，就是边界了 const BYTE *pi = pSrc; // 这个是图像数据,图像的灰度值 float *pii = pII; // 这是积分图 const BYTE *piw = pi + nImgWid; // 指向第一行的末尾位置，用于循环控制 float *pii_p = pii; // 指向积分图指针的指针 *pii++ = float(*pi) * pCosMtx[*pi]; // 灰度值乘以余弦系数 ++pi; while (pi < piw) // 第一行遍历 { *pii++ = (*pii_p++) + float(*pi) * pCosMtx[*pi]; ++pi; } const BYTE *piend = pSrc + nImgHei * nImgWid; while (pi < piend) { piw += nImgWid; pii_p = pii; *pii++ = float(*pi) * pCosMtx[*pi]; ++pi; while (pi < piw) { (*pii++) = (*pii_p++) + float(*pi) * pCosMtx[*pi]; ++pi; } float *piiw = pii; pii = pii - nImgWid; // 上一行起点 float *pii_p1 = pii - nImgWid; // 上上一行的起始点 while (pii < piiw) // 循环一行 { (*pii++) += (*pii_p1++); //其实就是在列的位置上加上上一行 } } float *pres = pDest + ri1*nImgWid; float *pend = pDest + (nImgHei - nPatchWid) * nImgWid; float *pii1 = NULL; float *pii2 = NULL; float *pii3 = NULL; float *pii4 = NULL; pii1 = pII + ri21 * nImgWid + rj21; // 积分图搜索区域最大位置的元素，堪称右下角的元素 pii2 = pII + ri21 * nImgWid; // 可以看成搜索区域中左下角的像素位置 pii3 = pII + rj21; // 搜索区域右上角像素位置 pii4 = pII; // 搜索区域左上角像素位置 pi = pSrc + ri1 * nImgWid; // 定位要处理的像素的位置的那一行 while (pres < pend) { float *pe = pres + nImgWid - nPatchHei; //限定了宽度的范围 pres = pres + rj1; //定位了要处理的像素点的归一化系数矩阵的位置 pi = pi + rj1; //定位了要处理的像素点 while (pres < pe) //遍历整个行 { (*pres++) = (*pres) + pCosDctcMtx[*pi++] * ( (*pii1++) - (*pii2++) - (*pii3++) + (*pii4++) ); //这个其实是计算整个像素块 } pres += nPatchHei; pi += nPatchHei; pii1 += rj21; pii2 += rj21; pii3 += rj21; pii4 += rj21; } }
