光流法

光流法

---

定义

光流法是空间运动物体在观察成像平面上的像素运动的瞬时速度，是利用图像序列中像素在时间域上的变化以及相邻帧之间的相关性来找到上一帧跟当前帧之间存在的对应关系，从而计算出相邻帧之间物体的运动信息的一种方法。一般而言，光流是由于场景中前景目标本身的移动、相机的运动，或者两者的共同运动所产生的。
简单来说，光流是空间运动物体在观测成像平面上的像素运动的“瞬时速度”。光流的研究是利用图像序列中的像素强度数据的时域变化和相关性来确定各自像素位置的“运动”。研究光流场的目的就是为了从图片序列中近似得到不能直接得到的运动场。

光流法的前提假设：
 （1）相邻帧之间的颜色恒定，对于灰度图来说，亮度恒定；
 （2）相邻视频帧的取帧时间连续，或者，相邻帧之间物体的运动比较“微小”；
 （3）保持空间一致性；即，同一子图像的像素点具有相同的运动

上面的三个假设对光流法的原理公式推导很重要，具体的推导公式不在具体阐述。同时也是判断使用光流法的使用条件，光流法主要分为两种，计算部分像素运动的：稀疏光流，以Lucas-Kanade为代表，计算所有像素运动的：稠密光流。下面将结合代码具体理解。

OpenCV中的LK光流：

在OpenCV库提供了一个完整的函数，cv2.calcOpticalFlowPyrLK()。
python的OpenCV 光流函数如下

该函数计算基于图像金字塔的稀疏光流

nextPts,status,err = cv.calcOpticalFlowPyrLK(   prevImg, nextImg, prevPts, nextPts[, status[, err[, winSize[, maxLevel[, criteria[, flags[, minEigThreshold]]]]]]])

返回值：

nextPtrs 输出一个二维点的向量，这个向量可以是用来作为光流算法的输入特征点，也是光流算法在当前帧找到特征点的新位置（浮点数）
status 标志，在当前帧当中发现的特征点标志status==1，否则为0
err 向量中的每个特征对应的错误率

输入值：

• prevImg 上一帧图片
• nextImg 当前帧图片
• prevPts 上一帧找到的特征点向量
• nextPts 与返回值中的nextPtrs相同
• status 与返回的status相同
• err 与返回的err相同
• winSize 在计算局部连续运动的窗口尺寸（在图像金字塔中）
• maxLevel 图像金字塔层数，0表示不使用金字塔
• criteria 寻找光流迭代终止的条件
• flags 有两个宏，表示两种计算方法，分别是OPTFLOW_USE_INITIAL_FLOW表示使用估计值作为寻找到的初始光流，OPTFLOW_LK_GET_MIN_EIGENVALS表示使用最小特征值作为误差测量
• minEigThreshold 该算法计算光流方程的2×2规范化矩阵的最小特征值，除以窗口中的像素数; 如果此值小于minEigThreshold，则会过滤掉相应的功能并且不会处理该光流，因此它允许删除坏点并获得性能提升。

大致流程就是，首先获取视频或者摄像头的第一帧图像。用goodFeaturesToTrack函数获取初始化的角点，然后开始无限循环获取视频图像帧，将新图像和上一帧图像放入calcOpticalFlowPyrLK函数当中，从而获取新图像的光流。

import numpy as np
import cv2


cap = cv2.VideoCapture(0)


# ShiTomasi 角点检测参数
feature_params = dict( maxCorners = 100,
                       qualityLevel = 0.3,
                       minDistance = 7,
                       blockSize = 7 )

# lucas kanade光流法参数
lk_params = dict( winSize  = (15,15),
                  maxLevel = 2,
                  criteria = (cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, 10, 0.03))

# 创建随机颜色
color = np.random.randint(0,255,(100,3))

# 获取第一帧，找到角点
ret, old_frame = cap.read()
#找到原始灰度图
old_gray = cv2.cvtColor(old_frame, cv2.COLOR_BGR2GRAY)

#获取图像中的角点，返回到p0中
p0 = cv2.goodFeaturesToTrack(old_gray, mask = None, **feature_params)

# 创建一个蒙版用来画轨迹
mask = np.zeros_like(old_frame)

while(1):
    ret,frame = cap.read()
    frame_gray = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)

    # 计算光流
    p1, st, err = cv2.calcOpticalFlowPyrLK(old_gray, frame_gray, p0, None, **lk_params)
    # 选取好的跟踪点
    good_new = p1[st==1]
    good_old = p0[st==1]

    # 画出轨迹
    for i,(new,old) in enumerate(zip(good_new,good_old)):
        a,b = new.ravel()
        c,d = old.ravel()
        mask = cv2.line(mask, (a,b),(c,d), color[i].tolist(), 2)
        frame = cv2.circle(frame,(a,b),5,color[i].tolist(),-1)
    img = cv2.add(frame,mask)

    cv2.imshow('frame',img)
    k = cv2.waitKey(30) & 0xff
    if k == 27:
        break

    # 更新上一帧的图像和追踪点
    old_gray = frame_gray.copy()
    p0 = good_new.reshape(-1,1,2)

cv2.destroyAllWindows()
cap.release()

源代码是使用一个视频，我这里面改成了摄像头了。

