一种简单并行的轨迹关联算法及其在小目标检测中的应用

    轨迹关联是利用目标在时空域中状态变化的连续性进行目标检测和噪声滤除的有效手段，常用于多目标检测和跟踪中，一个常见的应用是基于车载视觉的行人以及车辆的检测和跟踪。而对于大量的、缺乏丰富纹理特征的目标的情况，比如我在另一篇博文里提到的CLIF数据库，一副图像中存在上千个汽车，一般的多目标跟踪算法都将难以适用。比CLIF的情况更差的，是如下图所示的卫星影像图像中的运动目标检测：

    首先贴一下这里用到的几个视频的来源，都是公开的，广告视频，被压缩过，且有各种水印。

    视频一：http://www.berlin-space-tech.com/index.php_id=26.html

    视频二：http://www.firstimagery.skybox.com/hd-video/

    也可以去YouTube上分别搜索LAPAN-TUBSAT和Skybox，都能找到很多。

    上图是一张示意图，黄色圆圈中间一坨大概15个像素的亮团其实是一架飞机，红色那些点是我们做完图像匹配以后差分得到的，其中大部分是噪声和干扰，我们的轨迹关联算法能够检测出其中的飞机，做法大概就是在连续若干幅图像里搜索运动目标轨迹，由若干约束条件来去除不合理轨迹。

    我把我们的处理结果放到了csdn资源里，都是算法一遍跑出来的，没有任何修饰。不过，算法还做不到参数自适应学习，每个视频都有一个对应的参数配置文件，都是微调，不需要太大改动。

视频1、2：变化背景，低分辨率图像上检测运动飞机 (结果已经上传到CSDN资源，点击图像打开资源页)

 

视频3、4：静止背景，较高分辨率图像上检测运动汽车

 

    视频3、4中有较多运动目标，我们用中值背景建模的方法来分割运动目标，然后用轨迹关联来检测出真实目标。我们将检测到的目标依次编号。对以检测的目标进行跟踪，如果一个目标被跟丢，则其标号消失。如果一段时间以后该目标又被检测到，则它会得到一个新的编号。视频3、4和视频1、2的整体算法框架是一致的，区别在于视频1、2是运动背景，我们队不同图像进行匹配以统一坐标系，目标搜索和跟踪是在同一个参考坐标系下完成的。四个视频分别有四个不同的参数配置文件，其中1、2的参数设置基本相同，3、4的参数设置基本相同。

    然后是我们搜索算法的c++代码，一个普通的带剪枝策略的蛮搜索算法以及一个并行加速算法。这里说的并行不是指GPU或多线程，而是算法上由原来的逐条轨迹依次测试改成批处理，其结果是一样的，只是速度快了很多。

