65.9K
CodeProject 正在变化。 阅读更多。
Home

基于多项式基近似的密集光流扩展

pi19404
emptyStarIconemptyStarIconemptyStarIconemptyStarIconemptyStarIcon

0/5 (0投票)

2014年10月23日

CPOL

4分钟阅读

viewsIcon

22019

基于多项式展开的稠密运动估计引言在本文中，我们将探讨基于图像多项式表示的稠密运动估计。通过逼近图像局部邻域，可以获得图像的多项式基函数表示。

引言

在本文中，我们将探讨基于图像多项式表示的稠密运动估计。通过使用二次多项式基函数逼近图像的局部邻域，可以获得图像的多项式基函数表示。通过等同基函数的系数，可以得到相邻帧之间的位移。

引言

  • 本文描述了一种快速稠密光流计算算法，由[4]提出。在之前的文章中，我们已经看到可以使用多项式基函数来表示图像的局部邻域。利用这种表示方法，可以估计图像中每个点的稠密光流。

二维基函数

多项式基

  • 假设多项式在平移下如何转换，使用从当前帧和前一帧派生的多项式展开系数，并估计位移向量。
  • 多项式展开的思想是用一个多项式逼近二维函数中一个点的邻域。考虑二次多项式基函数 \(1,x^2,y^2,x,y,xy\)，图像邻域中的像素值由以下方式表示：
    $f(x) ~ x^T A x + b^T x + c$
    其中A是对称矩阵，b是向量，c是标量。
  • 正如在之前的文章中所见，可以通过对邻域像素值进行加权最小二乘估计来估计系数。与所有光流算法一样，我们假设亮度恒定。相邻帧中图像路径的亮度是恒定的。
  • 考虑在图像中的点 \(\mathbf{(x,y)}\) 处遇到的平移运动 \(d\)。
    $\begin{eqnarray*}\\ f_1(x) = x^T A_1 x + b_{1}^Tx +c_1 \\ f_2(x) = f_1(x-d) = (x-d)^T A_1 (x-d) + b_{1}^T(x-d) +c_1 \\ f_2(x) = f_1(x-d) = (x)^T A_1 (x) + (b_{1}-2A_1d)^T(x) \\ + d^TAd-b_1^Td+c_1 \\ \end{eqnarray*}$

    根据亮度恒定假设，等同于两个多项式中的系数。

    $\begin{eqnarray*}\\ A_2 = A_1 \\ b_2=(b_{1}-2A_1d)\\ c_2=d^TAd-b_1^Td+c_1 \end{eqnarray*}$

    假设 \(A\) 非奇异

    $\begin{eqnarray*}\\ d=-\frac{1}{2}A^-1(b_2-b_1)\\ \end{eqnarray*}$

    因此，通过等同多项式的系数，可以得到图像中每个点的位移向量，前提是相邻帧中的感兴趣区域（即图像邻域）之间存在重叠。

  • 假设我们有位移估计 \(\bar{d}\)。我们提取点 \(P(x,y)\) 和点 \(P(x+d.x,y+d.y)\) 附近的感兴趣区域（ROI）。提取多项式基，并执行上面所示的计算。
  • 总位移可以估计为：
    $\begin{eqnarray*} &\bar{b}_2=b_{1}-2A_1 (\bar{d}) \\ & \text { 我们知道 $\bar{d}$，$A_1$，$b_1$} \\ & d = -0.5*A^{-1} (b_2+\bar{b}_2 - b_1) \end{eqnarray*}$
  • 因此，可以使用迭代方案，其中在每次连续迭代中获得更好的位移向量估计。当连续迭代中位移向量的变化低于阈值或完成特定数量的迭代时，可以终止迭代。位移向量的初始估计假设为 \(\mathbf{(0,0)}\)。因此，当前帧和前一帧的感兴趣区域或图像块是相同的。
