光流法-基于梯度

. 基于正则化的光流法估计
 

Horn和Schunck得到了能量泛函：

由于对d做高阶梯度必须没有不连续点，并且此做法光流场的扩散机制是各项同性的，则运动场的边界保持的不好。
由于上述模型的缺陷：后来又有了运动边界保持的光流法正则化模型和图像驱动的光流法正则化模型。
运动边界保持的光流法正则化模型：

图像驱动的光流法正则化模型：

上面这两个屌模型反正现在还没看懂~~~~~~~~~~~~~。
 2.TV+L1变分模型
因为上面是采用亮哥的符号，这里采用TV-L1 Optical Flow Estimation的符号，重新定义光束方程符号：

I表示图像序列，u表示瞬时速度。
上面的模型基本上都是基于光束约束方程展开，光束方程对于图像在时间上有很好的连续性。事实上，存在着其它情况(大尺度运动，光照不均)，这里就需要改进光束方程。
非线性模型：

非线性模型的泰勒展开的线性模型：

u0和u非常接近。
得到能力方程：

采用L1范式作为保真项
为了解决此方程，引用其他变量，采用凸松弛求解：

求解：

 第一个子问题解决可以通过Rudin–Osher–Fatemi的全变分去噪模型解决，第二个子问题可以通过逐点阈值法解决。
解决上面第一个子问题，需要引入对偶变量p1和p2。

 算法：
 

结果：

算法过程：
1.求I1的方向导数：gradient(I1, I1x, I1y, nx, ny);//I1x为I1的x方向导数，I1y为I1的y方向导数，nx宽,ny高

void centered_gradient(float *input,float *dx,float *dy,int nx,int ny)

{

for (int i = 1; i < ny-1; i++)

{

for(int j = 1; j < nx-1; j++)

{

int k = i * nx + j;

dx[k] = float(0.5*(input[k+1] - input[k-1]));

dy[k] = float(0.5*(input[k+nx] - input[k-nx]));

}

}

for (int j = 1; j < nx-1; j++)

{

dx[j] = float(0.5*(input[j+1] - input[j-1]));

dy[j] = float(0.5*(input[j+nx] - input[j]));

int k = (ny - 1) * nx + j;

dx[k] = float(0.5*(input[k+1] - input[k-1]));

dy[k] = float(0.5*(input[k] - input[k-nx]));

}

for(int i = 1; i < ny-1; i++)

{

int p=i * nx;

dx[p] = float(0.5*(input[p+1] - input[p]));

dy[p] = float(0.5*(input[p+nx] - input[p-nx]));

int k = (i+1) * nx - 1;

dx[k] = float(0.5*(input[k] - input[k-1]));

dy[k] = float(0.5*(input[k+nx] - input[k-nx]));

}

// 计算四个角

dx[0] = float(0.5*(input[1] - input[0]));

dy[0] = float(0.5*(input[nx] - input[0]));

dx[nx-1] = float(0.5*(input[nx-1] - input[nx-2]));

dy[nx-1] = float(0.5*(input[2*nx-1] - input[nx-1]));

dx[(ny-1)*nx] = float(0.5*(input[(ny-1)*nx + 1] - input[(ny-1)*nx]));

dy[(ny-1)*nx] = float(0.5*(input[(ny-1)*nx] - input[(ny-2)*nx]));

dx[ny*nx-1] = float(0.5*(input[ny*nx-1] - input[ny*nx-1-1]));

dy[ny*nx-1] = float(0.5*(input[ny*nx-1] - input[(ny-1)*nx-1]));

}

for(w=0;w<Nwarps;w++)
{
2.双三线性插值得到流场u0对I1的变形
bicubic_interpolation_warp(I1,  u1, u2, I1w,  nx, ny, 1);//对I1插值
bicubic_interpolation_warp(I1x, u1, u2, I1wx, nx, ny, 1);//对I1的x方向导数插值
bicubic_interpolation_warp(I1y, u1, u2, I1wy, nx, ny, 1);//y方向导数插值
 3.求I1(x+u0(x))的梯度
grad[i]=I1wx[i] * I1wx[i]+I1wy[i] * I1wy[i]

{
4.求TH

5.计算p1和p2散度

divergence(p11, p12, div_p1, nx ,ny);

divergence(p21, p22, div_p2, nx ,ny); 
6.计算u的方向导数(文章说利用向前差分)，跟求I1的方向导数类似
 

forward_gradient(u1, u1x, u1y, nx ,ny);

forward_gradient(u2, u2x, u2y, nx ,ny);

7.求p

for (int i = 0; i < size; i++)

{

float taut = float(tau / theta);

float g1 = float(hypot(u1x[i], u1y[i]));

float g2 = float(hypot(u2x[i], u2y[i]));

float ng1 = float(1.0 + taut * g1);

float ng2  = float(1.0 + taut * g2);

p11[i] = (p11[i] + taut * u1x[i]) / ng1;

p12[i] = (p12[i] + taut * u1y[i]) / ng1;

p21[i] = (p21[i] + taut * u2x[i]) / ng2;

p22[i] = (p22[i] + taut * u2y[i]) / ng2;

} 

}
} 

//求散度
 

void divergence(float *v1,float *v2,float *div,int nx,int ny)

{

for (int i = 1; i < ny-1; i++)

{

for(int j = 1; j < nx-1; j++)

{

int p  = i * nx + j;

int p1 = p - 1;

int p2 = p - nx;

float v1x = v1[p] - v1[p1];

float v2y = v2[p] - v2[p2];

div[p] = v1x + v2y;

}

}

for (int j = 1; j < nx-1; j++)

{

int p = (ny-1) * nx + j;

div[j] = v1[j] - v1[j-1] + v2[j];

div[p] = v1[p] - v1[p-1] - v2[p-nx];

}

for (int i = 1; i < ny-1; i++)

{

int p1 = i * nx;

int p2 = (i+1) * nx - 1;

div[p1] =  v1[p1]   + v2[p1] - v2[p1 - nx];

div[p2] = -v1[p2-1] + v2[p2] - v2[p2 - nx];

}

div[0]  =  v1[0] + v2[0];

div[nx-1]  = -v1[nx - 2] + v2[nx - 1];

div[(ny-1)*nx] =  v1[(ny-1)*nx] - v2[(ny-2)*nx];

div[ny*nx-1] = -v1[ny*nx - 2] - v2[(ny-1)*nx - 1];

}