Code covered by the BSD License

# geom2d

### David Legland (view profile)

13 Jun 2005 (Updated )

Geometry library for matlab. Performs geometric computations on points, lines, circles, polygons...

minDistancePoints(p1, varargin)
```function varargout = minDistancePoints(p1, varargin)
%MINDISTANCEPOINTS Minimal distance between several points
%
%   DIST = minDistancePoints(PTS)
%   Returns the minimum distance between all couple 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
%
%   points2d, distancePoints, nndist, hausdorffDistance
%
%   ---------
%   author: David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 15/06/2004.
%

% HISTORY:
% 22/06/2005 compute sqrt only at the end (faster), and change behaviour
%   for 2 inputs: compute min distance for each point in PTS1.
%   Also add support for different norms.
% 15/08/2005 make difference when 1 array or 2 arrays of points
% 25/10/2006 also returns indices of closest points
% 30/10/2006 generalize to points of any dimension

%% Initialisations

% default norm (euclidean)
n = 2;

% flag for processing of all points
allPoints = false;

% process input variables
if isempty(varargin)
% specify only one array of points, not the norm
p2 = p1;

elseif length(varargin) == 1
var = varargin{1};
if length(var) > 1
% specify two arrays of points
p2  = var;
allPoints = true;
else
% specify array of points and the norm
n   = var;
p2  = p1;
end

else
% specify two array of points and the norm
p2  = varargin{1};
n   = varargin{2};
allPoints = true;
end

% number of points in each array
n1  = size(p1, 1);
n2  = size(p2, 1);

% dimension of points
d   = size(p1, 2);

%% Computation of distances

% allocate memory
dist = zeros(n1, n2);

% different behaviour depending on the norm used
if n == 2
% Compute euclidian distance (default case).
% Compute difference of coordinate for each pair of point and for each
% dimension. Result "dist" is a n1-by-n2 array.
% in 2D: dist = dx.*dx + dy.*dy;
for i = 1:d
dist = dist + (repmat(p1(:,i), [1 n2])-repmat(p2(:,i)', [n1 1])).^2;
end

% compute minimal distance:
if ~allPoints
% either on all couple of points
mat = repmat((1:n1)', [1 n1]);
ind = mat < mat';
[minSqDist, ind] = min(dist(ind));
else
% or for each point of P1
[minSqDist, ind] = min(dist, [], 2);
end

% convert squared distance to distance
minDist = sqrt(minSqDist);

elseif n == inf
% infinite norm corresponds to maximum absolute value of differences
% in 2D: dist = max(abs(dx) + max(abs(dy));
for i = 1:d
dist = max(dist, abs(p1(:,i)-p2(:,i)));
end

else
% compute distance using the specified norm.
% in 2D: dist = power(abs(dx), n) + power(abs(dy), n);
for i = 1:d
dist = dist + power((abs(repmat(p1(:,i), [1 n2])-repmat(p2(:,i)', [n1 1]))), n);
end

% compute minimal distance
if ~allPoints
% either on all couple of points
mat = repmat((1:n1)', [1 n1]);
ind = mat < mat';
[minSqDist, ind] = min(dist(ind));
else
% or for each point of P1
[minSqDist, ind] = min(dist, [], 2);
end

% convert squared distance to distance
minDist = power(minSqDist, 1/n);

end

if ~allPoints
% convert index in array to row and column subindices.
% This uses the fact that index are sorted in a triangular matrix,
% with the last index of each column being a so-called triangular
% number
ind2 = ceil((-1+sqrt(8*ind+1))/2);
ind1 = ind - ind2*(ind2-1)/2;
ind2 = ind2 + 1;
end

%% format output parameters

% format output depending on number of asked parameters
if nargout <= 1
varargout{1} = minDist;

elseif nargout == 2
% If two arrays are asked, 'ind' is an array of indices, one for each
% point in PTS1, corresponding to the result in minDist
varargout{1} = minDist;
varargout{2} = ind;

elseif nargout == 3
% If only one array is asked, minDist is a scalar, ind1 and ind2 are 2
% indices corresponding to the closest points.
varargout{1} = minDist;
varargout{2} = ind1;
varargout{3} = ind2;
end
```