01	/**
02	 * @brief updatePoly : The function computes the optical flow displacement vector
03	 * @param ptr1 : pointer to array containing the polynomial basis components
04	 * @param ptr2 : pointer to array containing the polynomial basis components
05	 * @param d    : estimate of displacement vector at the current point
06	 * @param flag : the flag indicating if the point is border pixels or not
07	 * @param M    : pointer to output array returning the computed coefficients
08	 * @param x    : co-ordinate at the current point where the coefficients is evaluated
09	 * @param y    :
10	 * @param s    : windows size for averaging
11	 */
12	void updatePoly(const float *ptr1,const float *ptr2,Point2f d,bool flag,
13	float *M,Point p,Size s)
14	{
15	    int x=p.x;
16	    int y=p.y;
17	    const int BORDER = 5;
18	    static const float border[BORDER] = {0.14f, 0.14f, 0.4472f, 0.4472f, 0.4472f};
19	 
20	    float r[5];
21	    float r2,r3,r4,r5,r6;
22	    if(flag==true)
23	    {
24	    //average A_1 and A_2
25	    r4=(ptr1[2]+ptr2[2])*0.5;
26	    r5=(ptr1[3]+ptr2[3])*0.5;
27	    r6=(ptr1[4]+ptr2[4])*0.25;
28	 
29	    }
30	    else
31	    {
32	        r[0]=0.f;
33	        r[1]=0.f;
34	        r[2]=ptr1[2];
35	        r[3]=ptr1[3];
36	        r[4]=ptr1[4]*0.5;
37	        r2=r[0];
38	        r3=r[1];
39	        r4=r[2];
40	        r5=r[3];
41	        r6=ptr1[4]*0.5;
42	 
43	 
44	    }
45	     
46	    //computing -(b1-b2)       
47	    r2=(ptr1[0]-ptr2[0])*0.5;
48	    r3=(ptr1[1]-ptr2[1])*0.5;
49	 
50	    //sum for iterative estimation b2=b2+\bar{b2}
51	    //r2 += r4*d.y + r6*d.x;
52	    ///r3 += r6*d.y + r5*d.x;
53	 
54	    r2 += r4*d.x + r6*d.y;
55	    r3 += r6*d.x + r5*d.y;
56	 
57	 
58	    if( (unsigned)(x - BORDER) >= (unsigned)(s.width - BORDER*2) ||
59	        (unsigned)(y - BORDER) >= (unsigned)(s.height - BORDER*2))
60	    {
61	        float scale = (x < BORDER ? border[x] : 1.f)*
62	            (x >= s.width - BORDER ? border[s.width - x - 1] : 1.f)*
63	            (y < BORDER ? border[y] : 1.f)*
64	            (y >= s.height - BORDER ? border[s.height - y - 1] : 1.f);
65	 
66	        r2 *= scale; r3 *= scale; r4 *= scale;
67	        r5 *= scale; r6 *= scale;
68	    }
69	 
70	 
71	   //computing final displacement d=A^-1(b2-b1)*0.5
72	    M[0]   = r4*r4 + r6*r6; // G(1,1)
73	    M[1] = (r4 + r5)*r6;  // G(1,2)=G(2,1)
74	    M[2] = r5*r5 + r6*r6; // G(2,2)
75	    M[3] = r4*r2 + r6*r3; // h(1)
76	    M[4] = r6*r2 + r5*r3; // h(2)
77	}
  • 方法 \(\textbf{EstimateFlow}\) 计算位移场计算所需的系数 \(A,b_2,b_1\)。\(\textbf{EstimateFlow}\) 函数为图像中的每个像素调用 \(\textbf{UpdatePoly}\) 方法。
  • 获得的位移场可能是不连续的，包含噪声和其他伪影。因为可以合理地假设，如果一个点遇到运动，那么其邻域像素也会遇到相同的运动。位移向量可以跨邻域平均，以获得更好的位移场估计。
  • 方法 \(\textbf{AverageFlow}\) 计算系数 \(A_1,b_2,b_1\) 的平均值，然后计算位移流场。
