function M = slmetric_pw(X1, X2, mtype, varargin)
%SLMETRIC_PW Compute the metric between column vectors pairwisely
%
% [ Syntax ]
% - M = slmetric_pw(X1, X2, mtype);
% - M = slmetric_pw(X1, X2, mtype, ...);
%
% [ Arguments ]
% - X1, X2: the sample matrices
% - mtype: the string indicating the type of metric
% - M: the resulting metric matrix
%
% [ Description ]
% - M = slmetric_pw(X1, X2, mtype) Computes the metrics between
% column vectors of X1 and X2 pairwisely, using the metric
% specified by mtype.
%
% Both X1 and X2 are matrices with each column representing a
% sample. X1 and X2 should have the same number of rows. Suppose
% the size of X1 is d x n1, and the size of X2 is d x n2. Then
% the output metric matrix M will be of size n1 x n2, in which
% M(i, j) is the metric value between X1(:,i) and X2(:,j).
%
% - M = slmetric_pw(X1, X2, mtype, ...) Some metric types requires
% extra parameters, which should be specified in params.
%
% The supported metrics of this function are listed as follows:
% \{:
% - eucdist: Euclidean distance:
% $ ||x - y|| $
%
% - sqdist: Square of Euclidean distance:
% $ ||x - y||^2 $
%
% - dotprod: Canonical dot product:
% $ <x,y> = x^T * y $
%
% - nrmcorr: Normalized correlation (cosine angle):
% $ (x^T * y ) / (||x|| * ||y||) $
%
% - corrdist: Normalized Correlation distance
% $ 1 - nrmcorr(x, y) $
%
% - angle: Angle between two vectors (in radian)
% $ arccos (nrmcorr(x, y)) $
% - quadfrm: Quadratic form:
% $ x^T * Q * y $
% Q is specified in the 1st extra parameter
%
% - quaddiff: Quadratic form of difference:
% $ (x - y)^T * Q * (x - y) $
% Q is specified in the 1st extra parameter
%
% - cityblk: City block distance (abssum of difference)
% $ sum_i |x_i - y_i| $
%
% - maxdiff: Maximum absolute difference
% $ max_i |x_i - y_i| $
%
% - mindiff: Minimum absolute difference
% $ min_i |x_i - y_i| $
%
% - minkowski: Minkowski distance
% $ (\sum_i |x_i - y_i|^p)^(1/p) $
% The order p is specified in the 1st extra parameter
%
% - wsqdist: Weighted square of Euclidean distance
% $ \sum_i w_i (x_i - y_i)^2 $
% the weights w is specified in 1st extra parameter
% as a d x 1 column vector
%
% - hamming: Hamming distance with threshold t
% \{
% ht1 = x > t
% ht2 = y > t
% d = sum(ht1 ~= ht2)
% \}
% use threshold t as the first extra param.
% (by default, t is set to zero).
%
% - hamming_nrm: Normalized hamming distance, which equals the
% ratio of the elements that differ.
% \{
% ht1 = x > t
% ht2 = y > t
% d = sum(ht1 ~= ht2) / length(ht1)
% \}
% use threshold t as the first extra param.
% (by default, t is set to zero).
%
% - intersect: Histogram Intersection
% $ d = sum min(x, y) / min(sum(x), sum(y))$
%
% - intersectdis: Histogram intersection distance
% $ d = 1 - sum min(x, y) / min(sum(x), sum(y)) $
%
% - chisq: Chi-Square Distance
% $ d = sum (x(i) - y(i))^2/(2 * (x(i)+y(i))) $
%
% - kldiv: Kull-back Leibler divergence
% $ d = sum x(i) log (x(i) / y(i)) $
%
% - jeffrey: Jeffrey divergence
% $ d = KL(h1, (h1+h2)/2) + KL(h2, (h1+h2)/2) $
% \:}
%
% [ Remarks ]
% - Both X1 and X2 should be a matrix of numeric values, except
% for case when metric type is 'hamming' or 'hamming_nrm'.
% For hamming or hamming_nrm metric, the input matrix can be logical.
%
% [ Examples ]
% - Compute different types of metrics in pairwise manner
% \{
% % prepare sample matrix
% X1 = rand(10, 100);
% X2 = rand(10, 150);
%
% % compute the euclidean distances (L2)
% % between the samples in X1 and X2
% M = slmetric_pw(X1, X2, 'eucdist');
%
% % compute the eucidean distances between the samples
% % in X1 in a pairwise manner
% M = slmetric_pw(X1, X1, 'eucdist');
%
% % compute the city block distances (L1)
% M = slmetric_pw(X1, X2, 'cityblk');
%
% % compute the normalize correlations
% M = slmetric_pw(X1, X2, 'nrmcorr');
%
% % compute hamming distances
% M = slmetric_pw(X1, X2, 'hamming', 0.5);
% M2 = slmetric_pw((X1 > 0.5), (X2 > 0.5), 'hamming');
% assert(isequal(M, M2));
% \}
%
% - Compute the parameterized metrics
% \{
% % compute weighted squared distances with user-supplied weights
% weights = rand(10, 1);
% M = slmetric_pw(X1, X2, 'wsqdist', weights);
%
% % compute quadratic distances (x-y)^T * Q (x-y)
% Q = rand(10, 10);
% M = slmetric_pw(X1, X2, 'quaddiff', Q);
%
% % compute Minkowski distance of order 3
% M = slmetric_pw(X1, X2, 'minkowski', 3);
% \}
%
% [ History ]
% - Created by Dahua Lin on Dec 06th, 2005
% - Modified by Dahua Lin on Apr 21st, 2005
% - regularize the error reporting
% - Modified by Dahua Lin on Sep 11st, 2005
% - completely rewrite the core codes based on new mex computation
% cores, and the runtime efficiency in both time and space is
% significantly increased.
% - Modified by Dahua Lin on Jul 02, 2007
% - rewrite the core computation based on the bsxfun introduced in
% MATLAB R2007a
% - rewrite the core-mex for cityblk, maxdiff, mindiff
% - introduce new metrics: corrdist, minkowski
% - Modified by Dahua Lin on Jul 30, 2007
% - Add the metric types for histograms, which are originally
% implemented in slhistmetric_pw in sltoolbox v1.
% - Modified by Dahua Lin on Aug 16, 2007
% - revise some of the help contents
%
%% parse and verify input arguments
error(nargchk(3, inf, nargin));
assert(ischar(mtype), 'sltoolbox:slmetric_pw:invalidarg', ...
'The metric type should be a string.');
if strcmp(mtype, 'hamming') || strcmp(mtype, 'hamming_nrm')
assert((isnumeric(X1) || islogical(X1)) && ndims(X1) == 2 && ...
(isnumeric(X2) || islogical(X2)) && ndims(X2) == 2, ...
'sltoolbox:slmetric_pw:invalidarg', 'X1 and X2 should be numeric or logical matrices.');
else
assert(isnumeric(X1) && ndims(X1) == 2 && isnumeric(X2) && ndims(X2) == 2, ...
'sltoolbox:slmetric_pw:invalidarg', 'X1 and X2 should be numeric matrices.');
end
assert(isa(X2, class(X1)), ...
'sltoolbox:slmetric_pw:invalidarg', 'X1 and X2 should be of the same class.');
if isempty(X1) || isempty(X2)
M = [];
return;
end
%% compute
switch mtype
case {'eucdist', 'sqdist'}
checkdim(X1, X2);
M = bsxfun(@plus, sum(X1 .* X1, 1)', (-2) * X1' * X2);
M = bsxfun(@plus, sum(X2 .* X2, 1), M);
M(M < 0) = 0;
if strcmp(mtype, 'eucdist')
M = sqrt(M);
end
case 'dotprod'
checkdim(X1, X2);
M = X1' * X2;
case {'nrmcorr', 'corrdist', 'angle'}
checkdim(X1, X2);
ns1 = sqrt(sum(X1 .* X1, 1));
ns2 = sqrt(sum(X2 .* X2, 1));
ns1(ns1 == 0) = 1;
ns2(ns2 == 0) = 1;
M = bsxfun(@times, X1' * X2, 1 ./ ns1');
M = bsxfun(@times, M, 1 ./ ns2);
switch mtype
case 'corrdist'
M = 1 - M;
case 'angle'
M = real(acos(M));
end
case 'quadfrm'
Q = varargin{1};
M = X1' * Q * X2;
case 'quaddiff'
checkdim(X1, X2);
Q = varargin{1};
M = X1' * (-(Q + Q')) * X2;
M = bsxfun(@plus, M, sum(X1 .* (Q * X1), 1)');
M = bsxfun(@plus, M, sum(X2 .* (Q * X2), 1));
case 'cityblk'
checkdim(X1, X2);
M = pwmetrics_cimp(X1, X2, int32(1));
case 'maxdiff'
checkdim(X1, X2);
M = pwmetrics_cimp(X1, X2, int32(3));
case 'mindiff'
checkdim(X1, X2);
M = pwmetrics_cimp(X1, X2, int32(2));
case 'minkowski'
checkdim(X1, X2);
pord = varargin{1};
if ~isscalar(pord)
error('sltoolbox:slmetric_pw:invalidparam', ...
'the mikowski order should be a scalar');
end
pord = cast(pord, class(X1));
M = pwmetrics_cimp(X1, X2, int32(4), pord);
case 'wsqdist'
d = checkdim(X1, X2);
w = varargin{1};
if ~isequal(size(w), [d, 1])
error('sltoolbox:slmetric_pw:invalidparam', ...
'the weights should be given as a d x 1 vector.');
end
wX2 = bsxfun(@times, X2, w);
M = bsxfun(@plus, (-2) * X1' * wX2, sum(wX2 .* X2, 1));
clear wX2;
wX1 = bsxfun(@times, X1, w);
M = bsxfun(@plus, M, sum(wX1 .* X1, 1)');
case {'hamming', 'hamming_nrm'}
checkdim(X1, X2);
if islogical(X1) && islogical(X2)
H1 = X1;
H2 = X2;
else
if isempty(varargin)
t = 0;
else
t = varargin{1};
assert(isnumeric(t) && isscalar(t), ...
'sltoolbox:slmetric_pw:invalidparam', 't should be a numeric scalar.');
end
H1 = X1 > t;
H2 = X2 > t;
end
M = pwhamming_cimp(H1, H2);
if strcmp(mtype, 'hamming_nrm')
M = M / size(H1, 1);
end
case 'intersect'
checkdim(X1, X2);
M = pwmetrics_cimp(X1, X2, int32(5));
case 'intersectdis'
checkdim(X1, X2);
M = 1 - pwmetrics_cimp(X1, X2, int32(5));
case 'chisq'
checkdim(X1, X2);
M = pwmetrics_cimp(X1, X2, int32(6));
case 'kldiv'
checkdim(X1, X2);
M = pwmetrics_cimp(X1, X2, int32(7));
case 'jeffrey'
checkdim(X1, X2);
M = pwmetrics_cimp(X1, X2, int32(8));
otherwise
error('sltoolbox:slmetric_pw:unknowntype', 'Unknown metric type %s', mtype);
end
%% Auxiliary function
function d = checkdim(X1, X2)
d = size(X1, 1);
if d ~= size(X2, 1)
error('sltoolbox:slmetric_pw:sizmismatch', ...
'X1 and X2 have different sample dimensions');
end
