opencv - Horn-Schunck 光流实现问题

Question

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

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

第 1 帧:

[[  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]]

第 2 帧:

[[  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]]

这只是在第 2 帧上移动 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)

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

Answer 1

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

要处理大位移，有下一种方法 Pyramidal Horn-Schunck