Horn-Schunck光流法
關於光流法請看我以前的博客Lukas-Kanade光流法。這裏介紹Horn-Schunck光流法。算法
Horn-Schunck光流法求得的是稠密光流，須要對每個像素計算光流值，計算量比較大。而Lucas-Kanade光流法只需計算若干點的光流，是一種稀疏光流。
函數
數學原理這裏就不介紹了，直接說算法步驟。性能
用uij與vij分別表示圖像像素點(i，j)處的水平方向光流值與垂直方向光流值，每次迭代後的更新方程爲
ui
n爲迭代次數，lamda反映了對圖像數據及平滑約束的可信度，當圖像數據自己含有較大噪聲時，此時須要加大lamda的值，相反，當輸入圖像含有較少的噪聲時，此時可減少lamda的值。this
表明u鄰域與v鄰域的平均值，通常採用相應4鄰域內的均值
spa
也能夠採用3*三、5*5的窗口用模板平滑，窗口不宜過大，過大會破壞光流假設。.net
Ix、Iy分別是圖像對x、y的偏導數。It是兩幀圖像間對時間的導數。code
固然你也能夠考慮相鄰像素及相鄰兩幀圖像的影響，Horn-Schunck 提出經過 4 次有限差分來獲得blog
這裏只考慮了先後兩幀圖像。考慮連續三幀圖像的話有以下方法：ci
一種性能更優的 3D-Sobel 算子 以下圖所示。該算子在x 、y 、t方向上分別使用不一樣的模板對連續3幀圖像進行卷積計算 得出中間幀的位於模板中心的像素在三個方向上的梯度 。
迭代必定次數後u、v收斂，光流計算中止。在實際的計算中迭代初值可取U(0) =0、V(0)=0。
算法改進
對於通常場景基本等式只有在圖像中灰度梯度值較大的點處才成立。所以爲了加強算法的穩定性和準確性 咱們僅在梯度較大的點處才使用亮度恆常性約束，而在梯度較小的點處只使用流場一致性約束。定義以下權函數
下面是個人實現，使用了圖像金字塔，關於圖像金字塔，請看Lukas-Kanade光流法。（寫代碼時傳錯一個參數，調了幾個小時）
#ifndef __HORNSCHUNCK__
#define __HORNSCHUNCK__
class HornSchunckTracker
{
private:
unsigned int max_pyramid_layer;
unsigned int original_imgH;
unsigned int original_imgW;
BYTE**pre_pyr;//the pyramid of previous frame image,img1_pyr[0] is of max size
BYTE**next_pyr;//the frame after img1_pyr
int*height;
int*width;
double*optical_field_U;
double*optical_field_V;
bool isusepyramid;
double lamda;//取20
const double precision = 1;
const int maxiteration=300;
double threshold;//最小的光流閾值
double scale_factor;//縮放因子
private:
void get_max_pyramid_layer();
void pyramid_down(BYTE*&src_gray_data, const int src_h,
const int src_w, BYTE*& dst, int&dst_h, int&dst_w);
void pyramid_up(double*src,int srcW,int srcH,double*dst,int dstW,int dstH);
void lowpass_filter(double*&src, const int H, const int W, double*&smoothed);
void get_fx_fy_ft(BYTE*img1, BYTE*img2, int w, int h, double*fx, double*fy, double*ft);
void build_pyramid(BYTE**&original_gray);
double get_average4(double*src, const int height, const int width, const int i, const int j);
void bilinear(double* lpSrc, double* lpDst, int nW, int nH, int H1, int W1);
void bilinear(BYTE* lpSrc, BYTE* lpDst, int nW, int nH, int H1, int W1);
public:
HornSchunckTracker();
~HornSchunckTracker();
void get_pre_frame(BYTE*&gray);//use only at the beginning
void discard_pre_frame();
//set the next frame as pre_frame,must dicard pre_pyr in advance
void get_pre_frame();
//use every time,must after using get_pre_frame(BYTE**pyr)
void get_next_frame(BYTE*&gray);
void get_info(const int nh, const int nw);
void set_paras(double lamda,double threshold,double scalefactor);
void run_single_frame();
void HornSchunck();
void get_optical_flow(double*&u, double*&v);
};
#endif
#include "stdafx.h"
#include "HornSchunckTracker.h"
HornSchunckTracker::HornSchunckTracker()
{
isusepyramid = true;
}
HornSchunckTracker::~HornSchunckTracker()
{
for (int i = 0; i < max_pyramid_layer; i++)
{
if (pre_pyr[i])
delete[]pre_pyr[i];
if (next_pyr[i])
delete[]next_pyr[i];
}
delete[]pre_pyr;
delete[]next_pyr;
if (height)
delete[]height;
if (width)
delete[]width;
}
void HornSchunckTracker::get_max_pyramid_layer()
{
double minsize = 80;
double temp = original_imgH > original_imgW ?
original_imgW : original_imgH;
double tt = log(temp / minsize) / log(scale_factor);
if (tt>4)
{
max_pyramid_layer = 5;
return;
}
max_pyramid_layer = tt;
}
void HornSchunckTracker::build_pyramid(BYTE**&pyramid)
{
for (int i = 1; i < max_pyramid_layer; i++)
{
pyramid_down(pyramid[i - 1], height[i - 1],
width[i - 1], pyramid[i], height[i], width[i]);
}
}
void HornSchunckTracker::pyramid_down(BYTE*&src_gray_data,
const int src_h, const int src_w, BYTE*& dst, int&dst_h, int&dst_w)
{
dst_h = src_h / scale_factor;
dst_w = src_w / scale_factor;
assert(dst_w > 3 && dst_h > 3);
//BYTE*smoothed = new BYTE[src_h*src_w];
dst = new BYTE[dst_h*dst_w];
//lowpass_filter(src_gray_data, src_h, src_w,smoothed);
bilinear(src_gray_data, dst, src_w, src_h, dst_h, dst_w);
/*for (int i = 0; i < dst_h - 1; i++)
for (int j = 0; j < dst_w - 1; j++)
{
int srcY = 2 * i + 1;
int srcX = 2 * j + 1;
double re = src_gray_data[srcY*src_w + srcX] * 0.25;
re += src_gray_data[(srcY - 1)*src_w + srcX] * 0.125;
re += src_gray_data[(srcY + 1)*src_w + srcX] * 0.125;
re += src_gray_data[srcY*src_w + srcX - 1] * 0.125;
re += src_gray_data[srcY*src_w + srcX + 1] * 0.125;
re += src_gray_data[(srcY - 1)*src_w + srcX + 1] * 0.0625;
re += src_gray_data[(srcY - 1)*src_w + srcX - 1] * 0.0625;
re += src_gray_data[(srcY + 1)*src_w + srcX - 1] * 0.0625;
re += src_gray_data[(srcY + 1)*src_w + srcX + 1] * 0.0625;
dst[i*dst_w + j] = re;
}
for (int i = 0; i < dst_h; i++)
dst[i*dst_w + dst_w - 1] = dst[i*dst_w + dst_w - 2];
for (int i = 0; i < dst_w; i++)
dst[(dst_h - 1)*dst_w + i] = dst[(dst_h - 2)*dst_w + i];*/
}
void HornSchunckTracker::get_pre_frame(BYTE*&gray)//use only at the beginning
{
pre_pyr[0] = gray;
build_pyramid(pre_pyr);
//save_gray("1.bmp", pre_pyr[1], height[1], width[1]);
}
void HornSchunckTracker::discard_pre_frame()
{
for (int i = 0; i < max_pyramid_layer; i++)
delete[]pre_pyr[i];
}
//set the next frame as pre_frame,must dicard pre_pyr in advance
void HornSchunckTracker::get_pre_frame()
{
for (int i = 0; i < max_pyramid_layer; i++)
pre_pyr[i] = next_pyr[i];
}
//use every time,must after using get_pre_frame(BYTE**pyr)
void HornSchunckTracker::get_next_frame(BYTE*&gray)
{
next_pyr[0] = gray;
build_pyramid(next_pyr);
//save_gray("1.bmp", next_pyr[1], height[1], width[1]);
}
void HornSchunckTracker::get_info(const int nh, const int nw)
{
original_imgH = nh;
original_imgW = nw;
if (isusepyramid)
get_max_pyramid_layer();
else
max_pyramid_layer = 1;
pre_pyr = new BYTE*[max_pyramid_layer];
next_pyr = new BYTE*[max_pyramid_layer];
height = new int[max_pyramid_layer];
width = new int[max_pyramid_layer];
height[0] = nh;
width[0] = nw;
}
//低通濾波
void HornSchunckTracker::lowpass_filter(double*&src, const int H, const int W, double*&smoothed)
{
//tackle with border
for (int i = 0; i < H; i++)
{
smoothed[i*W] = src[i*W];
smoothed[i*W + W - 1] = src[i*W + W - 1];
}
for (int i = 0; i < W; i++)
{
smoothed[i] = src[i];
smoothed[(H - 1)*W + i] = src[(H - 1)*W + i];
}
for (int i = 1; i < H - 1; i++)
for (int j = 1; j < W - 1; j++)
{
double re = 0;
re += src[i*W + j] * 0.25;
re += src[(i - 1)*W + j] * 0.125;
re += src[i*W + j + 1] * 0.125;
re += src[i*W + j - 1] * 0.125;
re += src[(i + 1)*W + j] * 0.125;
re += src[(i - 1)*W + j - 1] * 0.0625;
re += src[(i + 1)*W + j - 1] * 0.0625;
re += src[(i - 1)*W + j + 1] * 0.0625;
re += src[(i + 1)*W + j + 1] * 0.0625;
smoothed[i*W + j] = re;
}
}
void HornSchunckTracker::get_fx_fy_ft(BYTE*img1, BYTE*img2, int w, int h, double*fx, double*fy, double*ft)
{
for (int i = 0; i < h - 1; i++)
for (int j = 0; j < w - 1; j++)
{
fx[i*w + j] = 0.25*(img1[i*w + j + 1] - img1[i*w + j] + img1[(i + 1)*w + j + 1] - img1[(i + 1)*w + j]
+ img2[i*w + j + 1] - img2[i*w + j] + img2[(i + 1)*w + j + 1] - img2[(i + 1)*w + j]);
fy[i*w + j] = 0.25 * (img1[(i + 1)*w + j] - img1[i*w + j] +img1[(i + 1)*w + j + 1] - img1[i*w + j + 1]
+ img2[(i + 1)*w + j] - img2[i*w + j] + img2[(i + 1)*w + j + 1] - img2[i*w + j + 1]);
ft[i*w + j] = 0.25 * (img2[i*w + j] - img1[i*w + j] +img2[(i + 1)*w + j] - img1[(i + 1)*w + j] +
img2[(i + 1)*w + j + 1] - img1[(i + 1)*w + j + 1] + img2[i*w + j + 1] - img1[i*w + j + 1]);
}
for (int i = 0; i < h; i++)
{
//fx[i*w] = fx[i*w + w - 2];
fx[i*w + w - 1] = fx[i*w + w - 2];
//fy[i*w] = fy[i*w + w - 2];
fy[i*w + w - 1] = fy[i*w + w - 2];
//ft[i*w] = ft[i*w + w - 2];
ft[i*w + w - 1] = ft[i*w + w - 2];
}
for (int j = 0; j < w; j++)
{
//fx[j] = fx[h + j];
fx[(h - 1)*w + j] = fx[(h - 2)*w + j];
//fy[j] = fy[h + j];
fy[(h - 1)*w + j] = fy[(h - 2)*w + j];
//ft[j] = ft[h + j];
ft[(h - 1)*w + j] = ft[(h - 2)*w + j];
}
}
//取得計算獲得的光流值，u、v爲out型參數
void HornSchunckTracker::get_optical_flow(double*&u, double*&v)
{
assert(optical_field_U&&optical_field_V);
u = optical_field_U;
v = optical_field_V;
}
//int save_gray(char * savebmp_file, LPBYTE gray, int height, int width);
//返回求得的光流場，大小爲原始圖像大小
void HornSchunckTracker::HornSchunck()
{
//save_gray("22.bmp", pre_pyr[0], height[0], width[0]);
//初始化光流場爲0
if (optical_field_U)
delete[]optical_field_U;
if (optical_field_V)
delete[]optical_field_V;
optical_field_U = new double[width[max_pyramid_layer - 1]
* height[max_pyramid_layer - 1]];
optical_field_V = new double[width[max_pyramid_layer - 1]
* height[max_pyramid_layer - 1]];
memset(optical_field_U, 0, sizeof(double)*width[max_pyramid_layer - 1]
* height[max_pyramid_layer - 1]);
memset(optical_field_V, 0, sizeof(double)*width[max_pyramid_layer - 1]
* height[max_pyramid_layer - 1]);
//使用金字塔計算光流
for (int i = max_pyramid_layer - 1; i >= 0; i--)
{
double*Ix = new double[width[i] * height[i]];
double*Iy = new double[width[i] * height[i]];
double*It = new double[width[i] * height[i]];
//求偏導
get_fx_fy_ft(pre_pyr[i], next_pyr[i], width[i], height[i], Ix, Iy, It);
//將光流場平滑
double*smoothed_U = new double[width[i] * height[i]];
double*smoothed_V = new double[width[i] * height[i]];
if (i == max_pyramid_layer - 1)
{
memset(smoothed_U, 0, sizeof(double)*width[i] * height[i]);
memset(smoothed_V, 0, sizeof(double)*width[i] * height[i]);
}
else
{
lowpass_filter(optical_field_U, height[i], width[i], smoothed_U);
lowpass_filter(optical_field_V, height[i], width[i], smoothed_V);
}
double error = 1000000;
int iteration = 0;
//迭代計算每一個像素的光流，直到收斂或達到最大迭代次數
while (error > precision&&iteration < maxiteration)
{
iteration++;
error = 0;
//計算該層金字塔的光流
for (int j = 0; j < height[i]; j++)
for (int k = 0; k < width[i]; k++)
{
//採用改進方法，光流速度需大於閾值，這樣不只準確度增長，計算量也會減少
double w = Ix[j*width[i] + k] * Ix[j*width[i] + k]
+ Iy[j*width[i] + k] * Iy[j*width[i] + k] > threshold ? 1 : 0;
double u_pre = optical_field_U[j*width[i] + k];
double v_pre = optical_field_V[j*width[i] + k];
double utemp = smoothed_U[j*width[i] + k];//get_average4(optical_field_U, height[i], width[i], j, k);
double vtemp = smoothed_V[j*width[i] + k]; //get_average4(optical_field_V, height[i], width[i], j, k);
double denominator = lamda + w*(Ix[j*width[i] + k] * Ix[j*width[i] + k]
+ Iy[j*width[i] + k] * Iy[j*width[i] + k]);
double numerator = Ix[j*width[i] + k] * utemp + Iy[j*width[i] + k] *
vtemp + It[j*width[i] + k];
optical_field_U[j*width[i] + k] = utemp - Ix[j*width[i] + k] *w*numerator / denominator;
optical_field_V[j*width[i] + k] = vtemp - Iy[j*width[i] + k] *w*numerator / denominator;
error += pow(optical_field_U[j*width[i] + k] - u_pre,2) +
pow(optical_field_V[j*width[i] + k] - v_pre,2);
}
//下一次迭代前從新平滑光流場
if (error >exp(double(max_pyramid_layer-i))*precision&&iteration < maxiteration)
{
lowpass_filter(optical_field_U, height[i], width[i], smoothed_U);
lowpass_filter(optical_field_V, height[i], width[i], smoothed_V);
}
}
delete[]smoothed_U, smoothed_V,Ix,Iy,It;
if (i == 0)//獲得最終光流場
{
return;
}
//下一層的光流場
double*new_of_u = new double[width[i - 1] * height[i - 1]];
double*new_of_v = new double[width[i - 1] * height[i - 1]];
//上採樣
pyramid_up(optical_field_U, width[i], height[i], new_of_u, width[i - 1], height[i - 1]);
pyramid_up(optical_field_V, width[i], height[i], new_of_v, width[i - 1], height[i - 1]);
//將每一個像素的光流按縮放因子放大，獲得下一層的光流場的初值
//double scale = double(height[i - 1]) / height[i];
for (int j = 0; j < height[i - 1]; j++)
for (int k = 0; k < width[i - 1]; k++)
{
new_of_u[j*width[i - 1] + k] *= scale_factor;
new_of_v[j*width[i - 1] + k] *= scale_factor;
}
delete[]optical_field_U, optical_field_V;
optical_field_U = new_of_u;
optical_field_V = new_of_v;
}
}
//上採樣，採用雙線性插值，用雙立方插值應該更精確
void HornSchunckTracker::pyramid_up(double*src, int srcW, int srcH, double*dst, int dstW, int dstH)
{
bilinear(src, dst, srcW, srcH, dstH, dstW);
}
//雙線性插值
void HornSchunckTracker::bilinear(double* lpSrc, double* lpDst, int nW, int nH, int H1, int W1)
{
float fw = float(nW) / W1;
float fh = float(nH) / H1;
int y1, y2, x1, x2, x0, y0;
float fx1, fx2, fy1, fy2;
for (int i = 0; i < H1; i++)
{
y0 = i*fh;
y1 = int(y0);
if (y1 == nH - 1) y2 = y1;
else y2 = y1 + 1;
fy1 = y1 - y0;
fy2 = 1.0f - fy1;
for (int j = 0; j < W1; j++)
{
x0 = j*fw;
x1 = int(x0);
if (x1 == nW - 1) x2 = x1;
else x2 = x1 + 1;
fx1 = y1 - y0;
fx2 = 1.0f - fx1;
float s1 = fx1*fy1;
float s2 = fx2*fy1;
float s3 = fx2*fy2;
float s4 = fx1*fy2;
double c1, c2, c3, c4;
c1 = lpSrc[y1*nW + x1];
c2 = lpSrc[y1*nW + x2];
c3 = lpSrc[y2*nW + x1];
c4 = lpSrc[y2*nW + x2];
double r;
r = (c1*s3) + (c2*s4) + (c3*s2) + (c4*s1);
lpDst[i*W1 + j] = r;
}
}
}
//雙線性插值
void HornSchunckTracker::bilinear(BYTE* lpSrc, BYTE* lpDst, int nW, int nH, int H1, int W1)
{
float fw = float(nW) / W1;
float fh = float(nH) / H1;
int y1, y2, x1, x2, x0, y0;
float fx1, fx2, fy1, fy2;
for (int i = 0; i < H1; i++)
{
y0 = i*fh;
y1 = int(y0);
if (y1 == nH - 1) y2 = y1;
else y2 = y1 + 1;
fy1 = y1 - y0;
fy2 = 1.0f - fy1;
for (int j = 0; j < W1; j++)
{
x0 = j*fw;
x1 = int(x0);
if (x1 == nW - 1) x2 = x1;
else x2 = x1 + 1;
fx1 = y1 - y0;
fx2 = 1.0f - fx1;
float s1 = fx1*fy1;
float s2 = fx2*fy1;
float s3 = fx2*fy2;
float s4 = fx1*fy2;
double c1, c2, c3, c4;
c1 = lpSrc[y1*nW + x1];
c2 = lpSrc[y1*nW + x2];
c3 = lpSrc[y2*nW + x1];
c4 = lpSrc[y2*nW + x2];
double r;
r = (c1*s3) + (c2*s4) + (c3*s2) + (c4*s1);
lpDst[i*W1 + j] = BYTE(r);
}
}
}
void HornSchunckTracker::set_paras(double lamda, double threshold, double scalefactor)
{
this->lamda = lamda;
this->threshold = threshold;
scale_factor = scalefactor;
}
//double HornSchunckTracker::get_average4(double*src, const int height, const int width, const int i, const int j)
//{
// if (j < 0 || j >= width) return 0;
// if (i < 0 || i >= height) return 0;
//
// double val = 0.0;
// int tmp = 0;
// if ((i - 1) >= 0){
// ++tmp;
// val += src[(i - 1)*width + j];
// }
// if (i + 1<height){
// ++tmp;
// val += src[(i + 1)*width + j];
// }
// if (j - 1 >= 0){
// ++tmp;
// val += src[i*width + j - 1];
// }
// if (j + 1<width){
// ++tmp;
// val += src[i*width + j + 1];
// }
// return val / tmp;
//
//}
能夠看出對邊緣的光流檢測較好，內部點的光流檢測較難。
相關文章
- 1. 光流算法
- 2. KLT光流法
- 3. 光流法
- 4. 光流法—光流場估計概念
- 5. Lukas-Kanade光流法
- 6. 光流法（optical flow）
- 7. opencv python 光流法
- 8. 經典光流計算方法（HS光流法、Lucas-Kanada方法、Pyramidal LK方法）
- 9. LK光流法小結
- 10. 光流法簡單介紹
- 更多相關文章...
- • XSL-FO 流 - XSL-FO 教程
- • Lua 流程控制 - Lua 教程
- • Docker容器實戰(七) - 容器眼光下的文件系統
- • 算法總結-回溯法
相關標籤/搜索
每日一句
-
每一个你不满意的现在，都有一个你没有努力的曾经。
歡迎關注本站公眾號,獲取更多信息
相關文章
- 1. 光流算法
- 2. KLT光流法
- 3. 光流法
- 4. 光流法—光流場估計概念
- 5. Lukas-Kanade光流法
- 6. 光流法（optical flow）
- 7. opencv python 光流法
- 8. 經典光流計算方法（HS光流法、Lucas-Kanada方法、Pyramidal LK方法）
- 9. LK光流法小結
- 10. 光流法簡單介紹