LK金字塔光流法与简单实现

LK金字塔光流法与简单实现LK金字塔光流法

大家好,又见面了,我是你们的朋友全栈君。

闲谈时刻

不务正业预警
眼看着一个学期又告一段落,几个月来拢共还是没写几篇博客。不过手头上倒是还积累着不少资料值得一写,趁着新春得闲可以好好梳理梳理了。

介绍

开场依照惯例是得简单介绍什么是光流,来自 Wiki

光流(Optical flow or optic flow)是关于视域中的物体运动检测中的概念。用来描述相对于观察者的运动所造成的观测目标、表面或边缘的运动。

简单来讲,光流描述了场景中物体运动在视觉中的变化。光流的概念由Gibson在1950年提出,其通过相邻帧之间像素点的对应关系计算像素点的瞬时速度,从而描述物体信息。

为了应用光流计算物体运动法,相邻帧需要满足两个假设:

  • 亮度不变:两个间隔帧之间的像素点的亮度需要保持恒定,从而对两个相邻帧之间的像素点进行对应;
  • 运动较小:需要保证像素点并不会随着时间的变化而剧烈变化,从而通过相邻帧之间对应点位置变化引起的灰度值变化来近似为灰度对位置的偏导。

以二维的图像为例,以 I(x, y, t) 描述图像在 t 时刻的灰度值。其灰度值在两个间隔帧之间的变化值为 Δx,Δy,Δt。由于假设1,根据间隔帧的灰度值恒定,可以得出公式1:
公式1 (1)
由于假设2,可以通过泰勒级数对该方程进行展开,得到公式2:
公式2 (2)
H.O.T 为高阶小量,其在移动足够小的时候可以忽略,从而得到公式3:
公式3 (3)
公式3对左右两边除以 Δt 可以得到公式4:
公式4 (4)
也即是公式5:
公式5 (5)
其中,Vx,Vy 是像素点在 x 与 y 方向上的速度,也即是 I(x, y, t) 的光流 。用 Ix,Iy 和 It 分别表示像素点对于 x,y,t 的偏导后,可以得到最终的光流表达式:
公式6 (6)
抑或是:
公式7 (7)
了解了光流法的基本要素,那么应当如何求解光流呢?

Lucas–Kanade光流算法

首先介绍Lucas–Kanade光流算法(L-K算法),其本质是通过最小二乘法以不需要迭代的方法求解光流,是光流算法中最简单的一种。
为了简化光流计算的流程,L-K光流算法为其增加了一个额外的假设:

  • 空间一致:某个像素点领域内保持相同的瞬时速度。

对于给定大小的窗口,可以列出公式:
LK金字塔光流法与简单实现
在这里插入图片描述

在这里插入图片描述
其中,qn 代表窗口中的像素点,Ix(qn)是像素qn在x方向上的偏导。写成矩阵形式也就是公式8:
Av = b, 其中 公式8 (8)
该方程组为超定方程,等式个数多于未知数的个数,因此通过最小二乘法计算出近似解,如公式9所示。
公式9 (9)
也即是公式10:
公式10 (10)

L-K 金字塔光流算法

虽然光流法假设2本身并不容易满足,在实际应用中容易遇到间隔帧出现较大变化的情况,从而对计算结果造成较大的误差。为了减缓这一假设不成立情况所导致的问题,可以在计算中缩小图像的尺寸,从而使得像素点的运动减少(图像尺寸减半时,像素运动自然也同时减半)。
通过对图像尺寸进行缩放,建立图像金字塔,可以使得L-K金字塔算法应用于较大的运动。

算法原理

L-K金字塔算法的主要的步骤可以分为三步:建立金字塔,金字塔跟踪以及迭代过程。

建立金字塔

