Optical Flow 光流法概述

主要參考論文：Two-frame motion estimation based on polynomial expansion (Farnbäck), 2003

---

前言

OpenCV 就有提供不同光流法的函式供開發者使用，基本上沒什麼太大的實作困難，等於是直接 call function
除非有要在加寫優化或是直接去改 OpenCV 的源碼

其實光流法背後的數學模型有很多種，我也並沒有每個都深掘，這邊就大概描述下原理

概念

光流法目的是在一連串影像中找到相鄰影像在時間域上的變化，通常會應用在影片上，找到上一個 frame 與當前的 frame 存在的對應關係，然後計算出物體的運動資訊。把每個影像中每個像素的運動速度和運動方向找出來

假設要追蹤影像中的A像素的流向，我們必須要知道A像素在下一個 frame 的位置在哪，也就是說我們要有辦法辨認出 A 像素才能得知他的位置，通常會是一個近似解

光流法最基本的假設是：光流場幾乎是光滑的，基於此假設用數學去分析像素的運動方向

稀疏光流與稠密光流

光流法有很多種實現方式，以結果來分類的話一個是稀疏光流，一個是稠密光流

Lucas-Kanade 的方法是稀疏光流，只考慮部分像素的方向，通常會是特徵明顯的像素，像是角點或是邊界點

特徵點可以讓OpenCV幫你定義或是看需求自定義

我所採用的 Farnback 的方法是屬於稠密光流，考慮所有像素的方向

Lucas 和 Farnback 都有利用影像金字塔來解決無法找到動量太大的向量
就是把影像解析度降低，像素量減少來讓大動量能夠被偵測到，因為 Lucas 是稀疏光流，運算速度比 Farnback 的稠密光流快

左邊是 LK，右邊是 Farnback

Farnback optical flow

用二次多項式對於二維圖像建模，將影像的二維座標 $(x,y)$ 當成是二維信號的函數，轉換到以 $(1,x^2,y^2,x,y,xy)$ 做為基函數的空間

代表影像需要一個六維向量來帶入不同位置求出像素值，Farnback將每個像素周圍都設定一個 $(2n+1)*(2n+1)$ 的鄰近區域

利用區域內的像素點當作最小誤差平方法的樣本點，再利用高斯分布賦予樣本點不同的權重，求解後近似得到中心像素的六維向量

Code

我用了兩種視覺化的方式，一種就是箭頭方向，另一種是塗色的方式

使用箭頭方向視覺化 flow

void drawArrow(cv::Mat& img, cv::Point pStart, cv::Point pEnd, int len, int alpha, cv::Scalar& color, int thickness){
	const double PI = 3.1415926;
	Point arrow;
	 
	double angle = atan2((double)(pStart.y - pEnd.y), (double)(pStart.x - pEnd.x));
	line(img, pStart, pEnd, color, thickness);
	
	arrow.x = pEnd.x + len * cos(angle + PI * alpha / 180);
	arrow.y = pEnd.y + len * sin(angle + PI * alpha / 180);
	line(img, pEnd, arrow, color, thickness);
	arrow.x = pEnd.x + len * cos(angle - PI * alpha / 180);
	arrow.y = pEnd.y + len * sin(angle - PI * alpha / 180);
	line(img, pEnd, arrow, color, thickness);
}

void drawLines(Mat flow, Mat &target) {
    // By y += 5, x += 5 you can specify the grid 
    for (int y = 0; y < target.rows; y += 10) {
        for (int x = 0; x < target.cols; x += 10)
        {
            Vec2f flow_at_point = flow.at<Vec2f>(y, x);
            float fx = flow_at_point[0];
            float fy = flow_at_point[1];
            if ((fabs(fx) > UNKNOWN_FLOW_THRESH) || (fabs(fy) >  UNKNOWN_FLOW_THRESH))
                continue;
            // get the flow from y, x position * 10 for better visibility
            const Point2f flowatxy = flow.at<Point2f>(y, x) * 2;

            // draw initial point
            //circle(target, Point(x, y), 1, Scalar(255, 0, 0), -1);
            // draw line at flow direction
            Point startP = Point(x, y);
			Point endP = Point(cvRound(x + flowatxy.x), cvRound(y + flowatxy.y));
            if (startP == endP)
                continue;
                
            line(target, startP, endP, Scalar(255, 0, 0), 1);
            // draw arrow 
            drawArrow(target, startP, endP, 5, 30, Scalar(255, 0, 0), 1);
        }
    }
}

使用塗色的方式視覺化flow

#define UNKNOWN_FLOW_THRESH 1e9

