OpenCV源码分析（四）：中值滤波
2018-12-10 本文已影响0人
奔向火星005
注：为方便分析和解释，文章中的代码是OpenCV源码的简化版，会略去源码中一些如断言，OpenCL，硬件加速等代码。因此详情最好参看源码。
中值滤波对于去除与原图像差异较大的噪声有很好的效果，因为它是通过查找邻域中的所有像素的中值作为目标像素，而噪声被选中的概率很小。
先看下效果：
左边图像加了椒盐噪声，右图为中值滤波处理后的图像，可见效果不错。
OpenCV的中值滤波的接口为为：
void cv::medianBlur( InputArray _src0, OutputArray _dst, int ksize );
可以看出接口非常简单，_src0为输入，_dst为输出，ksize为滤波窗口大小。下面以ksize为5为例分析源码。
void cv::medianBlur( InputArray _src0, OutputArray _dst, int ksize )
{
Mat src0 = _src0.getMat();
_dst.create( src0.size(), src0.type() );
Mat dst = _dst.getMat();
bool useSortNet = ksize == 3 || (ksize == 5
#if !(CV_SIMD128)
&& ( src0.depth() > CV_8U || src0.channels() == 2 || src0.channels() > 4 )
#endif
);
Mat src;
if( useSortNet )
{
if( dst.data != src0.data )
src = src0;
else
src0.copyTo(src);
if( src.depth() == CV_8U ) //图像元素为字节
medianBlur_SortNet<MinMax8u, MinMaxVec8u>( src, dst, ksize );
else if( src.depth() == CV_16U )
medianBlur_SortNet<MinMax16u, MinMaxVec16u>( src, dst, ksize );
else if( src.depth() == CV_16S )
medianBlur_SortNet<MinMax16s, MinMaxVec16s>( src, dst, ksize );
else if( src.depth() == CV_32F )
medianBlur_SortNet<MinMax32f, MinMaxVec32f>( src, dst, ksize );
else
CV_Error(CV_StsUnsupportedFormat, "");
return;
}
else
{
//省略...
}
}
当输入ksize选5，图像每个元素为uchar类型时，会调用medianBlur_SortNet<MinMax8u, MinMaxVec8u>( src, dst, ksize )函数，下面我们看下它的实现：
template<class Op, class VecOp>
static void
medianBlur_SortNet( const Mat& _src, Mat& _dst, int m )
{
typedef typename Op::value_type T;
typedef typename Op::arg_type WT;
typedef typename VecOp::arg_type VT;
const T* src = _src.ptr<T>(); //const unsigned char*
T* dst = _dst.ptr<T>();
int sstep = (int)(_src.step/sizeof(T)); //2250
int dstep = (int)(_dst.step/sizeof(T));
Size size = _dst.size(); //750,1024
int i, j, k, cn = _src.channels(); //cn == 3
Op op;
VecOp vop;
volatile bool useSIMD = hasSIMD128();
if( m == 3 )
{
//省略...
}
else if( m == 5 )
{
if( size.width == 1 || size.height == 1 ) //跳过
{
int len = size.width + size.height - 1;
int sdelta = size.height == 1 ? cn : sstep;
int sdelta0 = size.height == 1 ? 0 : sstep - cn;
int ddelta = size.height == 1 ? cn : dstep;
for( i = 0; i < len; i++, src += sdelta0, dst += ddelta )
for( j = 0; j < cn; j++, src++ )
{
int i1 = i > 0 ? -sdelta : 0;
int i0 = i > 1 ? -sdelta*2 : i1;
int i3 = i < len-1 ? sdelta : 0;
int i4 = i < len-2 ? sdelta*2 : i3;
WT p0 = src[i0], p1 = src[i1], p2 = src[0], p3 = src[i3], p4 = src[i4];
op(p0, p1); op(p3, p4); op(p2, p3); op(p3, p4); op(p0, p2);
op(p2, p4); op(p1, p3); op(p1, p2);
dst[j] = (T)p2;
}
return;
}
size.width *= cn;
for( i = 0; i < size.height; i++, dst += dstep )
{
const T* row[5]; //row[5]的每个元素指向一行
row[0] = src + std::max(i - 2, 0)*sstep; //当i小于2时，说明是上边界，直接用src+0位置的像素
row[1] = src + std::max(i - 1, 0)*sstep; //当i小于2时，~
row[2] = src + i*sstep;
row[3] = src + std::min(i + 1, size.height-1)*sstep;
row[4] = src + std::min(i + 2, size.height-1)*sstep;
int limit = useSIMD ? cn*2 : size.width; //如果使用useSIMD，则需要先单组对比cn*2个数据，因为开头cn*2 个数据实际上是补边数据（它实际上是不存在的，是利用指针指向边界内的数据的），而vop操作（useSIMD）需要一次加载16个连续内存数据，不能在边界进行。
for(j = 0;; )
{
for( ; j < limit; j++ )
{
WT p[25]; //5x5的窗口，WT是
//没有j2？因为j就是j2
int j1 = j >= cn ? j - cn : j;
int j0 = j >= cn*2 ? j - cn*2 : j1;
int j3 = j < size.width - cn ? j + cn : j;
int j4 = j < size.width - cn*2 ? j + cn*2 : j3;
for( k = 0; k < 5; k++ ) //把5x5窗口内的数据赋给p[0]~p[24]
{
const T* rowk = row[k];
p[k*5] = rowk[j0]; p[k*5+1] = rowk[j1];
p[k*5+2] = rowk[j]; p[k*5+3] = rowk[j3];
p[k*5+4] = rowk[j4];
}
//op是优化过的比较函数，将较小的赋给第一个参数，较大的赋给第二个参数
op(p[1], p[2]); op(p[0], p[1]); op(p[1], p[2]); op(p[4], p[5]); op(p[3], p[4]);
op(p[4], p[5]); op(p[0], p[3]); op(p[2], p[5]); op(p[2], p[3]); op(p[1], p[4]);
op(p[1], p[2]); op(p[3], p[4]); op(p[7], p[8]); op(p[6], p[7]); op(p[7], p[8]);
op(p[10], p[11]); op(p[9], p[10]); op(p[10], p[11]); op(p[6], p[9]); op(p[8], p[11]);
op(p[8], p[9]); op(p[7], p[10]); op(p[7], p[8]); op(p[9], p[10]); op(p[0], p[6]);
op(p[4], p[10]); op(p[4], p[6]); op(p[2], p[8]); op(p[2], p[4]); op(p[6], p[8]);
op(p[1], p[7]); op(p[5], p[11]); op(p[5], p[7]); op(p[3], p[9]); op(p[3], p[5]);
op(p[7], p[9]); op(p[1], p[2]); op(p[3], p[4]); op(p[5], p[6]); op(p[7], p[8]);
op(p[9], p[10]); op(p[13], p[14]); op(p[12], p[13]); op(p[13], p[14]); op(p[16], p[17]);
op(p[15], p[16]); op(p[16], p[17]); op(p[12], p[15]); op(p[14], p[17]); op(p[14], p[15]);
op(p[13], p[16]); op(p[13], p[14]); op(p[15], p[16]); op(p[19], p[20]); op(p[18], p[19]);
op(p[19], p[20]); op(p[21], p[22]); op(p[23], p[24]); op(p[21], p[23]); op(p[22], p[24]);
op(p[22], p[23]); op(p[18], p[21]); op(p[20], p[23]); op(p[20], p[21]); op(p[19], p[22]);
op(p[22], p[24]); op(p[19], p[20]); op(p[21], p[22]); op(p[23], p[24]); op(p[12], p[18]);
op(p[16], p[22]); op(p[16], p[18]); op(p[14], p[20]); op(p[20], p[24]); op(p[14], p[16]);
op(p[18], p[20]); op(p[22], p[24]); op(p[13], p[19]); op(p[17], p[23]); op(p[17], p[19]);
op(p[15], p[21]); op(p[15], p[17]); op(p[19], p[21]); op(p[13], p[14]); op(p[15], p[16]);
op(p[17], p[18]); op(p[19], p[20]); op(p[21], p[22]); op(p[23], p[24]); op(p[0], p[12]);
op(p[8], p[20]); op(p[8], p[12]); op(p[4], p[16]); op(p[16], p[24]); op(p[12], p[16]);
op(p[2], p[14]); op(p[10], p[22]); op(p[10], p[14]); op(p[6], p[18]); op(p[6], p[10]);
op(p[10], p[12]); op(p[1], p[13]); op(p[9], p[21]); op(p[9], p[13]); op(p[5], p[17]);
op(p[13], p[17]); op(p[3], p[15]); op(p[11], p[23]); op(p[11], p[15]); op(p[7], p[19]);
op(p[7], p[11]); op(p[11], p[13]); op(p[11], p[12]);
dst[j] = (T)p[12];
}
if( limit == size.width )
break;
//下面的实现是个人猜测，没经过严格认证
//VecOp应该是可以一次性处理VecOp::SIZE对数值的比较，利用opencv的优化指令
for( ; j <= size.width - VecOp::SIZE - cn*2; j += VecOp::SIZE ) //每次步进VecOp::SIZE个数，断点看VecOp::SIZE是16
{
VT p[25];
for( k = 0; k < 5; k++ )
{
const T* rowk = row[k];
//猜测vop.load是一次性加载VecOp::SIZE(16)个内存数据
p[k*5] = vop.load(rowk+j-cn*2); p[k*5+1] = vop.load(rowk+j-cn);
p[k*5+2] = vop.load(rowk+j); p[k*5+3] = vop.load(rowk+j+cn);
p[k*5+4] = vop.load(rowk+j+cn*2);
}
//p[i]指向连续的16个内存数据，vop一次性对比16个数，并将较小的放到第一个参数，较大的放到第二个参数
vop(p[1], p[2]); vop(p[0], p[1]); vop(p[1], p[2]); vop(p[4], p[5]); vop(p[3], p[4]);
vop(p[4], p[5]); vop(p[0], p[3]); vop(p[2], p[5]); vop(p[2], p[3]); vop(p[1], p[4]);
vop(p[1], p[2]); vop(p[3], p[4]); vop(p[7], p[8]); vop(p[6], p[7]); vop(p[7], p[8]);
vop(p[10], p[11]); vop(p[9], p[10]); vop(p[10], p[11]); vop(p[6], p[9]); vop(p[8], p[11]);
vop(p[8], p[9]); vop(p[7], p[10]); vop(p[7], p[8]); vop(p[9], p[10]); vop(p[0], p[6]);
vop(p[4], p[10]); vop(p[4], p[6]); vop(p[2], p[8]); vop(p[2], p[4]); vop(p[6], p[8]);
vop(p[1], p[7]); vop(p[5], p[11]); vop(p[5], p[7]); vop(p[3], p[9]); vop(p[3], p[5]);
vop(p[7], p[9]); vop(p[1], p[2]); vop(p[3], p[4]); vop(p[5], p[6]); vop(p[7], p[8]);
vop(p[9], p[10]); vop(p[13], p[14]); vop(p[12], p[13]); vop(p[13], p[14]); vop(p[16], p[17]);
vop(p[15], p[16]); vop(p[16], p[17]); vop(p[12], p[15]); vop(p[14], p[17]); vop(p[14], p[15]);
vop(p[13], p[16]); vop(p[13], p[14]); vop(p[15], p[16]); vop(p[19], p[20]); vop(p[18], p[19]);
vop(p[19], p[20]); vop(p[21], p[22]); vop(p[23], p[24]); vop(p[21], p[23]); vop(p[22], p[24]);
vop(p[22], p[23]); vop(p[18], p[21]); vop(p[20], p[23]); vop(p[20], p[21]); vop(p[19], p[22]);
vop(p[22], p[24]); vop(p[19], p[20]); vop(p[21], p[22]); vop(p[23], p[24]); vop(p[12], p[18]);
vop(p[16], p[22]); vop(p[16], p[18]); vop(p[14], p[20]); vop(p[20], p[24]); vop(p[14], p[16]);
vop(p[18], p[20]); vop(p[22], p[24]); vop(p[13], p[19]); vop(p[17], p[23]); vop(p[17], p[19]);
vop(p[15], p[21]); vop(p[15], p[17]); vop(p[19], p[21]); vop(p[13], p[14]); vop(p[15], p[16]);
vop(p[17], p[18]); vop(p[19], p[20]); vop(p[21], p[22]); vop(p[23], p[24]); vop(p[0], p[12]);
vop(p[8], p[20]); vop(p[8], p[12]); vop(p[4], p[16]); vop(p[16], p[24]); vop(p[12], p[16]);
vop(p[2], p[14]); vop(p[10], p[22]); vop(p[10], p[14]); vop(p[6], p[18]); vop(p[6], p[10]);
vop(p[10], p[12]); vop(p[1], p[13]); vop(p[9], p[21]); vop(p[9], p[13]); vop(p[5], p[17]);
vop(p[13], p[17]); vop(p[3], p[15]); vop(p[11], p[23]); vop(p[11], p[15]); vop(p[7], p[19]);
vop(p[7], p[11]); vop(p[11], p[13]); vop(p[11], p[12]);
vop.store(dst+j, p[12]);
}
limit = size.width;
}
}
}
}
可以看到，medianBlur_SortNet函数有两个模板类Op和VecOp，也就是MinMax8u和MinMaxVec8u，它们的作用是可以快速的比较两个数，并把较小的数赋给第一个参数，较大的数赋给第二个参数，MinMax8u只对比一组数据，而MinMaxVec8u可以同时（并行）对比16个数据，它应该是采用了SIMD指令。他们的实现会根据不同硬件调用不同的优化指令接口。
函数中有一个大循环
int limit = useSIMD ? cn*2 : size.width; //根据是否使用SIMD来决定limit大小
for(j = 0;; )
{
//用op（一次比较一组数据）
for( ; j < limit; j++ )
{
//...
}
//用vop（一次比较VecOp::SIZE个数据）
for( ; j <= size.width - VecOp::SIZE - cn*2; j += VecOp::SIZE )
{
//...
}
}
我们一定会有个疑问，为什么非要先用op对比limit个数，后面才用vop？一开始就用vop不是更快吗？我个人的猜测，是因为需要补边的原因。如下图
如图，在求边界上的像素的中值时，如图中的红点，边界外的邻域实际上是不存在的，如图中虚线部分，代码中row[0]和row[1]代表边界外的像素，其实他们都是指向图像第一行，也就是row[2]。利用vop时，会一次性加载16个连续内存数据，因此在计算边界上的中值时，不能用vop。
再看下第二个循环，用vop，顾名思义，就是数组的操作（vector operation）。它可以一次性处理16个（VecOp::SIZE）个数据，可以加速处理，如下图：
