octave_packages/geometry-1.5.0/geom2d/minDistancePoints.m

   1 %% Copyright (c) 2011, INRA
   2 %% 2004-2011, David Legland <david.legland@grignon.inra.fr>
   3 %% 2011 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
   4 %%
   5 %% All rights reserved.
   6 %% (simplified BSD License)
   7 %%
   8 %% Redistribution and use in source and binary forms, with or without
   9 %% modification, are permitted provided that the following conditions are met:
  10 %%
  11 %% 1. Redistributions of source code must retain the above copyright notice, this
  12 %%    list of conditions and the following disclaimer.
  13 %%
  14 %% 2. Redistributions in binary form must reproduce the above copyright notice,
  15 %%    this list of conditions and the following disclaimer in the documentation
  16 %%    and/or other materials provided with the distribution.
  17 %%
  18 %% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  19 %% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  20 %% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  21 %% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
  22 %% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  23 %% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  24 %% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  25 %% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  26 %% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  27 %% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  28 %% POSSIBILITY OF SUCH DAMAGE.
  29 %%
  30 %% The views and conclusions contained in the software and documentation are
  31 %% those of the authors and should not be interpreted as representing official
  32 %% policies, either expressed or implied, of copyright holder.
  33
  34 %% -*- texinfo -*-
  35 %% @deftypefn {Function File} {@var{dist} = } minDistancePoints (@var{pts})
  36 %% @deftypefnx {Function File} {@var{dist} = } minDistancePoints (@var{pts1},@var{pts2})
  37 %% @deftypefnx {Function File} {@var{dist} = } minDistancePoints (@dots{},@var{norm})
  38 %% @deftypefnx {Function File} {[@var{dist} @var{i} @var{j}] = } minDistancePoints (@var{pts1}, @var{pts2}, @dots{})
  39 %% @deftypefnx {Function File} {[@var{dist} @var{j}] = } minDistancePoints (@var{pts1}, @var{pts2}, @dots{})
  40 %% Minimal distance between several points.
  41 %%
  42 %%   Returns the minimum distance between all couple of points in @var{pts}. @var{pts} is
  43 %%   an array of [NxND] values, N being the number of points and ND the
  44 %%   dimension of the points.
  45 %%
  46 %%   Computes for each point in @var{pts1} the minimal distance to every point of
  47 %%   @var{pts2}. @var{pts1} and @var{pts2} are [NxD] arrays, where N is the number of points,
  48 %%   and D is the dimension. Dimension must be the same for both arrays, but
  49 %%   number of points can be different.
  50 %%   The result is an array the same length as @var{pts1}.
  51 %%
  52 %%   When @var{norm} is provided, it uses a user-specified norm. @var{norm}=2 means euclidean norm (the default),
  53 %%   @var{norm}=1 is the Manhattan (or "taxi-driver") distance.
  54 %%   Increasing @var{norm} growing up reduces the minimal distance, with a limit
  55 %%   to the biggest coordinate difference among dimensions.
  56 %%
  57 %%
  58 %%   Returns indices @var{i} and @var{j} of the 2 points which are the closest. @var{dist}
  59 %%   verifies relation:
  60 %%   @var{dist} = distancePoints(@var{pts}(@var{i},:), @var{pts}(@var{j},:));
  61 %%
  62 %%   If only 2 output arguments are given, it returns the indices of points which are the closest. @var{j} has the
  63 %%   same size as @var{dist}. for each I It verifies the relation :
  64 %%   @var{dist}(I) = distancePoints(@var{pts1}(I,:), @var{pts2}(@var{J},:));
  65 %%
  66 %%
  67 %%   Examples:
  68 %%
  69 %% @example
  70 %%   % minimal distance between random planar points
  71 %%       points = rand(20,2)*100;
  72 %%       minDist = minDistancePoints(points);
  73 %%
  74 %%   % minimal distance between random space points
  75 %%       points = rand(30,3)*100;
  76 %%       [minDist ind1 ind2] = minDistancePoints(points);
  77 %%       minDist
  78 %%       distancePoints(points(ind1, :), points(ind2, :))
  79 %%   % results should be the same
  80 %%
  81 %%   % minimal distance between 2 sets of points
  82 %%       points1 = rand(30,2)*100;
  83 %%       points2 = rand(30,2)*100;
  84 %%       [minDists inds] = minDistancePoints(points1, points2);
  85 %%       minDists(10)
  86 %%       distancePoints(points1(10, :), points2(inds(10), :))
  87 %%   % results should be the same
  88 %% @end example
  89 %%
  90 %%   @seealso{points2d, distancePoints}
  91 %% @end deftypefn
  92
  93 function varargout = minDistancePoints(p1, varargin)
  94
  95   %% Initialisations
  96
  97   % default norm (euclidean)
  98   n = 2;
  99
 100   % flag for processing of all points
 101   allPoints = false;
 102
 103   % process input variables
 104   if isempty(varargin)
 105       % specify only one array of points, not the norm
 106       p2 = p1;
 107
 108   elseif length(varargin)==1
 109       var = varargin{1};
 110       if length(var)>1
 111           % specify two arrays of points
 112           p2  = var;
 113           allPoints = true;
 114       else
 115           % specify array of points and the norm
 116           n   = var;
 117           p2  = p1;
 118       end
 119
 120   else
 121       % specify two array of points and the norm
 122       p2  = varargin{1};
 123       n   = varargin{2};
 124       allPoints = true;
 125   end
 126
 127
 128   % number of points in each array
 129   n1  = size(p1, 1);
 130   n2  = size(p2, 1);
 131
 132   % dimension of points
 133   d   = size(p1, 2);
 134
 135
 136   %% Computation of distances
 137
 138   % allocate memory
 139   dist = zeros(n1, n2);
 140
 141   % different behaviour depending on the norm used
 142   if n==2
 143       % Compute euclidian distance. this is the default case
 144       % Compute difference of coordinate for each pair of point ([n1*n2] array)
 145       % and for each dimension. -> dist is a [n1*n2] array.
 146       % in 2D: dist = dx.*dx + dy.*dy;
 147       for i=1:d
 148           dist = dist + (repmat(p1(:,i), [1 n2])-repmat(p2(:,i)', [n1 1])).^2;
 149       end
 150
 151       % compute minimal distance:
 152       if ~allPoints
 153           % either on all couple of points
 154           mat = repmat((1:n1)', [1 n1]);
 155           ind = mat < mat';
 156           [minSqDist ind] = min(dist(ind));
 157       else
 158           % or for each point of P1
 159           [minSqDist ind] = min(dist, [], 2);
 160       end
 161
 162       % convert squared distance to distance
 163       minDist = sqrt(minSqDist);
 164   elseif n==inf
 165       % infinite norm corresponds to maximum absolute value of differences
 166       % in 2D: dist = max(abs(dx) + max(abs(dy));
 167       for i=1:d
 168           dist = max(dist, abs(p1(:,i)-p2(:,i)));
 169       end
 170   else
 171       % compute distance using the specified norm.
 172       % in 2D: dist = power(abs(dx), n) + power(abs(dy), n);
 173       for i=1:d
 174           dist = dist + power((abs(repmat(p1(:,i), [1 n2])-repmat(p2(:,i)', [n1 1]))), n);
 175       end
 176
 177       % compute minimal distance
 178       if ~allPoints
 179           % either on all couple of points
 180           mat = repmat((1:n1)', [1 n1]);
 181           ind = mat < mat';
 182           [minSqDist ind] = min(dist(ind));
 183       else
 184           % or for each point of P1
 185           [minSqDist ind] = min(dist, [], 2);
 186       end
 187
 188       % convert squared distance to distance
 189       minDist = power(minSqDist, 1/n);
 190   end
 191
 192
 193
 194   if ~allPoints
 195       % convert index in array to row ad column subindices.
 196       % This uses the fact that index are sorted in a triangular matrix,
 197       % with the last index of each column being a so-called triangular
 198       % number
 199       ind2 = ceil((-1+sqrt(8*ind+1))/2);
 200       ind1 = ind - ind2*(ind2-1)/2;
 201       ind2 = ind2 + 1;
 202   end
 203
 204
 205   %% format output parameters
 206
 207   % format output depending on number of asked parameters
 208   if nargout<=1
 209       varargout{1} = minDist;
 210   elseif nargout==2
 211       % If two arrays are asked, 'ind' is an array of indices, one for each
 212       % point in var{pts}1, corresponding to the result in minDist
 213       varargout{1} = minDist;
 214       varargout{2} = ind;
 215   elseif nargout==3
 216       % If only one array is asked, minDist is a scalar, ind1 and ind2 are 2
 217       % indices corresponding to the closest points.
 218       varargout{1} = minDist;
 219       varargout{2} = ind1;
 220       varargout{3} = ind2;
 221   end
 222
 223 endfunction
 224
 225 %!test
 226 %!  pts = [50 10;40 60;30 30;20 0;10 60;10 30;0 10];
 227 %!  assert (minDistancePoints(pts), 20);
 228
 229 %!test
 230 %!  pts = [10 10;25 5;20 20;30 20;10 30];
 231 %!  [dist ind1 ind2] = minDistancePoints(pts);
 232 %!  assert (10, dist, 1e-6);
 233 %!  assert (3, ind1, 1e-6);
 234 %!  assert (4, ind2, 1e-6);
 235
 236 %!test
 237 %!  pts = [0 80;10 60;20 40;30 20;40 0;0 0;100 0;0 100;0 -10;-10 -20];
 238 %!  assert (minDistancePoints([40 50], pts), 10*sqrt(5), 1e-6);
 239 %!  assert (minDistancePoints([25 30], pts), 5*sqrt(5), 1e-6);
 240 %!  assert (minDistancePoints([30 40], pts), 10, 1e-6);
 241 %!  assert (minDistancePoints([20 40], pts), 0, 1e-6);
 242
 243 %!test
 244 %!  pts1 = [40 50;25 30;40 20];
 245 %!  pts2 = [0 80;10 60;20 40;30 20;40 0;0 0;100 0;0 100;0 -10;-10 -20];
 246 %!  res = [10*sqrt(5);5*sqrt(5);10];
 247 %!  assert (minDistancePoints(pts1, pts2), res, 1e-6);
 248
 249 %!test
 250 %!  pts = [50 10;40 60;40 30;20 0;10 60;10 30;0 10];
 251 %!  assert (minDistancePoints(pts, 1), 30, 1e-6);
 252 %!  assert (minDistancePoints(pts, 100), 20, 1e-6);
 253
 254 %!test
 255 %!  pts = [0 80;10 60;20 40;30 20;40 0;0 0;100 0;0 100;0 -10;-10 -20];
 256 %!  assert (minDistancePoints([40 50], pts, 2), 10*sqrt(5), 1e-6);
 257 %!  assert (minDistancePoints([25 30], pts, 2), 5*sqrt(5), 1e-6);
 258 %!  assert (minDistancePoints([30 40], pts, 2), 10, 1e-6);
 259 %!  assert (minDistancePoints([20 40], pts, 2), 0, 1e-6);
 260 %!  assert (minDistancePoints([40 50], pts, 1), 30, 1e-6);
 261 %!  assert (minDistancePoints([25 30], pts, 1), 15, 1e-6);
 262 %!  assert (minDistancePoints([30 40], pts, 1), 10, 1e-6);
 263 %!  assert (minDistancePoints([20 40], pts, 1), 0, 1e-6);
 264
 265 %!test
 266 %!  pts1 = [40 50;25 30;40 20];
 267 %!  pts2 = [0 80;10 60;20 40;30 20;40 0;0 0;100 0;0 100;0 -10;-10 -20];
 268 %!  res1 = [10*sqrt(5);5*sqrt(5);10];
 269 %!  assert (minDistancePoints(pts1, pts2, 2), res1, 1e-6);
 270 %!  res2 = [30;15;10];
 271 %!  assert (minDistancePoints(pts1, pts2, 1), res2);
 272
 273 %!test
 274 %!  pts1    = [40 50;20 30;40 20];
 275 %!  pts2    = [0 80;10 60;20 40;30 20;40 0;0 0;100 0;0 100;0 -10;-10 -20];
 276 %!  dists0  = [10*sqrt(5);10;10];
 277 %!  inds1   = [3;3;4];
 278 %!  [minDists inds] = minDistancePoints(pts1, pts2);
 279 %!  assert (dists0, minDists);
 280 %!  assert (inds1, inds);