gusucode.com > 用粒子滤波算法进行跟踪的matlab代码 > gmm_utilities/gmm_distance_KLD.m
function K = gmm_distance_KLD(g1, g2, N) if nargin == 3 K = gmm_KLD_montecarlo(g1, g2, N); else K = gmm_KLD_unscented(g1, g2); end % % function K = gmm_KLD_montecarlo(g1, g2, N) % Monte Carlo approximation of KLD for Gaussian mixtures. s = gmm_samples(g1, N); w1 = gmm_evaluate(g1, s); w2 = gmm_evaluate(g2, s); %K = sum(log(w1) - log(w2)) / N; ii = find(w1 ~= 0 & w2 ~= 0); % warning: this step biases the result; we need a more precise way to compute w1 and w2 K = sum(log(w1(ii)) - log(w2(ii))) / length(ii); % % function K = gmm_KLD_unscented(g1, g2) % Unscented gmm KLD % Reference: Goldberger, An Efficient Image Similarity Measure Based on % Approximations of KL-Divergence Between Two Gaussian Mixtures, 2003. [D,N] = size(g1.x); Ds = sqrt(D); K = 0; for i=1:N Ps = Ds * matrix_square_root(g1.P(:,:,i), 1); x = repvec(g1.x(:,i), D); s = [x+Ps, x-Ps]; % unscented samples for i-th component of g1 if 0 % Golberger: int f log g w = gmm_evaluate(g2, s); K = K + g1.w(i)*sum(log(w)); else % KLD: int f (log f - log g) w1 = gmm_evaluate(g1, s); w2 = gmm_evaluate(g2, s); K = K + g1.w(i)*sum(log(w1)-log(w2)); end end K = K/(2*D); % % function R = matrix_square_root(P, type) switch type case 1 % svd decomposition, P = U*D*U' = R*R' (UDU form is also called modified Cholesky decomposition) [U,D,V] = svd(P); R = U*sqrt(D); %R = reprow(sqrt(diag(D))', size(U,1)) .* U; case 2 % cholesky decomposition (triangular), P = R*R' R = chol(P)'; case 3 % principal square root (symmetric), P = R*R R = sqrtm(P); end