TLD之L-K光流算法代码篇
2015-12-14 16:54
232 查看
代码是国外的大神Musawir Ali Shah 写的,我拿过来研究了一下,代码原理和步骤都明白了,代码里面都有我个人参考的资料名称或者链接,以及我个人的注释,这次就先上代码!下一次上原理!
运行结果:
这幅是通过高斯金字塔获得的图像(层数分别为:1,2,3,4)
主函数: main.m
clc; clear all; close all; file1='yos9.tif'; file2='yos10.tif'; density=10; %% Read Images %% img1 = im2double (imread (file1)); %% Take alternating rows and columns %% [odd1, even1] = split (img1); %在行方向分解成奇数行和偶数行 img2 = im2double (imread (file2)); [odd2, even2] = split (img2); %% Run Lucas Kanade %% [Dx, Dy] = Estimate (odd1, odd2); %% Plot %% figure; [maxI,maxJ]=size(Dx); Dx=Dx(1:density:maxI,1:density:maxJ); Dy=Dy(1:density:maxI,1:density:maxJ); quiver(1:density:maxJ,(maxI):(-density):1,Dx,-Dy,1); axis square;
Estimate.m
<span style="font-size:18px;">% Run Lucas Kanade on all levels and interpolate(插入) %% function [Dx, Dy] = Estimate (img1, img2) %[Dx, Dy] = Estimate (odd1, odd2); img1=odd1; img2=odd2; level = 4; half_window_size=2; [m, n] = size (img1); G00 = img1; G10 = img2; if (level>0) G01 = reduce (G00); G11 = reduce (G10); %进行的是高斯金字塔分解算法。 end if (level>1) G02 = reduce (G01); G12 = reduce (G11); end if (level>2) G03 = reduce (G02); G13 = reduce (G12); end if (level>3) G04 = reduce (G03); G14 = reduce (G13); end subplot(2,2,1); imshow(G00) subplot(2,2,2); imshow(G01) subplot(2,2,3); imshow(G02) subplot(2,2,4); imshow(G03) l = level; for i=level:-1:0, if (l == level) switch (l) case 4, Dx = zeros (size (G04)); Dy = zeros (size (G04)); case 3, Dx = zeros (size (G03)); Dy = zeros (size (G03)); case 2, Dx = zeros (size (G02)); Dy = zeros (size (G02)); case 1, Dx = zeros (size (G01)); Dy = zeros (size (G01)); case 0, Dx = zeros (size (G00)); Dy = zeros (size (G00)); end else Dx = expand (Dx); Dy = expand (Dy); Dx = Dx .* 2; Dy = Dy .* 2; end switch (l) case 4, W = warp (G04, Dx, Dy); [Vx, Vy] = EstimateMotion (W, G14, half_window_size); case 3, W = warp (G03, Dx, Dy); [Vx, Vy] = EstimateMotion (W, G13, half_window_size); case 2, W = warp (G02, Dx, Dy); [Vx, Vy] = EstimateMotion (W, G12, half_window_size); case 1, W = warp (G01, Dx, Dy); [Vx, Vy] = EstimateMotion (W, G11, half_window_size); case 0, W = warp (G00, Dx, Dy); [Vx, Vy] = EstimateMotion (W, G10, half_window_size); end [m, n] = size (W); Dx(1:m, 1:n) = Dx(1:m,1:n) + Vx; Dy(1:m, 1:n) = Dy(1:m, 1:n) + Vy; smooth (Dx); smooth (Dy); l = l - 1; end </span>
算法的核心函数:EstimateMotion.m
<span style="font-size:18px;">% Lucas Kanade on the image sequence at pyramid step %% function [Vx, Vy] = EstimateMotion (W, G1, half_window_size) % [Vx, Vy] = EstimateMotion (W, G14, half_window_size);G1=G14 [m, n] = size (W); Vx = zeros (size (W)); Vy = zeros (size (W));%初始化 N = zeros (2*half_window_size+1, 5);%为空间邻域5*5 此处参考的是:基于光流场的视频运动检测 for i = 1:m, l = 0; for j = 1-half_window_size:1+half_window_size, l = l + 1; N (l,:) = getSlice (W, G1, i, j, half_window_size);%求出想,方向的偏导数 end replace = 1; for j = 1:n, t = sum (N); [v, d] = eig ([t(1) t(2);t(2) t(3)]);%梯度矩阵 namda1 = d(1,1); namda2 = d(2,2); if (namda1 > namda2) tmp = namda1; namda1 = namda2; namda2 = tmp; tmp1 = v (:,1); v(:,1) = v(:,2); v(:,2) = tmp1; end if (namda2 < 0.001)%不满足 Vx (i, j) = 0; Vy (i, j) = 0; elseif (namda2 > 100 * namda1)%不满足 n2 = v(1,2) * t(4) + v(2,2) * t(5); Vx (i,j) = n2 * v(1,2) / namda2; Vy (i,j) = n2 * v(2,2) / namda2; else%执行 n1 = v(1,1) * t(4) + v(2,1) * t(5);%计算的是图像的光流(inv(G)*b) t(4)和t(5)表示的是图像的比匹配向量 n2 = v(1,2) * t(4) + v(2,2) * t(5); Vx (i,j) = n1 * v(1,1) / namda1 + n2 * v(1,2) / namda2;%就这个地方的公式仅仅针对半正定(类似于[a,b;b,c]的形式)的矩阵可以有如下的计算:这4行等价于inv(G)*b 即inv([t(1) t(2);t(2) t(3)])*[t(4,t(5))]' Vy (i,j) = n1 * v(2,1) / namda1 + n2 * v(2,2) / namda2;%代码的结尾后面有相关的举例说明。 end N (replace, :) = getSlice (W, G1, i, j+half_window_size+1, half_window_size); replace = replace + 1; if (replace == 2 * half_window_size + 2) replace = 1; end end end %%解释那个不清楚的公式:只有在半正定的情况下才能使用: % 计算方法一“ % a=[1,2;2,3]; % t=[2,3]'; % inv(a)*t % 计算方法二 % %仅仅只有在半正定的情况下才能使用 % [v,d]=eig(a); % namda1 = d(1,1); namda2 = d(2,2); % n1 = v(1,1) * t(1) + v(2,1) * t(2); % n2 = v(1,2) * t(1) + v(2,2) * t(2); % n1 * v(1,1) / namda1 + n2 * v(1,2) / namda2 % n1 * v(2,1) / namda1 + n2 * v(2,2) / namda2</span>
getSlice.m
<span style="font-size:18px;">%计算图像的不匹配向量 function = getSlice (W, G1, i, j, half_window_size) % N (l,:) = getSlice (W, G1, i, j, half_window_size) N = zeros (1, 5); [m, n] = size (W); for y=-half_window_size:half_window_size, Y1 = y +i; if (Y1 < 1) Y1 = Y1 + m; elseif (Y1 > m) Y1 = Y1 - m; end X1 = j; if (X1 < 1) X1 = X1 + n; elseif (X1 > n) X1 = X1 - n; end DeriX = Derivative (G1, X1, Y1, 'x'); DeriY = Derivative (G1, X1, Y1, 'y'); N = N + [ DeriX * DeriX, ... DeriX * DeriY, ... DeriY * DeriY, ... DeriX * (G1 (Y1, X1) - W (Y1, X1)), ... %G1 (Y1, X1) - W (Y1, X1)计算的是图像的像素差 DeriX * (G1 (Y1, X1) - W (Y1, X1))求和后就是算的是图像的不匹配向量 DeriY * (G1 (Y1, X1) - W (Y1, X1))]; % end </span>
Derivative.m
<span style="font-size:18px;">% Calculates the Fx Fy %% DeriX = Derivative (G1, X1, Y1, 'x'); function result = Derivative (img, x, y, direction) [m, n] = size (img); switch (direction) case 'x', if (x == 1) result = img (y, x+1) - img (y, x); elseif (x == n) result = img (y, x) - img (y, x-1); else result = 0.5 * (img (y, x+1) - img (y, x-1)); end case 'y', if (y == 1) result = img (y+1, x) - img (y, x); elseif (y == m) result = img (y, x) - img (y-1, x); else result = 0.5 * (img (y+1, x) - img (y-1, x)); end end</span>
smooth.m
<span style="font-size:18px;">function result = smooth (img) result = expand (reduce (img)); </span>
expand.m
<span style="font-size:18px;">% The Expansion Function for pyramid %% function result = expand (ori) [m,n] = size (ori); mid = zeros (m, n); m1 = m * 2; n1 = n * 2; result = zeros (m1, n1); w = generateFilter (0.4); for j=1:m, t = zeros (1, n1); t(1:2:n1-1) = ori (j,1:n); tmp = conv ([ori(j,n) 0 t ori(j,1) 0], w); %tmp长度为:n1+5+4-1 mid(j,1:n1) = 2 .* tmp (5:n1+4); end for i=1:n1, t = zeros (1, m1); t(1:2:m1-1) = mid (1:m,i)'; tmp = conv([mid(m,i) 0 t mid(1,i) 0], w); result(1:m1,i) = 2 .* tmp (5:m1+4)'; end </span>
Warp.m
<span style="font-size:18px;">% Interpolation %% function result = warp (img, Dx, Dy) % W = warp (G04, Dx, Dy); img=G04; [m, n] = size (img); [x,y] = meshgrid (1:n, 1:m); %x,y均是m*n的矩阵 x = x + Dx (1:m, 1:n); y = y + Dy (1:m,1:n); for i=1:m, for j=1:n, if x(i,j)>n x(i,j) = n; end if x(i,j)<1 x(i,j) = 1; end if y(i,j)>m y(i,j) = m; end if y(i,j)<1 y(i,j) = 1; end end end result = interp2 (img, x, y, 'linear');%看帮助文档 surf(x,y,result);</span>
reduce.m
<span style="font-size:18px;">% The Reduce Function for pyramid %% function result = reduce (ori) [m,n] = size (ori); mid = zeros (m, n); m1 = round (m/2); n1 = round (n/2); result = zeros (m1, n1); w = generateFilter (0.4); %产生的是一个长度为5的滤波器。窗口函数 for j=1:m, %行方向的线性卷积 tmp = conv([ori(j,n-1:n) ori(j,1:n) ori(j,1:2)], w); %计算的是两个向量的线性卷积(长度为2+n+2+4-1=n+7):http://blog.csdn.net/styxqdz/article/details/4943089 mid(j,1:n1) = tmp(5:2:n+4); end for i=1:n1, %列方向的线性卷积 tmp = conv([mid(m-1:m,i); mid(1:m,i); mid(1:2,i)]', w); result(1:m1,i) = tmp(5:2:m+4)'; end </span>
generateFilter .m
<span style="font-size:18px;">function filter = generateFilter (alpha) filter = [0.25-alpha/2 0.25 alpha 0.25 0.25-alpha/2]; %0.25=1/4(4表示的是滤波器的长度)</span>
split .m
<span style="font-size:18px;">function [odd, even] = split (img); [m, n] = size (img); odd = img (1:2:m, :); even = img (2:2:m, :); </span>
运行结果:
这幅是通过高斯金字塔获得的图像(层数分别为:1,2,3,4)
相关文章推荐
- Remove Element leetcode java
- java 将字符串数组变为字典顺序排序后的字符串数组
- Windows2008 64位系统asp连接access数据库解决方法
- php gd库
- php数组函数-array_reduce()
- Fragment与FragmentPagerAdapter的使用
- Git与Github折腾记-常用命令汇总
- Python基本内置数据类型有哪些?
- Git与Github折腾记-常用命令汇总
- 老李分享:qtp自动化测试框架赏析-关键字自动化测试框架 2
- c# string类型转变成Stream类型
- Page.java
- 老李分享:qtp自动化测试框架赏析-关键字自动化测试框架 1
- C语言利用Cairo图形库绘制太极图
- windows下设置QT程序的版本信息、程序图标和可执行文件图标
- Microsoft Visual Studio与Firefly 加载的项目已建议,更新源代码地位问题
- java socket编程
- 【Python】利用count函数求list中每个元素出现的次数,求众数的改进
- sping的SpringBootServletInitializer组件
- 【Java学习笔记】JDBC连接mySql数据库