立体匹配——Box Filter

基本原理

在给定的滑动窗口大小下，对每个窗口内的像素值进行快速相加求和。它可以将复杂度为0(MN)的求和、求方差等运算降低到O(1)或近似于O(1)的复杂度。

BoxFilter初始化过程

1）给定一张图像，大小为M*N，首先确定待求矩形模板的宽高(m,n),如图紫色区域。图中每个黑色方块代表一个像素。
2）开辟一段大小为M的数组，记为buff，用来存储计算过程的中间变量，用红色方块表示。
3）将矩形模板(紫色)从左上角（0，0)开始，逐像素向右滑动，到达行末时，矩形移动到下一行的开头(0,1)，如此反复，每移动到一个新位置时，计算矩形内的像素和，保存在数组A中。以(0,0)位置为例进行说明：首先将绿色矩形内的每一列像素求和，结果放在buff内（红色方块），再对蓝色矩形内的像素求和，结果即为紫色特征矩形内的像素和，把它存放到数组A中，如此便完成了第一次求和运算。
4）每次紫色矩形向右移动时，实际上就是求对应的蓝色矩形的像素和，此时只要把上一次的求和结果减去蓝色矩形内的第一个红色块，再加上它右面的一个红色块，就是当前位置的和了，用公式表示 sum[i] = sum[i-1] - buff[x-1] + buff[x+m-1]。
5）当紫色矩形移动到行末时，需要对buff进行更新。因为整个绿色矩形下移了一个像素，所以对于每个buff[i], 需要加上一个新进来的像素，再减去一个出去的像素，然后便开始新的一行的计算了。

源码

// 计算累计和
Mat CumSum( const Mat& src, const int d )
{
	int H = src.rows;
	int W = src.cols;
	Mat dest = Mat::zeros( H, W, src.type() );
    // 列 求 和
	if( d == 1 ) {
		// summation over column
		for( int y = 0; y < H; y ++ ) {
			double* curData = ( double* ) dest.ptr<double>( y );
			double* preData = ( double* ) dest.ptr<double>( y );
			if( y ) {
				// 不是第一行的话
				preData = ( double* ) dest.ptr<double>( y - 1 ); 
			}
			double* srcData = ( double* ) src.ptr<double>( y );
			for( int x = 0; x < W; x ++ ) {
				curData[ x ] = preData[ x ] + srcData[ x ];
			}
		}
	} else {
		// 行求和
		for( int y = 0; y < H; y ++ ) {
			double* curData = ( double* ) dest.ptr<double>( y );
			double* srcData = ( double* ) src.ptr<double>( y );
			for( int x = 0; x < W; x ++ ) {
				if( x ) {
					curData[ x ] = curData[ x - 1 ] + srcData[ x ];
				} else {
					curData[ x ] = srcData[ x ];
				}
			}
		}
	}
	return dest;
}


//Box Filter函数
Mat BoxFilter( const Mat& imSrc, const int r )
{
	int H = imSrc.rows;
	int W = imSrc.cols;
	// 图像大小应大于滤波窗口大小
	CV_Assert( W >= r && H >= r );
	Mat imDst = Mat::zeros( H, W, imSrc.type() );
	// Y轴累计和
	Mat imCum = CumSum( imSrc, 1);
	//( [ 0, r ], [r + 1, H - r - 1], [ H - r, H ] )沿Y方向的差
	for( int y = 0; y < r + 1; y ++ ) {
		double* dstData = ( double* ) imDst.ptr<double>( y );
		double* plusData = ( double* ) imCum.ptr<double>( y + r );
		for( int x = 0; x < W; x ++ ) {
			dstData[ x ] = plusData[ x ];
		}
	}
	for( int y = r + 1; y < H - r; y ++ ) {
		double* dstData = ( double* ) imDst.ptr<double>( y );
		double* minusData = ( double*  ) imCum.ptr<double>( y - r - 1);
		double* plusData = ( double* ) imCum.ptr<double>( y + r );
		for( int x = 0; x < W; x ++ ) {
			dstData[ x ] = plusData[ x ] - minusData[ x ];
		}
	}
	for( int y = H - r; y < H; y ++ ) {
		double* dstData = ( double* ) imDst.ptr<double>( y );
		double* minusData = ( double*  ) imCum.ptr<double>( y - r - 1);
		double* plusData = ( double* ) imCum.ptr<double>( H - 1 );
		for( int x = 0; x < W; x ++ ) {
			dstData[ x ] = plusData[ x ] - minusData[ x ];
		}
	}

	// cumulative sum over X axis
	imCum = CumSum( imDst, 2 );
	for( int y = 0; y < H; y ++ ) {
		double* dstData = ( double* ) imDst.ptr<double>( y );
		double* cumData = ( double* ) imCum.ptr<double>( y );
		for( int x = 0; x < r + 1; x ++ ) {
			dstData[ x ] = cumData[ x + r ];
		}
		for( int x = r + 1; x < W - r; x ++ ) {
			dstData[ x ] = cumData[ x + r ] - cumData[ x - r - 1 ];
		}
		for( int x = W - r; x < W; x ++ ) {
			dstData[ x ] = cumData[ W - 1 ] - cumData[ x - r - 1 ];
		}
	}
	return imDst;
}