distmat

 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/07

0001 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 % process inputs
0061 if nargin<1 || nargin>2;
0062    error('Call distmat with 1 or 2 args');
0063 end
0064 if ischar(xy) && strcmp(xy,'UNIT_TEST'); do_unit_test; return; end
0065 
0066 
0067 [n,dims] = size(xy);
0068 numels = n*n*dims;
0069 opt = 2; if numels > 5e4, opt = 3; elseif n < 20, opt = 1; end
0070 for var = varargin
0071     if length(var{1}) == 1
0072         opt = max(1, min(4, round(abs(var{1}))));
0073     else
0074         error('Invalid input argument.');
0075     end
0076 end
0077 
0078 % distance matrix calculation options
0079 switch opt
0080     case 1 % half as many computations (symmetric upper triangular property)
0081         [k,kk] = find(triu(ones(n),1));
0082         dmat = zeros(n);
0083         dmat(k+n*(kk-1)) = sqrt(sum((xy(k,:) - xy(kk,:)).^2,2));
0084         dmat(kk+n*(k-1)) = dmat(k+n*(kk-1));
0085     case 2 % fully vectorized calculation (very fast for medium inputs)
0086         a = reshape(xy,1,n,dims);
0087         b = reshape(xy,n,1,dims);
0088         dmat = sqrt(sum((a(ones(n,1),:,:) - b(:,ones(n,1),:)).^2,3));
0089     case 3 % partially vectorized (smaller memory requirement for large inputs)
0090         dmat = zeros(n,n);
0091         for k = 1:n
0092             dmat(k,:) = sqrt(sum((xy(k*ones(n,1),:) - xy).^2,2));
0093         end
0094     case 4 % another compact method, generally slower than the others
0095         a = (1:n);
0096         b = a(ones(n,1),:);
0097         dmat = reshape(sqrt(sum((xy(b,:) - xy(b',:)).^2,2)),n,n);
0098 end
0099 
0100 function do_unit_test
0101    a = 1:3; b=4:6;
0102    xy= [a;b];
0103    out = norm(a-b)*(1-eye(2));
0104    unit_test_cmp('distmat 0', distmat(xy), out, 1e-14); 
0105    unit_test_cmp('distmat 1', distmat(xy,1), out, 1e-14); 
0106    unit_test_cmp('distmat 2', distmat(xy,2), out, 1e-14); 
0107    unit_test_cmp('distmat 3', distmat(xy,3), out, 1e-14); 
0108    unit_test_cmp('distmat 4', distmat(xy,4), out, 1e-14); 
0109