基于粒子滤波的目标追踪

基于粒子滤波的目标追踪
基于粒子滤波的目标追踪
particle filter
object tracking
读"K. Nummiaro, E. Koller-Meier, L. Van Gool. An adaptive color-based particle filter[J], Image and Vision Computing, 2002"笔记

粒子滤波中最重要的一个过程就是重要性重采样， Sampling Importance Resampling (SIR).

这篇博客基于粒子滤波的物体跟踪讲的挺形象的 :) 。这里拿过来记录一下。

初始化阶段
基于粒子滤波的目标追踪方法是一种生成式跟踪方法，所以要有一个初始化的阶段。对于第一帧图像，人工标定出待检测的目标，对该目标区域提出特征。论文里将图像的HSV空间划分成不同的bins，然后统计直方图。另外

To increase the reliability of the color distribution when boundary pixels belong to the background or get occluded, smaller weights are assigned to the pixels that are further away from the region center by employing a weighting function

是到区域中心的距离，所以直方图分布可如下计算

是归一化因子，是区域半径，是示性函数，是将加权后的值添加到对应的区间内。

搜索阶段 - 放狗
现在已经知道了目标的特征，然后就在目标的周围撒点(particle),即放狗。放狗的方法有很多。如a)均匀的撒点;b)按高斯分布撒点，就是近的地方撒得多，远的地方撒的少。论文里使用的是后一种方法。每一个粒子都计算所在区域内的颜色直方图，如初始化提取特征一样，然后对所有的相似度进行归一化。文中相似性使用的是巴氏距离。

重要性重采样
根据第2阶段获得的相似度重新撒粒子，相似度高的粒子周围多撒，相似度低的地方少撒。这个过程其实可以重复几次，但为了检测速度，重采样一次也可以。

状态转移
对上一阶段多次重采样后获得的粒子计算下一时刻的位置。

是多元高斯分布变量。

观测阶段
i. 在出计算概率颜色直方图
ii. 计算各个粒子和目标的巴氏距离
iii. 更新每个粒子的权重

决策阶段
每个粒子都能够获得一个和目标的相似度，这个相似度体现了该区域是目标的置信度，可以把所有例子使用相似度加权后的结果作为目标可能的位置。

粒子滤波更新过程

iteration step of the color-based particle filter.jpg

matlab代码实现

% REFERENCE:
% An adaptive color-based particle filter [J]. Image and Vision
% Computing, 21 (2003) 99-110
%
% note:
% the structure of sample
% {1} x,
% {2} y,
% {3} vx,
% {4} vy,
% {5} hx,
% {6} hy,
% {7} a.
% and the selected region is an ellipse.
% These functions are programmed for RGB images.

function [sample,initHist]=ACPF(image,initSample,N,targetHist,bin)
% INPUTS:
% image - current frame
% initSample - init region to be detected
% N - the numbers of candidates
% targetHist - the histgram of the target
% bins -the element represents the number of bins in each channel
% OUTPUTS:
% sample - detected region

% predefined parameters
global curFrame; % current frame
global velocityTurb; % the turbulence of the velocity
global scaleTurb; % the turbulence of the scale
global bins; % the element represents the number of bins in each channel
global deltaT; % the time margin of consecutive frames
global SIGMA2; % the variance of the Gaussian specifing the weight in function: observe
global shiftTurb;

curFrame=image;
velocityTurb=4;
scaleTurb=4;
shiftTurb=40;
bins=bin;
deltaT=0.05;
SIGMA2=0.02; % 2*sigma^2
tolerance=initSample{3}*deltaT+scaleTurb;
A=eye(7);
beta=0.9; % the similarity ratio between new histgram and the origin histgram
alpha=0.1; % the learning factor of object histgram

[initHist,sampleSet,weight]=initialization(initSample,N);

tempHist=sqrt(initHist.*targetHist);
rho=sum(tempHist,1);
PIE=exp(-(1-rho)/SIGMA2);
fprintf('similarity:%4f ',PIE);
if PIE>beta
initHist=alpha*initHist+(1-alpha)*targetHist;
else
initHist=targetHist;
end