/* File: trajectory_search.cpp */#include <vector>using namespace std;struct Candidate // target candidate{float x, y, gray;};struct SearchTemporary // temporary for searchTraces(){float x, y, gray;float vx, vy, v2;};//---------------------------------------------------------------// Exhaustive search with pruningvoid  yuTrajectorySearch_V0(vector< vector<Candidate>* > &candidates, // input, each candidate is a 3-tuple <x,y,gray>vector< vector<int> > &validPaths, // output, each path is a series of candidate indexfloat thresh_dGray, // input, thresh for gray differencefloat thresh_minVelocity, // input, thresh for min velocityfloat thresh_maxVelocity, // input, thresh for max velocityfloat thresh_cosTheta // input, thresh for included angles of consecutive velocities){const int numFrames = candidates.size();vector<int> numCandidates( numFrames );for( int i=0; i<numFrames; i++ )numCandidates[i] = candidates[i]->size();vector<int> onePath( numFrames, 0 );vector<Candidate> pathInfo( numFrames );float dgray, dx, dy, dv2, DX, DY, DV2;const float thresh_minVelocity2 = thresh_minVelocity * thresh_minVelocity;const float thresh_maxVelocity2 = thresh_maxVelocity * thresh_maxVelocity;while(1){for( int i=0; i<numFrames; i++ ){int idxi = onePath[i];pathInfo[i] = (*candidates[i])[idxi];}//int ret = testPath( pathInfo );int ret = 1;do{// Test the validity of a path. If the path is valid, return num.// Else, if the path evaluation ends at i-th vertex, then return i. i ranges between [1,num-1]vector<Candidate>::iterator itr = pathInfo.begin();Candidate &a0 = *(itr++);Candidate &a1 = *(itr++);dgray = a1.gray - a0.gray;if( dgray<0 ) dgray = - dgray;if( dgray > thresh_dGray )break;dx = a1.x - a0.x;dy = a1.y - a0.y;dv2 = dx*dx + dy*dy;if( dv2<thresh_minVelocity2 || thresh_maxVelocity2<dv2 )break;Candidate ak1 = a1;for( ret=2; ret<numFrames; ret++ ){Candidate &ak = *(itr++);dgray = ak.gray - ak1.gray;if( dgray<0 ) dgray = - dgray;if( dgray > thresh_dGray )break;DX = ak.x - ak1.x;DY = ak.y - ak1.y;DV2 = DX*DX + DY*DY;if( DV2<thresh_minVelocity2 || thresh_maxVelocity2<DV2 )break;float cosTheta = (dx*DX + dy*DY) / sqrtf(dv2*DV2);if( cosTheta<thresh_cosTheta )break;ak1 = ak;dx = DX, dy = DY, dv2 = DV2;}} while(0);if( ret==numFrames ){ // valid pathvalidPaths.push_back( onePath );ret = numFrames - 1;}// go to next pathint i;for( i=ret; i>=0; i-- ){onePath[i]++;if( onePath[i]<numCandidates[i] ) // allow traceIdx[i] in [0,Xs[i].size()]break;onePath[i] = 0;}if( i<0 )break; // all paths already have been evaluated}}//---------------------------------------------------------------// Parallel search algorithm, an accelerated version of exhaustive searchvoid  yuTrajectorySearch_V1(vector< vector<Candidate>* > &candidates, // input, each candidate is a 3-tuple <x,y,gray>vector< vector<int> > &validPaths, // output, each path is a series of candidate indexfloat thresh_dGray, // input, thresh for gray differencefloat thresh_minVelocity, // input, thresh for min velocityfloat thresh_maxVelocity, // input, thresh for max velocityfloat thresh_cosTheta // input, thresh for included angles of consecutive velocities){validPaths.clear();const float thresh_minVelocity2 = thresh_minVelocity * thresh_minVelocity;const float thresh_maxVelocity2 = thresh_maxVelocity * thresh_maxVelocity;const int numFrames = candidates.size();vector<int> num(numFrames);bool unnecessary = false;for( int i=0; i<numFrames; i++ ){num[i] = candidates[i]->size();if( !num[i] )unnecessary = true;}if( unnecessary )return;typedef std::pair<int,int> tracklet;vector< vector<tracklet> > records(numFrames-1);vector<SearchTemporary> b, b2;float vx, vy, v2, x, y, gray;float dx, dy, dv2, dgray;SearchTemporary elem;// 第一步：寻找第0帧到第1帧的合理路径vector<Candidate>::iterator candidate_0 = candidates[0]->begin();for( int i=0; i<num[0]; i++ ){Candidate &a0i = *(candidate_0++);x = a0i.x, y = a0i.y, gray = a0i.gray;vector<Candidate>::iterator candidate_1 = candidates[1]->begin();for( int j=0; j<num[1]; j++ ){Candidate &a1j = *(candidate_1++);dgray = gray<a1j.gray ? a1j.gray-gray : gray-a1j.gray;if( thresh_dGray<dgray )continue;dx = a1j.x - x;dy = a1j.y - y;dv2 = dx*dx + dy*dy;if( dv2<thresh_minVelocity2 || thresh_maxVelocity2<dv2 )continue;elem.x = a1j.x; elem.y = a1j.y; elem.gray = a1j.gray;elem.vx = dx; elem.vy = dy; elem.v2 = dv2;b.push_back( elem );records[0].push_back( tracklet(i,j) );}}// 第二步：寻找到第k帧的合理路径for( int k=2; k<numFrames; k++ ){vector<SearchTemporary>::iterator candidate_track = b.begin();for( unsigned i=0; i<b.size(); i++ ){x = candidate_track->x;y = candidate_track->y;gray = candidate_track->gray;vx = candidate_track->vx;vy = candidate_track->vy;v2 = candidate_track->v2;candidate_track++;vector<Candidate>::iterator candidate_k = candidates[k]->begin();for( int j=0; j<num[k]; j++ ){Candidate &akj = *(candidate_k++);dgray = gray<akj.gray ? akj.gray-gray : gray-akj.gray;if( thresh_dGray<dgray )continue;dx = akj.x - x;dy = akj.y - y;dv2 = dx*dx + dy*dy;if( dv2<thresh_minVelocity2 || thresh_maxVelocity2<dv2 )continue;float cosTheta = (dx*vx + dy*vy) / sqrtf(dv2*v2);if( cosTheta<thresh_cosTheta )continue;elem.x = akj.x; elem.y = akj.y; elem.gray = akj.gray;elem.vx = dx; elem.vy = dy; elem.v2 = dv2;b2.push_back( elem );records[k-1].push_back( tracklet(i,j) );}}if( b2.empty() )return; // searched through 0 to k th frame, found no valid pathb = b2;b2.clear();//printf("searchTraces -- valid paths until %02dth frame : %d\n",k,b.size());}// 回溯法查找路径int num_valid_paths = records[numFrames-2].size();validPaths.assign( num_valid_paths, vector<int>(numFrames) );for( int m=0; m<num_valid_paths; m++ ){int ind = m;vector<int> &pt = validPaths[m];for( int k=numFrames-2; k>=0; k-- ){tracklet &t = records[k][ind];pt[k+1] = t.second;ind = t.first;}pt[0] = ind;}return;}

-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------