lucaskanade演算法

Ⅰ 光流（Optical Flow）

光流是由觀察者和場景之間的[相對運動]引起的視覺場景中物體、表面和邊緣的運動模式。一般而言，光流是由於場景中前景目標本身的移動、觀測者運動，或者兩者的共同運動所產生的。

光流在很多領域中都被用到，例如視頻中的運動目標檢測，視頻壓縮等等。

在分析光流時，需要用到兩個重要假設： 1.對象的像素強度在連續幀之間不會改變。2.相鄰像素具有相似的運動。 下面我們運用這兩個假設來推導光流公式。

光流法實際是通過檢測圖像像素點的強度隨時間的變化進而推斷出物體移動速度及方向的方法。假設該移動很小，假設該移動很小，那麼可以根據泰勒級數得出：

根據假設1，由於移動很小，此時這時對象位移不會改變對象的像素強度:

進一步得出:

最終得出光流公式：

這里的 , 是x, y方向上的速率，或稱為I(x, y, t)的光流。而 , , 則是圖像(x, y, t)在對應方向上的偏導數。可以看到光流公式中有兩個未知量 ， , 無法直接求解。因此我們需要更多的外部條件來求解方程。

Lucas-Kanade 方法

求解光流方程有很多方法，其中最著名的便是 Lucas-Kanade方法。它應用了之前提到的第二個假設，即所有相鄰像素都將具有相似的運動。對於每一個像素，Lucas-Kanade 方法選取與它相鄰的8個像素進行分析。根據假設，所有 9 個像素都有相同的運動。所以現在我們的問題求解 只有2個未知變數的9個方程組。

......

這樣的方程組沒有唯一解，這里我們使用最 最小二乘擬合 方法獲得一個最優近似解。

Lucas-Kanade方法的局限性

由於上述假設和分析都是針對較小的運動，光流演算法會受到突然移動的影響。如果兩幀圖像之間的運動太大，光流法可能會失效。對於這個問題我們可以通過圖像金字塔來解決，當我們使用更高層的金字塔圖像時，小的運動被去除，大的運動變成小運動。然後再應用 Lucas-Kanade，就可以得到正確的光流。

OpenCV 光流演示：https://docs.opencv.org/3.4/d4/dee/tutorial_optical_flow.html 

Ⅱ 請高人幫忙解讀一下下面這段MATLAB程序 主要是一些矩陣運算 但我看不懂 謝謝啦

還要嗎？100分挺誘惑我的哦～～

Ⅲ 光流計演算法Lucas_Kanade.m中的matlab函數的輸出矩陣是什麼怎麼理解而且輸入的第三個參數density是什麼

原函數比較長，我這里就不列了。根據函數頭部。
function [Dx, Dy] =Lucas_Kanade(img1,img2,density)
輸出的矩陣就是X（橫軸）方向和Y（縱軸）方向的光流場，你可以通過:
>>quiver(Dx,Dy)
查看光流場的圖形。
density表示輸出光流場的粒度（即每個點都輸出，那麼太多了，因此稀疏取一下輸出就可以了），例如：

a =
1 13 15 28
2 17 12 48
9 1 3 99

>> density=2;
>>a(1:density:size(a,1),1:density:size(a,2))
ans =
1 15
9 3
原來a有很多點，那麼只要間隔輸出一些代表點就可以了。

Ⅳ 光流(Optical flow)-視頻分析基礎概念

光流 是空間運動物體在 觀察成像平面 上的像素運動的 瞬時速度 ，是利用圖像序列中像素在時間域上的變化以及相鄰幀之間的 相關性 來找到上一幀跟當前幀之間存在的對應關系，從而計算出相鄰幀之間物體的運動信息的一種方法。一般而言，光流是由於場景中前景目標本身的移動、相機的運動，或者兩者的共同運動所產生的。考慮下圖(圖片來自 維基網路 )：

