ECG-Kit 1.0
(2,584 bytes)
%DISTM Compute square Euclidean distance matrix
%
% D = DISTM(A,B)
%
% INPUT
% A,B Datasets or matrices; B is optional, default B = A
%
% OUTPUT
% D Square Euclidean distance dataset or matrix
%
% DESCRIPTION
% Computation of the square Euclidean distance matrix D between two
% sets A and B. If A has M objects and B has N objects, then D is
% [M x N]. If A and B are datasets, then D is a dataset as well with
% the labels defined by the labels of A and the feature labels defined
% by the labels of B.
%
% Unlabeled objects in B are neglected, unless B is entirely unlabeled.
%
% If A is not a dataset, but a matrix of doubles then D is also not a
% dataset, but a set of doubles.
%
% NOTE
% DISTM(A,B) is equivalent to A*PROXM(B,'d',2)).
%
% SEE ALSO
% DATASETS, PROXM
% Copyright: R.P.W. Duin, r.p.w.duin@prtools.org
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands
% $Id: distm.m,v 1.9 2006/03/07 16:14:41 duin Exp $
function D = distm_new(A,B,DirectionalFeatures)
prtrace(mfilename);
if nargin < 2
B = A;
end
B = cdats(B,1);
[ma,ka] = size(A);
[mb,kb] = size(B);
if (ka ~= kb)
error('Feature sizes should be equal')
end
if( isempty(DirectionalFeatures) )
% The order of operations below is good for the accuracy.
D = ones(ma,1)*sum(B'.*B',1);
D = D + sum(A'.*A',1)'*ones(1,mb);
D = D - 2 .* (+A)*(+B)';
else
%simulo los dos wraps posibles
B1 = B(:,DirectionalFeatures) + 2*pi;
B2 = B(:,DirectionalFeatures) - 2*pi;
%mido las distancias con y sin los 2 wraps
D = ones(ma,1)*sum(B'.*B',1);
D = D + sum(A'.*A',1)'*ones(1,mb);
D = D - 2 .* (+A)*(+B)';
D1 = ones(ma,1)*sum(B1'.*B1',1);
D1 = D1 + sum(A'.*A',1)'*ones(1,mb);
D1 = D1 - 2 .* (+A)*(+B1)';
D2 = ones(ma,1)*sum(B2'.*B2',1);
D2 = D2 + sum(A'.*A',1)'*ones(1,mb);
D2 = D2 - 2 .* (+A)*(+B2)';
%me quedo con la menor distancia.
D = reshape(min([ D(:) D1(:) D2(:) ],[],2), ma, ma);
end
J = find(D<0); % Check for a numerical inaccuracy.
D(J) = zeros(size(J)); % D should be nonnegative.
if ((nargin < 2) & (ma == mb)) % take care of symmetric distance matrix
D = (D + D')/2;
D([1:ma+1:ma*ma]) = zeros(1,ma);
end
if isa(A,'prdataset') % set object and feature labels
if isa(B,'prdataset')
D = setdata(A,D,getlab(B));
else
D = setdata(A,D);
end
end
return