---

OpenCV中的稠密光流：

LK算法计算的是稀疏的特征点光流，如样例当中计算的是使用 Shi-Tomasi算法得到的特征点。opencv当总提供了查找稠密光流的方法。该方法计算一帧图像当中的所有点。该方法是基于Gunner Farneback提出的一篇论文Two-Frame Motion Estimation Based on Polynomial Expansion。

下面样例显示如何找到稠密光流，我们得到的一个两个通道的向量(u,v)。得到的该向量的大小和方向。用不同的颜色编码来使其可视化。
方向与Hue值相关，大小与Value值相关。

使用calcOpticalFlowFarneback函数得到

flow=cv.calcOpticalFlowFarneback(prev, next, flow, pyr_scale, levels, winsize, iterations, poly_n, poly_sigma, flags)

返回值是每个像素点的位移

参数

• prev 输入8位单通道图片
• next 下一帧图片，格式与prev相同
• flow 与返回值相同，得到一个CV_32FC2格式的光流图，与prev大小相同
• pyr_scale 构建图像金字塔尺度
• levels 图像金字塔层数
• winsize 窗口尺寸，值越大探测高速运动的物体越容易，但是越模糊，同时对噪声的容错性越强
• iterations 对每层金字塔的迭代次数
• poly_n 每个像素中找到多项式展开的邻域像素的大小。越大越光滑，也越稳定
• poly_sigma 高斯标准差，用来平滑倒数
• flags 光流的方式，有OPTFLOW_USE_INITIAL_FLOW 和OPTFLOW_FARNEBACK_GAUSSIAN 两种

import numpy as np
import cv2
import matplotlib.pyplot as plt

cap = cv2.VideoCapture('1/12311sc.mp4')
# params for ShiTomasi corner detection 特征点检测
feature_params = dict( maxCorners = 10,
                       qualityLevel = 0.1,
                       minDistance = 10,
                       blockSize = 3 )

# Parameters for lucas kanade optical flow光流法参数
lk_params = dict(winSize = (15,15),
                  maxLevel = 0,
                  criteria = (cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, 10, 0.03))

# Create some random colors 画轨迹
color = np.random.randint(0,255,(100,3))

# Take first frame and find corners in it
ret, old_frame = cap.read()
old_gray = cv2.cvtColor(old_frame, cv2.COLOR_BGR2GRAY)
roi = np.zeros_like(old_gray)
x,y,w,h = 266,143,150,150
roi[y:y+h, x:x+w] = 255
p0 = cv2.goodFeaturesToTrack(old_gray, mask = roi, **feature_params)

# Create a mask image for drawing purposes
mask = np.zeros_like(old_frame)

while(1):
    ret,frame = cap.read()
    if not ret:
      break
    frame_gray = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)

    # calculate optical flow
    p1, st, err = cv2.calcOpticalFlowPyrLK(old_gray, frame_gray, p0, None, **lk_params)

    # Select good points
    good_new = p1[st==1]
    good_old = p0[st==1]

    # draw the tracks
    for i,(new,old) in enumerate(zip(good_new,good_old)):
        a,b = new.ravel()
        c,d = old.ravel()
        mask = cv2.line(mask, (a,b),(c,d), color[i].tolist(), 2)
        frame = cv2.circle(frame,(a,b),3,color[i].tolist(),-1)
    img = cv2.add(frame,mask)

    cv2.imshow('frame',img)
    key = cv2.waitKey(60) & 0xff
    if key == 27:  # 按下ESC时，退出
        break
    elif key == ord(' '):  # 按下空格键时，暂停
        cv2.waitKey(0)

    # Now update the previous frame and previous points
    old_gray = frame_gray.copy()
    p0 = good_new.reshape(-1,1,2)

cv2.destroyAllWindows()
cap.release()

import numpy as np
import cv2


cap = cv2.VideoCapture(0)

#获取第一帧
ret, frame1 = cap.read()
prvs = cv2.cvtColor(frame1,cv2.COLOR_BGR2GRAY)
hsv = np.zeros_like(frame1)

#遍历每一行的第1列
hsv[...,1] = 255


while(1):
    ret, frame2 = cap.read()
    next = cv2.cvtColor(frame2,cv2.COLOR_BGR2GRAY)

    #返回一个两通道的光流向量，实际上是每个点的像素位移值
    flow = cv2.calcOpticalFlowFarneback(prvs,next, None, 0.5, 3, 15, 3, 5, 1.2, 0)

    #print(flow.shape)
    print(flow)

    #笛卡尔坐标转换为极坐标，获得极轴和极角
    mag, ang = cv2.cartToPolar(flow[...,0], flow[...,1])
    hsv[...,0] = ang*180/np.pi/2
    hsv[...,2] = cv2.normalize(mag,None,0,255,cv2.NORM_MINMAX)
    rgb = cv2.cvtColor(hsv,cv2.COLOR_HSV2BGR)

    cv2.imshow('frame2',rgb)
    k = cv2.waitKey(30) & 0xff
    if k == 27:
        break
    elif k == ord('s'):
        cv2.imwrite('opticalfb.png',frame2)
        cv2.imwrite('opticalhsv.png',rgb)
    prvs = next

cap.release()
cv2.destroyAllWindows()