本文介绍四种常见的图像滤波算法,并附上源码。图像滤波是一种非常重要的图像处理技术,现在大火的卷积神经网络其实也是滤波的一种,都是用卷积核去提取图像的特征模式。不过,传统的滤波,使用的卷积核是固定的参数,是由经验非常丰富的人去手动设计的,也称为手工特征。而卷积神经网络的卷积核参数初始时未知的,根据不同的任务由数据和神经网络反向传播算法去学习得到的参数,更能适应于不同的任务。
中值滤波器
中值滤波器是一种常用的非线性滤波器,其基本原理是:选择待处理像素的一个邻域中各像素值的中值来代替待处理的像素。主要功能使某像素的灰度值与周围领域内的像素比较接近,从而消除一些孤立的噪声点,所以中值滤波器能够很好的消除椒盐噪声。不仅如此,中值滤波器在消除噪声的同时,还能有效的保护图像的边界信息,不会对图像造成很大的模糊(相比于均值滤波器)。
中值滤波器的效果受滤波窗口尺寸的影响较大,在消除噪声和保护图像的细节存在着矛盾:滤波窗口较小,则能很好的保护图像中的某些细节,但对噪声的过滤效果就不是很好,因为实际中的噪声不可能只占一个像素位置;反之,窗口尺寸较大有较好的噪声过滤效果,但是会对图像造成一定的模糊。另外,根据中值滤波器原理,如果在滤波窗口内的噪声点的个数大于整个窗口内非噪声像素的个数,则中值滤波就不能很好的过滤掉噪声。
自适应中值滤波器
常规的中值滤波器,在噪声的密度不是很大的情况下,效果不错。但是当噪声出现的概率较高时,常规的中值滤波的效果就不是很好了。有一个选择就是增大滤波器的窗口大小,这虽然在一定程度上能解决上述的问题,但是会给图像造成较大的模糊。
常规的中值滤波器的窗口尺寸是固定大小不变的,就不能同时兼顾去噪和保护图像的细节。这时就要寻求一种改变,根据预先设定好的条件,在滤波的过程中,动态的改变滤波器的窗口尺寸大小,这就是
自适应中值滤波器 Adaptive Median Filter
。在滤波的过程中,自适应中值滤波器会根据预先设定好的条件,改变滤波窗口的尺寸大小,同时还会根据一定的条件判断当前像素是不是噪声,如果是则用邻域中值替换掉当前像素;不是,则不作改变。
自适应中值滤波器有三个目的:
-
-
-
尽可能的保护图像中细节信息,避免图像边缘的细化或者粗化。
自适应中值滤波算法描述
自适应滤波器不但能够滤除概率较大的椒盐噪声,而且能够更好的保护图像的细节,这是常规的中值滤波器做不到的。自适应的中值滤波器也需要一个矩形的窗口 ,和常规中值滤波器不同的是这个窗口的大小会在滤波处理的过程中进行改变(增大)。需要注意的是,滤波器的输出是一个像素值,该值用来替换点
处的像素值,点
是滤波窗口的中心位置。
在描述自适应中值滤波器时需要用到如下的符号:
自适应中值滤波器有两个处理过程,分别记为:
和
。
A
:
如果A1 > 0 且 A2 < 0,跳转到 B;
否则,增大窗口的尺寸 如果增大后窗口的尺寸
,则重复A过程。否则,输出 𝑍𝑚𝑒𝑑
B
:
如果B1 > 0 且 B2 <0则输出
,否则输出
自适应中值滤波原理说明
过程A的目的是确定当前窗口内得到中值 𝑍𝑚𝑒𝑑 是否是噪声。如果 𝑍𝑚𝑖𝑛
如果在过程A中,得到则 𝑍𝑚𝑒𝑑 不符合条件 𝑍𝑚𝑖𝑛
从上面分析可知,噪声出现的概率较低,自适应中值滤波器可以较快的得出结果,不需要去增加窗口的尺寸;反之,噪声的出现的概率较高,则需要增大滤波器的窗口尺寸,这也符合种中值滤波器的特点:噪声点比较多时,需要更大的滤波器窗口尺寸。
算法实现
有了算法的详细描述,借助于OpenCV对图像的读写,自适应中值滤波器实现起来也不是很困难。首先定义滤波器最小的窗口尺寸以及最大的窗口尺寸。要进行滤波处理,首先要扩展图像的边界,以便对图像的边界像素进行处理。copyMakeBorder根据选择的BorderTypes使用不同的值扩充图像的边界像素,具体可参考OpenCV的文档信息。下面就是遍历图像的像素,对每个像素进行滤波处理。需要注意一点,不论滤波器多么的复杂,其每次的滤波过程,都是值返回一个值,来替换掉当前窗口的中心的像素值。函数adpativeProcess就是对当前像素的滤波过程,其代码如下:
uchar adaptiveProcess(const Mat &im, int row,int col,int kernelSize,int maxSize)
{
vector pixels;
for (int a = -kernelSize / 2; a <= kernelSize / 2; a++)
for (int b = -kernelSize / 2; b <= kernelSize / 2; b++)
{
pixels.push_back(im.at(row + a, col + b));
}
sort(pixels.begin(), pixels.end());
auto min = pixels[0];
auto max = pixels[kernelSize * kernelSize - 1];
auto med = pixels[kernelSize * kernelSize / 2];
auto zxy = im.at(row, col);
if (med > min && med < max)
{
if (zxy > min && zxy < max)
return zxy;
else
return med;
}
else
{
kernelSize += 2;
if (kernelSize <= maxSize)
return adpativeProcess(im, row, col, kernelSize, maxSize);
else
return med;
}
}
有了上面这个函数,剩下的只需要对全部像素做一个遍历即可,更为完整的代码,请见我的Github地址:
https://github.com/zhangqizky/common-image-filteringgithub.com
高斯滤波
高斯滤波也是一种非常常见的滤波方法,其核的形式为:
其中
是图像中的点的坐标,
为标准差,高斯模板就是利用这个函数来计算的,x和y都是代表,以核中心点为坐标原点的坐标值。这里想说一下
的作用,当
比较小的时候,生成的高斯模板中心的系数比较大,而周围的系数比较小,这样对图像的平滑效果不明显。而当
比较大时,生成的模板的各个系数相差就不是很大,比较类似于均值模板,对图像的平滑效果比较明显。
高斯滤波没有特别多可说的,最主要的作用是滤除高斯噪声,即符合正态分布的噪声。
实现的方式有两种,第一种是按照公式暴力实现,代码如下:
void GaussianFilter(const Mat &src, Mat &dst, int ksize, double sigma)
{
CV_Assert(src.channels() || src.channels() == 3);
double **GaussianTemplate = new double *[ksize];
for(int i = 0; i < ksize; i++){
GaussianTemplate[i] = new double [ksize];
}
generateGaussianTemplate(GaussianTemplate, ksize, sigma);
int border = ksize / 2;
copyMakeBorder(src, dst, border, border, border, border, BORDER_CONSTANT);
int channels = dst.channels();
int rows = dst.rows - border;
int cols = dst.cols - border;
for(int i = border; i < rows; i++){
for(int j = border; j< cols; j++){
double sum[3] = {0};
for(int a = -border; a <= border; a++){
for(int b = -border; b <= border; b++){
if(channels == 1){
sum[0] += GaussianTemplate[border+a][border+b] * dst.at(i+a, j+b);
}else if(channels == 3){
Vec3b rgb = dst.at(i+a, j+b);
auto k = GaussianTemplate[border+a][border+b];
sum[0] += k * rgb[0];
sum[1] += k * rgb[1];
sum[2] += k * rgb[2];
}
}
}
for(int k = 0; k < channels; k++){
if(sum[k] < 0) sum[k] = 0;
else if(sum[k] > 255) sum[k] = 255;
}
if(channels == 1){
dst.at(i, j) = static_cast(sum[0]);
}else if(channels == 3){
Vec3b rgb = {static_cast(sum[0]), static_cast(sum[1]), static_cast(sum[2])};
dst.at(i, j) = rgb;
}
}
}
for(int i = 0; i < ksize; i++)
delete[] GaussianTemplate[i];
delete[] GaussianTemplate;
}
当核比较大时,高斯滤波会比较费时,此时可以使用分离X和Y通道的形式来实现可分离的高斯滤波。为什么可以这么做?因为高斯函数中没有
这样的耦合项,即x和y是相对独立的,此时就可以将两个维度分离开来。
void separateGaussianFilter(const Mat &src, Mat &dst, int ksize, double sigma){
assert(src.channels()==1 || src.channels() == 3);
double *matrix = new double[ksize];
double sum = 0;
int origin = ksize / 2;
for(int i = 0; i < ksize; i++){
double g = exp(-(i-origin) * (i-origin) / (2 * sigma * sigma));
sum += g;
matrix[i] = g;
}
for(int i = 0; i < ksize; i++) matrix[i] /= sum;
int border = ksize / 2;
copyMakeBorder(src, dst, border, border, border, border, BORDER_CONSTANT);
int channels = dst.channels();
int rows = dst.rows - border;
int cols = dst.cols - border;
for(int i = border; i < rows; i++){
for(int j = border; j < cols; j++){
double sum[3] = {0};
for(int k = -border; k<=border; k++){
if(channels == 1){
sum[0] += matrix[border + k] * dst.at(i, j+k);
}else if(channels == 3){
Vec3b rgb = dst.at(i, j+k);
sum[0] += matrix[border+k] * rgb[0];
sum[1] += matrix[border+k] * rgb[1];
sum[2] += matrix[border+k] * rgb[2];
}
}
for(int k = 0; k < channels; k++){
if(sum[k] < 0) sum[k] = 0;
else
if(sum[k] > 255) sum[k] = 255;
}
if(channels == 1)
dst.at(i, j) = static_cast(sum[0]);
else if(channels == 3){
Vec3b rgb = {static_cast(sum[0]), static_cast(sum[1]), static_cast(sum[2])};
dst.at(i, j) = rgb;
}
}
}
for(int i = border; i < rows; i++){
for(int j = border; j < cols; j++){
double sum[3] = {0};
for(int k = -border; k<=border; k++){
if(channels == 1){
sum[0] += matrix[border + k] * dst.at(i+k, j);
}else if(channels == 3){
Vec3b rgb = dst.at(i+k, j);
sum[0] += matrix[border+k] * rgb[0];
sum[1] += matrix[border+k] * rgb[1];
sum[2] += matrix[border+k] * rgb[2];
}
}
for(int k = 0; k < channels; k++){
if(sum[k] < 0) sum[k] = 0;
else if(sum[k] > 255) sum[k] = 255;
}
if(channels == 1)
dst.at(i, j) = static_cast(sum[0]);
else if(channels == 3){
Vec3b rgb = {static_cast(sum[0]), static_cast(sum[1]), static_cast(sum[2])};
dst.at(i, j) = rgb;
}
}
}
delete [] matrix;
}
同样,完整的代码请见:
https://github.com/zhangqizky/common-image-filtering/tree/maingithub.com
双边滤波
双边滤波是一种非线性滤波方法,是结合了图像的邻近度和像素值相似度的一种折中,在滤除噪声的同时可以保留原图的边缘信息。整个双边滤波是由两个函数构成:一个函数是由空间距离决定的滤波器系数,另外一个诗由像素差值决定的滤波器系数。整个双边滤波的公式如下:
其中权重系数
取决于定义域核:
和值域核
的乘积。其中定义域核影响的是空间位置,如果把图像看成一个二维函数,那么定义域就是图像的坐标,值域就是该坐标处对应的像素值。定义域核就是普通的高斯核,全局使用一个就可以。但值域核是需要对每个像素点滑动进行计算的。
那么如何理解双边滤波呢
高斯滤波的滤波核的意义是,滤波后的像素值等于窗口内的像素值的加权平均值,权值系数是符合高斯分布,距离该点越近,权值越大。但是没有考虑像素值与当前点的差距。现在加上值域核,意义就在,滤波后当前点的像素值还会受到领域内像素值与自身的像素值差异的影响,不仅仅是距离来决定。这样,在平缓的区域里,由于像素值差异非常小,则值域的权重趋向于1,所以双边滤波就近似为高斯滤波。而在边缘区域中,由于像素值的差异比较大,则值域核趋向于0,权重下降,即当前像素受到领域内像素影响比较小,从而保留了边缘信息。
双边滤波的代码
opencv中提供了bilateralFilter()函数来实现双边滤波操作,其原型如下:
void cv::bilateralFilter(InputArray src,
OutputArray dst,
int d,
double sigmaColor,
double sigmaSpace,
int borderType = BORDER_DEFAULT
)
-
InputArray src: 输入图像,可以是Mat类型,图像必须是8位整型或浮点型单通道、三通道的图像。
-
OutputArray dst: 输出图像,和原图像有相同的尺寸和类型。
-
int d: 表示在过滤过程中每个像素邻域的直径范围。如果这个值是非正数,则函数会从第五个参数sigmaSpace计算该值。
-
double sigmaColor: 颜色空间过滤器的
值,这个参数的值月大,表明该像素邻域内有越宽广的颜色会被混合到一起,产生较大的半相等颜色区域。(这个参数可以理解为值域核的
和
)
-
double sigmaSpace: 坐标空间中滤波器的sigma值,如果该值较大,则意味着越远的像素将相互影响,从而使更大的区域中足够相似的颜色获取相同的颜色。当d>0时,d指定了邻域大小且与sigmaSpace无关,否则d正比于sigmaSpace. (这个参数可以理解为空间域核的
和
)
-
int borderType=BORDER_DEFAULT: 用于推断图像外部像素的某种边界模式,有默认值BORDER_DEFAULT.
具体代码如下:
#include
#include
using namespace std;
using namespace cv;
const int g_ndMaxValue = 100;
const int g_nsigmaColorMaxValue = 200;
const int g_nsigmaSpaceMaxValue = 200;
int g_ndValue;
int g_nsigmaColorValue;
int g_nsigmaSpaceValue;
Mat g_srcImage;
Mat g_dstImage;
void on_bilateralFilterTrackbar(int, void*);
int main()
{
g_srcImage = imread("lena.jpg");
if(g_srcImage.empty())
{
cout << "图像加载失败!" << endl;
return -1;
}
else
cout << "图像加载成功!" << endl << endl;
namedWindow("原图像", WINDOW_AUTOSIZE);
imshow("原图像", g_srcImage);
namedWindow("双边滤波图像", WINDOW_AUTOSIZE);
g_ndValue = 10;
g_nsigmaColorValue = 10;
g_nsigmaSpaceValue = 10;
char dName[20];
sprintf(dName, "邻域直径 %d", g_ndMaxValue);
char sigmaColorName[20];
sprintf(sigmaColorName, "sigmaColor %d", g_nsigmaColorMaxValue);
char sigmaSpaceName[20];
sprintf(sigmaSpaceName, "sigmaSpace %d", g_nsigmaSpaceMaxValue);
createTrackbar(dName, "双边滤波图像", &g_ndValue, g_ndMaxValue, on_bilateralFilterTrackbar);
on_bilateralFilterTrackbar(g_ndValue, 0);
createTrackbar(sigmaColorName, "双边滤波图像", &g_nsigmaColorValue,
g_nsigmaColorMaxValue, on_bilateralFilterTrackbar);
on_bilateralFilterTrackbar(g_nsigmaColorValue, 0);
createTrackbar(sigmaSpaceName, "双边滤波图像", &g_nsigmaSpaceValue,
g_nsigmaSpaceMaxValue, on_bilateralFilterTrackbar);
on_bilateralFilterTrackbar(g_nsigmaSpaceValue, 0);
waitKey(0);
return 0;
}
void on_bilateralFilterTrackbar(int, void*)
{
bilateralFilter(g_srcImage, g_dstImage, g_ndValue, g_nsigmaColorValue, g_nsigmaSpaceValue);
imshow("双边滤波图像", g_dstImage);
}
需要有高斯滤波和双边滤波的相关知识背景才能更好的理解导向滤波。在导向滤波中,首先利用了局部线性模型。这个模型认为某函数上一点与其近邻部分的点成线性关系,一个复杂的函数就可以用很多局部的线性函数来表示,当需要求该函数上某一点的值时,只需要计算所有包含该点的线性函数的值并取平均值即可。这种模型,在表示非解析函数上,非常有用。
同理,我们可以认为图像是一个二维函数,并且假设该函数的输出与输入在一个二维窗口内满足线性关系,如下:
其中,
是输出像素的值,
是输入图像的值,
和
是像素索引,
和
是当窗口中心位于k时该线性函数的系数。其实,输入图像不一定是待滤波的图像本身,也可以是其他图像即引导图像,这也是为何称为引导滤波的原因。对上式两边取梯度,可以得到:
即当输入图像
有梯度时,输出
也有类似的梯度,现在可以解释为什么引导滤波有边缘保持特性了。下一步是求出线性函数的系数,也就是线性回归,即希望拟合函数的输出值与真实值
之间的差距最小,也就是让下式最小:
这里
只能是待滤波图像,并不像
那样可以是其他图像。同时,
之前的系数用于防止求得的
过大,也是调节滤波器滤波效果的重要参数(相当于L2正则化的权重惩罚)。接下来利用最小二乘法的原理令
和
得到2个二元一次方程,求解得到:
其中
是
在窗口
的平均值,
是
在窗口
的方差,
是窗口
中的像素个数,
是待滤波图像
在窗口
中的均值。在计算每个窗口的线性系数时,我们可以发现一个像素会被多个窗口包含,也就是说,每个像素都由多个线性函数所描述。因此,如之前所说,要具体求某一点的输出值时,只需将所有包含该点的线性函数值平均即可,如下: