물체 추적은 은닉 변수를 가진 Markov모델에 Bayesian 추론을 적용하는 것이다$^{주1}$. 시간 t 프레임에서 관찰 이미지(observation) 집합 $Y_t$={y1,y2,...,yt}가 주어진다고 할 때, 은닉 상태 변수 $x_t$는 recursive하게 추정할수 있다.
여기서, p(xt|xt-1)은 두 연속 상태 사이의 시스템 동적 모델이다. p(yt|xt)는 상태 xt에서 관찰(observation) 모델의 우도(likelihood)를 나타낸다.
이 식의 우변에서 시간 인덱스 $t-1$을 $t$로($t$는 $t+1$로) 바꾸고 좌변의 $p(x_t|Y_t)$를 다시 대입하면 시간 인덱스만 진전된 동일한 표현이 되므로 위 식은 반복해서 적용 가능한 재귀적(recursive) 표현이다.
Particle filter로 위 식을 표현하면, t프레임까지의 모든 관찰값이 주어질 때 추적 물체의 최적 상태는, 시간 t에서 N개의 샘플에서의 최대 사후확률을 추정하는 것이다.
이 식은 재귀적 형태를 포함하고 있다. $p(x_t^i|x_{t-1})$에 따라 $t-1$시간까지 추정한 상태 값 $x_{t-1}$의 근방에서 $x_t^i$를 확률적으로 발생시키고 각각의 파티클 $x_t^i$에서 관찰 값 $p(y_t^i|x_t^i)$를 얻는다. 시간 $t$에서의 새로운 상태는 두 확률 곱의 최대 값을 주는 $x_t^i$가 된다.
Particle filter의 simulation 예이 식의 우변에서 시간 인덱스 $t-1$을 $t$로($t$는 $t+1$로) 바꾸고 좌변의 $p(x_t|Y_t)$를 다시 대입하면 시간 인덱스만 진전된 동일한 표현이 되므로 위 식은 반복해서 적용 가능한 재귀적(recursive) 표현이다.
Particle filter로 위 식을 표현하면, t프레임까지의 모든 관찰값이 주어질 때 추적 물체의 최적 상태는, 시간 t에서 N개의 샘플에서의 최대 사후확률을 추정하는 것이다.
이 식은 재귀적 형태를 포함하고 있다. $p(x_t^i|x_{t-1})$에 따라 $t-1$시간까지 추정한 상태 값 $x_{t-1}$의 근방에서 $x_t^i$를 확률적으로 발생시키고 각각의 파티클 $x_t^i$에서 관찰 값 $p(y_t^i|x_t^i)$를 얻는다. 시간 $t$에서의 새로운 상태는 두 확률 곱의 최대 값을 주는 $x_t^i$가 된다.
아래 코드는 1차원 운동을 하는 물체를 tracking하는 pf를 보여준다.
움직이는 물체의 운동은 2개의 sine wave와 하나의 linear motion형태의 궤적이 결합되어 있으며 결합부에 모션 불연속부가 존재한다.
불연속이 있어도 pf는 물체를 놓치지 않고 성공적으로 대상을 추적하는 것을 보여준다.
-Green dots: particles (150개)
-Red line: true state
-Blue line: estimated state
-모션 불연속(red line)이 있는 위치가 2군데 있음(위치 데이터는 크게 삼등분 가능)
-특정 위치에서의 거동 설명을 위해 13frame, 53frame 위치가 blue 색의 수직 라인으로 표시되고, 아래 그림에서 설명.
위 그림은 100프레임 동안 물체를 추적한 결과 이다. 수직 위치를 추적하는 것이며, 약 25프레임 근방과 49프레임 근방에서 모션 불연속이 존재한다.
particle들은 운동물체의 궤적을 잘 따라 가고 있다.
13프레임의 푸른색 수직선을 보면 물체의 궤적이 연속곡선으로 부드럽게 이동한다. 따라서, particle들은 mean주위에 집중되어 분포한다.
추정된 현재 위치(녹색점) 근방에서 측정치(measurement) 값도 크므로 큰 변화없이 부드러운 추적이 발생한다.
%
% 1D particle filter tracking a height position
% Summarized by funMV, 12/2013
% [Ref.] T. Svoboda, J. Kybic, V. Hlavac, Image processing, Analysis,
% and Machine Vision, A Matlab companion, Thomson press.
% pp.232
%
% Original version is in http://visionbook.felk.cvut.cz
%
clear all;
clear classes;
% 파라메터 설정
conf.steps = 100; % number of steps(100프레임 동안)
conf.sigma_1 = 0.02; % standard deviation - system
conf.N = 150; % number of samples(파티클)
conf.sigma_2 = 0.5; % standard deviation for the measurement function
conf.Cov = 0.1; % covariance matrix - sampling (scalar for 1D)
conf.alpha = 0.3; % factor for exponential forgetting - motion model
% tracking할 대상 데이터의 생성(true values)
% 크게 3부분의 데이터 생성(첨부 자료 참조)
third = round(conf.steps/4);
X(1) = normrnd(0,conf.sigma_1);
for i = 1:(third-1)
X(i+1) = X(i) + (sin(2*pi*i/third)-sin(2*pi*(i-1)/third)) + ...
normrnd(0,conf.sigma_1);
end
X(third) = sin(3*pi/2) + normrnd(0,conf.sigma_1);
for i = third:(2*third-1)
X(i+1) = X(i) + (sin(2*pi*i/third+3*pi/2)-sin(2*pi*(i-1)/third+3*pi/2)) + ...
normrnd(0,conf.sigma_1);
end
increment = -2/(conf.steps-2*third);
X(2*third) = 1 + normrnd(0,conf.sigma_1);
for i = (2*third+1):conf.steps
X(i) = X(i-1) + increment + normrnd(0,conf.sigma_1);
end
% 파티클 값들의 설정
%
S.s = normrnd( 0, conf.Cov, 1, conf.N ); % initialize samples
% 평균을 0, 분산을 conf.Cov로 한 정규분포의 샘플을 conf.N개 발생
% S.s=1x150, mean(S.s)=-0.0011
S.pi = ones(1,conf.N) / conf.N; % set uniform weights
% 처음에는 가중치를 전부 1/N의 값을 가지게 함. size(S.pi)=1x150
S.v = 0; % initialize motion model
S.est = estimation( S, conf ); % set estimated value
% S.pi*S.s' = -0.001(평균이 0에서 샘플이 발생 했고,
% 모든 가중치는 동일하므로 기대값(가중평균)은 0)
data =[];
% tracking의 실행
for step = 2:conf.steps % for all generated states
S.x = X(step); % states of the system. 추적할 실제 값
S = particle_filtering( S, conf ); % run particle filtering
% 결과 확인을 위한 데이터 저장
data = [data; step X(step) S.est abs(S.est-X(step))]; % real(true), estimation, error
end
function S = particle_filtering(S,conf)
% PARTICLE_FILTERING Particle filtering (Condensation)
% CMP Vision Algorithms http://visionbook.felk.cvut.cz
% Petr Lhotsky, Tomas Svoboda, 2007
%
% Usage: S = particle_filtering(S,conf)
% Inputs:
% S struct Structure with samples and weights.
% .x [dim x 1] State of the system. This represent the true
% state of the system. It is used for demonstration purposes
% only. It is normally unknown. dim denotes dimensionality of
% the state space.
% .s [dim x N] Matrix containing N samples (particles).
% .pi [1 x N] Weights of the samples.
% .v [dim x 1] Predicted velocity - motion model.
% .est [dim x 1] Estimated value - estimated state of the observed system.
% conf struct Structure with configuration.
% .N 1x1 Number of samples.
% .Cov [dim x dim]Noise covariance matrix.
% .alpha 1x1 Factor for the exponential forgetting of the
% motion model;
% lower values put the accent on the history and
% higher values on the actual movement; 0 turns the motion model off.[-.75ex]
% Outputs:
% S struct The same structure as for input, but at time step k+1.
%
% The filtering procedure consists of the following steps:
% Resample data S_ using importance sampling. The sampling
% procedure draws samples from the set such that samples with higher
% weights are likely to be picked more times. Since the total
% number of samples is preserved, some samples with lower weights might
% not be
% selected at all. Note that we are still using the states at time t-1.
% The importance sampling step is performed by
% function importance_sampling.
S.s = importance_sampling( S.s, S.pi ); % S.s: 정규분포(0-mean), S.pi(weight, 1/N)
% Predict the state for time t, and add noise.
% Prediction is here only a simple motion drift, see drift_samples.
% See function noisify for details of adding noise.
S.s = drift_samples( S, conf );
S.s = noisify( S, conf );
% The correction (measurement) step evaluates the likelihood of how well the samples
% This function depends on the application, see a simple
% measurement1D or a more realistic measurement2D.
S.pi = measurement( S, conf );
% normalize weights
S.pi = S.pi ./ sum(S.pi);
% Compute a pointwise state estimate from the samples. It is worth mentioning
% that the computation of the state estimate may be meaningless in
% some situations, especially for
% multimodal distributions.
S.est_old = S.est;
S.est = estimation( S, conf );
% Update the state velocity which is needed for the simple
% drift. Simple exponential forgetting is applied. The velocity
% does not need to be used in some applications.
S.v = conf.alpha*(S.est-S.est_old) + (1-conf.alpha)*S.v;
return; % end of particle filtering
function s_noise = noisify(S,conf)
%NOISIFY Noisify the samples with normal noise
%CMP Vision Algorithms http://visionbook.felk.cvut.cz
% Usage: s_noise = noisify(S,conf)
% Inputs:
% S struct Structure with samples and weights.
% conf struct Structure with configuration.
% Outputs:
% s_noise [1 x N] Vector with samples.
% Add noise to samples for a given covariance matrix conf.Cov.
% This implementation is based on .
%
%
% Based on:
% About: Statistical Pattern Recognition Toolbox
% (C) 1999-2003, Written by Vojtech Franc and Vaclav Hlavac
% <a href="http://www.cvut.cz">Czech Technical University Prague</a>
% <a href="http://www.feld.cvut.cz">Faculty of Electrical Engineering</a>
% <a href="http://cmp.felk.cvut.cz">Center for Machine Perception</a>
[dim N] = size( S.s ); % get dimension
[U,L] = eig( conf.Cov ); % compute eigenvalues and eigenvectors
s_noise = S.s + inv(U')*sqrt(L)*randn(dim,N); % dewhitening transform
% S.s를 빼고 생각해 보면, Cov행렬을 고유값분해 해서 고유벡터 U와 그 기저에서의
% 좌표 값 L을 얻었다. 좌표 값 L은 Cov행렬의 크기이다. 여기에 곱해 평균 0과 분산 1을
% 가진 정규분포에서 샘플을 N개 생성한다면 평균 0근방에서 분산 L을 가진 분포의
% 샘플을 얻는 것이다.
% U'.a는 U기저로 어떤 벡터 a를 표현하는 것이다. U'a는 U기저에서의 벡터 a의 좌표
% 값이다. 따라서,
% U'a = sqrt(L)*randn(dim,N)에서 sqrt(L)*randn(dim,N)는 U기저에서의
% a벡터의 좌표 값이다.
% a = inv(U')*sqrt(L)*randn(dim,N)는 표준 기저에서의 좌표벡터를 얻는 것이다.
% (원래 상태 벡터의 값들은 표준 기저에서의 값들이었다)
% 정리하면, 샘플 분포는 기저 U에 대한 정규분포를 갖도록 발생시키고
% 표준 기저의 벡터로 다시 변환해 주는 것이다.
return
function [w] = measurement(S,conf)
%MEASUREMENT Evaluation of the samples
%CMP Vision Algorithms http://visionbook.felk.cvut.cz
% Usage: w = measurement(S,conf)
% S struct Structure with samples and weights.
% conf struct Structure with configuration.
% .sigma_2 [1] Standard deviation of the noise.
% w [1 x N] Vector with measurements for each sample.
% Function measurement evaluates each sample according to a measurement
% function and returns a likelihood of each sample.
% The measurement function is modeled by a Gaussian centered around the
% true value and a small constant is added to represent a non-zero
% response of a real measurement function.
% This is an approximation of an ideal which should be smooth and
% peak at the true value. The non-zero observation likelihood further away from the true state
% also prevents samples from dying prematurely.
% See demo2D/measurement for a more
% realistic measurement function.
%
w = normpdf((S.s-S.x),0,conf.sigma_2) + 0.1;
return
function [x_res,freq] = importance_sampling(x,p,varargin)
% IMPORTANCE_SAMPLING importance sampling
% CMP Vision Algorithms http://visionbook.felk.cvut.cz
% Tomas Svoboda, 2007
%
% Usage: [x_res,freq] = importance_sampling(x,p,Nnew,show,fid,save)
% Inputs:
% x [dim x N] Matrix of N samples of dimension .
% p [1 x N] Probabilities of samples.
% Nnew (default N) Number of new samples.
% show (default 0) If set to 1 the sampling steps are
% graphically visualized. It works only for 1D data.
% fid (default gcf) Figure handle for the graphical visualization.
% save (default 0) Save the graphically visualized steps to disk.
% Outputs:
% x_res [dim x Nnew] Selected (sampled) samples from x.
% freq [1 x N] Frequencies of selections, giving the number of
% times a particular sample from x was selected.
Nnew = length(x);
cum_prob = cumsum(p); % cumulative probability, cum_sum=1x150
x_res = zeros( size(x,1), Nnew); % space for resampled samples, x_res=1x150(전부 0)
freq = zeros( size(x) ); % selection frequency of the old samples, freq=1x150, 전부 0
for i = 1:Nnew
% selection sample from a uniform density 0...1
r = rand; % 0과 1사이의 난수 발생
j = find( cum_prob==min(cum_prob(cum_prob>r)), 1 );
% cum_prob>r이 되는 cum_prob요소중 최소 값을 주는 index계산
% take the j-sample to the new sample set
x_res(:,i) = x(:,j); % j인덱스에 있는 값을 저장
% x_res에는 더 많이 샘플링된 j인덱스가 주는 x값이 더 여러번 저장됨.
% increment the frequency table
freq(j) = freq(j) + 1; % j인덱스의 count를 증가
end
return; % end of importance_sampling, x_res만 리턴됨.
function est = estimation(S,conf);
%ESTIMATION The best guess for the state of the system - estimated value
%CMP Vision Algorithms http://visionbook.felk.cvut.cz
% Usage: est = estimation(S,conf)
% S struct Structure with samples and weights.
% conf struct Structure with configuration.
% est [dim x 1] Value of the predicted state.
% Computes the best guess for the state of the system
% as the weighted average of the samples
% x = sum_i^N,
% which is implemented as a scalar product of vectors
% and s.
%
est = [S.pi]*[S.s]';
% 평균을 구한다고 보면 됨. 샘플이 S.s에 들어 있고, 각 샘플의 가중치가 S.pi
% 이므로 각 샘플들의 가중 합이 계산. 즉, 기대값, 평균임.
return % end of estimation
function s_drift = drift_samples(S,conf)
%DRIFT_SAMPLES Drift samples according to the velocity of the motion model
%CMP Vision Algorithms http://visionbook.felk.cvut.cz
% Usage: s_drift = drift_samples(S,conf)
% S struct Structure with samples and weights.
% conf struct Structure with configuration.
% s_drift [dim x N] Drifted samples.
% Drift samples according to the velocity model.
% This is a simple demonstration function.
% A more general predictor is naturally possible.
%
s_drift = S.s + S.v;
return
(주1) 은닉변수(hidden variable)는 시스템에 숨겨져 있는 상태 변수(vector) $x$를 의미한다. 즉 관찰 값 $Y$를 이용하여 숨겨진 $x$를 추정하는 문제이다. Markov 모델이란 측정 값이 서로 독립적이지 않은 연속 값의 체인으로 발생한다는 것이다. 서로 인접한 값은 상관성이 크고 멀어질 수록 상관성이 작아진다. Bayesian 추론이란 수식 자체가 Bayesian theory를 이용하여 얻어지기 때문에 이렇게 표현한다.
