目标跟踪之Horn-Schunck光流法
关于光流法请看我之前的博客Lukas-Kanade光流法。这里介绍Horn-Schunck光流法。
Horn-Schunck光流法求得的是稠密光流，需要对每一个像素计算光流值，计算量比较大。而Lucas-Kanade光流法只需计算若干点的光流，是一种稀疏光流。
数学原理这里就不介绍了，直接说算法步骤。
用uij与vij分别表示图像像素点(i，j)处的水平方向光流值与垂直方向光流值，每次迭代后的更新方程为
n为迭代次数，lamda反映了对图像数据及平滑约束的可信度，当图像数据本身含有较大噪声时，此时需要加大lamda的值，相反，当输入图像含有较少的噪声时，此时可减小lamda的值。
代表u邻域与v邻域的平均值，一般采用相应4邻域内的均值
也可以采用33、55的窗口用模板平滑，窗口不宜过大，过大会破坏光流假设。
Ix、Iy分别是图像对x、y的偏导数。It是两帧图像间对时间的导数。
当然你也可以考虑相邻像素及相邻两帧图像的影响，Horn-Schunck提出通过4次有限差分来得到
这里只考虑了前后两帧图像。考虑连续三帧图像的话有如下方法：
一种性能更优的3D-Sobel算子如下图所示。该算子在x、y、t方向上分别使用不同的模板对连续3帧图像进行卷积计算得出中间帧的位于模板中心的像素在三个方向上的梯度。
迭代一定次数后u、v收敛，光流计算停止。在实际的计算中迭代初值可取U(0)=0、V(0)=0。
算法改进
对于一般场景基本等式只有在图像中灰度梯度值较大的点处才成立。因此为了增强算法的稳定性和准确性我们仅在梯度较大的点处才使用亮度恒常性约束，而在梯度较小的点处只使用流场一致性约束。定义如下权函数
下面是我的实现，使用了图像金字塔，关于图像金字塔，请看Lukas-Kanade光流法。（写代码时传错一个参数，调了几个小时哭）
[cpp]viewplaincopy
#ifndef__HORNSCHUNCK__
#define__HORNSCHUNCK__
classHornSchunckTracker
{
private:
unsignedintmax_pyramid_layer;
unsignedintoriginal_imgH;
unsignedintoriginal_imgW;
BYTEpre_pyr;//thepyramidofpreviousframeimage,img1_pyr[0]isofmaxsize
BYTEnext_pyr;//theframeafterimg1_pyr
intheight;
intwidth;
doubleoptical_field_U;
doubleoptical_field_V;
boolisusepyramid;
doublelamda;//取20
constdoubleprecision=1;
constintmaxiteration=300;
doublethreshold;//最小的光流阈值
doublescale_factor;//缩放因子
private:
voidget_max_pyramid_layer();
voidpyramid_down(BYTE&src_gray_data,constintsrc_h,
constintsrc_w,BYTE&dst,int&dst_h,int&dst_w);
voidpyramid_up(doublesrc,intsrcW,intsrcH,doubledst,intdstW,intdstH);
voidlowpass_filter(double&src,constintH,constintW,double&smoothed);
voidget_fx_fy_ft(BYTEimg1,BYTEimg2,intw,inth,doublefx,doublefy,doubleft);
voidbuild_pyramid(BYTE&original_gray);
doubleget_average4(doublesrc,constintheight,constintwidth,constinti,constintj);
voidbilinear(doublelpSrc,doublelpDst,intnW,intnH,intH1,intW1);
voidbilinear(BYTElpSrc,BYTElpDst,intnW,intnH,intH1,intW1);
public:
HornSchunckTracker();
~HornSchunckTracker();
voidget_pre_frame(BYTE&gray);//useonlyatthebeginning
voiddiscard_pre_frame();
//setthenextframeaspre_frame,mustdicardpre_pyrinadvance
voidget_pre_frame();
//useeverytime,mustafterusingget_pre_frame(BYTEpyr)
voidget_next_frame(BYTE&gray);
voidget_info(constintnh,constintnw);
voidset_paras(doublelamda,doublethreshold,doublescalefactor);
voidrun_single_frame();
voidHornwww.baiyuewang.netSchunck();
voidget_optical_flow(double&u,double&v);
};
#endif
[cpp]viewplaincopy
#include"stdafx.h"
#include"HornSchunckTracker.h"
HornSchunckTracker::HornSchunckTracker()
{
isusepyramid=true;
}
HornSchunckTracker::~HornSchunckTracker()
{
for(inti=0;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;
}
voidHornSchunckTracker::get_max_pyramid_layer()
{
doubleminsize=80;
doubletemp=original_imgH>original_imgW?
original_imgW:original_imgH;
doublett=log(temp/minsize)/log(scale_factor);
if(tt>4)
{
max_pyramid_layer=5;
return;
}
max_pyramid_layer=tt;
}
voidHornSchunckTracker::build_pyramid(BYTE&pyramid)
{
for(inti=1;i
{
pyramid_down(pyramid[i-1],height[i-1],
width[i-1],pyramid[i],height[i],width[i]);
}
}
voidHornSchunckTracker::pyramid_down(BYTE&src_gray_data,
constintsrc_h,constintsrc_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);
//BYTEsmoothed=newBYTE[src_hsrc_w];
dst=newBYTE[dst_hdst_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(inti=0;i
for(intj=0;j
{
intsrcY=2i+1;
intsrcX=2j+1;
doublere=src_gray_data[srcYsrc_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[srcYsrc_w+srcX-1]0.125;
re+=src_gray_data[srcYsrc_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[idst_w+j]=re;
}
for(inti=0;i
dst[idst_w+dst_w-1]=dst[idst_w+dst_w-2];
for(inti=0;i
dst[(dst_h-1)dst_w+i]=dst[(dst_h-2)dst_w+i];/
}
voidHornSchunckTracker::get_pre_frame(BYTE&gray)//useonlyatthebeginning
{
pre_pyr[0]=gray;
build_pyramid(pre_pyr);
//save_gray("1.bmp",pre_pyr[1],height[1],width[1]);
}
voidHornSchunckTracker::discard_pre_frame()
{
for(inti=0;i
delete[]pre_pyr[i];
}
//setthenextframeaspre_frame,mustdicardpre_pyrinadvance
voidHornSchunckTracker::get_pre_frame()
{
for(inti=0;i
pre_pyr[i]=next_pyr[i];
}
//useeverytime,mustafterusingget_pre_frame(BYTEpyr)
voidHornSchunckTracker::get_next_frame(BYTE&gray)
{
next_pyr[0]=gray;
build_pyramid(next_pyr);
//save_gray("1.bmp",next_pyr[1],height[1],width[1]);
}
voidHornSchunckTracker::get_info(constintnh,constintnw)
{
original_imgH=nh;
original_imgW=nw;
if(isusepyramid)
get_max_pyramid_layer();
else
max_pyramid_layer=1;
pre_pyr=newBYTE[max_pyramid_layer];
next_pyr=newBYTE[max_pyramid_layer];
height=newint[max_pyramid_layer];
width=newint[max_pyramid_layer];
height[0]=nh;
width[0]=nw;
}
//低通滤波
voidHornSchunckTracker::lowpass_filter(double&src,constintH,constintW,double&smoothed)
{
//tacklewithborder
for(inti=0;i
{
smoothed[iW]=src[iW];
smoothed[iW+W-1]=src[iW+W-1];
}
for(inti=0;i
{
smoothed[i]=src[i];
smoothed[(H-1)W+i]=src[(H-1)W+i];
}
for(inti=1;i
for(intj=1;j
{
doublere=0;
re+=src[iW+j]0.25;
re+=src[(i-1)W+j]0.125;
re+=src[iW+j+1]0.125;
re+=src[iW+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[iW+j]=re;
}
}
voidHornSchunckTracker::get_fx_fy_ft(BYTEimg1,BYTEimg2,intw,inth,doublefx,doublefy,doubleft)
{
for(inti=0;i
for(intj=0;j
{
fx[iw+j]=0.25(img1[iw+j+1]-img1[iw+j]+img1[(i+1)w+j+1]-img1[(i+1)w+j]
+img2[iw+j+1]-img2[iw+j]+img2[(i+1)w+j+1]-img2[(i+1)w+j]);
fy[iw+j]=0.25(img1[(i+1)w+j]-img1[iw+j]+img1[(i+1)w+j+1]-img1[iw+j+1]
+img2[(i+1)w+j]-img2[iw+j]+img2[(i+1)w+j+1]-img2[iw+j+1]);
ft[iw+j]=0.25(img2[iw+j]-img1[iw+j]+img2[(i+1)w+j]-img1[(i+1)w+j]+
img2[(i+1)w+j+1]-img1[(i+1)w+j+1]+img2[iw+j+1]-img1[iw+j+1]);
}
for(inti=0;i
{
//fx[iw]=fx[iw+w-2];
fx[iw+w-1]=fx[iw+w-2];
//fy[iw]=fy[iw+w-2];
fy[iw+w-1]=fy[iw+w-2];
//ft[iw]=ft[iw+w-2];
ft[iw+w-1]=ft[iw+w-2];
}
for(intj=0;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型参数
voidHornSchunckTracker::get_optical_flow(double&u,double&v)
{
assert(optical_field_U&&optical_field_V);
u=optical_field_U;
v=optical_field_V;
}
//intsave_gray(charsavebmp_file,LPBYTEgray,intheight,intwidth);
//返回求得的光流场，大小为原始图像大小
voidHornSchunckTracker::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=newdouble[width[max_pyramid_layer-1]
height[max_pyramid_layer-1]];
optical_field_V=newdouble[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(inti=max_pyramid_layer-1;i>=0;i--)
{
doubleIx=newdouble[width[i]height[i]];
doubleIy=newdouble[width[i]height[i]];
doubleIt=newdouble[width[i]height[i]];
//求偏导
get_fx_fy_ft(pre_pyr[i],next_pyr[i],width[i],height[i],Ix,Iy,It);
//将光流场平滑
doublesmoothed_U=newdouble[width[i]height[i]];
doublesmoothed_V=newdouble[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);
}
doubleerror=1000000;
intiteration=0;
//迭代计算每个像素的光流，直到收敛或达到最大迭代次数
while(error>precision&&iteration
{
iteration++;
error=0;
//计算该层金字塔的光流
for(intj=0;j
for(intk=0;k
{
//采用改进方法，光流速度需大于阈值，这样不仅准确度增加，计算量也会减小
doublew=Ix[jwidth[i]+k]Ix[jwidth[i]+k]
+Iy[jwidth[i]+k]Iy[jwidth[i]+k]>threshold?1:0;
doubleu_pre=optical_field_U[jwidth[i]+k];
doublev_pre=optical_field_V[jwidth[i]+k];
doubleutemp=smoothed_U[jwidth[i]+k];//get_average4(optical_field_U,height[i],width[i],j,k);
doublevtemp=smoowww.tt951.comthed_V[jwidth[i]+k];//get_average4(optical_field_V,height[i],width[i],j,k);
doubledenominator=lamda+w(Ix[jwidth[i]+k]Ix[jwidth[i]+k]
+Iy[jwidth[i]+k]Iy[jwidth[i]+k]);
doublenumerator=Ix[jwidth[i]+k]utemp+Iy[jwidth[i]+k]
vtemp+It[jwidth[i]+k];
optical_field_U[jwidth[i]+k]=utemp-Ix[jwidth[i]+k]wnumerator/denominator;
optical_field_V[jwidth[i]+k]=vtemp-Iy[jwidth[i]+k]wnumerator/denominator;
error+=pow(optical_field_U[jwidth[i]+k]-u_pre,2)+
pow(optical_field_V[jwidth[i]+k]-v_pre,2);
}
//下一次迭代前重新平滑光流场
if(error>exp(double(max_pyramid_layer-i))precision&&iteration
{
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;
}
//下一层的光流场
doublenew_of_u=newdouble[width[i-1]height[i-1]];
doublenew_of_v=newdouble[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]);
//将每个像素的光流按缩放因子放大，得到下一层的光流场的初值
//doublescale=double(height[i-1])/height[i];
for(intj=0;j
for(intk=0;k
{
new_of_u[jwidth[i-1]+k]=scale_factor;
new_of_v[jwidth[i-1]+k]=scale_factor;
}
delete[]optical_field_U,optical_field_V;
optical_field_U=new_of_u;
optical_field_V=new_of_v;
}
}
//上采样，采用双线性插值，用双立方插值应该更精确
voidHornSchunckTracker::pyramid_up(doublesrc,intsrcW,intsrcH,doubledst,intdstW,intdstH)
{
bilinear(src,dst,srcW,srcH,dstH,dstW);
}
//双线性插值
voidHornSchunckTracker::bilinear(doublelpSrc,doublelpDst,intnW,intnH,intH1,intW1)
{
floatfw=float(nW)/W1;
floatfh=float(nH)/H1;
inty1,y2,x1,x2,x0,y0;
floatfx1,fx2,fy1,fy2;
for(inti=0;i
{
y0=ifh;
y1=int(y0);
if(y1==nH-1)y2=y1;
elsey2=y1+1;
fy1=y1-y0;
fy2=1.0f-fy1;
for(intj=0;j
{
x0=jfw;
x1=int(x0);
if(x1==nW-1)x2=x1;
elsex2=x1+1;
fx1=y1-y0;
fx2=1.0f-fx1;
floats1=fx1fy1;
floats2=fx2fy1;
floats3=fx2fy2;
floats4=fx1fy2;
doublec1,c2,c3,c4;
c1=lpSrc[y1nW+x1];
c2=lpSrc[y1nW+x2];
c3=lpSrc[y2nW+x1];
c4=lpSrc[y2nW+x2];
doubler;
r=(c1s3)+(c2s4)+(c3s2)+(c4s1);
lpDst[iW1+j]=r;
}
}
}
//双线性插值
voidHornSchunckTracker::bilinear(BYTElpSrc,BYTElpDst,intnW,intnH,intH1,intW1)
{
floatfw=float(nW)/W1;
floatfh=float(nH)/H1;
inty1,y2,x1,x2,x0,y0;
floatfx1,fx2,fy1,fy2;
for(inti=0;i
{
y0=ifh;
y1=int(y0);
if(y1==nH-1)y2=y1;
elsey2=y1+1;
fy1=y1-y0;
fy2=1.0f-fy1;
for(intj=0;j
{
x0=jfw;
x1=int(x0);
if(x1==nW-1)x2=x1;
elsex2=x1+1;
fx1=y1-y0;
fx2=1.0f-fx1;
floats1=fx1fy1;
floats2=fx2fy1;
floats3=fx2fy2;
floats4=fx1fy2;
doublec1,c2,c3,c4;
c1=lpSrc[y1nW+x1];
c2=lpSrc[y1nW+x2];
c3=lpSrc[y2nW+x1];
c4=lpSrc[y2nW+x2];
doubler;
r=(c1s3)+(c2s4)+(c3s2)+(c4s1);
lpDst[iW1+j]=BYTE(r);
}
}
}
voidHornSchunckTracker::set_paras(doublelamda,doublethreshold,doublescalefactor)
{
this->lamda=lamda;
this->threshold=threshold;
scale_factor=scalefactor;
}
//doubleHornSchunckTracker::get_average4(doublesrc,constintheight,constintwidth,constinti,constintj)
//{
//if(j<0||j>=width)return0;
//if(i<0||i>=height)return0;
//
//doubleval=0.0;
//inttmp=0;
//if((i-1)>=0){
//++tmp;
//val+=src[(i-1)width+j];
//}
//if(i+1
//++tmp;
//val+=src[(i+1)width+j];
//}
//if(j-1>=0){
//++tmp;
//val+=src[iwidth+j-1];
//}
//if(j+1
//++tmp;
//val+=src[iwidth+j+1];
//}
//returnval/tmp;
//
//}
|
|
