光流演算法金字塔

① 如何使用opencv實現金字塔光流lk跟蹤演算法

#include <stdio.h>
#include <windows.h>
#include "cv.h"
#include "cxcore.h"
#include "highgui.h"
#include <opencv2\opencv.hpp>
using namespace cv;

static const double pi = 3.14159265358979323846;
inline static double square(int a)
{
return a * a;
}
/*該函數目的：給img分配內存空間，並設定format，如位深以及channel數*/
inline static void allocateOnDemand(IplImage **img, CvSize size, int depth, int channels)
{
if (*img != NULL) return;
*img = cvCreateImage(size, depth, channels);
if (*img == NULL)
{
fprintf(stderr, "Error: Couldn't allocate image. Out of memory?\n");
exit(-1);
}
}
/*主函數，原程序是讀取avi視頻文件，然後處理，我簡單改成從攝像頭直接讀取數據*/
int main(int argc, char *argv[])
{

//讀取攝像頭
VideoCapture cap(0);
//讀取視頻文件

//VideoCapture cap; cap.open("optical_flow_input.avi");
if (!cap.isOpened())
{
return -1;
}
Mat frame;

/*
bool stop = false;
while (!stop)
{
cap >> frame;
// cvtColor(frame, edges, CV_RGB2GRAY);
// GaussianBlur(edges, edges, Size(7, 7), 1.5, 1.5);
// Canny(edges, edges, 0, 30, 3);
// imshow("當前視頻", edges);
imshow("當前視頻", frame);
if (waitKey(30) >= 0)
stop = true;
}
*/

//CvCapture *input_video = cvCaptureFromFile( "optical_flow_input.avi" );
//cv::VideoCapture cap = *(cv::VideoCapture *) userdata;

//if (input_video == NULL)
// {
// fprintf(stderr, "Error: Can't open video device.\n");
// return -1;
// }

/*先讀取一幀，以便得到幀的屬性，如長、寬等*/
//cvQueryFrame(input_video);

/*讀取幀的屬性*/
CvSize frame_size;
frame_size.height = cap.get(CV_CAP_PROP_FRAME_HEIGHT);
frame_size.width = cap.get(CV_CAP_PROP_FRAME_WIDTH);

/*********************************************************/

/*用於把結果寫到文件中去,非必要
int frameW = frame_size.height; // 744 for firewire cameras
int frameH = frame_size.width; // 480 for firewire cameras
VideoWriter writer("VideoTest.avi", -1, 25.0, cvSize(frameW, frameH), true);

/*開始光流法*/
//VideoWriter writer("VideoTest.avi", CV_FOURCC('D', 'I', 'V', 'X'), 25.0, Size(640, 480), true);

while (true)
{
static IplImage *frame = NULL, *frame1 = NULL, *frame1_1C = NULL,
*frame2_1C = NULL, *eig_image = NULL, *temp_image = NULL,
*pyramid1 = NULL, *pyramid2 = NULL;

Mat framet;
/*獲取第一幀*/
// cap >> framet;
cap.read(framet);
Mat edges;
//黑白抽象濾鏡模式
// cvtColor(framet, edges, CV_RGB2GRAY);
// GaussianBlur(edges, edges, Size(7, 7), 1.5, 1.5);
// Canny(edges, edges, 0, 30, 3);

//轉換mat格式到lpiimage格式
frame = &IplImage(framet);
if (frame == NULL)
{
fprintf(stderr, "Error: Hmm. The end came sooner than we thought.\n");
return -1;
}

/*由於opencv的光流函數處理的是8位的灰度圖，所以需要創建一個同樣格式的
IplImage的對象*/
allocateOnDemand(&frame1_1C, frame_size, IPL_DEPTH_8U, 1);

/* 把攝像頭圖像格式轉換成OpenCV慣常處理的圖像格式*/
cvConvertImage(frame, frame1_1C, 0);

/* 我們需要把具有全部顏色信息的原幀保存，以備最後在屏幕上顯示用*/
allocateOnDemand(&frame1, frame_size, IPL_DEPTH_8U, 3);
cvConvertImage(frame, frame1, 0);

/* 獲取第二幀 */
//cap >> framet;
cap.read(framet);
// cvtColor(framet, edges, CV_RGB2GRAY);
// GaussianBlur(edges, edges, Size(7, 7), 1.5, 1.5);
// Canny(edges, edges, 0, 30, 3);
frame = &IplImage(framet);
if (frame == NULL)
{
fprintf(stderr, "Error: Hmm. The end came sooner than we thought.\n");
return -1;
}

/*原理同上*/
allocateOnDemand(&frame2_1C, frame_size, IPL_DEPTH_8U, 1);
cvConvertImage(frame, frame2_1C, 0);

/*********************************************************
開始shi-Tomasi演算法，該演算法主要用於feature selection,即一張圖中哪些是我
們感興趣需要跟蹤的點(interest point)
input:
* "frame1_1C" 輸入圖像.
* "eig_image" and "temp_image" 只是給該演算法提供可操作的內存區域.
* 第一個".01" 規定了特徵值的最小質量，因為該演算法要得到好的特徵點，哪就
需要一個選擇的閾值
* 第二個".01" 規定了像素之間最小的距離，用於減少運算復雜度，當然也一定
程度降低了跟蹤精度
* "NULL" 意味著處理整張圖片，當然你也可以指定一塊區域
output:
* "frame1_features" 將會包含fram1的特徵值
* "number_of_features" 將在該函數中自動填充上所找到特徵值的真實數目,
該值<= 400
**********************************************************/

/*開始准備該演算法需要的輸入*/

/* 給eig_image,temp_image分配空間*/
allocateOnDemand(&eig_image, frame_size, IPL_DEPTH_32F, 1);
allocateOnDemand(&temp_image, frame_size, IPL_DEPTH_32F, 1);

/* 定義存放frame1特徵值的數組，400隻是定義一個上限 */
CvPoint2D32f frame1_features[400];
int number_of_features = 400;

/*開始跑shi-tomasi函數*/
cvGoodFeaturesToTrack(frame1_1C, eig_image, temp_image,
frame1_features, &number_of_features, .01, .01, NULL);

/**********************************************************
開始金字塔Lucas Kanade光流法，該演算法主要用於feature tracking,即是算出
光流，並跟蹤目標。
input:
* "frame1_1C" 輸入圖像，即8位灰色的第一幀
* "frame2_1C" 第二幀，我們要在其上找出第一幀我們發現的特徵點在第二幀
的什麼位置
* "pyramid1" and "pyramid2" 是提供給該演算法可操作的內存區域，計算中間
數據
* "frame1_features" 由shi-tomasi演算法得到的第一幀的特徵點.
* "number_of_features" 第一幀特徵點的數目
* "optical_flow_termination_criteria" 該演算法中迭代終止的判別，這里是
epsilon<0.3，epsilon是兩幀中對應特徵窗口的光度之差的平方，這個以後的文
章會講
* "0" 這個我不知道啥意思，反正改成1就出不來光流了，就用作者原話解釋把
means disable enhancements. (For example, the second array isn't
pre-initialized with guesses.)
output:
* "frame2_features" 根據第一幀的特徵點，在第二幀上所找到的對應點
* "optical_flow_window" lucas-kanade光流演算法的運算窗口,具體lucas-kanade
會在下一篇詳述
* "5" 指示最大的金字塔層數，0表示只有一層，那就是沒用金字塔演算法
* "optical_flow_found_feature" 用於指示在第二幀中是否找到對應特徵值，
若找到，其值為非零
* "optical_flow_feature_error" 用於存放光流誤差
**********************************************************/

/*開始為pyramid lucas kanade光流演算法輸入做准備*/
CvPoint2D32f frame2_features[400];

/* 該數組相應位置的值為非零，如果frame1中的特徵值在frame2中找到 */
char optical_flow_found_feature[400];

/* 數組第i個元素表對應點光流誤差*/
float optical_flow_feature_error[400];

/*lucas-kanade光流法運算窗口,這里取3*3的窗口,可以嘗試下5*5,區別就是5*5
出現aperture problem的幾率較小,3*3運算量小，對於feature selection即shi-tomasi演算法來說足夠了*/
CvSize optical_flow_window = cvSize(5, 5);
// CvSize optical_flow_window = cvSize(5, 5);
/* 終止規則，當完成20次迭代或者當epsilon<=0.3，迭代終止，可以嘗試下別的值*/
CvTermCriteria optical_flow_termination_criteria= cvTermCriteria(CV_TERMCRIT_ITER | CV_TERMCRIT_EPS, 20, .3);

/*分配工作區域*/
allocateOnDemand(&pyramid1, frame_size, IPL_DEPTH_8U, 1);
allocateOnDemand(&pyramid2, frame_size, IPL_DEPTH_8U, 1);

/*開始跑該演算法*/
cvCalcOpticalFlowPyrLK(frame1_1C, frame2_1C, pyramid1, pyramid2,frame1_features, frame2_features, number_of_features,
optical_flow_window, 5, optical_flow_found_feature,optical_flow_feature_error, optical_flow_termination_criteria, 0);

/*畫光流場，畫圖是依據兩幀對應的特徵值，
這個特徵值就是圖像上我們感興趣的點，如邊緣上的點P(x,y)*/
for (int i = 0; i< number_of_features; i++)
{
/* 如果沒找到對應特徵點 */
if (optical_flow_found_feature[i] == 0)
continue;
int line_thickness;
line_thickness = 1;

/* CV_RGB(red, green, blue) is the red, green, and blue components
* of the color you want, each out of 255.
*/
CvScalar line_color;
line_color = CV_RGB(255, 0, 0);

/*畫箭頭,因為幀間的運動很小，所以需要縮放，不然看不見箭頭，縮放因子為3*/
CvPoint p, q;
p.x = (int)frame1_features[i].x;
p.y = (int)frame1_features[i].y;
q.x = (int)frame2_features[i].x;
q.y = (int)frame2_features[i].y;

double angle;
angle = atan2((double)p.y - q.y, (double)p.x - q.x);
double hypotenuse;
hypotenuse = sqrt(square(p.y - q.y) + square(p.x - q.x));

/*執行縮放*/
q.x = (int)(p.x - 5 * hypotenuse * cos(angle));
q.y = (int)(p.y - 5 * hypotenuse * sin(angle));

/*畫箭頭主線*/
/* "frame1"要在frame1上作畫.
* "p" 線的開始點.
* "q" 線的終止點.
* "CV_AA" 反鋸齒.
* "0" 沒有小數位.
*/
cvLine(frame1, p, q, line_color, line_thickness, CV_AA, 0);

/* 畫箭的頭部*/
p.x = (int)(q.x + 9 * cos(angle + pi / 4));
p.y = (int)(q.y + 9 * sin(angle + pi / 4));
cvLine(frame1, p, q, line_color, line_thickness, CV_AA, 0);
p.x = (int)(q.x + 9 * cos(angle - pi / 4));
p.y = (int)(q.y + 9 * sin(angle - pi / 4));
cvLine(frame1, p, q, line_color, line_thickness, CV_AA, 0);
}
/*顯示圖像*/

/*創建一個名為optical flow的窗口，大小自動改變*/
cvNamedWindow("Optical Flow", CV_WINDOW_NORMAL);
cvFlip(frame1, NULL, 2);
cvShowImage("Optical Flow", frame1);

/*延時，要不放不了*/
cvWaitKey(33);

/*寫入到文件中去*/

// cv::Mat m = cv::cvarrToMat(frame1);//轉換lpimgae到mat格式
// writer << m;//opencv3.0 version writer

}
cap.release();
cvWaitKey(33);
system("pause");
}

