Minimal distance between several points.
DIST = minDistancePoints(PTS)
Returns the minimum distance between all pairs of points in PTS. PTS
is a N-by-D array of values, N being the number of points and D the
dimension of the points.
DIST = minDistancePoints(PTS1, PTS2)
Computes for each point in PTS1 the minimal distance to every point of
PTS2. PTS1 and PTS2 are N-by-D arrays, where N is the number of points,
and D is the dimension. Dimension must be the same for both arrays, but
number of points can be different.
The result is an array the same length as PTS1.
DIST = minDistancePoints(..., NORM)
Uses a user-specified norm. NORM=2 means euclidean norm (the default),
NORM=1 is the Manhattan (or "taxi-driver") distance.
Increasing NORM growing up reduces the minimal distance, with a limit
to the biggest coordinate difference among dimensions.
[DIST, I, J] = minDistancePoints(PTS)
Returns indices I and J of the 2 points which are the closest. DIST
verifies relation:
DIST = distancePoints(PTS(I,:), PTS(J,:));
[DIST, J] = minDistancePoints(PTS1, PTS2, ...)
Also returns the indices of points which are the closest. J has the
same size as DIST. It verifies relation:
DIST(I) = distancePoints(PTS1(I,:), PTS2(J,:));
for I comprised between 1 and the number of rows in PTS1.
Examples:
% minimal distance between random planar points
points = rand(20,2)*100;
minDist = minDistancePoints(points);
% minimal distance between random space points
points = rand(30,3)*100;
[minDist ind1 ind2] = minDistancePoints(points);
minDist
distancePoints(points(ind1, :), points(ind2, :))
% results should be the same
% minimal distance between 2 sets of points
points1 = rand(30,2)*100;
points2 = rand(30,2)*100;
[minDists inds] = minDistancePoints(points1, points2);
minDists(10)
distancePoints(points1(10, :), points2(inds(10), :))
% results should be the same
% Find the (approximated) orthogonal projection onto an ellipse
elli = [50 50 40 20 30];
poly = ellipseToPolygon(elli, 200);
figure; axis equal; axis([0 100 0 100]); hold on;
drawPolygon(poly, 'k')
pts = [20 20; 50 20; 80 30];
[dists, inds] = minDistancePoints(pts, poly);
drawPoint(pts, 'bo');
drawPoint(poly(inds,:), 'ko');
drawEdge([pts poly(inds,:)], 'k')
See Also
points2d, distancePoints, nndist, findClosestPoint, hausdorffDistance0001 function varargout = minDistancePoints(p1, varargin) 0002 % Minimal distance between several points. 0003 % 0004 % DIST = minDistancePoints(PTS) 0005 % Returns the minimum distance between all pairs of points in PTS. PTS 0006 % is a N-by-D array of values, N being the number of points and D the 0007 % dimension of the points. 0008 % 0009 % DIST = minDistancePoints(PTS1, PTS2) 0010 % Computes for each point in PTS1 the minimal distance to every point of 0011 % PTS2. PTS1 and PTS2 are N-by-D arrays, where N is the number of points, 0012 % and D is the dimension. Dimension must be the same for both arrays, but 0013 % number of points can be different. 0014 % The result is an array the same length as PTS1. 0015 % 0016 % 0017 % DIST = minDistancePoints(..., NORM) 0018 % Uses a user-specified norm. NORM=2 means euclidean norm (the default), 0019 % NORM=1 is the Manhattan (or "taxi-driver") distance. 0020 % Increasing NORM growing up reduces the minimal distance, with a limit 0021 % to the biggest coordinate difference among dimensions. 0022 % 0023 % 0024 % [DIST, I, J] = minDistancePoints(PTS) 0025 % Returns indices I and J of the 2 points which are the closest. DIST 0026 % verifies relation: 0027 % DIST = distancePoints(PTS(I,:), PTS(J,:)); 0028 % 0029 % [DIST, J] = minDistancePoints(PTS1, PTS2, ...) 0030 % Also returns the indices of points which are the closest. J has the 0031 % same size as DIST. It verifies relation: 0032 % DIST(I) = distancePoints(PTS1(I,:), PTS2(J,:)); 0033 % for I comprised between 1 and the number of rows in PTS1. 0034 % 0035 % 0036 % Examples: 0037 % % minimal distance between random planar points 0038 % points = rand(20,2)*100; 0039 % minDist = minDistancePoints(points); 0040 % 0041 % % minimal distance between random space points 0042 % points = rand(30,3)*100; 0043 % [minDist ind1 ind2] = minDistancePoints(points); 0044 % minDist 0045 % distancePoints(points(ind1, :), points(ind2, :)) 0046 % % results should be the same 0047 % 0048 % % minimal distance between 2 sets of points 0049 % points1 = rand(30,2)*100; 0050 % points2 = rand(30,2)*100; 0051 % [minDists inds] = minDistancePoints(points1, points2); 0052 % minDists(10) 0053 % distancePoints(points1(10, :), points2(inds(10), :)) 0054 % % results should be the same 0055 % 0056 % % Find the (approximated) orthogonal projection onto an ellipse 0057 % elli = [50 50 40 20 30]; 0058 % poly = ellipseToPolygon(elli, 200); 0059 % figure; axis equal; axis([0 100 0 100]); hold on; 0060 % drawPolygon(poly, 'k') 0061 % pts = [20 20; 50 20; 80 30]; 0062 % [dists, inds] = minDistancePoints(pts, poly); 0063 % drawPoint(pts, 'bo'); 0064 % drawPoint(poly(inds,:), 'ko'); 0065 % drawEdge([pts poly(inds,:)], 'k') 0066 % 0067 % 0068 % See Also 0069 % points2d, distancePoints, nndist, findClosestPoint, hausdorffDistance 0070 % 0071 0072 % ------ 0073 % Author: David Legland 0074 % e-mail: david.legland@inrae.fr 0075 % Copyright 2009 INRAE - Cepia Software Platform. 0076 % created the 15/06/2004. 0077 0078 0079 % HISTORY: 0080 % 22/06/2005 compute sqrt only at the end (faster), and change behaviour 0081 % for 2 inputs: compute min distance for each point in PTS1. 0082 % Also add support for different norms. 0083 % 15/08/2005 make difference when 1 array or 2 arrays of points 0084 % 25/10/2006 also returns indices of closest points 0085 % 30/10/2006 generalize to points of any dimension 0086 % 28/08/2007 code cleanup, add comments and help 0087 % 01/02/2017 complete re-write by JuanPi Carbajal 0088 % 01/02/2017 fix bugs, update code, fix MLInt Warning (D. Legland) 0089 % 01/02/2017 remove use of vech 0090 0091 0092 %% Initialisations 0093 0094 % default norm (euclidean) 0095 n = 2; 0096 0097 % a single array is given 0098 one_array = true; 0099 0100 % process input variables 0101 if nargin == 1 0102 % specify only one array of points, not the norm 0103 p2 = p1; 0104 elseif nargin == 2 0105 if isscalar (varargin{1}) 0106 % specify array of points and the norm 0107 n = varargin{1}; 0108 p2 = p1; 0109 else 0110 % specify two arrays of points 0111 p2 = varargin{1}; 0112 one_array = false; 0113 end 0114 elseif nargin == 3 0115 % specify two array of points and the norm 0116 p2 = varargin{1}; 0117 n = varargin{2}; 0118 one_array = false; 0119 else 0120 error ('Wrong number of input arguments'); 0121 end 0122 0123 % number of points in each array 0124 n1 = size (p1, 1); 0125 n2 = size (p2, 1); 0126 0127 % dimensionality of points 0128 d = size (p1, 2); 0129 0130 0131 %% Computation of distances 0132 0133 % allocate memory 0134 dist = zeros (n1, n2); 0135 0136 % Compute difference of coordinate for each pair of point (n1-by-n2 array) 0137 % and for each dimension. -> dist is a n1-by-n2 array. 0138 % in 2D: dist = dx.*dx + dy.*dy; 0139 if n == inf 0140 % infinite norm corresponds to maximum absolute value of differences 0141 % in 2D: dist = max(abs(dx) + max(abs(dy)); 0142 for i = 1:d 0143 dist = max (dist, abs(bsxfun (@minus, p1(:,i), p2(:,i).'))); 0144 end 0145 else 0146 for i = 1:d 0147 dist = dist + abs (bsxfun (@minus, p1(:,i), p2(:,i).')).^n; 0148 end 0149 end 0150 % TODO the previous could be optimized when a single array is given (maybe!) 0151 0152 if ~one_array 0153 % If two array of points where given 0154 [minSqDist, ind] = min (dist, [], 2); 0155 minDist = power (minSqDist, 1/n); 0156 [ind2, ind1] = ind2sub ([n1 n2], ind); 0157 0158 else 0159 % A single array was given 0160 dist = dist + diag (inf (n1,1)); % remove zeros from diagonal 0161 dist = dist (tril(true(n1, n1))); 0162 [minSqDist, ind] = min (dist); % index on packed lower triangular matrix 0163 minDist = power (minSqDist, 1/n); 0164 0165 [ind2, ind1] = ind2sub_tril (n1, ind); 0166 ind2 = ind2(1); 0167 ind1 = ind1(1); 0168 ind = sub2ind ([n1 n1], ind2, ind1); 0169 end 0170 0171 0172 %% format output parameters 0173 0174 % format output depending on number of asked parameters 0175 if nargout <= 1 0176 varargout{1} = minDist; 0177 0178 elseif nargout == 2 0179 % If two arrays are asked, 'ind' is an array of indices of p2, one for each 0180 % point in p1, corresponding to the result in minDist 0181 varargout{1} = minDist; 0182 varargout{2} = ind; 0183 0184 elseif nargout == 3 0185 % If only one array is asked, minDist is a scalar, ind1 and ind2 are 2 0186 % indices corresponding to the closest points. 0187 varargout{1} = minDist; 0188 varargout{2} = ind1; 0189 varargout{3} = ind2; 0190 end 0191 0192 end 0193 0194 function [r, c] = ind2sub_tril (N, idx) 0195 % [r, c] = ind2sub_tril (N, idx) 0196 % Convert a linear index to subscripts of a trinagular matrix. 0197 % 0198 % An example of triangular matrix linearly indexed follows 0199 % 0200 % N = 4; 0201 % A = -repmat (1:N,N,1); 0202 % A = [A repmat (diagind, N,1) - A.']; 0203 % A = tril(A) 0204 % => A = 0205 % 1 0 0 0 0206 % 2 5 0 0 0207 % 3 6 8 0 0208 % 4 7 9 10 0209 % 0210 % The following example shows how to convert the linear index `6' in 0211 % the 4-by-4 matrix of the example into a subscript. 0212 % 0213 % [r, c] = ind2sub_tril (4, 6) 0214 % => r = 3 0215 % c = 2 0216 % 0217 % when idx is a row or column matrix of linear indeces then r and 0218 % c have the same shape as idx. 0219 % 0220 % See also 0221 % ind2sub 0222 0223 endofrow = 0.5 * (1:N) .* (2*N:-1:N + 1); 0224 c = zeros(size(endofrow)); 0225 for i = 1:length(endofrow) 0226 ind = find(endofrow <= idx - 1, 1, 'last') + 1; 0227 if isempty(ind) 0228 ind = 1; 0229 end 0230 c(i) = ind; 0231 % c = lookup (endofrow, idx - 1) + 1; 0232 end 0233 r = N - endofrow(c) + idx ; 0234 0235 end