基于粒子滤波的目标追踪

particle filter

object tracking

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

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

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

参考文献

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

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