void makecolorwheel(vector<Scalar> &colorwheel)
{
	int RY = 15;
	int YG = 6;
	int GC = 4;
	int CB = 11;
	int BM = 13;
	int MR = 6;

	int i;

	for (i = 0; i < RY; i++) colorwheel.push_back(Scalar(255, 255 * i / RY, 0));
	for (i = 0; i < YG; i++) colorwheel.push_back(Scalar(255 - 255 * i / YG, 255, 0));
	for (i = 0; i < GC; i++) colorwheel.push_back(Scalar(0, 255, 255 * i / GC));
	for (i = 0; i < CB; i++) colorwheel.push_back(Scalar(0, 255 - 255 * i / CB, 255));
	for (i = 0; i < BM; i++) colorwheel.push_back(Scalar(255 * i / BM, 0, 255));
	for (i = 0; i < MR; i++) colorwheel.push_back(Scalar(255, 0, 255 - 255 * i / MR));
}

void motionToColor(Mat flow, Mat &color){
    if (color.empty())
        color.create(flow.rows, flow.cols, CV_8UC3);
        
    static vector<Scalar> colorwheel; //Scalar r,g,b   
    
    if (colorwheel.empty())
        makecolorwheel(colorwheel);
        
    // determine motion range:   
    float maxrad = -1;
    
    // Find max flow to normalize fx and fy   
    for (int i = 0; i < flow.rows; ++i){
        for (int j = 0; j < flow.cols; ++j){
            Vec2f flow_at_point = flow.at<Vec2f>(i, j);
            float fx = flow_at_point[0];
            float fy = flow_at_point[1];
            if ((fabs(fx) > UNKNOWN_FLOW_THRESH) || (fabs(fy) >  UNKNOWN_FLOW_THRESH))
                continue;
            float rad = sqrt(fx * fx + fy * fy);
            maxrad = maxrad > rad ? maxrad : rad;
        }
    }
    for (int i = 0; i < flow.rows; ++i){
        for (int j = 0; j < flow.cols; ++j){
            uchar *data = color.data + color.step[0] * i + color.step[1] * j;
            Vec2f flow_at_point = flow.at<Vec2f>(i, j);
            
            float fx = flow_at_point[0] / maxrad;
            float fy = flow_at_point[1] / maxrad;
            if ((fabs(fx) > UNKNOWN_FLOW_THRESH) || (fabs(fy) >  UNKNOWN_FLOW_THRESH)){
                data[0] = data[1] = data[2] = 0;
                continue;
            }
            float rad = sqrt(fx * fx + fy * fy);
            float angle = atan2(-fy, -fx) / CV_PI;
            float fk = (angle + 1.0) / 2.0 * (colorwheel.size() - 1);
            int k0 = (int)fk;
            int k1 = (k0 + 1) % colorwheel.size();
            float f = fk - k0;
            
            for (int b = 0; b < 3; b++){
                float col0 = colorwheel[k0][b] / 255.0;
                float col1 = colorwheel[k1][b] / 255.0;
                float col = (1 - f) * col0 + f * col1;
                if (rad <= 1)
                    col = 1 - rad * (1 - col); // increase saturation with radius   
                else
                    col *= .75; // out of range   
                data[2 - b] = (int)(255.0 * col);
            }
        }
    }
}

主程式

int main(){
    
    VideoCapture cap("video.mp4");  // load your video

    if (!cap.isOpened()){
        cout << "Video Capture Fail" << endl;
        return 0;
    }
    
    double frameNum = cap.get(CV_CAP_PROP_FRAME_COUNT);
    cout << frameNum << endl;

	Mat flow;
    
    // some faster than mat image container
    UMat  flowUmat, prevgray;
    for (int i = 0; i < 500; i++) {   // loop number depends on you (<= frames)
        Mat img;
        Mat original;
        // capture frame from video file
        cap >> img;
        imshow("original", img);
        
        /* save original for later */
        img.copyTo(original);
        /* just make current frame gray */
        cvtColor(img, img, COLOR_BGR2GRAY);
        
        /* not a fisrt frame, do optical flow */
        if (prevgray.empty() == false) {
            // calculate optical flow 
            calcOpticalFlowFarneback(prevgray, img, flowUmat, 0.5, 3, 15, 3, 5, 1.2, 0);
            // copy Umat container to standard Mat
            flowUmat.copyTo(flow);
            
            /* -------------- visulize flow ----------------*/
            Mat motion2color;
            motionToColor(flow, motion2color);
            imshow("visualize flow 1", motion2color);
            Mat blackbg(original.rows, original.cols, CV_8UC3, cv::Scalar(0, 0, 0));
            
            drawLines(flow, original);
            imshow("optical flow", original);
            
            // fill previous image again
            img.copyTo(prevgray);
        }
        else {  // first frame, fill previous image
            img.copyTo(prevgray);
        }
        waitKey(1);
    }
}