C++计算中值图像，增量算法

这篇文章最初发表在http://blog.csdn.net/j56754gefge/article/details/44598381，均是我个人原创，转载请注明出处！

% input: im1, im2, im3, im4, im5% output: background image FF = im1;for i=1:3    A = cat(3,im1(:,:,i),im2(:,:,i),im3(:,:,i),im4(:,:,i),im5(:,:,i));    F(:,:,i) = median(A,3);endfigure,imshow(F)

        在工程应用中，常常需要对固定相机拍摄的n幅图像做一个背景建模，然后用图像减去背景就得到运动变化部分。背景建模的方法很多，常用的有GM，GMM，这些经典算法在openCV里都集成了，但是速度上还不是特别快。对于类似的情况其实有更简单的方法建背景--中值背景就是在工程中用得最多，效果最理想的方法之一。注意，这里说的中值，就是median，不是mean（均值）。

       计算中值的过程：举个例子，对5个数[20 30 32 23 90]计算中值，首选需要对其排序，用冒泡排序算法得到的结果应该是[20 23 30 32 90]，然后取中间一个数，即30。

       计算中值图像( F )的时候，就是要对n幅图像( A1, A2, ..., An )上每个像素位置分别求中值，大概就是：

F = median_background( A1, A2, A3, ..., An )

s.t., F(x,y) = median( [ A1(x,y), A2(x,y), ..., An(x,y) ] )

        另外，我们希望算法是增量式的：比如刚才我已经对A1到An共n幅图像计算出了一个中值图像F，那么现在我又得到了第n+1幅图像，我希望更新这个背景：

F_new = median_background( A2, A3, A4, ..., A(n+1) )

       如果从头开始计算F_new，感觉之前的结果没用上，有点可惜。为了减少计算量，我希望这个算法是增量的，即F_new可以在F的基础上得到。那么我想到了在opencv和matlab里面都有一个index_filter之类的函数，就是对数组元素进行排序，只保存排序的索引，而不改变数组元素本身。基于这个思想，我写出了本文代码。

下面贴代码：

// 文件 yuMedianBackground.h

#ifndef YU_MEDIAN_BACKGROUND_H#define YU_MEDIAN_BACKGROUND_H#include "cv.h"using namespace cv;// **********************************// calculating a median image of a set of images.// require the image set to be CV_8UC1 type.// by : Yu Xianguo, Oct.29.2014// **********************************class yuMedianBackground{public:yuMedianBackground(){};Mat  calc( const Mat *frames, int num_frames );Mat  calcOn( const Mat &frame );Mat compounded;Mat indexes;Mat background;};#endif

// 文件 yuMedianBackground.cpp

