DISTMAT Distance matrix for a set of points
Returns the point-to-point distance between all pairs of points in XY
(similar to PDIST in the Statistics Toolbox)
DMAT = DISTMAT(XY) Calculates the distance matrix using an automatic option
DMAT = DISTMAT(XY,OPT) Uses the specified option to compute the distance matrix
[DMAT,OPT] = DISTMAT(XY) Also returns the automatic option used by the function
Inputs:
XY is an NxP matrix of coordinates for N points in P dimensions
OPT (optional) is an integer between 1 and 4 representing the chosen
method for computing the distance matrix (see note below)
Outputs:
DMAT is an NxN matrix, where the value of DMAT(i,j) corresponds to
the distance from XY(i,:) to XY(j,:)
OPT (optional) is an integer between 1 and 4 representing the method
used to compute the distance matrix (see note below)
Note:
DISTMAT contains 4 methods for computing the distance matrix
OPT=1 Usually fastest for small inputs. Takes advantage of the symmetric
property of distance matrices to perform half as many calculations
OPT=2 Usually fastest for medium inputs. Uses a fully vectorized method
OPT=3 Usually fastest for large inputs. Uses a partially vectorized
method with relatively small memory requirement
OPT=4 Another compact calculation, but usually slower than the others
Example:
% Test computation times for the options
n = [10 100 1000];
dmat = distmat(10*rand(10,3),1); % First call is always really slow
for k=1:3
for opt=1:4
tic; [dmat,opt] = distmat(10*rand(n(k),3),opt); t=toc;
disp(sprintf('n=%d, opt=%d, t=%0.6f', n(k), opt, t))
end
end
Example:
xy = 10*rand(25,2); % 25 points in 2D
dmat = distmat(xy);
figure; plot(xy(:,1),xy(:,2),'.');
for k=1:25, text(xy(k,1),xy(k,2),[' ' num2str(k)]); end
figure; imagesc(dmat); set(gca,'XTick',1:25,'YTick',1:25); colorbar
Example:
xyz = 10*rand(20,3); % 20 points in 3D
dmat = distmat(xyz);
figure; plot3(xyz(:,1),xyz(:,2),xyz(:,3),'.');
for k=1:20, text(xyz(k,1),xyz(k,2),xyz(k,3),[' ' num2str(k)]); end
figure; imagesc(dmat); set(gca,'XTick',1:20,'YTick',1:20); colorbar
Author: Joseph Kirk
Email: jdkirk630 at gmail dot com
Release: 1.0
Release Date: 5/29/070001 function [dmat,opt] = distmat(xy,varargin) 0002 % DISTMAT Distance matrix for a set of points 0003 % Returns the point-to-point distance between all pairs of points in XY 0004 % (similar to PDIST in the Statistics Toolbox) 0005 % 0006 % DMAT = DISTMAT(XY) Calculates the distance matrix using an automatic option 0007 % DMAT = DISTMAT(XY,OPT) Uses the specified option to compute the distance matrix 0008 % [DMAT,OPT] = DISTMAT(XY) Also returns the automatic option used by the function 0009 % 0010 % Inputs: 0011 % XY is an NxP matrix of coordinates for N points in P dimensions 0012 % OPT (optional) is an integer between 1 and 4 representing the chosen 0013 % method for computing the distance matrix (see note below) 0014 % 0015 % Outputs: 0016 % DMAT is an NxN matrix, where the value of DMAT(i,j) corresponds to 0017 % the distance from XY(i,:) to XY(j,:) 0018 % OPT (optional) is an integer between 1 and 4 representing the method 0019 % used to compute the distance matrix (see note below) 0020 % 0021 % Note: 0022 % DISTMAT contains 4 methods for computing the distance matrix 0023 % OPT=1 Usually fastest for small inputs. Takes advantage of the symmetric 0024 % property of distance matrices to perform half as many calculations 0025 % OPT=2 Usually fastest for medium inputs. Uses a fully vectorized method 0026 % OPT=3 Usually fastest for large inputs. Uses a partially vectorized 0027 % method with relatively small memory requirement 0028 % OPT=4 Another compact calculation, but usually slower than the others 0029 % 0030 % Example: 0031 % % Test computation times for the options 0032 % n = [10 100 1000]; 0033 % dmat = distmat(10*rand(10,3),1); % First call is always really slow 0034 % for k=1:3 0035 % for opt=1:4 0036 % tic; [dmat,opt] = distmat(10*rand(n(k),3),opt); t=toc; 0037 % disp(sprintf('n=%d, opt=%d, t=%0.6f', n(k), opt, t)) 0038 % end 0039 % end 0040 % 0041 % Example: 0042 % xy = 10*rand(25,2); % 25 points in 2D 0043 % dmat = distmat(xy); 0044 % figure; plot(xy(:,1),xy(:,2),'.'); 0045 % for k=1:25, text(xy(k,1),xy(k,2),[' ' num2str(k)]); end 0046 % figure; imagesc(dmat); set(gca,'XTick',1:25,'YTick',1:25); colorbar 0047 % 0048 % Example: 0049 % xyz = 10*rand(20,3); % 20 points in 3D 0050 % dmat = distmat(xyz); 0051 % figure; plot3(xyz(:,1),xyz(:,2),xyz(:,3),'.'); 0052 % for k=1:20, text(xyz(k,1),xyz(k,2),xyz(k,3),[' ' num2str(k)]); end 0053 % figure; imagesc(dmat); set(gca,'XTick',1:20,'YTick',1:20); colorbar 0054 % 0055 % Author: Joseph Kirk 0056 % Email: jdkirk630 at gmail dot com 0057 % Release: 1.0 0058 % Release Date: 5/29/07 0059 0060 % We need this because pdist and pdist2 are in the stats toolbox in matlab 0061 0062 % process inputs 0063 if nargin<1 || nargin>2; 0064 error('Call distmat with 1 or 2 args'); 0065 end 0066 if ischar(xy) && strcmp(xy,'UNIT_TEST'); do_unit_test; return; end 0067 0068 0069 [n,dims] = size(xy); 0070 numels = n*n*dims; 0071 opt = 2; if numels > 5e4, opt = 3; elseif n < 20, opt = 1; end 0072 for var = varargin 0073 if length(var{1}) == 1 0074 opt = max(1, min(4, round(abs(var{1})))); 0075 else 0076 error('Invalid input argument.'); 0077 end 0078 end 0079 0080 % distance matrix calculation options 0081 switch opt 0082 case 1 % half as many computations (symmetric upper triangular property) 0083 [k,kk] = find(triu(ones(n),1)); 0084 dmat = zeros(n); 0085 dmat(k+n*(kk-1)) = sqrt(sum((xy(k,:) - xy(kk,:)).^2,2)); 0086 dmat(kk+n*(k-1)) = dmat(k+n*(kk-1)); 0087 case 2 % fully vectorized calculation (very fast for medium inputs) 0088 a = reshape(xy,1,n,dims); 0089 b = reshape(xy,n,1,dims); 0090 dmat = sqrt(sum((a(ones(n,1),:,:) - b(:,ones(n,1),:)).^2,3)); 0091 case 3 % partially vectorized (smaller memory requirement for large inputs) 0092 dmat = zeros(n,n); 0093 for k = 1:n 0094 dmat(k,:) = sqrt(sum((xy(k*ones(n,1),:) - xy).^2,2)); 0095 end 0096 case 4 % another compact method, generally slower than the others 0097 a = (1:n); 0098 b = a(ones(n,1),:); 0099 dmat = reshape(sqrt(sum((xy(b,:) - xy(b',:)).^2,2)),n,n); 0100 end 0101 0102 function do_unit_test 0103 a = 1:3; b=4:6; 0104 xy= [a;b]; 0105 out = norm(a-b)*(1-eye(2)); 0106 unit_test_cmp('distmat 0', distmat(xy), out, 1e-14); 0107 unit_test_cmp('distmat 1', distmat(xy,1), out, 1e-14); 0108 unit_test_cmp('distmat 2', distmat(xy,2), out, 1e-14); 0109 unit_test_cmp('distmat 3', distmat(xy,3), out, 1e-14); 0110 unit_test_cmp('distmat 4', distmat(xy,4), out, 1e-14); 0111