浅谈OpenCV均值滤波的实现原理

1.滤波、核和卷积基础知识滤波简单来说就是按照某种规则用周围的数据来重新计算当前数据从而改变信号或图像的特征。在图像处理中图像可以看成一张二维的数字表格每个格子是一个像素。滤波就是拿一个小窗口在图像上滑动对窗口里的像素做某种运算算出一个新值代替窗口中心的像素这个小窗口就是滤波器也称核。这种用核在图像上滑动、对局部像素进行加权求和的操作在数学上称为卷积Convolution。严格来说核不翻转、直接加权求和的操作在数学上称为相关Correlation卷积则需要先把核旋转。OpenCV实现的其实是相关运算。由于本文讨论的盒状核与高斯核都中心对称旋转后不变因此文中沿用了业界习惯的叫法“卷积”。日常图像处理中我们所说的“做一次滤波”通常就是指对图像进行一次卷积运算。图1. 5x5盒状核与5x5高斯核“锚点”由粗体表示图1中展示了图像处理中常用的两种5x5核用于对图像进行平滑或模糊处理。以盒状核为例矩阵中有一个点的值加粗且带有下划线这个点就称作“核的锚点anchor”。锚点用于定义核与源图像的对齐关系在计算输出图像的某个像素时将该像素的位置与锚点对齐再用核覆盖源图像中对应的邻域像素。注意锚点不一定在核的中心图中它位于中心只是为了展示一种常见情形。图1中矩阵旁边的1/25代表在计算时将矩阵内的数字与图像对应像素相乘后求和最后再乘以1/25。实际上这就是5x5的均值滤波核它会对周围25个像素求平均值起到平滑图像、去除噪声的作用但会让图像边缘变得模糊。再来说明一下图中的高斯核与盒状核相同的是它也是5x5的矩阵不同的是核里的数字呈现“中间大、四周小”的对称分布。旁边的分数1/273中273是矩阵里所有25个权重数字相加的总和用来对加权求和的结果进行归一化即除以273求平均从而保证图像整体亮度不变矩阵里的数字如1、4、7、26、41是赋予像素的实际权重而并非实际的像素值。这一点与盒状核原理一样其他核同理只不过盒状核赋予所有像素的权重都是平等的全是1而高斯核则根据距离中心的远近赋予了不同的权重。2.滤波核实际工作流程1.非边界值滤波核处理流程在上一章中我们了解了滤波的基础原理这一章中我们将学习均值滤波的完整流程并对照OpenCV中均值滤波的实现原理展开。图2. 3x3均值滤波核的工作流程图2中展示了3x3均值滤波核的工作流程src为源图像dst为输出图像src中边缘加粗的方框代表3x3滤波核核中心灰色的元素是锚点当滤波核在源图像上滑动到某一位置时会将核内元素的像素值进行加权求和然后进行归一化处理即除以核内元素的总数3x3核即除以9计算的结果四舍五入最终计算的结果会存入dst中与锚点对应的元素内一般不会对原图造成任何影响。但是如果把同一个Mat同时当src和dst传入(blur(img, img, ...))两个名字会指向同一块内存语义上就变成原地写入某些实现下结果不正确所以尽量不要把src和dst传同一个Mat虽然在源码中确实有路径专门支持这么做但仍建议src与dst使用不同的Mat。2.边界值滤波核处理流程然而当滤波核滑动到图像边缘时核的一部分会超出图像边界出现“无像素可取”的情况例如计算图像第一行第一列的像素时核的左上角会悬空。为了解决这个问题我们需要在滤波前对图像边界进行边界值处理。常见的处理方式有以下几种边框类型含义效果BORDER_CONSTANT用指定常量默认 0填充边缘易出现黑边BORDER_WRAP复制对边的像素上下/左右强行缝合极不自然BORDER_REPLICATE复制最边缘的像素无黑边但边缘有拉伸拖尾感BORDER_REFLECT镜像反射包含边界元素最边缘像素被重复一次过渡不够平滑BORDER_REFLECT_101镜像反射不包含边界元素对称轴在边缘外侧过渡平滑且不引入新颜色BORDER_DEFAULTBORDER_REFLECT_101的别名默认值图3. 3x3滤波核在镜像反射BORDER_REFLECT_101处理后的矩阵中进行计算在边界处理中镜像反射既是OpenCV默认的边界处理类型也是实际使用中最常见的类型图3展示了3x3滤波核在进行镜像反射处理后的矩阵中进行计算的过程。边界需要向外扩展的宽度由核的大小决定当锚点位于核中心时扩展宽度等于核的半径图3以3x3核为例半径是1因此在原矩阵外围扩展了1圈若换成5x5核就需要向外扩展2圈7x7核则扩展3圈依此类推。在表格中的几种边界处理类型中BORDER_REFLECT_101对于普通图像的处理一般是最平滑自然的因此OpenCV将其作为默认的边界类型所以在实际编写代码时如果不手动指定borderType边界类型处理OpenCV底层自动调用的就是BORDER_REFLECT_101。图4.通过不同的边界处理方法处理图像的结果Original作为对比的源图像3.浅谈opencv中实现的均值滤波1.分离式滤波在之前的章节中我们学习了以3x3滤波核为例的均值滤波工作流程可以发现如果直接使用二维核9个元素每个像素都需要进行九次加法运算。当核变大如31x31时计算量会呈平方级爆炸式增长。而分离式滤波就是解决这个问题的核心技巧。由于均值滤波核是一个二维的全1矩阵在数学上它可以被分解为两个一维向量的乘积。以3x3均值核为例这意味着原本的二维卷积可以拆分为两步一维卷积。1.行滤波先对图像的每一行做一维均值滤波如图5所示1x3的核在水平方向滑动计算水平方向3个像素的平均值。2.列滤波对行滤波的结果每一列做一维均值滤波如图6所示3x1的核在垂直方向滑动将行滤波的结果再做垂直方向的平均得到最终输出。顺序问题先做行滤波再做列滤波或者先列后行最终结果理论上是一样的OpenCV底层通常先做行再做列因为这样对内存缓存的利用更友好。边界问题分离式滤波同样需要进行边界处理按照BORDER_REFLECT_101等规则对边界进行填充。图5.行滤波工作流程图6.列滤波工作流程为什么分离式滤波这么快因为对比直接二维滤波计算复杂度从O(N²)降到了O(N)。我们具体对比一下计算量。当N3时直接二维滤波每个像素需要N²次运算也就是3x39次加法运算而分离式滤波需要NN2N次计算,也就是336次加法运算差距尚不明显当N11时直接二维滤波需要进行11x11121次加法运算而分离式滤波只需要111122次加法运算当N31时直接二维滤波需要进行31x31961次加法运算而分离式滤波只需要313162次加法运算。也就是说核越大二者差距越大分离式滤波对于性能的提升越明显这就是OpenCV均值滤波底层默认采用分离式滤波的原因。但是OpenCV实际上还会再做一次优化下一节我们再详细说明。2.滑动窗口分离式滤波把二维核拆成了两个一维核但一维滤波内部仍然是朴素求和——每滑动一格都要把窗口内的N个数重新计算一遍。图7.核宽N3的滑块计算过程图8.核宽N5的滑块计算过程观察图7和图8中滑动的过程可以发现相邻两次的窗口大部分是重叠的。图中红色元素为离开的元素绿色元素为新进入的元素蓝色元素则是两次窗口共有、无需重新计算的元素。窗口每移动一格只有一端有一个像素离开、另一端有一个像素进入中间的N-1个元素不动。因此不必重新求和只需在上一次的结果上做一次修正新窗口和 旧窗口和 - 离开的像素 进入的像素。除第一次计算时需要进行N-1次加法剩下的每一步只需要进行两次计算。于是对N3、11、31这些尺寸每输出一个像素所需的运算次数都是4次行滤波中“加进一个新像素、减去一个离开的像素”是两次运算列滤波中“加进一行、减去一行”同样为两次合计为4次不会随核增大而增加运算次数。但是有一个例外在N5时详细原因在下一章中结合源代码来分析。4.opencv中均值滤波的源代码浅析在上面的章节中我们讨论了滤波的实现原理在这一章中我们将结合OpenCV中滤波的源代码来分析实现过程。1.均值滤波的实现blur 是 OpenCV 对外提供的均值滤波接口而它的实现只有一行void blur(InputArray src, OutputArray dst, Size ksize, Point anchor, int borderType) { CV_INSTRUMENT_REGION(); boxFilter( src, dst, -1, ksize, anchor, true, borderType ); }均值滤波本身不是一个独立的算法它只是“归一化的方框滤波”。前面讲的盒状核、分离式滤波、滑动窗口落到代码里全都发生在boxFilter里面。所以接下来我们只需读懂boxFilter就可以理解均值滤波。2.方框滤波的实现2.1.整体实现方框滤波整体结构如下void boxFilter(InputArray _src, OutputArray _dst, int ddepth, Size ksize, Point anchor, bool normalize, int borderType) { CV_INSTRUMENT_REGION(); CV_Assert(!_src.empty()); CV_OCL_RUN(_dst.isUMat() (borderType BORDER_REPLICATE || borderType BORDER_CONSTANT || borderType BORDER_REFLECT || borderType BORDER_REFLECT_101), ocl_boxFilter3x3_8UC1(_src, _dst, ddepth, ksize, anchor, borderType, normalize)) CV_OCL_RUN(_dst.isUMat(), ocl_boxFilter(_src, _dst, ddepth, ksize, anchor, borderType, normalize)) Mat src _src.getMat(); int stype src.type(), sdepth CV_MAT_DEPTH(stype), cn CV_MAT_CN(stype); if( ddepth 0 ) ddepth sdepth; _dst.create( src.size(), CV_MAKETYPE(ddepth, cn) ); Mat dst _dst.getMat(); if( borderType ! BORDER_CONSTANT normalize (borderType BORDER_ISOLATED) ! 0 ) { if( src.rows 1 ) ksize.height 1; if( src.cols 1 ) ksize.width 1; } Point ofs; Size wsz(src.cols, src.rows); if(!(borderTypeBORDER_ISOLATED)) src.locateROI( wsz, ofs ); CALL_HAL(boxFilter, cv_hal_boxFilter, src.ptr(), src.step, dst.ptr(), dst.step, src.cols, src.rows, sdepth, ddepth, cn, ofs.x, ofs.y, wsz.width - src.cols - ofs.x, wsz.height - src.rows - ofs.y, ksize.width, ksize.height, anchor.x, anchor.y, normalize, borderType~BORDER_ISOLATED); //CV_IPP_RUN_FAST(ipp_boxfilter(src, dst, ksize, anchor, normalize, borderType)); borderType (borderType~BORDER_ISOLATED); if(sdepth CV_32F src.type() dst.type() (ksize.height 5 ksize.width 5)) { CV_CPU_DISPATCH(blockSum, (src, dst, ksize, anchor, wsz, ofs, normalize, borderType), CV_CPU_DISPATCH_MODES_ALL); return; } PtrFilterEngine f createBoxFilter( src.type(), dst.type(), ksize, anchor, normalize, borderType ); f-apply( src, dst, wsz, ofs ); }boxFilter内部按优先级依次判断选择一条执行路径①尝试OpenCL(GPU加速)②尝试HAL厂商提供的加速库③浮点数据且核不超过5x5 → blockSum④以上均不适用 → FilterEngine通用实现2.2.具体实现2.2.1.函数签名void boxFilter(InputArray _src, OutputArray _dst, int ddepth, Size ksize, Point anchor, bool normalize, int borderType)各参数在前文均已出现ksize图1中的核尺寸。anchor图1中的锚点传-1-1表示取核中心。normalize图1中的1/25和1/273。borderType第二章的边界处理类型。有一处命名约定需要说明_src与_dst前面的下划线表示这是尚未确定具体类型的通用容器——它既可以接收Mat也可以接受UMatGPU上的矩阵。因此函数内的第一步是将其转换为确定的类型Mat src _src.getMat(); // 转换为 CPU 上的 Mat2.2.2.前置检查CV_INSTRUMENT_REGION(); CV_Assert(!_src.empty());CV_INSTRUMENT_REGION()是性能分析工具的插桩标记与算法本身无关。CV_Assert则是一项输入检查源图不能为空。第二章讨论的是“核部分超出图像”由边界填充解决而整幅图像都没有像素属于另一种情形此时直接断言失败并报错。2.2.3.尝试OpenCLCV_OCL_RUN(_dst.isUMat() (borderType BORDER_REPLICATE || borderType BORDER_CONSTANT || borderType BORDER_REFLECT || borderType BORDER_REFLECT_101), ocl_boxFilter3x3_8UC1(_src, _dst, ddepth, ksize, anchor, borderType, normalize)) CV_OCL_RUN(_dst.isUMat(), ocl_boxFilter(_src, _dst, ddepth, ksize, anchor, borderType, normalize))CV_OCL_RUN条件表达式的含义是若条件成立则执行该表达式若执行成功便立即从boxFilter返回。因此这两行合起来的效果是“能交给GPU就交给GPU交不出去就继续往下走”。两行的分工不同第一行调用的是专用内核它只服务于3x3核、8位单通道等一小组特定情形并且要求边界类型为四种常见类型之一。该内核自身还会检查设备型号与数据对齐方式条件不满足时会返回false于是控制权返回到boxFilter继续执行下一行。第二行调用的是通用内核适用范围更宽但不享受专用内核的优化。此外要注意两行的第一个条件_dst.isUMat()它表示目标位于GPU这是走OpenCL路径的前提。若目标是普通的Mat这两行都会被跳过。2.2.4.退回CPU若GPU路径未走通则退回CPU执行Mat src _src.getMat(); int stype src.type(), sdepth CV_MAT_DEPTH(stype), cn CV_MAT_CN(stype); if( ddepth 0 ) ddepth sdepth; _dst.create( src.size(), CV_MAKETYPE(ddepth, cn) ); Mat dst _dst.getMat(); if( borderType ! BORDER_CONSTANT normalize (borderType BORDER_ISOLATED) ! 0 ) { if( src.rows 1 ) ksize.height 1; if( src.cols 1 ) ksize.width 1; }前五行的作用如下拆包把通用容器_src转换为确定的Mat_dst同样如此。拆分类型CV_MAT_DEPTH取出数据深度如CV_8UCV_MAT_CN取出通道数。后续计算要分别用到这两项信息因此先把合成类型拆开。分配输出_dst.create(...)只在该矩阵的尺寸或类型与要求不符时才真正重新分配内存这是OpenCV的惰性分配约定。ddepth 0表示与源图保持相同深度1.1节中blur传入的-1正是此意。后面的if处理的是图像退化为单行或单列的情形。它只在三个条件同时成立时触发边界类型不是常量填充、需要归一化、且处于隔离模式。此时核在另一维度上找不到足够的真实像素而归一化又要除以核面积即核内元素总数若不作修正除数就会包含并未真正参与运算的像素导致结果整体偏暗。把该维度的核尺寸压成1等价于只除以真实参与运算的像素数。2.2.5.区域位置与隔离标志Point ofs; Size wsz(src.cols, src.rows); if(!(borderTypeBORDER_ISOLATED)) src.locateROI( wsz, ofs );src未必是一幅完整的图像它可能只是从大图中截取的一块区域例如img(Rect(...))。locateROI的作用是查询该区域在整幅图中的位置ofs该区域左上角在整幅图中的坐标。wsz整幅图的尺寸。滤波需要这些信息的原因是如果区域之外存在真实像素就没有必要人为构造边界像素。例如从大图中截取一块区域做滤波时它的边界之外在大图中确实存在真实像素直接使用它们显然比镜像复制的结果更准确。BORDER_ISOLATED标志正是用来区分这两种情形的标志已设置时该区域被视为一幅独立完整的图像边界像素按borderType构造标志未设置时则承认它是大图的一部分优先使用大图中真实存在的像素。这也解释了代码中那个if的作用——隔离模式下不需要查询位置因为区域之外的像素本来就不该被使用。2.2.6.尝试HALCALL_HAL(boxFilter, cv_hal_boxFilter, src.ptr(), src.step, dst.ptr(), dst.step, src.cols, src.rows, sdepth, ddepth, cn, ofs.x, ofs.y, wsz.width - src.cols - ofs.x, wsz.height - src.rows - ofs.y, ksize.width, ksize.height, anchor.x, anchor.y, normalize, borderType~BORDER_ISOLATED);HAL是硬件加速层的缩写。CALL_HAL的含义是若当前平台提供了cv_hal_boxFilter这一实现则调用它调用成功即直接从boxFilter返回。参数虽多但可以分组来看源src.ptr()、src.step。目标dst.ptr()、dst.step。尺寸与类型src.cols、src.rows、sdepth、ddepth、cn。区域位置ofs.x、ofs.y以及右、下两侧的距离。滤波参数ksize、anchor、normalize、borderType。在普通计算机上通常没有该实现因此这一行一般会被跳过。紧接着还有一行borderType (borderType~BORDER_ISOLATED);BORDER_ISOLATED不是一种边界填充方式而是一个标志位其含义已在2.2.5说明。到这里它的作用已经完成——区域的位置信息已由ofs与wsz显式表达后续调用不再需要这个标志所以先把它去掉以免干扰下面的判断。2.2.7.blockSum分支if(sdepth CV_32F src.type() dst.type() (ksize.height 5 ksize.width 5)) { CV_CPU_DISPATCH(blockSum, (src, dst, ksize, anchor, wsz, ofs, normalize, borderType), CV_CPU_DISPATCH_MODES_ALL); return; }触发条件是三项浮点数据、输入与输出类型相同、核不超过5x5。满足则提前返回。单独为它开一条路径是因为2.2.8中的FilterEngine是一套通用框架要管理行缓冲、列缓冲、边界、锚点与缩放系数框架本身的开销对小核而言并不划算。blockSum采用更直接的做法并省去了一次中间结果的内存往返因此在小核情形下更快。核增大之后这一优势不再明显所以限制在5x5以内。2.2.8.通用实现PtrFilterEngine f createBoxFilter( src.type(), dst.type(), ksize, anchor, normalize, borderType ); f-apply( src, dst, wsz, ofs );这是最后一条路径也是任何情况下都可用的实现。createBoxFilter负责装配它先构造一个行滤波器与一个列滤波器再把两者组装成一个滤波引擎。这正是第三章所讲的分离式滤波落地的地方——第三章从数学上说明了二维核可以拆成两个一维核而这里就是把两个一维核对应的滤波器实际创建出来。f-apply(src,dst,wsz,ofs)随后执行滤波wsz与ofs用于告知引擎这块区域在整幅图中的位置。至此boxFilter的工作结束。它完成的事情是判断采用哪条执行路径、确定中间结果使用的数据类型、装配行滤波器与列滤波器。至于两个滤波器具体如何计算boxFilter并未涉及那部分在box_filter.simd.hpp里。3.box_filter.simd.hpp中行滤波和列滤波的实现3.1.均值滤波会用到哪些代码2.2.8中提到createBoxFilter会构造一个行滤波器与一个列滤波器。这两个滤波器就定义在box_filter.simd.hpp中分别是RowSum与ColumnSum两个结构体。均值滤波真正用到的就是这两个。文件中还有另外两个结构体BlockSum服务于2.2.7的blockSum分支做法与上面两者不同SqrRowSum服务于sqrBoxFilter与均值滤波无关。在展开这两个结构体之前先明确两个文件的关系。它们并非并列而是同一份代码的两个层次box_filter.dispatch.cpp决定采用哪条执行路径并装配各组成部分box_filter.simd.hpp负责具体的计算逻辑。box_filter.simd.hpp并不是一个普通的头文件它会被反复包含多次每次针对不同的CPU指令集从而使同一套算法生成多份实现opt_BASELINE 中的 RowSum、ColumnSum…… opt_SSE4_2 中的 RowSum、ColumnSum…… opt_AVX2 中的 RowSum、ColumnSum…… opt_AVX512 中的 RowSum、ColumnSum……而2.2.7中出现的CV_CPU_DISPATCH(blockSum,...)其作用正是在运行时选择一份当前CPU支持的实现来调用。概括而言算法只编写一次机器码生成多份运行时择优选用。因此同一份OpenCV代码在不同计算机上执行的实际上是不同的机器码。3.2.行滤波RowSumRowSum的完整实现如下templatetypename T, typename ST struct RowSum : public BaseRowFilter { RowSum( int _ksize, int _anchor ) : BaseRowFilter() { ksize _ksize; anchor _anchor; } virtual void operator()(const uchar* src, uchar* dst, int width, int cn) CV_OVERRIDE { CV_INSTRUMENT_REGION(); const T* S (const T*)src; ST* D (ST*)dst; int i 0, k, ksz_cn ksize*cn; width (width - 1)*cn; if( ksize 3 ) { for( i 0; i width cn; i ) { D[i] (ST)S[i] (ST)S[icn] (ST)S[icn*2]; } } else if( ksize 5 ) { for( i 0; i width cn; i ) { D[i] (ST)S[i] (ST)S[icn] (ST)S[icn*2] (ST)S[i cn*3] (ST)S[i cn*4]; } } else if( cn 1 ) { ST s 0; for( i 0; i ksz_cn; i ) s (ST)S[i]; D[0] s; for( i 0; i width; i ) { s (ST)S[i ksz_cn] - (ST)S[i]; D[i1] s; } } else if( cn 3 ) { ST s0 0, s1 0, s2 0; for( i 0; i ksz_cn; i 3 ) { s0 (ST)S[i]; s1 (ST)S[i1]; s2 (ST)S[i2]; } D[0] s0; D[1] s1; D[2] s2; for( i 0; i width; i 3 ) { s0 (ST)S[i ksz_cn] - (ST)S[i]; s1 (ST)S[i ksz_cn 1] - (ST)S[i 1]; s2 (ST)S[i ksz_cn 2] - (ST)S[i 2]; D[i3] s0; D[i4] s1; D[i5] s2; } } else if( cn 4 ) { ST s0 0, s1 0, s2 0, s3 0; for( i 0; i ksz_cn; i 4 ) { s0 (ST)S[i]; s1 (ST)S[i1]; s2 (ST)S[i2]; s3 (ST)S[i3]; } D[0] s0; D[1] s1; D[2] s2; D[3] s3; for( i 0; i width; i 4 ) { s0 (ST)S[i ksz_cn] - (ST)S[i]; s1 (ST)S[i ksz_cn 1] - (ST)S[i 1]; s2 (ST)S[i ksz_cn 2] - (ST)S[i 2]; s3 (ST)S[i ksz_cn 3] - (ST)S[i 3]; D[i4] s0; D[i5] s1; D[i6] s2; D[i7] s3; } } else for( k 0; k cn; k, S, D ) { ST s 0; for( i 0; i ksz_cn; i cn ) s (ST)S[i]; D[0] s; for( i 0; i width; i cn ) { s (ST)S[i ksz_cn] - (ST)S[i]; D[icn] s; } } } };逐段说明构造函数只保存核尺寸ksize与锚点anchor。这里有一处细节值得注意anchor被保存下来但在后面的operator()中从头到尾没有出现。原因是滤波器的职责只有计算对准关系由FilterEngine在调用之前安排妥当。operator()接收一行数据输出这一行的行和。开头先把宽度换算成元素个数width (width - 1) * cn同时算出ksz_cn ksize * cn也就是一个窗口跨越的元素个数。接下来的if-else按核尺寸与通道数区分ksize 3把窗口内三个元素直接相加每个输出两次加法。ksize 5把窗口内五个元素直接相加每个输出四次加法。cn 1其余核尺寸这是滑动写法。先逐个累加出第一个窗口的和此后每前进一格只做一次减法与一次加法——减去离开的S[i]加上进入的S[i ksz_cn]。cn 3、cn 4以及最后的通用分支写法与cn 1相同只是为每个通道各维护一个累加器从而避免通道的交错存储带来的重复寻址。3.3.N5的例外第三章末尾留下了一个问题为什么N5时为例外从上面对RowSum的说明就可以看出答案N3与N5走的是展开相加并未滑动。因此各尺寸的运算量如下N行滤波列滤波合计3展开2次滑动2次45展开4次滑动2次611滑动2次滑动2次431滑动2次滑动2次4N5时行滤波多出两次加法合计即为6次。那么为何不始终采用滑动因为滑动存在一项代价每一步都必须等待上一步完成。s (ST)S[i ksz_cn] - (ST)S[i]; // 本句必须等待上一句的 s 计算完毕 D[i1] s;这构成了串行依赖无法并行执行。而展开相加不存在这一问题D[i] (ST)S[i] (ST)S[icn] (ST)S[icn*2];相邻的两个输出之间互不依赖处理器可以同时执行多个。因此核较小时展开相加更优——运算次数本就不多并行执行的收益更大。核较大时滑动更优——滑动恒为两次运算而展开的运算次数随N线性增长。两者的交叉点约在N6~7附近。因此OpenCV只为3和5编写了展开版本这两个恰好是最常用的核尺寸其余尺寸统一采用滑动。需要指出的是源码中并不存在“N≥7则采用滑动”这类判断只有ksize 3与ksize 5两条硬编码分支。因此更准确的表述是OpenCV只特化了两个最常用的尺寸其余尺寸统一使用滑动。3.4.列滤波ColumnSumColumnSum的通用版本如下此外还有若干针对具体数据类型的特化版本见本章节末尾说明templatetypename ST, typename T struct ColumnSum : public BaseColumnFilter { ColumnSum( int _ksize, int _anchor, double _scale ) : BaseColumnFilter() { ksize _ksize; anchor _anchor; scale _scale; sumCount 0; } virtual void reset() CV_OVERRIDE { sumCount 0; } virtual void operator()(const uchar** src, uchar* dst, int dststep, int count, int width) CV_OVERRIDE { CV_INSTRUMENT_REGION(); int i; ST* SUM; bool haveScale scale ! 1; double _scale scale; if( width ! (int)sum.size() ) { sum.resize(width); sumCount 0; } SUM sum[0]; if( sumCount 0 ) { memset((void*)SUM, 0, width*sizeof(ST)); for( ; sumCount ksize - 1; sumCount, src ) { const ST* Sp (const ST*)src[0]; for( i 0; i width; i ) SUM[i] Sp[i]; } } else { CV_Assert( sumCount ksize-1 ); src ksize-1; } for( ; count--; src ) { const ST* Sp (const ST*)src[0]; const ST* Sm (const ST*)src[1-ksize]; T* D (T*)dst; if( haveScale ) { for( i 0; i width - 2; i 2 ) { ST s0 SUM[i] Sp[i], s1 SUM[i1] Sp[i1]; D[i] saturate_castT(s0*_scale); D[i1] saturate_castT(s1*_scale); s0 - Sm[i]; s1 - Sm[i1]; SUM[i] s0; SUM[i1] s1; } for( ; i width; i ) { ST s0 SUM[i] Sp[i]; D[i] saturate_castT(s0*_scale); SUM[i] s0 - Sm[i]; } } else { for( i 0; i width - 2; i 2 ) { ST s0 SUM[i] Sp[i], s1 SUM[i1] Sp[i1]; D[i] saturate_castT(s0); D[i1] saturate_castT(s1); s0 - Sm[i]; s1 - Sm[i1]; SUM[i] s0; SUM[i1] s1; } for( ; i width; i ) { ST s0 SUM[i] Sp[i]; D[i] saturate_castT(s0); SUM[i] s0 - Sm[i]; } } dst dststep; } } double scale; int sumCount; std::vectorST sum; };逐段说明构造函数除了核尺寸与锚点还接收一个缩放系数scale。第三章讲过归一化系数只放在列滤波这一步就是通过这个参数传进来的。同样anchor在这里也没有被使用。成员变量中有一个vectorSTsum它是整段代码的关键累加器不是一个标量而是一个数组长度等于图像的宽度。这与RowSum只用一个标量s形成对照。operator()的参数与前两个不同src不是一整行而是若干个指针const uchar**。dststep输出每次前进的字节数。count本次要输出多少行。width每行多少列。函数体分三段。第一段是初始化if( sumCount 0 ) { memset((void*)SUM, 0, width*sizeof(ST)); for( ; sumCount ksize - 1; sumCount, src ) { const ST* Sp (const ST*)src[0]; for( i 0; i width; i ) SUM[i] Sp[i]; } }它先把ksize-1行累加进数组。之所以是ksize-1行是因为接下来每读入一行就会输出一行那么第一行输出所对应的窗口正好由这ksize-1行加上新读入的一行凑齐。第二段是主体循环外层推进一行for( ; count--; src ) { const ST* Sp (const ST*)src[0]; // 新进入的一整行 const ST* Sm (const ST*)src[1-ksize]; // 离开的一整行 ... }Sp指向新进入窗口的那一行Sm指向离开窗口的那一行。第三段是内层循环遍历这一行的所有列ST s0 SUM[i] Sp[i]; // 加上新进入的一行 D[i] saturate_castT(s0*_scale); // 输出需要时乘以缩放系数 SUM[i] s0 - Sm[i]; // 减去离开的一行因此列滤波的推进单位是一整行而不是单个像素内层虽然逐列进行但它们在同一个外层轮次中全部完成。为什么列滤波的累加器是数组行滤波的窗口是横向的N列滑动方向沿着一行向右每前进一格只影响一个位置因此一个标量s就足以记住“当前窗口的和”。列滤波的窗口是纵向的N行、同一列。要对某一列做滑动就需要记住“这一列当前的和”而图像有W列每一列都要各自记住于是累加器必然是W个数也就是一个数组。上面给出的是通用版本。文件中还有一批针对具体数据类型的特化版本例如ColumnSumint,uchar、ColumnSumushort,uchar、ColumnSumint,short等具体使用哪一个由2.2.8中createBoxFilter选定的中间类型决定。这些特化版本沿用的是同一套算法只是在此基础上增加了一些优化例如用整数乘法加移位代替浮点除法、以及对循环做SIMD向量化。这些内容超出本章范围此处不再展开。5.简单实现均值滤波下面来根据原理来简单实现一下opencv中的均值滤波。#include opencv2/opencv.hpp #include iostream #include windows.h using namespace cv; using namespace std; // 边界镜像把越界下标折回图内 // 实现 BORDER_REFLECT_101 —— 反射时不重复最外圈像素所以 -1→1 而非 0size→size-2 而非 size-1 static inline int mirrorIndex(int index, int size) { if (index 0) return -index; if (index size) return 2 * size - index - 2; return index; } // 分离式 滑块的均值滤波 // 朴素每像素 ksize² 次加法 → 分离式两趟一维2×ksize 次 → 滑块恒定 4 次与核大小无关 Mat myblur(const Mat src, int ksize) { const int r ksize / 2, h src.rows, w src.cols, count ksize * ksize; // 第一步行滤波横向滑动 // 中间图存和不存平均值用 32 位整数是因为和会溢出 uchar // ksize31 时单通道和最大 255×317905远超 255 Mat rowSum(src.size(), CV_32SC3); for (int row 0; row h; row) { // 累加器必须放在行循环内、列循环外值才能跨列传递滑动才成立 int b 0, g 0, rr 0; for (int col 0; col w; col) { if (col 0) { // 第一列没有上一格没有可减的像素只能把整个窗口加一遍 for (int k -r; k r; k) { Vec3b p src.atVec3b(row, mirrorIndex(col k, w)); b p[0]; g p[1]; rr p[2]; } } else { // 窗口右移一格左端出去一个、右端进来一个 // 输出 col-1 时窗口是 [col-1-r, col-1r]输出 col 时是 [col-r, colr] // 所以离开的是 col-r-1进入的是 colr Vec3b leave src.atVec3b(row, mirrorIndex(col - r - 1, w)); Vec3b enter src.atVec3b(row, mirrorIndex(col r, w)); b enter[0] - leave[0]; g enter[1] - leave[1]; rr enter[2] - leave[2]; } // 只存和除法留到最后做一次整条链路上只取整一次 rowSum.atVec3i(row, col) Vec3i(b, g, rr); } } // 第二步列滤波纵向滑动 Mat dst Mat::zeros(src.size(), CV_8UC3); // 循环顺序反过来了col 在外、row 在内 // 因为要沿纵向滑动必须把一列走完累加器才能跨行传递 for (int col 0; col w; col) { int b 0, g 0, rr 0; for (int row 0; row h; row) { if (row 0) { for (int k -r; k r; k) { Vec3i p rowSum.atVec3i(mirrorIndex(row k, h), col); b p[0]; g p[1]; rr p[2]; } } else { Vec3i leave rowSum.atVec3i(mirrorIndex(row - r - 1, h), col); Vec3i enter rowSum.atVec3i(mirrorIndex(row r, h), col); b enter[0] - leave[0]; g enter[1] - leave[1]; rr enter[2] - leave[2]; } // 此处 b/g/rr 是 ksize² 个像素的总和除以 count 才是平均值 // 整数除法是向下取整cv::blur 是四舍五入所以两者最多差 1 dst.atVec3b(row, col) Vec3b(b / count, g / count, rr / count); } } return dst; } // 对比与输出 int main() { SetConsoleOutputCP(CP_UTF8); Mat img imread(C:/Users/Administrator/Desktop/sucai/milkfrog2.jpg); // imread 失败不抛异常只返回空 Mat必须显式检查 if (img.empty()) { cout img not found endl; return -1; } if (img.channels() 1) cvtColor(img, img, COLOR_GRAY2BGR); // 灰度图按 Vec3b 读会崩溃 int ksize 11; int64 t getTickCount(); Mat mine myblur(img, ksize); double ms1 (getTickCount() - t) * 1000.0 / getTickFrequency(); // 周期差 ÷ 每秒周期数 秒×1000 得毫秒写 1000.0 防止整数除法截断小数 Mat ref; t getTickCount(); blur(img, ref, Size(ksize, ksize)); // 第 4 个参数省略 BORDER_DEFAULT BORDER_REFLECT_101与自己的实现一致才有可比性 double ms2 (getTickCount() - t) * 1000.0 / getTickFrequency(); cout myblur : ms1 ms endl; cout cv::blur : ms2 ms endl; cout ratio : ms1 / ms2 x endl; // 差异图差值最大只有 1乘 20 才看得出来 // 全黑完全一致均匀暗灰全是取整误差局部发亮那里真的算错了 Mat diff; absdiff(mine, ref, diff); diff.convertTo(diff, -1, 20); // 必须 clone()因为 putText 会直接修改传入的图像 Mat a img.clone(), b mine.clone(), c ref.clone(), d diff.clone(); putText(a, src, Point(10, 35), FONT_HERSHEY_SIMPLEX, 1.0, Scalar(0, 0, 255), 2); putText(b, myblur, Point(10, 35), FONT_HERSHEY_SIMPLEX, 1.0, Scalar(0, 0, 255), 2); putText(c, cv::blur, Point(10, 35), FONT_HERSHEY_SIMPLEX, 1.0, Scalar(0, 0, 255), 2); putText(d, diff x20, Point(10, 35), FONT_HERSHEY_SIMPLEX, 1.0, Scalar(0, 0, 255), 2); Mat t1, t2, canvas; hconcat(a, b, t1); hconcat(c, d, t2); vconcat(t1, t2, canvas); // 按宽、按高各算一个缩放比取小的那个 —— 只按宽度算的话 2×2 拼图会竖向超出屏幕 double sw 0.6 * GetSystemMetrics(SM_CXSCREEN); double sh 0.6 * GetSystemMetrics(SM_CYSCREEN); double s (sw / canvas.cols sh / canvas.rows) ? sw / canvas.cols : sh / canvas.rows; if (s 1.0) resize(canvas, canvas, Size(), s, s, INTER_AREA); namedWindow(comparison, WINDOW_NORMAL); imshow(comparison, canvas); imwrite(C:/Users/Administrator/Desktop/sucai/myblur.png, canvas); waitKey(0); return 0; }ksize11时结果如下图9.对比结果图速度对比如下图10.Debug下速度对比图11.Release下速度对比可以看出还是比官方慢且两种配置下的差距并不相同Debug下官方快了约15倍Release下官方快了约5倍Debug下的代码是未经优化过的而OpenCV中的代码是已经优化过的两者不在同一起跑线上所以真正需要探讨的是Release下的对比cv::blur使用SIMD指令一条指令同时处理8到16个像素列滤波按整行推进全程顺序访问内存中间图在核面积不超过256时只占16位边界也只对首尾若干像素单独处理。而两者的运算次数其实完全相同——都是每像素4次与核尺寸无关。换言之加法次数已经追平剩下的差距在实现细节。6.小结本文围绕均值滤波这一个算法走完了从原理到实现的全过程。均值滤波看起来是最简单的滤波——一个全1的核一次求和一次除法。但就是这样一个小算法背后涉及可分离性、滑动窗口、边界策略、缓存布局与指令集分发。图像处理中的大多数算法也是这样一层一层叠加起来的。本文为学习OpenCV时的心得如有遗漏或错误欢迎在评论区补充或指正。参考书籍《数字图像处理第三版》《学习OpenCV3中文版》《CPrimer Plus第6版中文版》