01	/**
02	  * @brief AverageFlow
03	  * @param _R0   : Polynomial basis of prev frame
04	  * @param _R1   : Polynomial basis coefficients of current frame
05	  * @param _flow :estimate of current displacement field
06	  * @param matM :containing coefficients of polynomial basis for computing displacemnt field
07	  */
08	 void AverageFlow( const Mat& _R0, const Mat& _R1,Mat& _flow, Mat& matM)
09	{
10	    int x, y, i, width = _flow.cols, height = _flow.rows;
11	 
12	//computing the average flow field
13	cv::GaussianBlur(matM,matM,Size(15,15),sigma);
14	     
15	     
16	    for( y = 0; y < height; y++ )
17	    {
18	        double g11, g12, g22, h1, h2;
19	        float* flow = (float*)(_flow.data + _flow.step*y);
20	        float *coeff=(float *)(matM.data+matM.step*y);
21	        // update borders
22	        for( x = 0; x < width; x++ )
23	        {
24	            g11 = coeff[x*5];
25	            g12 = coeff[x*5+1];
26	            g22 = coeff[x*5+2];
27	            h1 = coeff[x*5+3];
28	            h2 = coeff[x*5+4];
29	    //computing determinant for inverse
30	            double idet = 1./(g11*g22 - g12*g12 + 1e-3);
31	    //computing displacement flow field
32	            flow[x*2] = (float)((g11*h2-g12*h1)*idet);
33	            flow[x*2+1] = (float)((g22*h1-g12*h2)*idet);
34	        }
35	 
36	 
37	    }
38	 //calling EstimateFlow for updading coefficients
39	 //for computation of next iteration
40	     EstimateFlow(_R0, _R1, _flow, matM,0,flow.rows);
41	}
  • 这种方法在位移较大的情况下可能无效。因此，执行多尺度估计。流场估计在最小分辨率下进行。在较低分辨率下计算的位移被用作估计值，用于在更高分辨率下执行位移场计算。
  • 图 2b 显示了对两帧计算出的稠密光流。

参考系 1

参考系 2

光流

结论

本文描述了基于[4]论文的稠密光流算法的理论和实现细节。该算法的代码可以在github存储库 https://github.com/pi19404/OpenVision 的 DenseOf.cpp 和 DenseOf.hpp 文件中找到。在未来的文章中，我们将探讨使用 SSE、NEON 和 OpenCL 优化代码，以实现稠密光流场的实时计算。

参考文献

  1. Kenneth Andersson 和 Hans Knutsson。Continuous normalized convolution。在：ICME (1)。IEEE，2002，第 725-728 页。 isbn: 0-7803-7304-9。 url: http://dblp.uni-trier.de/db/conf/icmcs/icme2002-1.html
  2. Kenneth Andersson、Carl-Fredrik Westin 和 Hans Knutsson。Prediction from o -grid samples using continuous normalized convolution。在：Signal Processing 87.3 (2007 年 3 月 22 日)，第 353-365 页。 url: http://dblp.uni- trier.de/db/journals/sigpro/sigpro87.html
  3. Gunnar Farnebäck。Motion-based Segmentation of Image Sequences 。LiTH-ISY-EX-1596。MA 论文。SE-581 83 林雪平，瑞典：林雪平大学，1996。
  4. Gunnar Farnebäck。Polynomial Expansion for Orientation and Motion Estimation 。论文号 790，ISBN 91-7373-475-6。博士论文。SE-581 83 林雪平，瑞典：林雪平大学，瑞典，2002。
  5. Gunnar Farneback。Two-Frame Motion Estimation Based on Polynomial Expan-sion 。在：SCIA。LNCS 2749。哥德堡，瑞典，2003 年，第 363-370 页。
© . All rights reserved.