<span style="font-family: FangSong_GB2312; background-color: rgb(255, 255, 255);">这篇文章给出的代码将能解决如下问题：</span>
<span style="font-family: Arial; font-size: 14px; color: rgb(51, 204, 0);">% 我们有一系列灰度图像 images % 现在想要计算一幅中值图像作为背景 </span><span style="font-family:Arial;font-size:14px;color:#999999;">im1 = randi（500600); im2 = randi(500,600); im3 = randi(500,600); images = cat(3,im1,im2,im3); median_background = median(images,[],3);</span>
在工程应用中，常常需要对固定相机拍摄的n幅图像做一个背景建模，然后用图像减去背景就得到运动变化部分。背景建模的方法很多，常用的有GM，GMM，这些经典算法在openCV里都集成了，但是速度上还不是特别快。对于类似的情况其实有更简单的方法建背景--中值背景就是在工程中用得最多，效果最理想的方法之一。
注意，这里说的中值，就是median，不是mean（均值）。
计算中值还是比较麻烦的，比如说对5个数[20 30 32 23 90]计算中值，首选需要对其排序，用冒泡排序算法得到的结果应该是[20 23 30 32 90]，然后取中间一个数，即30。
计算中值图像的时候呢，就是要对n幅图像（相同尺寸）上每个像素位置分别求中值，大概就是：
<span style="font-family:Arial Black;font-size:18px;color:#33cc00;">F = median_background( A1, A2, A3, ..., An )</span>
<span style="font-family:Arial Black;font-size:18px;color:#33cc00;"></span><pre name="code" class="html" style="text-align: center;"><span style="font-family:Arial Black;font-size:18px;color:#33cc00;">s.t., F(x,y) = median( [ A1(x,y), A2(x,y), ..., An(x,y) ] )</span>
为了达到这个目的，我尝试过对每个位置进行冒泡排序，分别算，但是很慢；然后，我又将矩阵编程的手段用上，快了很多，但是面临这样一个难题：比如刚才我已经对A1到An共n幅图像计算出了一个中值图像，那么现在我又得到了第n+1幅图像，我希望更新这个背景：
<span style="font-family:Arial Black;font-size:18px;color:#33cc00;">F_new = median_background( A2, A3, A4, ..., A(n+1) )</span>如果从头开始计算F_new，感觉之前的结果没用上，有点可惜。为了减少计算量，我希望这个算法是增量的，即F_new可以在F的基础上得到。那么我想到了在opencv和matlab里面都有一个index_filter之类的函数，就是对数组元素进行排序，只保存排序的索引，而不改变数组元素本身。基于这个思想，我写出了本文代码。
下面贴代码：
<span style="font-size:18px;">// 文件 yuMedianBackground.h</span>
<span style="font-size:18px;">#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</span><span style="font-size:18px;">// 文件 yuMedianBackground.cpp</span>
<span style="font-size:18px;">#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
//}
}
}
}
}
// extract
background.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 beforehand
assert( frame.size()==background.size() && frame.type()==CV_8UC1 );
// --------- define variables
int 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 than
Mat 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++);
// backward
for( 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)
}
else break;
}
if( pt4!=pt5 )
continue;
// forward
pt4 = pt5;
for( pt5=pt2+num1; pt4!=pt5; ){
uchar &a = *(pt4+1);
if( pt1[a]<tmp ){
*(pt4++) = a; a = num1;
}
else
break;
}
}
}
// --------- extract
int 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;
}</span>里面主要接口是calc()和calcOn()两个函数，分别是计算F以及在F的基础上计算F_new的。
看这个代码可能会让你们感到头疼，实际上我现在去看也觉得头疼，主要是逻辑方面稍微麻烦点，在调试代码的时候让我花了不少时间，现在用着没法现bug就再没看过了。
速度上我自信这是你们见过的最快的或接近最快的计算中值图像的程序了，如果谁有更快的函数，希望不要吝啬传上来给我们看看，可以发csdn资源或博客，然后记得把网址贴在回复里，也可以给我邮件
如果有人用着本文代码发现了bug，请一定告诉我！
原文：http://blog.csdn.net/j56754gefge/article/details/44598381