Lucas Kanade光流实现无法正常工作

Question

我正在使用python和numpy上的lucas kanade method编写自己的光流脚本。但是我用这个算法（This is testing video）的opencv实现得到的流结果与我自己的不同。在

这是我的完整代码：

import cv2
import numpy as np

def drawlines (image, pts0, pts1, color):
    for pt0, pt1 in zip(pts0, pts1):
        cv2.line(image,(pt0[0],pt0[1]),(pt1[0],pt1[1]),color,thickness=2);

def calcFlowLib (prev, curr, points):
    #OPTICAL FLOW USING OPENCV IMPLEMENTATION
    currpts = np.reshape(points,(-1,1,2));
    lk_params = dict(winSize=(41,41), maxLevel=0, criteria=(cv2.TERM_CRITERIA_EPS|cv2.TERM_CRITERIA_COUNT,10,0.03));
    nextpts, stat, err = cv2.calcOpticalFlowPyrLK(prev,curr,currpts,None,**lk_params);
    nextpts = np.reshape(nextpts,(-1,2));
    return nextpts;

def calcFlow (im1, im2, points, halfwin=(20,20)):
    #MY OWN IMPLEMENTATION
    pts = np.copy(points)
    npts = points.shape[0]

    i_x = np.zeros(im1.shape)
    i_y = np.zeros(im1.shape)
    i_t = np.zeros(im1.shape)
    i_x[1:-1, 1:-1] = (im1[1:-1, 2:] - im1[1:-1, :-2]) / 2
    i_y[1:-1, 1:-1] = (im1[2:, 1:-1] - im1[:-2, 1:-1]) / 2
    i_t[1:-1, 1:-1] = im2[1:-1, 1:-1] - im1[1:-1, 1:-1]


    for i in xrange(0,npts,1):
        p = np.reshape(pts[i],(2,1))

        ix = i_x[p[1]-halfwin[1]:p[1]+halfwin[1]+1,p[0]-halfwin[0]:p[0]+halfwin[0]+1]
        iy = i_y[p[1]-halfwin[1]:p[1]+halfwin[1]+1,p[0]-halfwin[0]:p[0]+halfwin[0]+1]            
        it = i_t[p[1]-halfwin[1]:p[1]+halfwin[1]+1,p[0]-halfwin[0]:p[0]+halfwin[0]+1]

        ixx = ix*ix
        iyy = iy*iy
        ixy = ix*iy
        ixt = ix*it
        iyt = iy*it

        G = np.array([[np.sum(ixx),np.sum(ixy)],[np.sum(ixy),np.sum(iyy)]],dtype=np.float32)
        B = np.array([[np.sum(ixy)],[np.sum(iyt)]],dtype=np.float32)
        vel = np.linalg.lstsq(G,B)[0]
        pts[i] = pts[i] + np.ravel(vel)        

    return pts


cap = cv2.VideoCapture("OF.mp4")    
prevgray = None
currgray = None
prevpts0 = np.array([[320,180],[100,100]],dtype=np.float32)
prevpts1 = np.array([[320,180],[100,100]],dtype=np.float32)
while(cap.isOpened()):
    ret, frame = cap.read()
    if not ret:
        break
    frame = cv2.resize(frame,(640,360))
    currgray = cv2.cvtColor(frame,cv2.COLOR_BGR2GRAY)

    if(prevgray is not None):
        currpts0 = calcFlow(prevgray,currgray,prevpts0)
        drawlines(frame,prevpts0,currpts0,(0,255,0)) 
        currpts1 = calcFlowLib(prevgray,currgray,prevpts1)
        drawlines(frame,prevpts1,currpts1,(0,0,255))
        prevpts0 = np.copy(currpts0) 
        prevpts1 = np.copy(currpts1)

    prevgray = np.copy(currgray)
    cv2.imshow("Frame",frame)
    if cv2.waitKey(1) & 0xFF == ord('1'):
        break
cap.release()
cv2.destroyAllWindows()

我的实现的错误部分在哪里？