目标跟踪之Horn-Schunck光流法

缺乏、安全感 2022-08-09 14:59 226阅读 0赞

关于光流法请看我之前的博客Lukas-Kanade光流法。这里介绍Horn-Schunck光流法。

Horn-Schunck光流法求得的是稠密光流,需要对每一个像素计算光流值,计算量比较大。而Lucas-Kanade光流法只需计算若干点的光流,是一种稀疏光流。

数学原理这里就不介绍了,直接说算法步骤。

用uij与vij分别表示图像像素点(i,j)处的水平方向光流值与垂直方向光流值,每次迭代后的更新方程为

20151030145919176

n为迭代次数,lamda反映了对图像数据及平滑约束的可信度,当图像数据本身含有较大噪声时,此时需要加大lamda的值,相反,当输入图像含有较少的噪声时,此时可减小lamda的值。

代表u邻域与v邻域的平均值,一般采用相应4邻域内的均值

20151029233433165

也可以采用3*3、5*5的窗口用模板平滑,窗口不宜过大,过大会破坏光流假设。

Ix、Iy分别是图像对x、y的偏导数。It是两帧图像间对时间的导数。

20151030104632199

当然你也可以考虑相邻像素及相邻两帧图像的影响,Horn-Schunck 提出通过 4 次有限差分来得到

20151030113610321

这里只考虑了前后两帧图像。考虑连续三帧图像的话有如下方法:

一种性能更优的 3D-Sobel 算子 如下图所示。该算子在x 、y 、t方向上分别使用不同的模板对连续3帧图像进行卷积计算 得出中间帧的位于模板中心的像素在三个方向上的梯度 。

20151030114019786

迭代一定次数后u、v收敛,光流计算停止。在实际的计算中迭代初值可取U(0) =0、V(0)=0。

算法改进

对于一般场景基本等式只有在图像中灰度梯度值较大的点处才成立。因此为了增强算法的稳定性和准确性 我们仅在梯度较大的点处才使用亮度恒常性约束,而在梯度较小的点处只使用流场一致性约束。定义如下权函数

20151030114707342