② 運動目標檢測——光流法與opencv代碼實現

運動目標的檢測的其主要目的是 獲取目標對象的運動參數（位置、速度、加速度等）及運動軌跡 ，通過進一步分析處理，實現對目標行為更高層級上的理解。
運動目標檢測技術目的是 從序列圖像中將變化區域從背景圖像中提取出來 ，常用於視頻監視、圖像壓縮、三維重構、異常檢測等。

運動目標檢測主流方法有幀差法、背景差法、光流法等。光流法源於 仿生學 思想，更貼近於直覺，大量昆蟲的視覺機理便是基於光流法。
二十世紀五十年代心理學家Gibson在他的著作「The Perception of Visual World」中首次提出了以心理學實驗為基礎的光流法基本概念，而直到八十年代才由Horn、Kanade、Lucash和Schunck創造性地將灰度與二維速度場相聯系，引入光流約束方程的演算法，對光流計算做了奠基性的工作。

光流（optical flow）：由於目標對象或者攝像機的移動造成的圖像對象在連續兩幀圖像中的移動。

通俗說，對於一個圖片序列，把每張圖像每個像素在連續幀之間的運動速度和方向（ 某像素點在連續兩幀上的位移矢量 ）找出來就是光流場。

第t幀的時A點的位置是(x1, y1)，第t+1幀時A點位置是(x2,y2)，則像素點A的位移矢量：(ux, vy) = (x2, y2) - (x1,y1)