#include "yuMedianBackground.h"using namespace std;Mat  yuMedianBackground::calc( const Mat *frames, int num_frames ){assert( num_frames>0 && num_frames<32 && frames->type()==CV_8UC1 ); cv::merge( frames, num_frames, compounded );int rows = frames[0].rows, cols = frames[0].cols;indexes.create( rows, cols, CV_8UC(num_frames) );std::vector<uchar> v(num_frames);uchar num = num_frames, num1 = num_frames - 1;for( uchar i=0; i<num; i++ ) v[i]=i; indexes.setTo(v);uchar *pt1 = compounded.data, *pt2 = indexes.data;int gap1 = compounded.step-cols*num_frames;int gap2 = indexes.step-cols*num_frames;uchar tmp, *pt3;for( int y=0; y++<rows; pt1+=gap1, pt2+=gap2 ){for( int x=0; x++<cols; pt1+=num, pt2+=num ){for( uchar i=0; i<num1; i++ ){pt3 = pt2;for( uchar j=0, nd=num1-i; j<nd; j++ ){uchar &a = *(pt3++); // pt2[j]if( pt1[*pt3]<pt1[a] ){tmp = a; a = *pt3; *pt3 = tmp;}//uchar &a = pt2[j], &b = pt2[j+1];//if( pt1[b]<pt1[a] ){//tmp = a; a = b; b = tmp; // exchange a & b//}}}}}// extractbackground.create( rows, cols, CV_8UC1 );int mid = (int)( (num_frames-1)/2.f+.5f );pt1 = compounded.data, pt2 = indexes.data+mid,pt3 = background.data;int gap3 = background.step-cols;for( int y=0; y++<rows; pt1+=gap1, pt2+=gap2, pt3+=gap3 )for( int x=0; x++<cols; pt1+=num, pt2+=num )*(pt3++) = pt1[*pt2];return background;}Mat  yuMedianBackground::calcOn( const Mat &frame ){// --------- need calc() to be called beforehandassert( frame.size()==background.size() && frame.type()==CV_8UC1 );// --------- define variablesint rows = frame.rows, cols = frame.cols, num_frames = compounded.channels();uchar *pt1 = compounded.data, *pt2 = indexes.data, *pt3 = frame.data;int gap1 = compounded.step-cols*num_frames, gap2 = indexes.step-cols*num_frames, gap3 = frame.step-cols;uchar i, num = num_frames, num1 = num_frames-1, tmp;// --------- update "compounded" with "frame"// code 1.//y = 0; while(y++<rows){//x = 0; while(x++<cols){//i = 1; while(i++<num){//*pt1 = *(pt1+1); // overwrite frame[0]//pt1++;//}//*(pt1++) = *(pt3++);//}//pt1 += gap1, pt3 += gap3;//}//pt1 = compounded.data, pt3 = frame.data; // recover// code 2. alternative to code(1), but much faster thanMat from[] = { compounded, frame };Mat to[] = { compounded };int *from_to = new int [2*num_frames];int *pt4 = from_to;for( i=0; i<num; i++ ){*(pt4++) = i+1;*(pt4++) = i;}cv::mixChannels( from, 2, to, 1, from_to, num_frames );delete [] from_to;// --------- update "indexes"for( int y=0; y++<rows; pt1+=gap1, pt2+=gap2, pt3+=gap3 ){for( int x=0; x++<cols; pt1+=num, pt2+=num ){uchar *pt4 = pt2, *pt5;i = 0; while(i++<num){if( *pt4 )(*(pt4++))--;else*(pt5=pt4++) = num1;}tmp = *(pt3++);// backwardfor( pt4=pt5; pt4!=pt2; ){uchar &a = *(pt4-1); // corresponding index to position (i-1)if( tmp<pt1[a] ){ // corresponding gray value to index (a)*(pt4--)= a; a = num1; // exchange index in positions (i-1) & (i)}elsebreak;}if( pt4!=pt5 )continue;// forwardpt4 = pt5;for( pt5=pt2+num1; pt4!=pt5; ){uchar &a = *(pt4+1);if( pt1[a]<tmp ){*(pt4++) = a; a = num1;}elsebreak;}}}// --------- extractint mid = (int)( (num_frames-1)/2.f+.5f );pt1 = compounded.data, pt2 = indexes.data+mid;pt3 = background.data, gap3 = background.step-cols;for( int y=0; y++<rows; pt1+=gap1, pt2+=gap2, pt3+=gap3 )for( int x=0; x++<cols; pt1+=num, pt2+=num )*(pt3++) = pt1[*pt2];return background;}

1.里面主要接口是calc()和calcOn()两个函数，分别是计算F以及在F的基础上计算F_new的。

2.代码是我花了一些时间写的，在工程中也用过，还没发现bug。如果有谁发现了代码中存在bug，请一定在回复里告诉我，我会修改本文代码，以惠及后来者。

3.速度上我自信这是你们见过的最快的或接近最快的计算中值图像的程序了，如果谁有更快的代码，希望不要吝啬传上来给我们看看，可以发csdn资源或博客，然后记得把网址贴在回复里，让我学习一下