首先需要对原始图像建立金字塔,其中原始图像位于底层,其上每一层均基于上一层进行计算。金字塔中每一层均是上一层的下采样,其计算公式为:
公式11 (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 (12)
其中 Ix 为上x方向的偏导,Iy 为y方向上的偏导。
迭代计算间隔帧之间对应点的灰度误差,如公式13所示:
公式13 (13)
并计算对应的误差矩阵,如公式14所示:
公式14 (14)
通过L-K算法可以计算出当前的仿射光流,并对当前层的光流进行更新。当该仿射光流的模小于阈值时,可以认为迭代计算结果已经收敛,并结束迭代过程。

算法流程

上述的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++ 罢。

参考资料

  1. 光流Optical Flow介绍与OpenCV实现
  2. wiki 光流法
  3. Lucas–Kanade光流算法
  4. 光流(三)–LK算法改进(金字塔LK)
  5. 《Pyramidal Implementation of the Lucas Kanade Feature TrackerDescription of the algorithm》
版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 举报,一经查实,本站将立刻删除。

发布者:全栈程序员-用户IM,转载请注明出处:https://javaforall.cn/163527.html原文链接:https://javaforall.cn

【正版授权,激活自己账号】: Jetbrains全家桶Ide使用,1年售后保障,每天仅需1毛

【官方授权 正版激活】: 官方授权 正版激活 支持Jetbrains家族下所有IDE 使用个人JB账号...

(0)


相关推荐

  • CC2541蓝牙学习——ADC

    CC2541蓝牙学习——ADCCC2541的ADC支持多达14位的模拟数字转换与高达12位的有效位数。它包括一个模拟多路转换器,具有多达8个各自可独立配置的通道,一个参考电压发生器。转换结果通过DMA写入存储器。还具有若干运行模式

  • php属于前端还是后端_php实时推送到前端

    php属于前端还是后端_php实时推送到前端功能说明使用第三方平台goeasy实现服务端向前端推送数据基本原理WebSocket使用准备申请goeasy账号并创建应用官网http://www.goeasy.io安装并开启goeasy插件(注意清除缓存)在插件配置中填写应用的Appkeys等配置项使用说明使用插件集成的事件插件在前台(index模块)和后台(admin模块)各集成了两个默认的事件订阅,可以在js中通过监听top来处理,例:也…

  • Python之json文件

    json简介json是一种轻量级的数据交换格式完全独立于编程语言的文本格式来存储和表示数据简单和清晰的层次结构使得json成为理想的数据交换语言。易于阅读和编写,易于机器解析和生成,并有效地提升

    2021年12月19日
  • UFT12.02 LICENSE延期

    UFT12.02 LICENSE延期1.以win7系统为例,安装完成后,修改图片目录下红色目录,仅第一次需修改2.删除红色目录,每次延期都需删除或修改此目录3.运行instdemo.exe,无报错,应该就延期成功了

  • tcp为什么是三次握手不是两次握手_tcp四次挥手

    tcp为什么是三次握手不是两次握手_tcp四次挥手一、为什么握手是三次,而不是两次或者四次?答:两次不安全,四次没必要。tcp通信需要确保双方都具有数据收发的能力,因此双方都要发送SYN确保对方具有通信的能力二、为什么挥手是四次而不是三次?答:发送FIN包只能表示对方不再发送数据了,不代表对方不再接收数据,因此被动关闭方进行ACK回复之后有可能还会继续发送数据,等到不再发送数据了才会发送下一个FIN包,因此FIN包和ACK包是分开的…

  • 第八章 软件项目团队管理

    第八章 软件项目团队管理本章内容提纲8.1软件项目团队管理概述8.2项目组织的规划8.3团队人员获取8.4团队建设和日常管理8.5沟通管理8.6软件专业人员的非技术素养8.1软件项目团队管理概述什么是软件项目团队?   软件项目团队是由软件项目的不同干系人所组成的,具有共同目标、紧密协作的集体。软件项目团队包括所有项目干系人:项目发起人、资助者、项目组(开发团队)、供应商、客户等。有时,软件项目团队特指项…

发表回复

您的电子邮箱地址不会被公开。

关注全栈程序员社区公众号