基於粒子濾波的目標追蹤

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\n',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

參考文獻

免責聲明！

本站轉載的文章為個人學習借鑒使用，本站對版權不負任何法律責任。如果侵犯了您的隱私權益，請聯系本站郵箱yoyou2525@163.com刪除。

猜您在找 室內定位系列（六）——目標跟蹤（粒子濾波）粒子濾波粒子濾波卡爾曼濾波+單目標追蹤+python-opencv 粒子濾波的退化問題粒子濾波算法理解粒子濾波(particle filter) 二維粒子濾波純代碼 OpenCV--目標追蹤從貝葉斯到粒子濾波——Round 1