Horn-Schunck 光流实现问题

Posted

技术标签:

【中文标题】Horn-Schunck 光流实现问题【英文标题】:Horn-Schunck optical flow implementation issue 【发布时间】:2015-01-12 14:26:12 【问题描述】:

我正在尝试通过 NumPy 和 OpenCV 实现 Horn-Schunck 光流算法 我使用Horn-Schunck method on wiki 和original paper

但我的实现在以下简单示例中失败

第一帧:

[[  0   0   0   0   0   0   0   0   0   0]
 [  0 255 255   0   0   0   0   0   0   0]
 [  0 255 255   0   0   0   0   0   0   0]
 [  0   0   0   0   0   0   0   0   0   0]
 [  0   0   0   0   0   0   0   0   0   0]]

第二帧:

[[  0   0   0   0   0   0   0   0   0   0]
 [  0   0   0 255 255   0   0   0   0   0]
 [  0   0   0 255 255   0   0   0   0   0]
 [  0   0   0   0   0   0   0   0   0   0]
 [  0   0   0   0   0   0   0   0   0   0]]

这只是一个白色的小矩形,在 frame2 上移动了 2 个像素 我的实现产生以下流程 U 部分流程(我将 np.round 应用于流程的每个部分。原始值几乎相同):

[[ 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.]
 [ 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.]
 [ 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.]
 [ 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.]
 [ 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.]]

V 部分流:

