34运动与跟踪-运动模板

34运动与跟踪-运动模板

OpenCV中运动与跟踪这一章节中，在前面的介绍中，主要给出了LK光流法，以及基于概率统计，窗口搜索的meanshif算法以及meanshif算法的改进版-camshift算法，这两天主要在看运动模板的跟踪运动方法，下面就简要的介绍下。
运动模板的方法是美国的MIT实验室提出来的，是一种有效的跟踪普通运动的方法，尤其可应用在姿态识别中。运动模板的方法首先需要的是知道物体的轮廓，而轮廓的获取可以有很多的方法，关于轮廓的获取的方法在此不做多余的赘述。
运动模板的方法程序的大致编程思路应该是是这样的：获得当前帧与上一帧的差，然后对差图像进行二值化；更新运动历史图像；然后计算运动历史图像的梯度方向，并将整个运动分割为独立的运动部分，在用一个结构序列标记每一个运动分割，最后计算出选择区域的全局运动方向，从而获得运动目标的质心位置与运动方向。
关于上面的思路，给出几个主要的函数，对其参数进行介绍。

---

void cvUpdateMotionHistory(
    const CvArr* silhouette,//非0像素代表前景物体的最新的分割轮廓
    CvArr* mhi,//运动模板，也即运动历史图像,记录的是时间信息
    double timestamp,//当前系统时间
    double duration//运动跟踪的最大持续时间
);

函数 cvUpdateMotionHistory 用下面方式更新运动历史图像：
函数只是更新 像素点的运动历史。也就说更新的不是图像，而是对图像中像素点运动情况的更新。
　　silhouette（x,y） !=0时，即该像素点发生运动，所以要对其进行更新，即mhi(x,y) = timestamp 表示运动发生的时刻
　　silhouette（x,y） =0时，即该像素点未发生运动，但还需检测对该点的跟踪时间是否超过了预设最大跟踪时间，即判断mhi(x,y)与timestamp -duration的大小。此时mhi(x,y)即为该点最近一次发生运动的时刻值，如其小于timestamp-duration,表示该点运动时刻已 经超出跟踪时间，故可以舍弃。
　　而当mhi(x,y)大于或者等于timestamp-duration时，表示该点此刻虽未发生运动，但还在跟踪时间内，所以不对该点发生运动的时间标记进行操作。

---

一旦运动模板记录了不同时间的物体轮廓，就快成用计算mhi图像的梯度来获取全局运动信息。即有函数：

void cvCalcMotionGradient(
    const CvArrs* mhi,//运动历史图像
    CvArr* mask,//非0值代表此点处的梯度有效
    CvArr* orientation,//每一点梯度方向的角度
    double delta1,
    double delta2,//delta1,2分别是允许的最小和最大梯度值
    int aperture_size=3//设置梯度算子的宽和高
);
注：sobel梯度是带方向的：

 1. 由亮到暗为正方向（大-->小）
 2. 由暗到亮为负方向（小-->大）
 3. 大而无效的梯度值是由背景与运动轮廓相交处产生的，即没有相对运动的像素点
 4. 期望的梯度值是连续调用cvUpdateMotionHistory()的每个轮廓间当前时间标记的平均数；将deltal1设置为小于此平均值的一半，deltal2设置为大于此平均值的一半比较适合。

---

梯度方向求出后，即可以计算全局运动方向，给出如下函数：

double cvCalcGlobalOrientation(
    const CvArr* orientation,
    const CvArr* mask,
    const CvArr* mhi,//前三个参数都是根据上面两个函数得出来的
    double timestamp,
    double duration
);

---

关于运动模板后面的数学原理，我其实还没找到相关的文献阅读，有些数学原理地方不是特别明白。下面给出完整的运动模板示例程序。

#include "opencv2/video/tracking.hpp"
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc_c.h"
#include <time.h>
#include <stdio.h>
#include <ctype.h>
//不同的跟踪参数  
const double MHI_DURATION = 0.5;
const double MAX_TIME_DELTA = 0.5;
const double MIN_TIME_DELTA = 0.05;
// 用于运动检测的循环帧数，与机器速度及FPS设置有关  
const int N = 2;
IplImage **buf = 0;
int last = 0;//上一帧图像序号
// 临时图像  
IplImage *mhi = 0; // MHI: 运动历史图像  
IplImage *orient = 0; // 方向  
IplImage *mask = 0; // 有效的运动掩码  
IplImage *segmask = 0; // 运动分割映射  
CvMemStorage* storage = 0; // 临时存储区  // parameters:  
                           //  img - input video frame  
                           //  dst - resultant motion picture  
                           //  args - optional parameters  
