首页 > 代码库 > 基于粒子滤波的目标追踪

基于粒子滤波的目标追踪

2024-08-13 14:30:52   225人阅读

基于粒子滤波的目标追踪

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).

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

  1. 初始化阶段
    基于粒子滤波的目标追踪方法是一种生成式跟踪方法，所以要有一个初始化的阶段。对于第一帧图像，人工标定出待检测的目标，对该目标区域提出特征。论文里将图像的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

技术分享

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

技术分享

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

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

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

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

技术分享

技术分享是多元高斯分布变量。

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

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

粒子滤波更新过程

技术分享

iteration step of the color-based particle filter.jpg

matlab代码实现

  1. % REFERENCE: 
  2. % An adaptive color-based particle filter [J]. Image and Vision 
  3. % Computing, 21 (2003) 99-110 
  4. % 
  5. % note: 
  6. % the structure of sample 
  7. % {1} x,  
  8. % {2} y, 
  9. % {3} vx, 
  10. % {4} vy, 
  11. % {5} hx, 
  12. % {6} hy, 
  13. % {7} a. 
  14. % and the selected region is an ellipse. 
  15. % These functions are programmed for RGB images. 
  16.  
  17. function [sample,initHist]=ACPF(image,initSample,N,targetHist,bin) 
  18. % INPUTS: 
  19. % image - current frame 
  20. % initSample - init region to be detected 
  21. % N - the numbers of candidates 
  22. % targetHist - the histgram of the target 
  23. % bins -the element represents the number of bins in each channel 
  24. % OUTPUTS: 
  25. % sample - detected region 
  26.  
  27. % predefined parameters  
  28. global curFrame; % current frame 
  29. global velocityTurb; % the turbulence of the velocity 
  30. global scaleTurb; % the turbulence of the scale 
  31. global bins; % the element represents the number of bins in each channel 
  32. global deltaT; % the time margin of consecutive frames 
  33. global SIGMA2; % the variance of the Gaussian specifing the weight in function: observe 
  34. global shiftTurb; 
  35.  
  36. curFrame=image; 
  37. velocityTurb=4; 
  38. scaleTurb=4; 
  39. shiftTurb=40; 
  40. bins=bin; 
  41. deltaT=0.05; 
  42. SIGMA2=0.02; % 2*sigma^2 
  43. tolerance=initSample{3}*deltaT+scaleTurb; 
  44. A=eye(7); 
  45. beta=0.9; % the similarity ratio between new histgram and the origin histgram 
  46. alpha=0.1; % the learning factor of object histgram 
  47.  
  48. [initHist,sampleSet,weight]=initialization(initSample,N); 
  49.  
  50. tempHist=sqrt(initHist.*targetHist); 
  51. rho=sum(tempHist,1); 
  52. PIE=exp(-(1-rho)/SIGMA2); 
  53. fprintf(‘similarity:%4f\n‘,PIE); 
  54. if PIE>beta 
  55. initHist=alpha*initHist+(1-alpha)*targetHist; 
  56. else 
  57. initHist=targetHist; 
  58. end 
  59.  
  60. MAXITER=2; 
  61. iter=0; 
  62. oldsample=initSample; 
  63. while iter<MAXITER 
  64. sampleSet=reselect(sampleSet,weight); 
  65. sampleSet=propagate(sampleSet,A); 
  66. weight=observe(sampleSet,initHist); 
  67. sample=estimate(sampleSet,weight); 
  68. if different(sample,oldsample)<tolerance 
  69. break; 
  70. end 
  71. oldsample=sample; 
  72. iter=iter+1; 
  73. end 
  74. % draw the scatter dots ------------------------ 
  75. for i=1:N 
  76. tempSample=sampleSet{i}; 
  77. plot(tempSample{1},tempSample{2},‘o‘,‘MarkerSize‘,5,‘MarkerFaceColor‘,‘g‘); 
  78. hold on 
  79. end 
  80. plot(sample{1},sample{2},‘o‘,‘MarkerSize‘,5,‘MarkerFaceColor‘,‘r‘); 
  81. hold on 
  82. hold off 
  83. % ------------------------------------------------- 
  84. end 
  85.  
  86. function err=different(sample1,sample2) 
  87. % calculate the difference between sample1 and sample2 
  88. err=0; 
  89. for i=1:7 
  90. err=err+abs(sample1{i}-sample2{i}); 
  91. end 
  92. end 
  93.  
  94. function [initHist, sampleSet,weight]=initialization(initSample,N) 
  95. % Initializing the target in each frame at the beginning of 
  96. % tracking process. 
  97. % INPUTS: 
  98. % initSample -the sample which will be tracked in current frame 
  99. % N -the length of the sampleSets 
  100. % channel, respectively. 
  101. % OUTPUTS: 
  102. % initHist -the histgram of the target model 
  103. % sampleSet -the set of candidate samples 
  104. % weight -the initial weight of each sample in the sampleSet  
  105. global velocityTurb; 
  106. global scaleTurb; 
  107. global shiftTurb; 
  108.  
  109. sampleSet=cell(N,1); 
  110. weight=zeros(N,1); 
  111. numSamples=0; % the number of initialized samples 
  112. while numSamples<N 
  113. randoms=randn(7,1); 
  114. sampleTemp{1}=round(initSample{1}+randoms(1)*shiftTurb); 
  115. sampleTemp{2}=round(initSample{2}+randoms(2)*shiftTurb); 
  116. sampleTemp{3}=initSample{3}+randoms(3)*velocityTurb; 
  117. sampleTemp{4}=initSample{4}+randoms(4)*velocityTurb; 
  118. sampleTemp{5}=round(initSample{5}+randoms(5)*scaleTurb); 
  119. sampleTemp{6}=round(initSample{6}+randoms(6)*scaleTurb); 
  120. sampleTemp{7}=scaleTurb; 
  121. if isValidate(sampleTemp)==1 
  122. numSamples=numSamples+1; 
  123. sampleSet{numSamples}=sampleTemp; 
  124. weight(numSamples)=1/N; 
  125. end 
  126. end 
  127. initHist=calColorHist(initSample); 
  128.  
  129. end 
  130.  
  131. function colorHist=calColorHist(sample) 
  132. % calculating the color histgram of the region selected by the sample  
  133. % note: 
  134. % the weight function if r<1 then k(r)=1-r^2 else k(r)=0 
  135. global curFrame; 
  136. global bins; 
  137. x=sample{1}-sample{5}; 
  138. y=sample{2}-sample{6}; 
  139. region=curFrame(y:y+2*sample{6},x:x+2*sample{5},:); 
  140. [m,n,~]=size(region); 
  141. tempBins=zeros(bins); 
  142. margin=1.1 ./bins; 
  143. for i=1:m 
  144. for j=1:n 
  145. r=double(region(i,j,1)); 
  146. g=double(region(i,j,2)); 
  147. b=double(region(i,j,3)); 
  148. tempxy=(i-double(sample{6}))^2/(double(sample{6}))^2+1.0*(j-double(sample{5}))^2/(double(sample{5}))^2; 
  149. if tempxy<=1 
  150. tempBins(fix(r/margin(1)+1), fix(g/margin(2)+1), fix(b/margin(3)+1))=... 
  151. tempBins(fix(r/margin(1)+1), fix(g/margin(2)+1), fix(b/margin(3)+1))+... 
  152. 1-tempxy;  
  153. end 
  154. end 
  155. end 
  156. colorHist=tempBins(:); 
  157. colorHist=colorHist./sum(colorHist,1); 
  158. end 
  159.  
  160. function sampleSet=reselect(sampleSet,weight) 
  161. % reselect process. 
  162. % new samples are selected according the corresponding weight 
  163.  
  164. % following the steps in the referred paper. 
  165. N=size(weight,1); 
  166. cumulativeW=zeros(N,1); 
  167. cumulativeW(1)=weight(1); 
  168. % calculate the normalized cumulative probabilities 
  169. for i=2:N 
  170. cumulativeW(i)=cumulativeW(i-1)+weight(i); 
  171. end 
  172. cumulativeW=cumulativeW./cumulativeW(N); 
  173. index=zeros(N,1); 
  174. for i=1:N 
  175. % generate a uniformly distributed random number r belong [0,1] 
  176. r=rand(1); 
  177. [~,ind]=find(cumulativeW‘>=r); 
  178. if isempty(ind) 
  179. ind=N; 
  180. end 
  181. index(i)=ind(1); 
  182. end 
  183. tempSampleSet=cell(N,1); 
  184. for i=1:N 
  185. tempSampleSet{i}=sampleSet{index(i)}; 
  186. end 
  187. sampleSet=tempSampleSet; 
  188. end 
  189.  
  190. function sampleSet=propagate(sampleSet, A ) 
  191. global deltaT; 
  192. global velocityTurb; 
  193. global scaleTurb; 
  194. global shiftTurb; 
  195.  
  196. MINWIDTH=10; 
  197. MINHEIGHT=15; 
  198. N=length(sampleSet); 
  199. numSamples=0; 
  200. while numSamples<N 
  201. tempSample=sampleSet{numSamples+1}; 
  202. % convernient for expand to high order model A 
  203. tempVector=zeros(7,1); 
  204. for j=1:7 
  205. tempVector(j)=tempSample{j}; 
  206. end 
  207. tempVector=A*tempVector; 
  208. randoms=0.7*randn(7,1); 
  209. tempSample{1}=round(tempVector(1)+deltaT*tempSample{3}+randoms(1)*shiftTurb); 
  210. tempSample{2}=round(tempVector(2)+deltaT*tempSample{4}+randoms(2)*shiftTurb); 
  211. tempSample{3}=tempVector(3)+randoms(3)*velocityTurb; 
  212. tempSample{4}=tempVector(4)+randoms(4)*velocityTurb; 
  213. tempSample{5}=max(round(tempVector(5)+randoms(5)*scaleTurb),MINWIDTH); 
  214. tempSample{6}=max(round(tempVector(6)+randoms(6)*scaleTurb),MINHEIGHT); 
  215. tempSample{7}=scaleTurb; 
  216. if isValidate(tempSample)==1 
  217. numSamples=numSamples+1; 
  218. sampleSet{numSamples}=tempSample; 
  219. end 
  220. end 
  221. end 
  222.  
  223. function weight=observe(sampleSet,colorHist) 
  224. % calculate the color distribution of each sample in sampleSet 
  225. % return the weight 
  226. global SIGMA2; 
  227. N=length(sampleSet); 
  228. weight=zeros(N,1); 
  229. for i=1:N 
  230. sampleHist=calColorHist(sampleSet{i}); 
  231. tempHist=sqrt(sampleHist.*colorHist); 
  232. rho=sum(tempHist,1);% Bhattacharyya cofficient 
  233. weight(i)=exp(-(1-rho)/SIGMA2); 
  234. end 
  235. weight=weight./sum(weight,1); 
  236. end 
  237.  
  238. function sample=estimate(sampleSet,weight) 
  239. % estimate the mean state of the set 
  240. N=length(sampleSet); 
  241. sample=cell(7,1); 
  242. for i=1:7 
  243. sample{i}=0; 
  244. for j=1:N 
  245. tempSample=sampleSet{j}; 
  246. sample{i}=sample{i}+weight(j)*tempSample{i}; 
  247. end 
  248. end 
  249. sample{1}=round(sample{1}); 
  250. sample{2}=round(sample{2}); 
  251. sample{5}=round(sample{5}); 
  252. sample{6}=round(sample{6}); 
  253. end 
  254.  
  255. function validate=isValidate(sample) 
  256. % change the validation of the sample 
  257. % if the coordinate of the corner exceed the image edge,then the sample is not rightful 
  258. global curFrame; 
  259. [m,n,~]=size(curFrame); 
  260. x=sample{1}; 
  261. y=sample{2}; 
  262. hx=sample{5}; 
  263. hy=sample{6}; 
  264. validate=0; %false 
  265. 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 
  266. validate=1; 
  267. end 
  268.  
  269. end 
  270.  

参考文献

  • 基于粒子滤波的物体跟踪

  • 基于粒子滤波器的目标跟踪算法及实现

基于粒子滤波的目标追踪


联系
我们