[[ 0.  1.  0. -1. -0.  0.  0.  0.  0.  0.]
 [-0. -0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [-0. -1. -0.  1.  0.  0.  0.  0.  0.  0.]
 [-0. -0. -0.  0.  0.  0.  0.  0.  0.  0.]
 [-0. -0. -0.  0.  0.  0.  0.  0.  0.  0.]]

看起来这个流程是不正确的(因为如果我将 frame2 的每个像素都朝相应的流程组件的方向移动,我永远不会得到 frame1) 我的实现在真实图像上也失败了

但是,如果我将矩形向右(或向左或向上或向下)移动 1 个像素,我的实现会产生: U部分流量:

[[1 1 1 .....]
[1 1 1 .....]
......
[1 1 1 .....]]

V 部分流:

[[0 0 0 .....]
[0 0 0 .....]
......
[0 0 0 .....]]

我认为这个流程是正确的,因为我可以按照以下步骤重建第 1 帧

def translateBrute(img, u, v):
    res = np.zeros_like(img)
    u = np.round(u).astype(np.int)
    v = np.round(v).astype(np.int)
    for i in xrange(img.shape[0]):
        for j in xrange(img.shape[1]):
            res[i, j] = takePixel(img, i + v[i, j], j + u[i, j])
    return res

其中 takePixel 是一个简单的函数,如果输入坐标位于图像内部,则返回像素强度,否则返回图像边框上的强度

这是我的实现

import cv2
import sys
import numpy as np

def takePixel(img, i, j):
    i = i if i >= 0 else 0
    j = j if j >= 0 else 0    
    i = i if i < img.shape[0] else img.shape[0] - 1
    j = j if j < img.shape[1] else img.shape[1] - 1
    return img[i, j]

#Numerical derivatives from original paper: http://people.csail.mit.edu/bkph/papers/Optical_Flow_OPT.pdf
def xDer(img1, img2):
    res = np.zeros_like(img1)
    for i in xrange(res.shape[0]):
        for j in xrange(res.shape[1]):
            sm = 0
            sm += takePixel(img1, i,     j + 1) - takePixel(img1, i,     j)
            sm += takePixel(img1, i + 1, j + 1) - takePixel(img1, i + 1, j)
            sm += takePixel(img2, i,     j + 1) - takePixel(img2, i,     j)
            sm += takePixel(img2, i + 1, j + 1) - takePixel(img2, i + 1, j)
            sm /= 4.0
            res[i, j] = sm
    return res

def yDer(img1, img2):
    res = np.zeros_like(img1)
    for i in xrange(res.shape[0]):
        for j in xrange(res.shape[1]):
            sm = 0
            sm += takePixel(img1, i + 1, j    ) - takePixel(img1, i, j    )
            sm += takePixel(img1, i + 1, j + 1) - takePixel(img1, i, j + 1)
            sm += takePixel(img2, i + 1, j    ) - takePixel(img2, i, j    )
            sm += takePixel(img2, i + 1, j + 1) - takePixel(img2, i, j + 1)
            sm /= 4.0
            res[i, j] = sm
    return res

def tDer(img, img2):
    res = np.zeros_like(img)
    for i in xrange(res.shape[0]):
        for j in xrange(res.shape[1]):
            sm = 0
            for ii in xrange(i, i + 2):
                for jj in xrange(j, j + 2):
                    sm += takePixel(img2, ii, jj) - takePixel(img, ii, jj)
            sm /= 4.0
            res[i, j] = sm
    return res

averageKernel = np.array([[ 0.08333333,  0.16666667,  0.08333333],
                          [ 0.16666667,  0.        ,  0.16666667],
                          [ 0.08333333,  0.16666667,  0.08333333]], dtype=np.float32)
#average intensity around flow in point i,j. I use filter2D to improve performance. 
def average(img):
    return cv2.filter2D(img.astype(np.float32), -1, averageKernel)

def translateBrute(img, u, v):
    res = np.zeros_like(img)
    u = np.round(u).astype(np.int)
    v = np.round(v).astype(np.int)
    for i in xrange(img.shape[0]):
        for j in xrange(img.shape[1]):
            res[i, j] = takePixel(img, i + v[i, j], j + u[i, j])
    return res

#Core of algorithm. Iterative scheme from wiki: https://en.wikipedia.org/wiki/Horn%E2%80%93Schunck_method#Mathematical_details
def hornShunckFlow(img1, img2, alpha):
    img1 = img1.astype(np.float32)
    img2 = img2.astype(np.float32)

    Idx = xDer(img1, img2)
    Idy = yDer(img1, img2)
    Idt = tDer(img1, img2)

    u = np.zeros_like(img1)
    v = np.zeros_like(img1)

    #100 iterations enough for small example
    for iteration in xrange(100):
        u0 = np.copy(u)
        v0 = np.copy(v)

        uAvg = average(u0)
        vAvg = average(v0)
        # '*', '+', '/' operations in numpy works component-wise
        u = uAvg - 1.0/(alpha**2 + Idx**2 + Idy**2) * Idx * (Idx * uAvg + Idy * vAvg + Idt)
        v = vAvg - 1.0/(alpha**2 + Idx**2 + Idy**2) * Idy * (Idx * uAvg + Idy * vAvg + Idt)
        if  iteration % 10 == 0:
            print 'iteration', iteration, np.linalg.norm(u - u0) + np.linalg.norm(v - v0)

    return u, v

if __name__ == '__main__':
    img1c = cv2.imread(sys.argv[1])
    img2c = cv2.imread(sys.argv[2])
    img1g = cv2.cvtColor(img1c, cv2.COLOR_BGR2GRAY)
    img2g = cv2.cvtColor(img2c, cv2.COLOR_BGR2GRAY)

    u, v = hornShunckFlow(img1g, img2g, 0.1)
    imgRes = translateBrute(img2g, u, v)
    cv2.imwrite('res.png', imgRes)
    print img1g
    print translateBrute(img2g, u, v)

优化方案取自wikipedia,数值导数取自原始论文。

有人知道为什么我的实现会产生不正确的流程吗? 如有必要,我可以提供任何其他信息

PS 对不起我的英语不好

更新: 我实现了 Horn-Schunck 成本函数

def grad(img):
    Idx = cv2.filter2D(img, -1, np.array([
            [-1, -2, -1], 
            [ 0,  0,  0], 
            [ 1,  2,  1]], dtype=np.float32))
    Idy = cv2.filter2D(img, -1, np.array([
            [-1, 0, 1], 
            [-2, 0, 2], 
            [-1, 0, 1]], dtype=np.float32))
    return Idx, Idy

def hornShunckCost(Idx, Idy, Idt, u, v, alpha):
    #return sum(sum(It**2))
    udx, udy = grad(u)
    vdx, vdy = grad(v)
    return (sum(sum((Idx*u + Idy*v + Idt)**2)) + 
            (alpha**2)*(sum(sum(udx**2)) +
                        sum(sum(udy**2)) +
                        sum(sum(vdx**2)) +
                        sum(sum(vdy**2))
                ))

并在迭代中检查此函数的值

if  iteration % 10 == 0:
            print 'iter', iteration, np.linalg.norm(u - u0) + np.linalg.norm(v - v0)
            print hornShunckCost(Idx, Idy, Idt, u, v, alpha)

如果我使用一个简单的例子,矩形移动了一个像素,一切都很好:成本函数的值每一步都会减少。 但是在示例中,矩形移动了两个像素,成本函数的值每一步都会增加。 算法的这种行为真的很奇怪 也许我选择了不正确的方法来计算成本函数。

【问题讨论】:

【参考方案1】:

我失去了一个事实,即经典的 Horn-Schunck 方案使用线性化数据项 (I1(x, y) - I2(x + u(x, y), y + v(x, y)))。这种线性化使优化变得容易,但不允许大位移

为了处理大位移,有下一个方法Pyramidal Horn-Schunck

【讨论】:

以上是关于Horn-Schunck 光流实现问题的主要内容,如果未能解决你的问题,请参考以下文章

目标跟踪之Horn-Schunck光流法

图像配准基于Horn-Schunck和Lucas-Kanade等光流场实现图像配准matlab源码含GUI界面

Javascript 光流实现

神经光流网络——用卷积网络实现光流预测(FlowNet: Learning Optical Flow with Convolutional Networks)

opencv 怎样提取运动物体的轮廓?

OpenCV4中DIS光流算法与应用