LK金字塔光流法与简单实现
LK金字塔光流法与简单实现
- 闲谈时刻
- 介绍
- Lucas–Kanade光流算法
- L-K 金字塔光流算法
- 算法原理
- 建立金字塔
- 金字塔迭代
- 迭代过程
- 算法流程
- 算法实现
- 总结
- 参考资料
闲谈时刻
不务正业预警
眼看着一个学期又告一段落，几个月来拢共还是没写几篇博客。不过手头上倒是还积累着不少资料值得一写，趁着新春得闲可以好好梳理梳理了。
介绍
开场依照惯例是得简单介绍什么是光流，来自 Wiki：
光流(Optical flow or optic flow)是关于视域中的物体运动检测中的概念。用来描述相对于观察者的运动所造成的观测目标、表面或边缘的运动。
简单来讲，光流描述了场景中物体运动在视觉中的变化。光流的概念由Gibson在1950年提出，其通过相邻帧之间像素点的对应关系计算像素点的瞬时速度，从而描述物体信息。
为了应用光流计算物体运动法，相邻帧需要满足两个假设：
- 亮度不变：两个间隔帧之间的像素点的亮度需要保持恒定，从而对两个相邻帧之间的像素点进行对应；
- 运动较小：需要保证像素点并不会随着时间的变化而剧烈变化，从而通过相邻帧之间对应点位置变化引起的灰度值变化来近似为灰度对位置的偏导。
以二维的图像为例，以 I(x, y, t) 描述图像在 t 时刻的灰度值。其灰度值在两个间隔帧之间的变化值为 Δx，Δy，Δt。由于假设1，根据间隔帧的灰度值恒定，可以得出公式1：
（1）
由于假设2，可以通过泰勒级数对该方程进行展开，得到公式2：
（2）
H.O.T 为高阶小量，其在移动足够小的时候可以忽略，从而得到公式3：
（3）
公式3对左右两边除以 Δt 可以得到公式4：
（4）
也即是公式5：
（5）
其中，Vx，Vy 是像素点在 x 与 y 方向上的速度，也即是 I(x, y, t) 的光流 。用 Ix，Iy 和 It 分别表示像素点对于 x，y，t 的偏导后，可以得到最终的光流表达式：
（6）
抑或是：
（7）
了解了光流法的基本要素，那么应当如何求解光流呢？
Lucas–Kanade光流算法
首先介绍Lucas–Kanade光流算法（L-K算法），其本质是通过最小二乘法以不需要迭代的方法求解光流，是光流算法中最简单的一种。
为了简化光流计算的流程，L-K光流算法为其增加了一个额外的假设：
- 空间一致：某个像素点领域内保持相同的瞬时速度。
对于给定大小的窗口，可以列出公式：
…
其中，qn 代表窗口中的像素点，Ix(qn)是像素qn在x方向上的偏导。写成矩阵形式也就是公式8：
Av = b， 其中 （8）
该方程组为超定方程，等式个数多于未知数的个数，因此通过最小二乘法计算出近似解，如公式9所示。
（9）
也即是公式10：
（10）
L-K 金字塔光流算法
虽然光流法假设2本身并不容易满足，在实际应用中容易遇到间隔帧出现较大变化的情况，从而对计算结果造成较大的误差。为了减缓这一假设不成立情况所导致的问题，可以在计算中缩小图像的尺寸，从而使得像素点的运动减少（图像尺寸减半时，像素运动自然也同时减半）。
通过对图像尺寸进行缩放，建立图像金字塔，可以使得L-K金字塔算法应用于较大的运动。
算法原理
L-K金字塔算法的主要的步骤可以分为三步：建立金字塔，金字塔跟踪以及迭代过程。
建立金字塔
首先需要对原始图像建立金字塔，其中原始图像位于底层，其上每一层均基于上一层进行计算。金字塔中每一层均是上一层的下采样，其计算公式为：
（11）
其中 IL(x, y) 为 L 层图像中 x, y 位置所在像素点的灰度值，其相当于通过[0.25 0.5 0.25]的低通滤波器进行迭代计算。
金字塔迭代
金字塔迭代过程为：
- 建立金字塔后，首先对于最高层的图像计算出其上的光流；
- 通过上一层（L+1层）的计算结果对下一层（L层）的的图像进预平移，并在L层的基础上计算出该层的残余光流向量dL。由于金字塔中每一层的尺寸均是上一层的一半，因此其每一层的光流均是其上一层的二分之一。
通过L+1层计算出的光流作为初值计算L层的光流，可以保证每一层的残余光流向量较小，从而应用L-K光流算法； - 对于每一层迭代光流计算，最终得到的光流也即是所有层光流的叠加。
由于顶层图像尺寸较小，其初始的光流估计量可以设置为0，即 gL = [0, 0] 。
迭代过程
介绍了对于整个金字塔迭代过程，还需要提供对于每一层的残余光流计算方法。
在 L 层上计算残余光流，首先需要对于光流进行更新，即 v = 2 * v，并通过其某个像素点的偏导值计算出其空间梯度矩阵，如公式12所示：
（12）
其中 Ix 为上x方向的偏导，Iy 为y方向上的偏导。
迭代计算间隔帧之间对应点的灰度误差，如公式13所示：
（13）
并计算对应的误差矩阵，如公式14所示：
（14）
通过L-K算法可以计算出当前的仿射光流，并对当前层的光流进行更新。当该仿射光流的模小于阈值时，可以认为迭代计算结果已经收敛，并结束迭代过程。
算法流程
上述的L-K 金字塔光流算法原理较为清晰，具体的计算结果可以参照论文中所提供的伪代码：
算法实现
在OpenCV中对L-K金字塔光流算法提供了实现，其API为calcOpticalFlowPyrLK。这里给出一个简单的 C++ 版本的实现。
首先是头文件，lk_tracker.h：
#ifndef LK_TRACKER
#define LK_TRACKER#include <iostream>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/core/core.hpp>using namespace std;
using namespace cv;void trace(string out) {cout << out << endl;
};struct TrackedPoint {Point point;Mat opticalFlow;TrackedPoint() {};TrackedPoint(const Point p, const Mat of):point(p), opticalFlow(of) {};
};class LKTracker {private:int maxPyramidLayer = 3;int wx = 2;int wy = 2;int maxIteration = 50;double opticalflowResidual = 0.0001;private:static int getMatInt(Mat mat, int row, int col);static double getMatDouble(Mat mat, int row, int col);static double getMatDouble(Mat mat, double row, double col);void lowpassFilter(InputArray src, OutputArray dst);Mat calcGradientMatrix(InputArray frame, Point2f p);Mat calcMismatchVector(InputArray preFrame, InputArray curFrame, Point2f p, Mat g, Mat v);vector<Mat> buildPyramid(InputArray src);double calcResidual(Mat mat);bool isOpticalFlowValid(Mat of);double calcHarrisResponse(Mat gradient, double alpha);public:LKTracker();~LKTracker();vector<TrackedPoint> track(InputArray preFrame, InputArray curFrame, vector<Point> keyPoints);vector<TrackedPoint> trackAll(InputArray preFrame, InputArray curFrame, double qualityLevel);
};#endif
实现部分 lk_tracker.cc：
#include "lk_tracker.h"#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/core/core.hpp>using namespace cv;LKTracker::LKTracker() {}LKTracker::~LKTracker() {}int LKTracker::getMatInt(Mat mat, int row, int col) {Size size = mat.size();if (col < 0 || col >= size.width || row < 0 || row >= size.height) {return 0;}return mat.at<uchar>(row, col);
}double LKTracker::getMatDouble(Mat mat, int row, int col) {Size size = mat.size();if (col < 0 || col >= size.width || row < 0 || row >= size.height) {return 0;}if (mat.type() == CV_64F) {return mat.at<double>(row, col);} else {return (double)mat.at<uchar>(row, col);}
}double LKTracker::getMatDouble(Mat mat, double row, double col) {Size size = mat.size();if (col < 0 || col >= size.width || row < 0 || row >= size.height) {return 0;}int floorRow = floor(row);int floorCol = floor(col);double fracRow = row - floorRow;double fracCol = col - floorCol;int ceilRow = floorRow + 1;int ceilCol = floorCol + 1;return ((1.0 - fracRow) * (1.0 - fracCol) * getMatDouble(mat, floorRow, floorCol) + (fracRow * (1.0 - fracCol) * getMatDouble(mat, ceilRow, floorCol)) + ((1.0 - fracRow) * fracCol * getMatDouble(mat, floorRow, ceilCol)) +(fracRow * fracCol * getMatDouble(mat, ceilRow, ceilCol)));
}void LKTracker::lowpassFilter(InputArray src, OutputArray dst) {Mat srcMat = src.getMat();int dstWidth = srcMat.size().width / 2;int dstHeight = srcMat.size().height / 2;// dst.create(Size(dstWidth, dstHeight), srcMat.type());dst.create(Size(dstWidth, dstHeight), CV_64F);Mat dstMat = dst.getMat();for (int x = 0; x < dstHeight; x++) {for (int y = 0; y < dstWidth; y++) {double val = 0;val += getMatInt(srcMat, 2*x, 2*y) * 0.25;val += getMatInt(srcMat, 2*x-1, 2*y) * 0.125;val += getMatInt(srcMat, 2*x+1, 2*y) * 0.125;val += getMatInt(srcMat, 2*x, 2*y-1) * 0.125;val += getMatInt(srcMat, 2*x, 2*y+1) * 0.125;val += getMatInt(srcMat, 2*x-1, 2*y-1) * 0.0625;val += getMatInt(srcMat, 2*x+1, 2*y-1) * 0.0625;val += getMatInt(srcMat, 2*x-1, 2*y+1) * 0.0625;val += getMatInt(srcMat, 2*x+1, 2*y+1) * 0.0625;// dstMat.at<uchar>(x, y) = (uchar)val;dstMat.at<double>(x, y) = val;}}
}vector<Mat> LKTracker::buildPyramid(InputArray src) {vector<Mat> layers;Mat currLayer = src.getMat();layers.push_back(currLayer);for (int i = 0; i < this->maxPyramidLayer; i++) {Mat nextLayer;this->lowpassFilter(currLayer, nextLayer);layers.push_back(nextLayer);currLayer = nextLayer;}return layers;
}Mat LKTracker::calcGradientMatrix(InputArray frame, Point2f p) {Mat frameMat = frame.getMat();Mat mat(2, 2, CV_64F, Scalar(0));int pRow = p.y, pCol = p.x;for (int x = pRow - wx; x <= pRow + wx; x++) {for (int y = pCol - wy; y <= pCol + wy; y++) {double ix = (getMatDouble(frameMat, x + 1, y) - getMatDouble(frameMat, x - 1, y)) / 2.0;double iy = (getMatDouble(frameMat, x, y + 1) - getMatDouble(frameMat, x, y - 1)) / 2.0;mat.at<double>(0, 0) += ix * ix;mat.at<double>(0, 1) += ix * iy;mat.at<double>(1, 0) += ix * iy;mat.at<double>(1, 1) += iy * iy;}}return mat;
}Mat LKTracker::calcMismatchVector(InputArray preFrame, InputArray curFrame, Point2f p, Mat g, Mat v) {Mat preMat = preFrame.getMat();Mat curMat = curFrame.getMat();Mat mismatchVector(2, 1, CV_64F, Scalar(0));int pRow = p.y, pCol = p.x;for (int x = pRow - wx; x <= pRow + wx; x++) {for (int y = pCol - wy; y <= pCol + wy; y++) {double curX = 1.0 * x + g.at<double>(0, 0) + v.at<double>(0, 0);double curY = 1.0 * y + g.at<double>(1, 0) + v.at<double>(1, 0);double ix = (getMatDouble(preMat, x + 1, y) - getMatDouble(preMat, x - 1, y)) / 2.0;double iy = (getMatDouble(preMat, x, y + 1) - getMatDouble(preMat, x, y - 1)) / 2.0;double diff = getMatDouble(preMat, x, y) - getMatDouble(curMat, curX, curY);mismatchVector.at<double>(0, 0) += diff * ix;mismatchVector.at<double>(1, 0) += diff * iy;}}return mismatchVector;
}double LKTracker::calcResidual(Mat mat) {return sqrt(pow(mat.at<double>(0, 0), 2) + pow(mat.at<double>(1, 0), 2));
}bool LKTracker::isOpticalFlowValid(Mat of) {return (abs(of.at<double>(0, 0)) <= 50 && abs(of.at<double>(0, 0)) >= 1 && abs(of.at<double>(1, 0)) <= 50 && abs(of.at<double>(1, 0)) >= 1);
}vector<TrackedPoint> LKTracker::track(InputArray preFrame, InputArray curFrame, vector<Point> keyPoints) {vector<Mat> prePyramid = this->buildPyramid(preFrame);vector<Mat> curPyramid = this->buildPyramid(curFrame);vector<TrackedPoint> tPoints;for (unsigned int i = 0; i < keyPoints.size(); i++) {Mat g(2, 1, CV_64F, Scalar(0));Point srcPoint = keyPoints[i];bool isValid = true;for (int j = this->maxPyramidLayer - 1; j > 0; j--) {Mat preMat = prePyramid[j];Mat curMat = curPyramid[j];Point2f prePoint;prePoint.x = 1.0 * srcPoint.x / pow(2.0, j);prePoint.y = 1.0 * srcPoint.y / pow(2.0, j);Mat gradient = calcGradientMatrix(preMat, prePoint);Mat gradientInv = gradient.inv();Mat v(2, 1, CV_64F, Scalar(0));int iterCount = 0;double residual = 1;while(iterCount < this->maxIteration && residual > this->opticalflowResidual) {iterCount++;Mat mismatch = calcMismatchVector(preMat, curMat, prePoint, g, v);Mat delta = gradientInv * mismatch;v += delta;residual = calcResidual(delta);}if (iterCount >= this->maxIteration) {isValid = false;break;}if (j == 0) {g = g + v;} else {g = 2 * (g + v);}}if (isValid && isOpticalFlowValid(g)) {Point dstPoint;dstPoint.x = (int)(srcPoint.x + g.at<double>(1, 0));dstPoint.y = (int)(srcPoint.y + g.at<double>(0, 0));TrackedPoint tPoint(dstPoint, g);tPoints.push_back(tPoint);}}return tPoints;
}double LKTracker::calcHarrisResponse(Mat gradient, double alpha) {double A = gradient.at<double>(0, 0);double B = gradient.at<double>(0, 1);double C = gradient.at<double>(1, 1);double det = A * C - B * B;double trace = A + C;return det - trace * trace * alpha;
}vector<TrackedPoint> LKTracker::trackAll(InputArray preFrame, InputArray curFrame, double qualityLevel) {Mat preMat = preFrame.getMat();vector<Point> keyPoints;goodFeaturesToTrack(preFrame, keyPoints, 1000, 0.01, 10, Mat());return track(preFrame, curFrame, keyPoints);
}
总结
好久好久没写 C++ 了，突然布置作业拿出 C++ 开始干活之后发现自己连指针的概念都已经拎不清，结果只能全用 vector 一把梭。实现部分代码丑陋就不说了，性能也非常低下。等什么时候得空了就把 JS 先放下，好好复习 C++ 罢。
参考资料
- 光流Optical Flow介绍与OpenCV实现
- wiki 光流法
- Lucas–Kanade光流算法
- 光流（三）–LK算法改进（金字塔LK）
- 《Pyramidal Implementation of the Lucas Kanade Feature TrackerDescription of the algorithm》
LK金字塔光流法与简单实现相关推荐
- 【OpenCV 笔记】金字塔光流法追踪运动目标
原理部分见https://blog.csdn.net/zouxy09/article/details/8683859 环境:OpenCV3.3+VS2017 光流简言之就是像素在某时刻运动的瞬时速度, ...
- 基于金字塔LK的光流法实现—根据论文自己实现的c++代码
自己实现代码与opencv的接口函数对比 在VINS-Mono上效果对比 在两张图片上效果对比 PyrLKTracking()源码 总结 最近自己在研究VINS-Mono源码的特征提取与追踪部分的原理 ...
- 图像金字塔LK光流法原理分析
图像金字塔LK光流法原理分析 1.LK光流法原理分析 2.基于图像金字塔的LK光流法原理分析 本篇博客只讲述原理,c++代码实现请参考博客< 基于金字塔LK的光流法实现-根据论文自己实现的c++ ...
- 传统目标跟踪——光流法
目录 一.光流法 二.LK光流法 2.1 实现原理 2.2 API 三.代码 四.总结 基于特征点的光流跟踪,在目标上提取一些特征点,然后在下一帧计算这些特征点的光流匹配点,统计得到目标的位置.在跟踪 ...
- 学习MSCKF笔记——前端、图像金字塔光流、Two Point Ransac
学习MSCKF笔记--前端.图像金字塔光流.Two Point Ransac 学习MSCKF笔记--前端.图像金字塔光流.Two Point Ransac 1. 图像金字塔光流 2. Two Poin ...
- OpenCV Using Python——基于SURF特征提取和金字塔LK光流法的单目视觉三维重建 (光流、场景流)...
https://blog.csdn.net/shadow_guo/article/details/44312691 基于SURF特征提取和金字塔LK光流法的单目视觉三维重建 1. 单目视觉三维重建问题 ...
- OpenCV：金字塔LK光流法
金字塔LK光流法的三个假设 亮度恒定,即图像场景中目标的像素在帧间运动时外观上保持不变: 时间连续或者运动是"小运动",即图像的运动随时间的变化比较缓慢: 空间一致,即一个场景中同 ...
- 详解LK光流法（含金字塔多层光流），反向光流法（附代码）
LK光流法可用来跟踪特征点的位置. 比如在img1中的特征点,由于相机或物体的运动,在img2中来到了不同的位置.后面会称img1为Template(T),img2为I. 光流法有个假设: 灰度不变假 ...
- 【老生谈算法】matlab实现金字塔LK光流法源码——金字塔LK光流法
基于金字塔LK光流法的MATLAB代码 1.原文下载: 本算法原文如下,有需要的朋友可以点击进行下载 序号 原文(点击下载) 本项目原文 [老生谈算法]基于金字塔LK光流法的MATLAB代码.docx ...
最新文章
- string.Format 方法拼入{}
- 深度学习时出现的一些安装问题+ubuntu apt的一些问题+github release文件加速
- php7 mysql json 小程序_微信小程序JSON数组递交PHP服务端解析处理
- comps电磁场模拟软件_|Mentor Graphics IE3D(电磁场仿真软件)下载v15.0官方版 - 欧普软件下载...
- python画名侦探柯南_基于flask的可视化动漫分析网站【python入门必学】
- Duilib使用wke显示echarts
- 快手、携程等公司转战到 ClickHouse，ES 难道不行了？
- 一个树莓派集群 （VAX）
- mac pro 系统升级带来的问题
- 二进制 八进制 十进制 十六进制 之间进制转换(图解篇)
- 数据库 mysql 删除一列数据
- NIT 股市风云 按位与运算 F. 休赛季的引援#2
- crm客户管理系统的功能有哪些？
- Python实战案例，PyQt5模块，实现疫情信息快速查看工具(附源码)
- 如何在uni-app 平台快速实现一对一音视频通话应用
- 小程序js赋值在wxhl中作用不到
- 父元素和子元素外边距重合
- 微盟遭员工“删库跑路”：SaaS服务暂停，或涉及300万商户
- android采用MVP漫画APP、适配刘海屏、小黄车主界面、录音波浪动画、综合APP等源码
- 2020 非常火的 11 个微前端框架