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");
}