MAXITER=2;
iter=0;
oldsample=initSample;
while iter<MAXITER
sampleSet=reselect(sampleSet,weight);
sampleSet=propagate(sampleSet,A);
weight=observe(sampleSet,initHist);
sample=estimate(sampleSet,weight);
if different(sample,oldsample)<tolerance
break;
end
oldsample=sample;
iter=iter+1;
end
% draw the scatter dots ------------------------
for i=1:N
tempSample=sampleSet{i};
plot(tempSample{1},tempSample{2},'o','MarkerSize',5,'MarkerFaceColor','g');
hold on
end
plot(sample{1},sample{2},'o','MarkerSize',5,'MarkerFaceColor','r');
hold on
hold off
% -------------------------------------------------
end

function err=different(sample1,sample2)
% calculate the difference between sample1 and sample2
err=0;
for i=1:7
err=err+abs(sample1{i}-sample2{i});
end
end

function [initHist, sampleSet,weight]=initialization(initSample,N)
% Initializing the target in each frame at the beginning of
% tracking process.
% INPUTS:
% initSample -the sample which will be tracked in current frame
% N -the length of the sampleSets
% channel, respectively.
% OUTPUTS:
% initHist -the histgram of the target model
% sampleSet -the set of candidate samples
% weight -the initial weight of each sample in the sampleSet
global velocityTurb;
global scaleTurb;
global shiftTurb;

sampleSet=cell(N,1);
weight=zeros(N,1);
numSamples=0; % the number of initialized samples
while numSamples<N
randoms=randn(7,1);
sampleTemp{1}=round(initSample{1}+randoms(1)*shiftTurb);
sampleTemp{2}=round(initSample{2}+randoms(2)*shiftTurb);
sampleTemp{3}=initSample{3}+randoms(3)*velocityTurb;
sampleTemp{4}=initSample{4}+randoms(4)*velocityTurb;
sampleTemp{5}=round(initSample{5}+randoms(5)*scaleTurb);
sampleTemp{6}=round(initSample{6}+randoms(6)*scaleTurb);
sampleTemp{7}=scaleTurb;
if isValidate(sampleTemp)==1
numSamples=numSamples+1;
sampleSet{numSamples}=sampleTemp;
weight(numSamples)=1/N;
end
end
initHist=calColorHist(initSample);

end

function colorHist=calColorHist(sample)
% calculating the color histgram of the region selected by the sample
% note:
% the weight function if r<1 then k(r)=1-r^2 else k(r)=0
global curFrame;
global bins;
x=sample{1}-sample{5};
y=sample{2}-sample{6};
region=curFrame(y:y+2*sample{6},x:x+2*sample{5},:);
[m,n,~]=size(region);
tempBins=zeros(bins);
margin=1.1 ./bins;
for i=1:m
for j=1:n
r=double(region(i,j,1));
g=double(region(i,j,2));
b=double(region(i,j,3));
tempxy=(i-double(sample{6}))^2/(double(sample{6}))^2+1.0*(j-double(sample{5}))^2/(double(sample{5}))^2;
if tempxy<=1
tempBins(fix(r/margin(1)+1), fix(g/margin(2)+1), fix(b/margin(3)+1))=...
tempBins(fix(r/margin(1)+1), fix(g/margin(2)+1), fix(b/margin(3)+1))+...
1-tempxy;
end
end
end
colorHist=tempBins(:);
colorHist=colorHist./sum(colorHist,1);
end

function sampleSet=reselect(sampleSet,weight)
% reselect process.
% new samples are selected according the corresponding weight

