コンピュータサイエンス系勉強ノート

計算機科学に限らず日々学んだことを色々まとめていきます

オプティカルフローの実装(Python)

オプティカルフローを求めるための方法の1つであるLucas Kanade法を実装してみました.

勾配法

Lucal Kanade法を実装する前にまず勾配法について学ぶ必要があります.勾配法とは,連続する2枚の画像において追跡対象の物体の移動量は微小であると仮定し,オプティカルフローを求める方法です.まず時刻tにおける画像上の輝度値をI(x,y,t)としましょう.そして時間Δtの間にその画素が(Δx,Δy)動いたとしましょう.すると,

 I(x,y,t) = I(x + Δx, y + Δy, t + Δt) (1)

そしてIの微小変化を表現するするために,(1式)をテイラー展開で一次近似します.x,y,tが変化した分だけIも変化しその増分が加えられることを意味しています

 I(x,y,t) ≒ I(x,y,t) + \frac{∂I}{∂x}Δx + \frac{∂I}{∂y}Δy + \frac{∂I}{∂t}Δt (2)

ここで(2)について両辺のI(x,y,t)はほぼ同じとみなせるので移行して打ち消せます.そしてΔtで割ると,

\frac{Δx}{Δt}\frac{∂I}{∂x}+\frac{Δy}{Δt}\frac{∂I}{∂y}+\frac{∂I}{∂t}=0 (3)

そしてΔtが限りなく0に近いと仮定すると微分の形にもっていけます.つまり,

\frac{∂x}{∂t}\frac{∂I}{∂x}+\frac{∂y}{∂t}\frac{∂I}{∂y}+\frac{∂I}{∂t}=0 (4)

\frac{∂x}{∂t}=u\frac{∂y}{∂t}=vが時間変化に対する画素の移動方向を表しており,これがオプティカルフローとなります.また,\frac{∂I}{∂t}=I_x\frac{∂I}{∂y}=I_yが時間変化に対する輝度値の変化となります.そして最終的には,


I_xu+I_yv=-I_t (5)

と表され,時間変化に対するIの変化は,画素値の勾配と画素のオプティカルフローとの内積で表されることになります.

ここでオプティカルフロー(u,v)を求めたいわけですが,未知変数2つに対して式が一つしかありません.なのでこのままでは解けません.しかし,Lucal Kanade法ではある仮定を置くことで式の数を増やして(u,v)を求めることができます.


Lucal Kanade法

ではその仮定とは何かというと,注目画素の周辺の局所領域では(u,v)の値は同じだという仮定です.


f:id:clientver2:20160601001357p:plain:w500

こっからの数式は気力がなくなったのでそのうち書きます.今のところは他のサイトを参照・・

ソースコード

汚いコードですんません.なんか間違ってる気がしますが大体こんな感じです.

from scipy import ndimage
import cv2
import numpy as np

def opticalFlow(pre, next):
    image1 = cv2.cvtColor(pre, cv2.COLOR_RGB2GRAY).astype(np.float)
    image2 = cv2.cvtColor(next, cv2.COLOR_RGB2GRAY).astype(np.float)
    
    #輝度勾配
    dx = ndimage.gaussian_filter(image1, sigma=1.0, order=(0,1))
    dy = ndimage.gaussian_filter(image1, sigma=1.0, order=(1,0))
    
    #Iの時間変化 つまり次の時間の画像と前時間の画像との差分が時間変化
    de = image2 - image1

    height, width = image1.shape

    for i in range(1, height-1, 20):
        for j in range(1, width-1, 20):
            #局所領域(3*3)
            patchDx = dx[i-1:i+1+1,j-1:j+1+1]
            patchDy = dy[i-1:i+1+1,j-1:j+1+1]
            patchDe = de[i-1:i+1+1,j-1:j+1+1]

            A = np.zeros((2,2))
            A[0,0] = np.sum(patchDx**2)
            A[0,1] = np.sum(patchDx * patchDy)
            A[1,0] = A[0,1]
            A[1,1] = np.sum(patchDy**2)

            B = np.zeros((2))
            B[0] = -np.sum(patchDx * patchDe)
            B[1] = -np.sum(patchDy * patchDe)
            
            #ランクが2じゃないと逆行列もとめられん
            if np.linalg.matrix_rank(A) == 2:
                u, v = np.linalg.solve(A, B)
                #フローを描画
                next = cv2.line(next, (j,i), (int(j+u), int(i+v)), (255,0,0), 1)

    return next


def main():
    cap = cv2.VideoCapture(0)
    cv2.namedWindow("frame")
    preFrame = None

    while(True):
        ret, frame = cap.read()

        if preFrame != None:
            flow = opticalFlow(preFrame, frame)
            cv2.imshow("frame", flow)

        if cv2.waitKey(20) & 0xFF == 27:
            break

        preFrame = frame
    cap.release()
    cv2.destroyAllWindows()

main()