下面是我的实现,使用了图像金字塔,关于图像金字塔,请看Lukas-Kanade光流法。(写代码时传错一个参数,调了几个小时哭

  1. #ifndef __HORNSCHUNCK__
  2. #define __HORNSCHUNCK__
  3. class HornSchunckTracker
  4. {
  5. private:
  6. unsigned int max_pyramid_layer;
  7. unsigned int original_imgH;
  8. unsigned int original_imgW;
  9. BYTE**pre_pyr;//the pyramid of previous frame image,img1_pyr[0] is of max size
  10. BYTE**next_pyr;//the frame after img1_pyr
  11. int*height;
  12. int*width;
  13. double*optical_field_U;
  14. double*optical_field_V;
  15. bool isusepyramid;
  16. double lamda;//取20
  17. const double precision = 1;
  18. const int maxiteration=300;
  19. double threshold;//最小的光流阈值
  20. double scale_factor;//缩放因子
  21. private:
  22. void get_max_pyramid_layer();
  23. void pyramid_down(BYTE*&src_gray_data, const int src_h,
  24. const int src_w, BYTE*& dst, int&dst_h, int&dst_w);
  25. void pyramid_up(double*src,int srcW,int srcH,double*dst,int dstW,int dstH);
  26. void lowpass_filter(double*&src, const int H, const int W, double*&smoothed);
  27. void get_fx_fy_ft(BYTE*img1, BYTE*img2, int w, int h, double*fx, double*fy, double*ft);
  28. void build_pyramid(BYTE**&original_gray);
  29. double get_average4(double*src, const int height, const int width, const int i, const int j);
  30. void bilinear(double* lpSrc, double* lpDst, int nW, int nH, int H1, int W1);
  31. void bilinear(BYTE* lpSrc, BYTE* lpDst, int nW, int nH, int H1, int W1);
  32. public:
  33. HornSchunckTracker();
  34. ~HornSchunckTracker();
  35. void get_pre_frame(BYTE*&gray);//use only at the beginning
  36. void discard_pre_frame();
  37. //set the next frame as pre_frame,must dicard pre_pyr in advance
  38. void get_pre_frame();
  39. //use every time,must after using get_pre_frame(BYTE**pyr)
  40. void get_next_frame(BYTE*&gray);
  41. void get_info(const int nh, const int nw);
  42. void set_paras(double lamda,double threshold,double scalefactor);
  43. void run_single_frame();
  44. void HornSchunck();
  45. void get_optical_flow(double*&u, double*&v);
  46. };
  47. #endif
  48. #include "stdafx.h"
  49. #include "HornSchunckTracker.h"
  50. HornSchunckTracker::HornSchunckTracker()
  51. {
  52. isusepyramid = true;
  53. }
  54. HornSchunckTracker::~HornSchunckTracker()
  55. {
  56. for (int i = 0; i < max_pyramid_layer; i++)
  57. {
  58. if (pre_pyr[i])
  59. delete[]pre_pyr[i];
  60. if (next_pyr[i])
  61. delete[]next_pyr[i];
  62. }
  63. delete[]pre_pyr;
  64. delete[]next_pyr;
  65. if (height)
  66. delete[]height;
  67. if (width)
  68. delete[]width;
  69. }
  70. void HornSchunckTracker::get_max_pyramid_layer()
  71. {
  72. double minsize = 80;
  73. double temp = original_imgH > original_imgW ?
  74. original_imgW : original_imgH;
  75. double tt = log(temp / minsize) / log(scale_factor);
  76. if (tt>4)
  77. {
  78. max_pyramid_layer = 5;
  79. return;
  80. }
  81. max_pyramid_layer = tt;
  82. }
  83. void HornSchunckTracker::build_pyramid(BYTE**&pyramid)
  84. {
  85. for (int i = 1; i < max_pyramid_layer; i++)
  86. {
  87. pyramid_down(pyramid[i - 1], height[i - 1],
  88. width[i - 1], pyramid[i], height[i], width[i]);
  89. }
  90. }
  91. void HornSchunckTracker::pyramid_down(BYTE*&src_gray_data,
  92. const int src_h, const int src_w, BYTE*& dst, int&dst_h, int&dst_w)
  93. {
  94. dst_h = src_h / scale_factor;
  95. dst_w = src_w / scale_factor;
  96. assert(dst_w > 3 && dst_h > 3);
  97. //BYTE*smoothed = new BYTE[src_h*src_w];
  98. dst = new BYTE[dst_h*dst_w];
  99. //lowpass_filter(src_gray_data, src_h, src_w,smoothed);
  100. bilinear(src_gray_data, dst, src_w, src_h, dst_h, dst_w);
  101. /*for (int i = 0; i < dst_h - 1; i++)
  102. for (int j = 0; j < dst_w - 1; j++)
  103. {
  104. int srcY = 2 * i + 1;
  105. int srcX = 2 * j + 1;
  106. double re = src_gray_data[srcY*src_w + srcX] * 0.25;
  107. re += src_gray_data[(srcY - 1)*src_w + srcX] * 0.125;
  108. re += src_gray_data[(srcY + 1)*src_w + srcX] * 0.125;
  109. re += src_gray_data[srcY*src_w + srcX - 1] * 0.125;
  110. re += src_gray_data[srcY*src_w + srcX + 1] * 0.125;
  111. re += src_gray_data[(srcY - 1)*src_w + srcX + 1] * 0.0625;
  112. re += src_gray_data[(srcY - 1)*src_w + srcX - 1] * 0.0625;
  113. re += src_gray_data[(srcY + 1)*src_w + srcX - 1] * 0.0625;
  114. re += src_gray_data[(srcY + 1)*src_w + srcX + 1] * 0.0625;
  115. dst[i*dst_w + j] = re;
  116. }
  117. for (int i = 0; i < dst_h; i++)
  118. dst[i*dst_w + dst_w - 1] = dst[i*dst_w + dst_w - 2];
  119. for (int i = 0; i < dst_w; i++)
  120. dst[(dst_h - 1)*dst_w + i] = dst[(dst_h - 2)*dst_w + i];*/
  121. }
  122. void HornSchunckTracker::get_pre_frame(BYTE*&gray)//use only at the beginning
  123. {
  124. pre_pyr[0] = gray;
  125. build_pyramid(pre_pyr);
  126. //save_gray("1.bmp", pre_pyr[1], height[1], width[1]);
  127. }
  128. void HornSchunckTracker::discard_pre_frame()
  129. {
  130. for (int i = 0; i < max_pyramid_layer; i++)
  131. delete[]pre_pyr[i];
  132. }
  133. //set the next frame as pre_frame,must dicard pre_pyr in advance
  134. void HornSchunckTracker::get_pre_frame()
  135. {
  136. for (int i = 0; i < max_pyramid_layer; i++)
  137. pre_pyr[i] = next_pyr[i];
  138. }
  139. //use every time,must after using get_pre_frame(BYTE**pyr)
  140. void HornSchunckTracker::get_next_frame(BYTE*&gray)
  141. {
  142. next_pyr[0] = gray;
  143. build_pyramid(next_pyr);
  144. //save_gray("1.bmp", next_pyr[1], height[1], width[1]);
  145. }
  146. void HornSchunckTracker::get_info(const int nh, const int nw)
  147. {
  148. original_imgH = nh;
  149. original_imgW = nw;
  150. if (isusepyramid)
  151. get_max_pyramid_layer();
  152. else
  153. max_pyramid_layer = 1;
  154. pre_pyr = new BYTE*[max_pyramid_layer];
  155. next_pyr = new BYTE*[max_pyramid_layer];
  156. height = new int[max_pyramid_layer];
  157. width = new int[max_pyramid_layer];
  158. height[0] = nh;
  159. width[0] = nw;
  160. }
  161. //低通滤波
  162. void HornSchunckTracker::lowpass_filter(double*&src, const int H, const int W, double*&smoothed)
  163. {
  164. //tackle with border
  165. for (int i = 0; i < H; i++)
  166. {
  167. smoothed[i*W] = src[i*W];
  168. smoothed[i*W + W - 1] = src[i*W + W - 1];
  169. }
  170. for (int i = 0; i < W; i++)
  171. {
  172. smoothed[i] = src[i];
  173. smoothed[(H - 1)*W + i] = src[(H - 1)*W + i];
  174. }
  175. for (int i = 1; i < H - 1; i++)
  176. for (int j = 1; j < W - 1; j++)
  177. {
  178. double re = 0;
  179. re += src[i*W + j] * 0.25;
  180. re += src[(i - 1)*W + j] * 0.125;
  181. re += src[i*W + j + 1] * 0.125;
  182. re += src[i*W + j - 1] * 0.125;
  183. re += src[(i + 1)*W + j] * 0.125;
  184. re += src[(i - 1)*W + j - 1] * 0.0625;
  185. re += src[(i + 1)*W + j - 1] * 0.0625;
  186. re += src[(i - 1)*W + j + 1] * 0.0625;
  187. re += src[(i + 1)*W + j + 1] * 0.0625;
  188. smoothed[i*W + j] = re;
  189. }
  190. }
  191. void HornSchunckTracker::get_fx_fy_ft(BYTE*img1, BYTE*img2, int w, int h, double*fx, double*fy, double*ft)
  192. {
  193. for (int i = 0; i < h - 1; i++)
  194. for (int j = 0; j < w - 1; j++)
  195. {
  196. 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]
  197. + img2[i*w + j + 1] - img2[i*w + j] + img2[(i + 1)*w + j + 1] - img2[(i + 1)*w + j]);
  198. 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]
  199. + img2[(i + 1)*w + j] - img2[i*w + j] + img2[(i + 1)*w + j + 1] - img2[i*w + j + 1]);
  200. ft[i*w + j] = 0.25 * (img2[i*w + j] - img1[i*w + j] +img2[(i + 1)*w + j] - img1[(i + 1)*w + j] +
  201. img2[(i + 1)*w + j + 1] - img1[(i + 1)*w + j + 1] + img2[i*w + j + 1] - img1[i*w + j + 1]);
  202. }
  203. for (int i = 0; i < h; i++)
  204. {
  205. //fx[i*w] = fx[i*w + w - 2];
  206. fx[i*w + w - 1] = fx[i*w + w - 2];
  207. //fy[i*w] = fy[i*w + w - 2];
  208. fy[i*w + w - 1] = fy[i*w + w - 2];
  209. //ft[i*w] = ft[i*w + w - 2];
  210. ft[i*w + w - 1] = ft[i*w + w - 2];
  211. }
  212. for (int j = 0; j < w; j++)
  213. {
  214. //fx[j] = fx[h + j];
  215. fx[(h - 1)*w + j] = fx[(h - 2)*w + j];
  216. //fy[j] = fy[h + j];
  217. fy[(h - 1)*w + j] = fy[(h - 2)*w + j];
  218. //ft[j] = ft[h + j];
  219. ft[(h - 1)*w + j] = ft[(h - 2)*w + j];
  220. }
  221. }
  222. //取得计算得到的光流值,u、v为out型参数
  223. void HornSchunckTracker::get_optical_flow(double*&u, double*&v)
  224. {
  225. assert(optical_field_U&&optical_field_V);
  226. u = optical_field_U;
  227. v = optical_field_V;
  228. }
  229. //int save_gray(char * savebmp_file, LPBYTE gray, int height, int width);
  230. //返回求得的光流场,大小为原始图像大小
  231. void HornSchunckTracker::HornSchunck()
  232. {
  233. //save_gray("22.bmp", pre_pyr[0], height[0], width[0]);
  234. //初始化光流场为0
  235. if (optical_field_U)
  236. delete[]optical_field_U;
  237. if (optical_field_V)
  238. delete[]optical_field_V;
  239. optical_field_U = new double[width[max_pyramid_layer - 1]
  240. * height[max_pyramid_layer - 1]];
  241. optical_field_V = new double[width[max_pyramid_layer - 1]
  242. * height[max_pyramid_layer - 1]];
  243. memset(optical_field_U, 0, sizeof(double)*width[max_pyramid_layer - 1]
  244. * height[max_pyramid_layer - 1]);
  245. memset(optical_field_V, 0, sizeof(double)*width[max_pyramid_layer - 1]
  246. * height[max_pyramid_layer - 1]);
  247. //使用金字塔计算光流
  248. for (int i = max_pyramid_layer - 1; i >= 0; i--)
  249. {
  250. double*Ix = new double[width[i] * height[i]];
  251. double*Iy = new double[width[i] * height[i]];
  252. double*It = new double[width[i] * height[i]];
  253. //求偏导
  254. get_fx_fy_ft(pre_pyr[i], next_pyr[i], width[i], height[i], Ix, Iy, It);
  255. //将光流场平滑
  256. double*smoothed_U = new double[width[i] * height[i]];
  257. double*smoothed_V = new double[width[i] * height[i]];
  258. if (i == max_pyramid_layer - 1)
  259. {
  260. memset(smoothed_U, 0, sizeof(double)*width[i] * height[i]);
  261. memset(smoothed_V, 0, sizeof(double)*width[i] * height[i]);
  262. }
  263. else
  264. {
  265. lowpass_filter(optical_field_U, height[i], width[i], smoothed_U);
  266. lowpass_filter(optical_field_V, height[i], width[i], smoothed_V);
  267. }
  268. double error = 1000000;
  269. int iteration = 0;
  270. //迭代计算每个像素的光流,直到收敛或达到最大迭代次数
  271. while (error > precision&&iteration < maxiteration)
  272. {
  273. iteration++;
  274. error = 0;
  275. //计算该层金字塔的光流
  276. for (int j = 0; j < height[i]; j++)
  277. for (int k = 0; k < width[i]; k++)
  278. {
  279. //采用改进方法,光流速度需大于阈值,这样不仅准确度增加,计算量也会减小
  280. double w = Ix[j*width[i] + k] * Ix[j*width[i] + k]
  281. + Iy[j*width[i] + k] * Iy[j*width[i] + k] > threshold ? 1 : 0;
  282. double u_pre = optical_field_U[j*width[i] + k];
  283. double v_pre = optical_field_V[j*width[i] + k];
  284. double utemp = smoothed_U[j*width[i] + k];//get_average4(optical_field_U, height[i], width[i], j, k);
  285. double vtemp = smoothed_V[j*width[i] + k]; //get_average4(optical_field_V, height[i], width[i], j, k);
  286. double denominator = lamda + w*(Ix[j*width[i] + k] * Ix[j*width[i] + k]
  287. + Iy[j*width[i] + k] * Iy[j*width[i] + k]);
  288. double numerator = Ix[j*width[i] + k] * utemp + Iy[j*width[i] + k] *
  289. vtemp + It[j*width[i] + k];
  290. optical_field_U[j*width[i] + k] = utemp - Ix[j*width[i] + k] *w*numerator / denominator;
  291. optical_field_V[j*width[i] + k] = vtemp - Iy[j*width[i] + k] *w*numerator / denominator;
  292. error += pow(optical_field_U[j*width[i] + k] - u_pre,2) +
  293. pow(optical_field_V[j*width[i] + k] - v_pre,2);
  294. }
  295. //下一次迭代前重新平滑光流场
  296. if (error >exp(double(max_pyramid_layer-i))*precision&&iteration < maxiteration)
  297. {
  298. lowpass_filter(optical_field_U, height[i], width[i], smoothed_U);
  299. lowpass_filter(optical_field_V, height[i], width[i], smoothed_V);
  300. }
  301. }
  302. delete[]smoothed_U, smoothed_V,Ix,Iy,It;
  303. if (i == 0)//得到最终光流场
  304. {
  305. return;
  306. }
  307. //下一层的光流场
  308. double*new_of_u = new double[width[i - 1] * height[i - 1]];
  309. double*new_of_v = new double[width[i - 1] * height[i - 1]];
  310. //上采样
  311. pyramid_up(optical_field_U, width[i], height[i], new_of_u, width[i - 1], height[i - 1]);
  312. pyramid_up(optical_field_V, width[i], height[i], new_of_v, width[i - 1], height[i - 1]);
  313. //将每个像素的光流按缩放因子放大,得到下一层的光流场的初值
  314. //double scale = double(height[i - 1]) / height[i];
  315. for (int j = 0; j < height[i - 1]; j++)
  316. for (int k = 0; k < width[i - 1]; k++)
  317. {
  318. new_of_u[j*width[i - 1] + k] *= scale_factor;
  319. new_of_v[j*width[i - 1] + k] *= scale_factor;
  320. }
  321. delete[]optical_field_U, optical_field_V;
  322. optical_field_U = new_of_u;
  323. optical_field_V = new_of_v;
  324. }
  325. }
  326. //上采样,采用双线性插值,用双立方插值应该更精确
  327. void HornSchunckTracker::pyramid_up(double*src, int srcW, int srcH, double*dst, int dstW, int dstH)
  328. {
  329. bilinear(src, dst, srcW, srcH, dstH, dstW);
  330. }
  331. //双线性插值
  332. void HornSchunckTracker::bilinear(double* lpSrc, double* lpDst, int nW, int nH, int H1, int W1)
  333. {
  334. float fw = float(nW) / W1;
  335. float fh = float(nH) / H1;
  336. int y1, y2, x1, x2, x0, y0;
  337. float fx1, fx2, fy1, fy2;
  338. for (int i = 0; i < H1; i++)
  339. {
  340. y0 = i*fh;
  341. y1 = int(y0);
  342. if (y1 == nH - 1) y2 = y1;
  343. else y2 = y1 + 1;
  344. fy1 = y1 - y0;
  345. fy2 = 1.0f - fy1;
  346. for (int j = 0; j < W1; j++)
  347. {
  348. x0 = j*fw;
  349. x1 = int(x0);
  350. if (x1 == nW - 1) x2 = x1;
  351. else x2 = x1 + 1;
  352. fx1 = y1 - y0;
  353. fx2 = 1.0f - fx1;
  354. float s1 = fx1*fy1;
  355. float s2 = fx2*fy1;
  356. float s3 = fx2*fy2;
  357. float s4 = fx1*fy2;
  358. double c1, c2, c3, c4;
  359. c1 = lpSrc[y1*nW + x1];
  360. c2 = lpSrc[y1*nW + x2];
  361. c3 = lpSrc[y2*nW + x1];
  362. c4 = lpSrc[y2*nW + x2];
  363. double r;
  364. r = (c1*s3) + (c2*s4) + (c3*s2) + (c4*s1);
  365. lpDst[i*W1 + j] = r;
  366. }
  367. }
  368. }
  369. //双线性插值
  370. void HornSchunckTracker::bilinear(BYTE* lpSrc, BYTE* lpDst, int nW, int nH, int H1, int W1)
  371. {
  372. float fw = float(nW) / W1;
  373. float fh = float(nH) / H1;
  374. int y1, y2, x1, x2, x0, y0;
  375. float fx1, fx2, fy1, fy2;
  376. for (int i = 0; i < H1; i++)
  377. {
  378. y0 = i*fh;
  379. y1 = int(y0);
  380. if (y1 == nH - 1) y2 = y1;
  381. else y2 = y1 + 1;
  382. fy1 = y1 - y0;
  383. fy2 = 1.0f - fy1;
  384. for (int j = 0; j < W1; j++)
  385. {
  386. x0 = j*fw;
  387. x1 = int(x0);
  388. if (x1 == nW - 1) x2 = x1;
  389. else x2 = x1 + 1;
  390. fx1 = y1 - y0;
  391. fx2 = 1.0f - fx1;
  392. float s1 = fx1*fy1;
  393. float s2 = fx2*fy1;
  394. float s3 = fx2*fy2;
  395. float s4 = fx1*fy2;
  396. double c1, c2, c3, c4;
  397. c1 = lpSrc[y1*nW + x1];
  398. c2 = lpSrc[y1*nW + x2];
  399. c3 = lpSrc[y2*nW + x1];
  400. c4 = lpSrc[y2*nW + x2];
  401. double r;
  402. r = (c1*s3) + (c2*s4) + (c3*s2) + (c4*s1);
  403. lpDst[i*W1 + j] = BYTE(r);
  404. }
  405. }
  406. }
  407. void HornSchunckTracker::set_paras(double lamda, double threshold, double scalefactor)
  408. {
  409. this->lamda = lamda;
  410. this->threshold = threshold;
  411. scale_factor = scalefactor;
  412. }
  413. //double HornSchunckTracker::get_average4(double*src, const int height, const int width, const int i, const int j)
  414. //{
  415. // if (j < 0 || j >= width) return 0;
  416. // if (i < 0 || i >= height) return 0;
  417. //
  418. // double val = 0.0;
  419. // int tmp = 0;
  420. // if ((i - 1) >= 0){
  421. // ++tmp;
  422. // val += src[(i - 1)*width + j];
  423. // }
  424. // if (i + 1<height){
  425. // ++tmp;
  426. // val += src[(i + 1)*width + j];
  427. // }
  428. // if (j - 1 >= 0){
  429. // ++tmp;
  430. // val += src[i*width + j - 1];
  431. // }
  432. // if (j + 1<width){
  433. // ++tmp;
  434. // val += src[i*width + j + 1];
  435. // }
  436. // return val / tmp;
  437. //
  438. //}

下面是两帧图像和检测结果

20151030224744051 20151030224757826

20151030224949095

可以看出对边缘的光流检测较好,内部点的光流检测较难。

发表评论

表情:
评论列表 (有 0 条评论,226人围观)

还没有评论,来说两句吧...

相关阅读

    相关 基本原理

    摘自百度百科 光流是一种简单实用的图像运动的表达方式\[1\],通常定义为一个图像序列中的图像亮度模式的表观运动,即空间物体表面上的点的运动速度在视觉传感器的成像平面上的表达

    相关 Lucas-Kanade

    Lucas-Kanade光流法是通过先在前后两帧图像里分别建立一个固定大小窗口,然后找到让两个窗口间像素强度差的平方和最小的位移。然后将窗口内像素的移动近似为这样的位移向量,然