% following the steps in the referred paper.
N=size(weight,1);
cumulativeW=zeros(N,1);
cumulativeW(1)=weight(1);
% calculate the normalized cumulative probabilities
for i=2:N
cumulativeW(i)=cumulativeW(i-1)+weight(i);
end
cumulativeW=cumulativeW./cumulativeW(N);
index=zeros(N,1);
for i=1:N
% generate a uniformly distributed random number r belong [0,1]
r=rand(1);
[~,ind]=find(cumulativeW'>=r);
if isempty(ind)
ind=N;
end
index(i)=ind(1);
end
tempSampleSet=cell(N,1);
for i=1:N
tempSampleSet{i}=sampleSet{index(i)};
end
sampleSet=tempSampleSet;
end

function sampleSet=propagate(sampleSet, A )
global deltaT;
global velocityTurb;
global scaleTurb;
global shiftTurb;

MINWIDTH=10;
MINHEIGHT=15;
N=length(sampleSet);
numSamples=0;
while numSamples<N
tempSample=sampleSet{numSamples+1};
% convernient for expand to high order model A
tempVector=zeros(7,1);
for j=1:7
tempVector(j)=tempSample{j};
end
tempVector=A*tempVector;
randoms=0.7*randn(7,1);
tempSample{1}=round(tempVector(1)+deltaT*tempSample{3}+randoms(1)*shiftTurb);
tempSample{2}=round(tempVector(2)+deltaT*tempSample{4}+randoms(2)*shiftTurb);
tempSample{3}=tempVector(3)+randoms(3)*velocityTurb;
tempSample{4}=tempVector(4)+randoms(4)*velocityTurb;
tempSample{5}=max(round(tempVector(5)+randoms(5)*scaleTurb),MINWIDTH);
tempSample{6}=max(round(tempVector(6)+randoms(6)*scaleTurb),MINHEIGHT);
tempSample{7}=scaleTurb;
if isValidate(tempSample)==1
numSamples=numSamples+1;
sampleSet{numSamples}=tempSample;
end
end
end

function weight=observe(sampleSet,colorHist)
% calculate the color distribution of each sample in sampleSet
% return the weight
global SIGMA2;
N=length(sampleSet);
weight=zeros(N,1);
for i=1:N
sampleHist=calColorHist(sampleSet{i});
tempHist=sqrt(sampleHist.*colorHist);
rho=sum(tempHist,1);% Bhattacharyya cofficient
weight(i)=exp(-(1-rho)/SIGMA2);
end
weight=weight./sum(weight,1);
end

function sample=estimate(sampleSet,weight)
% estimate the mean state of the set
N=length(sampleSet);
sample=cell(7,1);
for i=1:7
sample{i}=0;
for j=1:N
tempSample=sampleSet{j};
sample{i}=sample{i}+weight(j)*tempSample{i};
end
end
sample{1}=round(sample{1});
sample{2}=round(sample{2});
sample{5}=round(sample{5});
sample{6}=round(sample{6});
end

function validate=isValidate(sample)
% change the validation of the sample
% if the coordinate of the corner exceed the image edge,then the sample is not rightful
global curFrame;
[m,n,~]=size(curFrame);
x=sample{1};
y=sample{2};
hx=sample{5};
hy=sample{6};
validate=0; %false
if fix(x-hx-0.5)>=1 && ceil(x+hx+0.5)<=n && fix(y-hy-0.5)>=1 && ceil(y+hy+0.5)<=m
validate=1;
end

end

参考文献

基于粒子滤波的物体跟踪

基于粒子滤波器的目标跟踪算法及实现
相关阅读:
uva694 The Collatz Sequence
WindowsPhone7开发简单豆瓣网应用程序之主页面功能实现
 使用codeplex管理WP7项目开发版本
 WindowsPhone7开发简单豆瓣网应用程序之界面设计
 WindowsPhone操作SkyDrive之获取共享文件
 Windows Phone7监测网络接口及状态变化
 《设计模式》杂记之里氏替换原则
 Silverlight杂记之HTTP通信WebClient介绍
 WindowsPhone获取是否进行拨电话信息[使用PhoneCallTask]
《设计模式》杂记之单一职责原则
原文地址：https://www.cnblogs.com/YiXiaoZhou/p/5877398.html

基于粒子滤波的目标追踪

基于粒子滤波的目标追踪

读"K. Nummiaro, E. Koller-Meier, L. Van Gool. An adaptive color-based particle filter[J], Image and Vision Computing, 2002"笔记

粒子滤波更新过程

matlab代码实现

参考文献