如何知道第t+1幀的時候A點的位置涉及到不同的光流計算方法，主要有四種：基於梯度的方法、基於匹配的方法、基於能量的方法、基於相位的方法。

光流法依賴於三個假設：

根據所形成的光流場中 二維矢量的疏密程度 ，光流法可分為稠密光流與稀疏光流。

稀疏光流只對有 明顯特徵的組點 （如角點）進行跟蹤，計算開銷小。

http://www.opencv.org.cn/opencvdoc/2.3.2/html/moles/video/doc/motion_analysis_and_object_tracking.html#calcopticalflowfarneback
（1）calcOpticalFlowPyrLK
基於金字塔LK光流演算法，計算某些點集的稀疏光流。
參考論文《Pyramidal Implementation of the Lucas Kanade Feature TrackerDescription of the algorithm》
（2）calcOpticalFlowFarneback
基於Gunnar Farneback 的演算法計算稠密光流。
參考論文《Two-Frame Motion Estimation Based on PolynomialExpansion》
（3）CalcOpticalFlowBM
通過塊匹配的方法來計算光流
（4）CalcOpticalFlowHS
基於Horn-Schunck 的演算法計算稠密光流。
參考論文《Determining Optical Flow》
（5）calcOpticalFlowSF
論文《SimpleFlow: A Non-iterative, Sublinear Optical FlowAlgo》的實現

③ 134+135+136+137+136+135+134怎麼用金字塔演算法

金字塔演算法是一種簡單的加法運算方法，可以將數字從下往上排列成金字塔形狀，從而便於悔族計算。對於這道題目，可以按照以下步驟進行金字塔演算法的計算：
Step 1：將題目中的數字從下往上排列成金字塔形狀，如下所示：
134 135 136 137 136 135 134
Step 2：從倒數第二行開始，每相鄰的兩個數字進行相加，並將結果寫在它們上面的位置上，如下所示：
134
135 136
137 136 135
134
134
135 136
273 271 270
Step 3：重復Step 2，直攜液到計算碧隱弊到金字塔的頂端，最終結果即為頂端數字的值，如下所示：
134
135 136
273 271 270
679 542 406
因此，134+135+136+137+136+135+134的結果為679。