圖中表示一個小球在連續5幀圖像中的移動，箭頭則表示小球的 位移矢量 。簡單來說，光流是空間運動物體在觀測成像平面上的像素運動的「瞬時速度」，光流的研究是利用圖像序列中的像素強度數據的時域變化和相關性來確定各自像素位置的「運動」，研究光流場的目的就是為了從圖片序列中近似得到不能直接得到的現實中的 運動場 。光流應用於諸多領域：

在介紹光流的計演算法之前有必要了解：光流之所以生效是依賴於這幾個假設：

其實也很好理解，如果不滿足以上條件，那麼也找不到該像素在下一幀的位置，自然也無法計算出它的運動。

假設第一幀圖像中的像素 I(x, y, t) 在時間 dt 後移動到第二幀圖像的 (x+dx, y+dy) 處。根據上述第一條假設：灰度值不變，我們可以得到：

對等號右側進行泰勒級數展開，消去相同項，兩邊都除以 dt ，得到如下方程：

其中：

上面的等式叫做光流方程。其中 f x 和 f y 的梯度，同樣 f t 是時間方向的梯度。但 (u, v) 是不知道的。我們不能在一個等式中求解兩個未知數。有幾個方法可以幫我們解決這個問題，其中的一個是 Lucas-Kanade 法 。

這里就要用到上面提到的第二個假設條件，領域內的所有像素點具有相同的運動。Lucas-Kanade法就是利用一個3x3的領域中的9個像素點具有相同的運動，就可以得到9個點的光流方程(即上述公式)，用這些方程來求得 (u, v) 這兩個未知數，顯然這是個約束條件過多的方程組，不能解得精確解，一個好的解決方法就是使用最小二乘來擬合。求解過程：

這樣我們跟蹤一些點就能得到這些點的光流向量，但是這里還存在 尺度空間 的問題，簡單來講，直到現在我們還只是處理一些很小的運動，如果是大的運動那該怎麼辦？ 圖像金字塔 。在圖像金字塔頂層，小的運動被移除，大的運動轉換成了小的運動，這樣就能跟蹤到了原本大的運動，重復計算圖像金字塔不同層的圖像的光流，我們就得到了在不同尺度空間上的光流。

Ⅳ 對三種圖像配准方法的說明

圖像配准，英文稱為image alignment。本文將分別對四種圖像配準的方法進行說明，即前向累加法（forward additive）、前向合成法（forward compositional）和逆向合成法（invers compositional）。

forward additive method又稱為Lucas-Kanade algorithm，它的目標是將一個模板圖像T配准（align）到輸入圖像I上，T表示圖像上的提取出來的一個小patch，它的目標函數如圖1-1所示：

forward compositional的迭代方式圖2-1所示：

泰勒展開線性化求解過程如圖2-2所示，W(x;0)相當於沒有對點進行變化，因此W(W(x;0);p) = W(x;p)：

inverse compositional方法將模板T和輸入圖像I的角色做了反轉，迭代方式如圖3-1所示：

泰勒展開線性化求解非線性最小二乘問題如圖3-2所示，注意到其中的Hessian矩陣與之前相比有所不同，用模板T的梯度代替了原來輸入圖像I的梯度：

由於inverse compositional的Hessian矩陣和待求參數無關，所以Hessian矩陣可以預先計算出來，而Hessian矩陣是整個演算法中最耗時的部分，不用在每一次迭代過程都計算Hessian矩陣就大大提高了演算法的效率。

以上就是對三種圖像配准方法的全部說明。

Ⅵ 非線性最小二乘法

一.梯度下降法以及Jacobian矩陣計算

在2010年的關於L-K和AAM的博客里提到，模板匹配公式的一階泰勒展開ΔT=J*Δp，J是用於梯度下降的Jacobian矩陣，是高維矢量函數值T=f(p)相對與參數矢量p變化時的增量(導數)。如果p是n維矢量，T是M維矢量，則J是一個[m*n]的矩陣。J在(i,j)處的元素值是(əTi/əpj)。

Ⅶ 如何使用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");
}