void  update_mhi(IplImage* img, IplImage* dst, int diff_threshold) {
    double timestamp = (double)clock() / CLOCKS_PER_SEC; // 获取当前时间,以秒为单位  
    CvSize size = cvSize(img->width, img->height); // 获取当前帧尺寸  
    int i, idx1 = last, idx2;
    IplImage* silh;
    CvSeq* seq;
    CvRect comp_rect;
    double count;
    double angle;
    CvPoint center;
    double magnitude;
    CvScalar color;

    // 开始时为图像分配内存 or 帧尺寸改变时重新分配内存  
    if (!mhi || mhi->width != size.width || mhi->height != size.height) {
        if (buf == 0) {//第一次分配内存

            buf = (IplImage**)malloc(N * sizeof(buf[0]));//分配内存
            memset(buf, 0, N * sizeof(buf[0]));
        }

        for (i = 0; i < N; i++) {//每次进入都做一次清除操作

            cvReleaseImage(&buf[i]);
            buf[i] = cvCreateImage(size, IPL_DEPTH_8U, 1);
            cvZero(buf[i]);
        }
        //执行算法前先清空图片内存
        cvReleaseImage(&mhi);
        cvReleaseImage(&orient);
        cvReleaseImage(&segmask);
        cvReleaseImage(&mask);

        mhi = cvCreateImage(size, IPL_DEPTH_32F, 1);//创建运动历史图像
        cvZero(mhi); // clear MHI at the beginning  
        orient = cvCreateImage(size, IPL_DEPTH_32F, 1);// 方向
        segmask = cvCreateImage(size, IPL_DEPTH_32F, 1);// 运动分割映射
        mask = cvCreateImage(size, IPL_DEPTH_8U, 1);// 有效的运动掩码  
    }

    cvCvtColor(img, buf[last], CV_BGR2GRAY); //RGB帧图像格式转换为gray  

    idx2 = (last + 1) % N; // index of (last - (N-1))th frame  last为上一帧index
    last = idx2;

    silh = buf[idx2];
    // 相邻两帧的差  
    cvAbsDiff(buf[idx1], buf[idx2], silh);//把运动的物体检测出来

    cvThreshold(silh, silh, diff_threshold, 1, CV_THRESH_BINARY); // 对差图像做二值化  
    cvUpdateMotionHistory(silh, mhi, timestamp, MHI_DURATION); // 更新运动历史  

                                                               // convert MHI to blue 8u image  
                                                               // cvCvtScale的第四个参数 shift = (MHI_DURATION - timestamp)*255./MHI_DURATION  
                                                               // 控制帧差的消失速率  
    cvCvtScale(mhi, mask, 255. / MHI_DURATION,
        (MHI_DURATION - timestamp)*255. / MHI_DURATION);//将时间信息转换为8 bit的图像,时间越长数值越小。
    /****************************************for test***********************************************
    cvNamedWindow("Motion1", 1);
    cvShowImage("Motion1", mask);
    while (true)
    {
        if (cvWaitKey(10) >= 0)
            break;
    }
    cvDestroyWindow("Motion");
    printf("%f\n%f\n", 255. / MHI_DURATION, (MHI_DURATION - timestamp)*255. / MHI_DURATION);
    ****************************************for test***********************************************/
    cvZero(dst);
    cvMerge(mask, 0, 0, 0, dst);//将mask复制到蓝色通道上去，用于显示

    // B,G,R,0 convert to BLUE image  

    // 计算运动的梯度方向以及正确的方向掩码  
    // Filter size = 3  
    cvCalcMotionGradient(mhi, mask, orient,
        MAX_TIME_DELTA, MIN_TIME_DELTA, 3);

    if (!storage)
        storage = cvCreateMemStorage(0);
    else
        cvClearMemStorage(storage);

    // 运动分割： 获得运动部件的连续序列  
    seq = cvSegmentMotion(mhi, segmask, storage, timestamp, MAX_TIME_DELTA);
    for (i = -1; i < seq->total; i++) {

        if (i < 0) {        // 对整幅图像操作  
            comp_rect = cvRect(0, 0, size.width, size.height);
            color = CV_RGB(255, 255, 255);//全局的运动方向为白色
            magnitude = 100;  // 画线长度以及圆半径的大小控制  
        }
        else {          // 第i个运动组件  
            comp_rect = ((CvConnectedComp*)cvGetSeqElem(seq, i))->rect;
            // 去掉小区域的部分  
            if (comp_rect.width + comp_rect.height < 100)
                continue;
            color = CV_RGB(255, 0, 0);//局部的运动方向为红色
            magnitude = 30;
            //if(seq->total > 0) MessageBox(NULL,"Motion Detected",NULL,0);  
        }
        // 选择组件ROI  
        cvSetImageROI(silh, comp_rect);//设置为每个区域的感兴趣区域
        cvSetImageROI(mhi, comp_rect);
        cvSetImageROI(orient, comp_rect);
        cvSetImageROI(mask, comp_rect);

        // 在选择的区域内，计算运动方向  
        angle = cvCalcGlobalOrientation(orient, mask, mhi, timestamp, MHI_DURATION);
        angle = 360.0 - angle;  // adjust for images with top-left origin  
                                // 在轮廓内计算点数  
                                // Norm(L1) = 所有像素值的和  
        count = cvNorm(silh, 0, CV_L1, 0);//silh帧差图像，count:帧差图像中像素点的个数
        /*如果 arr2 为空（NULL），函数 cvNorm 计算 arr1 的绝对范数：
            norm = ||arr1||C = maxI abs(arr1(I)), 如果 normType = CV_C
            norm = ||arr1||L1 = sumI abs(arr1(I)), 如果 normType = CV_L1
            norm = ||arr1||L2 = sqrt( sumI arr1(I)2), 如果 normType = CV_L2*/

        cvResetImageROI(mhi);//重置感兴趣区域
        cvResetImageROI(orient);
        cvResetImageROI(mask);
        cvResetImageROI(silh);

        // 检查小运动的情形  
        if (count < comp_rect.width*comp_rect.height * 0.05)  //  像素的5%  如果当前移动太小，则不重新绘制运动方向
            continue;

        // 画一个带箭头的记录以表示方向  
        center = cvPoint((comp_rect.x + comp_rect.width / 2), (comp_rect.y + comp_rect.height / 2));

        cvCircle(dst, center, cvRound(magnitude*1.2), color, 3, CV_AA, 0);//绘制圆形，由magnitude控制圆形大小
        cvLine(dst, center, cvPoint(cvRound(center.x + magnitude * cos(angle*CV_PI / 180)),
            cvRound(center.y - magnitude * sin(angle*CV_PI / 180))), color, 3, CV_AA, 0);//对一个double型的数进行四舍五入，并返回一个整型数
            //center起点，cvPoint终点，
    }
}

int main(int argc, char** argv) {

    IplImage* motion = 0;
    CvCapture* capture = 0;

    if (argc == 1 || (argc == 2 && strlen(argv[1]) == 1 && isdigit(argv[1][0])))//从摄像头中获取视频信息
        capture = cvCaptureFromCAM(argc == 2 ? argv[1][0] - '0' : 0);
    else if (argc == 2)
        capture = cvCaptureFromAVI(argv[1]);//从文件中获取视频信息

    if (capture) {

        cvNamedWindow("Motion", 1);
        for (;;) {

            IplImage* image = cvQueryFrame(capture);//获取图像帧
            if (!image)
                break;

            if (!motion)
            {
                motion = cvCreateImage(cvSize(image->width, image->height), 8, 3);
                cvZero(motion);
                motion->origin = image->origin;
            }

            update_mhi(image, motion, 30);//调用图像模板接口
            cvShowImage("Motion", motion);

            if (cvWaitKey(10) >= 0)
                break;
        }
        cvReleaseCapture(&capture);
        cvDestroyWindow("Motion");
    }

    return 0;
}
#ifdef _EiC
main(1, "motempl.c");
#endif

结果图如下：