--------滤波简介--------
滤波的是图像处理之中必备的手段,也是必经之路。模式识别,深度学习都会用到滤波的相关内容,没有好的图片你怎么识别?
个人理解的“滤波”二字:
滤波没有什么神奇的地方,就是通过旁边像素对比来做判断(其他的下面再说)举个例子:
假如上面的图形1是瑕疵,通过滤波给过滤掉
上面是一个两张简单的图片,两个问题: 1.把左边的'形状1'清除(过滤掉)。
2.把右边的'形状1'清除(过滤掉)。
对于第一个问题:那很简单了,前面的博文已经知道,图像是一个个像素点构成的,那对比蓝色和白色的像素点就可以知道,直接把蓝色排除了。
对于第二个问题:那就复杂了,因为图形通过旁边的像素来做判断根本无法识别,这就的用到一些算法知识了。
感觉说的不着边际,简单的意思就是滤波没什么神秘的,接下来举几个滤波实际例子就能说明:
--------均值滤波--------
从另一篇博文可以知道,对于“核”这个概念的知识,还有怎么操作定义核,这里就不啰嗦了。不知道的直接去看另一篇博文:掩膜操作实例分析
均值就是平均的意思,对邻近的像素做个平均:x[i][j]=1/9(y[i-1][j-1]+y[i-1][j]+y[i-1][j+1]+y[i][j-1]+y[i][j]+y[i][j+1]+y[i+1][j-1]+y[i+1][j]+y[i+1][j+1])
手写实例:
其实这些例子都不难,看懂理论基础和opencv源代码,自己再去写很简单,有时间这部分会补上,现在为了想快点把opencv学一遍。。。。所以偷懒了
API:
--------中值滤波--------
中值滤波就是旁边的值进行排序取中间值x[i][j]=Middle(y[i-1][j-1],y[i-1][j],y[i-1][j+1],y[i][j-1],y[i][j],y[i][j+1],y[i+1][j-1],y[i+1][j],y[i+1][j+1]);
手写实例:
API:
--------高斯滤波--------
这个得重点说明一下,刚开始学的时候真的感觉好难,数学公式也看不懂,源代码看着就烦。。。
公式大家自己看书本或者百度百科吧,这里只说一下我自己的理解
二维图像
三维图像
不说专业名词,看着这幅图就是对称的,并且两遍递减,想象一下这个函数去表示“核”是什么样子?
这是随便画的一个图像,从内到外逐渐减小,如果运用到图像“卷积”之后,那就是距离处理点越近那就权值越大了。
卷积:在数字图像上就是,简单理解就是图像和内核的乘积!
至于这个权值是怎么计算的,那就是直接把坐标带入高斯公式就行了,想深入了解的就去看一些数学的书籍。
原公式
简化公式
具体公式不介绍了,看百度、或者其他资料。。。
opencv源码分析:
1 /****************************************************************************************\ 2 Gaussian Blur 3 \****************************************************************************************/ 4 5 cv::Mat cv::getGaussianKernel( int n, double sigma, int ktype ) 6 { 7 //----自定义了几个简单的内核数据-------// 8 const int SMALL_GAUSSIAN_SIZE = 7; 9 static const float small_gaussian_tab[][SMALL_GAUSSIAN_SIZE] = 10 { 11 {1.0f}, 12 {0.25f, 0.5f, 0.25f}, 13 {0.0625f, 0.25f, 0.375f, 0.25f, 0.0625f}, 14 {0.03125f, 0.109375f, 0.21875f, 0.28125f, 0.21875f, 0.109375f, 0.03125f} 15 }; 16 //-------当 n是奇数 && n<7 && sigma<=0 --->>>直接取small_gaussian_tab里的数据-----// 17 //-------当不满足上面的情况时候,就自己按照公式计算,其实上面的值也是按照公式计算的----// 18 const float* fixed_kernel = n % 2 == 1 && n <= SMALL_GAUSSIAN_SIZE && sigma <= 0 ? 19 small_gaussian_tab[n>>1] : 0; 20 21 CV_Assert( ktype == CV_32F || ktype == CV_64F ); 22 Mat kernel(n, 1, ktype);//计算的核是n行1列,后面有程序会合并成二维数组 23 //----保证程序可以执行float、double两种类型都可以----// 24 float* cf = (float*)kernel.data; 25 double* cd = (double*)kernel.data; 26 //-------当sigma取值,这个是高斯函数规定的------// 27 double sigmaX = sigma > 0 ? sigma : ((n-1)*0.5 - 1)*0.3 + 0.8; 28 double scale2X = -0.5/(sigmaX*sigmaX);//高斯公式里面的一个系数 29 double sum = 0; 30 31 int i; 32 for( i = 0; i < n; i++ ) 33 { 34 double x = i - (n-1)*0.5; 35 //------就是上面说的看自己定义还是直接读取------// 36 double t = fixed_kernel ? (double)fixed_kernel[i] : std::exp(scale2X*x*x); 37 if( ktype == CV_32F ) 38 { 39 cf[i] = (float)t; 40 sum += cf[i]; 41 } 42 else 43 { 44 cd[i] = t; 45 sum += cd[i]; 46 } 47 } 48 49 sum = 1./sum;//系数归一化 50 //---归一化系数之后再把系数放在内核里面---// 51 for( i = 0; i < n; i++ ) 52 { 53 if( ktype == CV_32F ) 54 cf[i] = (float)(cf[i]*sum); 55 else 56 cd[i] *= sum; 57 } 58 59 return kernel; 60 } 61 62 //-----函数是处理x核、y核、size核,其中之间的关系和数据变换--------// 63 cv::Ptr<cv::FilterEngine> cv::createGaussianFilter( int type, Size ksize, 64 double sigma1, double sigma2, 65 int borderType ) 66 { 67 int depth = CV_MAT_DEPTH(type); 68 if( sigma2 <= 0 ) 69 sigma2 = sigma1; 70 71 // automatic detection of kernel size from sigma 72 if( ksize.width <= 0 && sigma1 > 0 ) 73 ksize.width = cvRound(sigma1*(depth == CV_8U ? 3 : 4)*2 + 1)|1; 74 if( ksize.height <= 0 && sigma2 > 0 ) 75 ksize.height = cvRound(sigma2*(depth == CV_8U ? 3 : 4)*2 + 1)|1; 76 77 CV_Assert( ksize.width > 0 && ksize.width % 2 == 1 && 78 ksize.height > 0 && ksize.height % 2 == 1 ); 79 80 sigma1 = std::max( sigma1, 0. ); 81 sigma2 = std::max( sigma2, 0. ); 82 83 Mat kx = getGaussianKernel( ksize.width, sigma1, std::max(depth, CV_32F) ); 84 Mat ky; 85 if( ksize.height == ksize.width && std::abs(sigma1 - sigma2) < DBL_EPSILON ) 86 ky = kx; 87 else 88 ky = getGaussianKernel( ksize.height, sigma2, std::max(depth, CV_32F) ); 89 90 return createSeparableLinearFilter( type, type, kx, ky, Point(-1,-1), 0, borderType ); 91 } 92 93 94 void cv::GaussianBlur( InputArray _src, OutputArray _dst, Size ksize, 95 double sigma1, double sigma2, 96 int borderType ) 97 { 98 Mat src = _src.getMat(); 99 _dst.create( src.size(), src.type() ); 100 Mat dst = _dst.getMat(); 101 //-----处理边界问题----// 102 if( borderType != BORDER_CONSTANT ) 103 { 104 if( src.rows == 1 ) 105 ksize.height = 1; 106 if( src.cols == 1 ) 107 ksize.width = 1; 108 } 109 110 if( ksize.width == 1 && ksize.height == 1 ) 111 { 112 src.copyTo(dst); 113 return; 114 } 115 116 #ifdef HAVE_TEGRA_OPTIMIZATION 117 if(sigma1 == 0 && sigma2 == 0 && tegra::gaussian(src, dst, ksize, borderType)) 118 return; 119 #endif 120 121 Ptr<FilterEngine> f = createGaussianFilter( src.type(), ksize, sigma1, sigma2, borderType ); 122 f->apply( src, dst ); 123 }
总结的注意点:
A.高斯核的计算是两种方法,第一种是经过计算放在一个数组里面,第二种是直接用高斯公式计算。
1 static const float small_gaussian_tab[][SMALL_GAUSSIAN_SIZE] = 2 10 { 3 11 {1.0f}, 4 12 {0.25f, 0.5f, 0.25f}, 5 13 {0.0625f, 0.25f, 0.375f, 0.25f, 0.0625f}, 6 14 {0.03125f, 0.109375f, 0.21875f, 0.28125f, 0.21875f, 0.109375f, 0.03125f} 7 15 };
B.高斯的内核是合成的,而不是一次型的计算,具体合成过程我没有细看,我们自己手写API的时候基本 都不会去合成,直接计算便可。
1 Mat kernel(n, 1, ktype);//计算的核是n行1列,后面有程序会合并成二维数组
1 Mat kx = getGaussianKernel( ksize.width, sigma1, std::max(depth, CV_32F) ); 2 84 Mat ky; 3 85 if( ksize.height == ksize.width && std::abs(sigma1 - sigma2) < DBL_EPSILON ) 4 86 ky = kx; 5 87 else 6 88 ky = getGaussianKernel( ksize.height, sigma2, std::max(depth, CV_32F) ); 7 89 8 90 return createSeparableLinearFilter( type, type, kx, ky, Point(-1,-1), 0, borderType );
C.滤波器的引擎里面有很多函数,通过一些参数的匹配来分别解决哪个滤波。就是一个封装多个滤波的函 数,看源代码就会发现。源代码在filter.cpp文件,想看的可以自己去看
1 Ptr<FilterEngine> f = createGaussianFilter( src.type(), ksize, sigma1, sigma2, borderType ); 2 f->apply( src, dst );
--------双边滤波--------
公式百度吧。。。。。
记得“高斯滤波”是考虑了空间问题,使图像的滤波提升了一个档次(相对于均值滤波),但是高斯滤波对比中值滤波不一定哪个更好,看下面例子:
这图片有两个突出的点,一个是0,一个是255,黑白点,其它的像素点都是110左右,这一看就知道是噪声,下面用滤波函数操作这个图像。
均值:0->80+
中值:0->110
高斯:0->90
这个例子举的不明显,但是意图就是高斯函数只能过滤空间的部分,但是距离的部分不能过滤,就是像素差别特别大的时候无法过滤,看网上有个a-截断滤波
高斯+截断=双边 这是简单表示的。
双边滤波就是把那些怪异的点去除之后在进行高斯滤波:
opencv源代码分析:
1 /****************************************************************************************\ 2 Bilateral Filtering 3 \****************************************************************************************/ 4 5 namespace cv 6 { 7 8 static void 9 bilateralFilter_8u( const Mat& src, Mat& dst, int d, 10 double sigma_color, double sigma_space, 11 int borderType ) 12 { 13 int cn = src.channels();//原图通道数 14 int i, j, k, maxk, radius; 15 Size size = src.size();//原图大小 16 17 CV_Assert( (src.type() == CV_8UC1 || src.type() == CV_8UC3) && 18 src.type() == dst.type() && src.size() == dst.size() && 19 src.data != dst.data );//只能处理单通道和三通道的uchar 20 //-------限定sigma的取值--------// 21 if( sigma_color <= 0 ) 22 sigma_color = 1; 23 if( sigma_space <= 0 ) 24 sigma_space = 1; 25 //------高斯因子,计算方便-------// 26 double gauss_color_coeff = -0.5/(sigma_color*sigma_color); 27 double gauss_space_coeff = -0.5/(sigma_space*sigma_space); 28 //-----d 表示滤波时像素邻域直径,d为负时由 sigaColor计算得到;d>5时不能实时处理--// 29 if( d <= 0 ) 30 radius = cvRound(sigma_space*1.5);//float->int 31 else 32 radius = d/2; 33 //---保证直径是基整数---// 34 radius = MAX(radius, 1); 35 d = radius*2 + 1; 36 37 Mat temp; 38 copyMakeBorder( src, temp, radius, radius, radius, radius, borderType ); 39 //---直接用公式计算权值,只有用到直接调用----// 40 vector<float> _color_weight(cn*256);//距离权值范围为channels*256 41 vector<float> _space_weight(d*d);//空间权值存储 42 vector<int> _space_ofs(d*d);//位置便宜存储->用于计算点的位置 43 float* color_weight = &_color_weight[0]; 44 float* space_weight = &_space_weight[0]; 45 int* space_ofs = &_space_ofs[0]; 46 47 // initialize color-related bilateral filter coefficients 48 for( i = 0; i < 256*cn; i++ )//把空间可能的权值权值全部存储起来,用的时候直接找 49 color_weight[i] = (float)std::exp(i*i*gauss_color_coeff);//省的到时候计算 50 51 // ------计算核的值和偏移值----// 52 for( i = -radius, maxk = 0; i <= radius; i++ ) 53 for( j = -radius; j <= radius; j++ ) 54 { 55 double r = std::sqrt((double)i*i + (double)j*j);//高斯公式系数之一 56 if( r > radius ) 57 continue; 58 space_weight[maxk] = (float)std::exp(r*r*gauss_space_coeff);//空间权值 59 space_ofs[maxk++] = (int)(i*temp.step + j*cn);//偏移值计算 60 } 61 //------------------进行滤波操作----------// 62 for( i = 0; i < size.height; i++ ) 63 { 64 //-----sptr核中心的像素地址,同时是输入像素变量----// 65 //-----dptr表示第i行的第一个像素地址,同时是输出像素变量-----// 66 const uchar* sptr = temp.data + (i+radius)*temp.step + radius*cn; 67 uchar* dptr = dst.data + i*dst.step; 68 69 if( cn == 1 ) 70 { 71 for( j = 0; j < size.width; j++ ) 72 { 73 float sum = 0, wsum = 0; 74 int val0 = sptr[j]; 75 for( k = 0; k < maxk; k++ )//maxk=d就是直径的大小 76 { 77 int val = sptr[j + space_ofs[k]];//核对应的像素值 78 float w = space_weight[k]*color_weight[std::abs(val - val0)];//总权值 79 sum += val*w;//权值*像素值总和 80 wsum += w;//权值和,为了下面归一化 81 } 82 // overflow is not possible here => there is no need to use CV_CAST_8U 83 dptr[j] = (uchar)cvRound(sum/wsum);//权值归一化 84 } 85 } 86 else 87 { 88 assert( cn == 3 ); 89 for( j = 0; j < size.width*3; j += 3 )//三个通道(一个像素)同时操作 90 { 91 float sum_b = 0, sum_g = 0, sum_r = 0, wsum = 0; 92 int b0 = sptr[j], g0 = sptr[j+1], r0 = sptr[j+2];//三个通道代表一个像素 93 for( k = 0; k < maxk; k++ )//直径的大小 94 { 95 const uchar* sptr_k = sptr + j + space_ofs[k];//核对应的像素地址 96 int b = sptr_k[0], g = sptr_k[1], r = sptr_k[2]; 97 //-----三个通道的差值和作为计算、所以三个通道的权重是相同的------// 98 float w = space_weight[k]*color_weight[std::abs(b - b0) + 99 std::abs(g - g0) + std::abs(r - r0)]; 100 sum_b += b*w; sum_g += g*w; sum_r += r*w; 101 wsum += w; 102 } 103 //----下面就是BGR三个通道的调整了-------// 104 wsum = 1.f/wsum; 105 b0 = cvRound(sum_b*wsum); 106 g0 = cvRound(sum_g*wsum); 107 r0 = cvRound(sum_r*wsum); 108 dptr[j] = (uchar)b0; dptr[j+1] = (uchar)g0; dptr[j+2] = (uchar)r0; 109 } 110 } 111 } 112 }
总结的注意点:
A.双边滤波如果是三通道的BGA的话,那是三个通道同时进行的,三个通道的差值求和,所以三个通道是一个权重。
为对于其他其他的滤波方式:
B.这个滤波函数使用了指针访问像素,但是这个指针是:*(M.data+i*M.step[0]+j*M.channels()),记得上次我们使用的指针是 M.ptr<Vec3b>(i)[j]访问的。下面就详细说明一下这个函数:
(a).其中M.data是M这个Mat类的首地址,也就是一幅图像的首地址。
(b).M.step[],这个函数比较绕人,三维空间咱们不说,就说我们使用的二维数据。
step[0] :表示图片一行数据字节。
step[1] :表示图片一个像素的字节。
扩展几个知识点:
elemSize() :一个像素的字节 == step[1]
elemSize1():一个通道的字节
比如:我们图片是16X16,CV_8UC3-->>>elemSize1()=1,elemSize()=step[1]=1*3,step[0]=1*3*16
我们图片是25*25,CV_16UC3-->>>elemSize1()=2,elemSize()=step[1]=2*3,step[0]=2*3*25
8位的uchar是一个字节,16位的uchar是2个字节。。。
16X16,CV_8UC3的图片 *(M.data+i*M.step[0]+j*M.channels()):*(0+1*3*16*I+3*J),这里省略step[1]=1
我们之前写的Mat.ptr和Mat.at内部已经帮我们做了这些事,现在直接写的话只有通过字节的地址操作,
就像char a[100];char*p=a;p++就可以操作a[100],int a[100];int* p=a;p+=4;
参考:
http://blog.csdn.net/qianqing13579/article/details/45318279#comments
http://blog.csdn.net/xiaowei_cqu/article/details/7785365
http://blog.csdn.net/zhaocj/article/details/39520187
还有其它参考,没来及记录,入您发现侵犯您的权利,请告之,立马清除!