MINDISTANCEPOINTS 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 %MINDISTANCEPOINTS 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 % Examples: 0036 % % minimal distance between random planar points 0037 % points = rand(20,2)*100; 0038 % minDist = minDistancePoints(points); 0039 % 0040 % % minimal distance between random space points 0041 % points = rand(30,3)*100; 0042 % [minDist ind1 ind2] = minDistancePoints(points); 0043 % minDist 0044 % distancePoints(points(ind1, :), points(ind2, :)) 0045 % % results should be the same 0046 % 0047 % % minimal distance between 2 sets of points 0048 % points1 = rand(30,2)*100; 0049 % points2 = rand(30,2)*100; 0050 % [minDists inds] = minDistancePoints(points1, points2); 0051 % minDists(10) 0052 % distancePoints(points1(10, :), points2(inds(10), :)) 0053 % % results should be the same 0054 % 0055 % % Find the (approximated) orthogonal projection onto an ellipse 0056 % elli = [50 50 40 20 30]; 0057 % poly = ellipseToPolygon(elli, 200); 0058 % figure; axis equal; axis([0 100 0 100]); hold on; 0059 % drawPolygon(poly, 'k') 0060 % pts = [20 20; 50 20; 80 30]; 0061 % [dists, inds] = minDistancePoints(pts, poly); 0062 % drawPoint(pts, 'bo'); 0063 % drawPoint(poly(inds,:), 'ko'); 0064 % drawEdge([pts poly(inds,:)], 'k') 0065 % 0066 % 0067 % See also 0068 % points2d, distancePoints, nndist, findClosestPoint, hausdorffDistance 0069 % 0070 0071 % ------ 0072 % Author: David Legland 0073 % E-mail: david.legland@inrae.fr 0074 % Created: 2004-06-15 0075 % Copyright 2004-2024 INRAE - Cepia Software Platform 0076 0077 %% Initialisations 0078 0079 % default norm (euclidean) 0080 n = 2; 0081 0082 % a single array is given 0083 one_array = true; 0084 0085 % process input variables 0086 if nargin == 1 0087 % specify only one array of points, not the norm 0088 p2 = p1; 0089 elseif nargin == 2 0090 if isscalar(varargin{1}) 0091 % specify array of points and the norm 0092 n = varargin{1}; 0093 p2 = p1; 0094 else 0095 % specify two arrays of points 0096 p2 = varargin{1}; 0097 one_array = false; 0098 end 0099 elseif nargin == 3 0100 % specify two array of points and the norm 0101 p2 = varargin{1}; 0102 n = varargin{2}; 0103 one_array = false; 0104 else 0105 error ('Wrong number of input arguments'); 0106 end 0107 0108 % number of points in each array 0109 n1 = size(p1, 1); 0110 n2 = size(p2, 1); 0111 0112 % dimensionality of points 0113 d = size(p1, 2); 0114 0115 0116 %% Computation of distances 0117 0118 % allocate memory 0119 dist = zeros(n1, n2); 0120 0121 % Compute difference of coordinate for each pair of point (n1-by-n2 array) 0122 % and for each dimension. -> dist is a n1-by-n2 array. 0123 % in 2D: dist = dx.*dx + dy.*dy; 0124 if n == inf 0125 % infinite norm corresponds to maximum absolute value of differences 0126 % in 2D: dist = max(abs(dx) + max(abs(dy)); 0127 for i = 1:d 0128 dist = max(dist, abs(bsxfun(@minus, p1(:,i), p2(:,i)'))); 0129 end 0130 else 0131 for i = 1:d 0132 dist = dist + abs(bsxfun(@minus, p1(:,i), p2(:,i)')).^n; 0133 end 0134 end 0135 0136 % compute minimum distance, and indices 0137 if ~one_array 0138 % If two array of points where given 0139 [minSqDist, ind] = min(dist, [], 2); 0140 minDist = power(minSqDist, 1/n); 0141 [ind2, ind1] = ind2sub([n1 n2], ind); 0142 0143 else 0144 % A single array was given 0145 dist = dist + diag(inf(n1,1)); % remove zeros from diagonal 0146 dist = dist(tril(true(n1, n1))); 0147 [minSqDist, ind] = min(dist); % index on packed lower triangular matrix 0148 minDist = power(minSqDist, 1/n); 0149 0150 [ind2, ind1] = ind2sub_tril(n1, ind); 0151 ind2 = ind2(1); 0152 ind1 = ind1(1); 0153 ind = sub2ind([n1 n1], ind2, ind1); 0154 end 0155 0156 0157 %% format output parameters 0158 0159 % format output depending on number of asked parameters 0160 if nargout <= 1 0161 varargout{1} = minDist; 0162 0163 elseif nargout == 2 0164 % If two arrays are asked, 'ind' is an array of indices of p2, one for each 0165 % point in p1, corresponding to the result in minDist 0166 varargout{1} = minDist; 0167 varargout{2} = ind; 0168 0169 elseif nargout == 3 0170 % If only one array is asked, minDist is a scalar, ind1 and ind2 are 2 0171 % indices corresponding to the closest points. 0172 varargout{1} = minDist; 0173 varargout{2} = ind1; 0174 varargout{3} = ind2; 0175 end 0176 0177 end 0178 0179 function [r, c] = ind2sub_tril (N, idx) 0180 % Convert a linear index to subscripts of a triangular matrix. 0181 % 0182 % [r, c] = ind2sub_tril (N, idx) 0183 % 0184 % An example of triangular matrix linearly indexed follows 0185 % 0186 % N = 4; 0187 % A = -repmat (1:N,N,1); 0188 % A = [A repmat (diagind, N,1) - A.']; 0189 % A = tril(A) 0190 % => A = 0191 % 1 0 0 0 0192 % 2 5 0 0 0193 % 3 6 8 0 0194 % 4 7 9 10 0195 % 0196 % The following example shows how to convert the linear index `6' in 0197 % the 4-by-4 matrix of the example into a subscript. 0198 % 0199 % [r, c] = ind2sub_tril (4, 6) 0200 % => r = 3 0201 % c = 2 0202 % 0203 % when idx is a row or column matrix of linear indeces then r and 0204 % c have the same shape as idx. 0205 % 0206 0207 endOfRow = 0.5 * (1:N) .* (2*N:-1:N + 1); 0208 c = zeros(size(endOfRow)); 0209 for i = 1:length(endOfRow) 0210 ind = find(endOfRow <= idx - 1, 1, 'last') + 1; 0211 if isempty(ind) 0212 ind = 1; 0213 end 0214 c(i) = ind; 0215 end 0216 0217 r = N - endOfRow(c) + idx ; 0218 0219 end