在 numpy 中加速矢量化眼球追踪算法

我正在尝试实现 Fabian Timm 的眼球追踪算法 [http://www.inb.uni-luebeck.de/publikationen/pdfs/TiBa11b.pdf] http://www.inb.uni-luebeck.de/publikationen/pdfs/TiBa11b.pdf%5D（在这里找到：[http://thume.ca/projects/2012/11/04/simple-accurate-eye-center-tracking-in-opencv/] http://thume.ca/projects/2012/11/04/simple-accurate-eye-center-tracking-in-opencv/%5D）在 numpy 和 OpenCV 中我遇到了障碍。我认为我已经对我的实现进行了足够的矢量化，但它仍然不够快，无法实时运行，而且它检测学生的准确度没有我希望的那么高。这是我第一次使用 numpy，所以我不确定我做错了什么。

def find_pupil(eye):
    eye_len = np.arange(eye.shape[0])
    xx,yy = np.meshgrid(eye_len,eye_len) #coordinates
    XX,YY = np.meshgrid(xx.ravel(),yy.ravel()) #all distance vectors
    Dx,Dy = [YY-XX, YY-XX] #y2-y1, x2-x1 -- simpler this way because YY = XXT
    Dlen = np.sqrt(Dx**2+Dy**2)
    Dx,Dy = [Dx/Dlen, Dy/Dlen] #normalized

    Gx,Gy = np.gradient(eye)
    Gmagn = np.sqrt(Gx**2+Gy**2)

    Gx,Gy = [Gx/Gmagn,Gy/Gmagn] #normalized
    GX,GY = np.meshgrid(Gx.ravel(),Gy.ravel())

    X = (GX*Dx+GY*Dy)**2
    eye = cv2.bitwise_not(cv2.GaussianBlur(eye,(5,5),0.005*eye.shape[1])) #inverting and blurring eye for use as w
    eyem = np.repeat(eye.ravel()[np.newaxis,:],eye.size,0)
    C = (np.nansum(eyem*X, axis=0)/eye.size).reshape(eye.shape)

    return np.unravel_index(C.argmax(), C.shape)

以及其余代码：

def find_eyes(face):
    left_x, left_y = [int(floor(0.5 * face.shape[0])), int(floor(0.2 * face.shape[1]))]
    right_x, right_y = [int(floor(0.1 * face.shape[0])), int(floor(0.2 * face.shape[1]))]
    area = int(floor(0.2 * face.shape[0]))
    left_eye = (left_x, left_y, area, area)
    right_eye = (right_x, right_y, area, area)

    return [left_eye,right_eye]



faceCascade = cv2.CascadeClassifier("haarcascade_frontalface_default.xml")
video_capture = cv2.VideoCapture(0)

while True:
    # Capture frame-by-frame
    ret, frame = video_capture.read()

    gray = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)

    faces = faceCascade.detectMultiScale(
        gray,
        scaleFactor=1.1,
        minNeighbors=5,
        minSize=(30, 30),
        flags=cv2.CASCADE_SCALE_IMAGE
    )

    # Draw a rectangle around the faces
    for (x, y, w, h) in faces:
        cv2.rectangle(frame, (x, y), (x+w, y+h), (0, 255, 0), 2)
        roi_gray = gray[y:y+h, x:x+w]
        roi_color = frame[y:y+h, x:x+w]
        eyes = find_eyes(roi_gray)
        for (ex,ey,ew,eh) in eyes:
            eye_gray = roi_gray[ey:ey+eh,ex:ex+ew]
            eye_color = roi_color[ey:ey+eh,ex:ex+ew]
            cv2.rectangle(roi_color,(ex,ey),(ex+ew,ey+eh),(255,0,0),2)
            px,py = find_pupil(eye_gray)
            cv2.rectangle(eye_color,(px,py),(px+1,py+1),(255,0,0),2)

    # Display the resulting frame
    cv2.imshow('Video', frame)

    if cv2.waitKey(1) & 0xFF == ord('q'):
        break

# When everything is done, release the capture
video_capture.release()
cv2.destroyAllWindows()

您可以执行许多保存复制元素的操作，然后通过在创建单例维度后直接执行数学运算来执行一些数学运算，这将允许NumPy broadcasting http://docs.scipy.org/doc/numpy-1.10.1/user/basics.broadcasting.html。因此，有两个好处：即时操作可以节省工作空间内存并提高性能。另外，最后我们可以替换nansum使用简化版本进行计算。因此，考虑到所有这些哲学，这里有一种修改后的方法 -

def find_pupil_v2(face, x, y, w, h):    
    eye = face[x:x+w,y:y+h]
    eye_len = np.arange(eye.shape[0])

    N = eye_len.size**2
    eye_len_diff = eye_len[:,None] - eye_len
    Dlen = np.sqrt(2*((eye_len_diff)**2))
    Dxy0 = eye_len_diff/Dlen 

    Gx0,Gy0 = np.gradient(eye)
    Gmagn = np.sqrt(Gx0**2+Gy0**2)
    Gx,Gy = [Gx0/Gmagn,Gy0/Gmagn] #normalized

    B0 = Gy[:,:,None]*Dxy0[:,None,:]
    C0 = Gx[:,None,:]*Dxy0
    X = ((C0.transpose(1,0,2)[:,None,:,:]+B0[:,:,None,:]).reshape(N,N))**2

    eye1 = cv2.bitwise_not(cv2.GaussianBlur(eye,(5,5),0.005*eye.shape[1]))
    C = (np.nansum(X,0)*eye1.ravel()/eye1.size).reshape(eye1.shape)

    return np.unravel_index(C.argmax(), C.shape)

There's one repeat still left in it at Dxy. It might be possible to avoid that step and Dxy0 could be fed directly into the step that uses Dxy to give us X, but I haven't worked through it. Everything's converted to broadcasting based!

运行时测试和输出验证 -

In [539]: # Inputs with random elements
     ...: face = np.random.randint(0,10,(256,256)).astype('uint8')
     ...: x = 40
     ...: y = 60
     ...: w = 64
     ...: h = 64
     ...: 

In [540]: find_pupil(face,x,y,w,h)
Out[540]: (32, 63)

In [541]: find_pupil_v2(face,x,y,w,h)
Out[541]: (32, 63)

In [542]: %timeit find_pupil(face,x,y,w,h)
1 loops, best of 3: 4.15 s per loop

In [543]: %timeit find_pupil_v2(face,x,y,w,h)
1 loops, best of 3: 529 ms per loop

看来我们正在接近8x加速！