octave: Identifying Points in Triangulation

 
 30.1.2 Identifying Points in Triangulation
 ------------------------------------------
 
 It is often necessary to identify whether a particular point in the
 N-dimensional space is within the Delaunay tessellation of a set of
 points in this N-dimensional space, and if so which N-simplex contains
 the point and which point in the tessellation is closest to the desired
 point.  The functions ‘tsearch’ and ‘dsearch’ perform this function in a
 triangulation, and ‘tsearchn’ and ‘dsearchn’ in an N-dimensional
 tessellation.
 
    To identify whether a particular point represented by a vector P
 falls within one of the simplices of an N-simplex, we can write the
 Cartesian coordinates of the point in a parametric form with respect to
 the N-simplex.  This parametric form is called the Barycentric
 Coordinates of the point.  If the points defining the N-simplex are
 given by N + 1 vectors ‘T(I,:)’, then the Barycentric coordinates
 defining the point P are given by
 
      P = BETA * T
 
 where BETA contains N + 1 values that together as a vector represent the
 Barycentric coordinates of the point P.  To ensure a unique solution for
 the values of BETA an additional criteria of
 
      sum (BETA) == 1
 
 is imposed, and we can therefore write the above as
 
      P - T(end, :) = BETA(1:end-1) * (T(1:end-1, :)
                      - ones (N, 1) * T(end, :)
 
 Solving for BETA we can then write
 
      BETA(1:end-1) = (P - T(end, :)) /
                      (T(1:end-1, :) - ones (N, 1) * T(end, :))
      BETA(end) = sum (BETA(1:end-1))
 
 which gives the formula for the conversion of the Cartesian coordinates
 of the point P to the Barycentric coordinates BETA.  An important
 property of the Barycentric coordinates is that for all points in the
 N-simplex
 
      0 <= BETA(I) <= 1
 
 Therefore, the test in ‘tsearch’ and ‘tsearchn’ essentially only needs
 to express each point in terms of the Barycentric coordinates of each of
 the simplices of the N-simplex and test the values of BETA.  This is
 exactly the implementation used in ‘tsearchn’.  ‘tsearch’ is optimized
 for 2-dimensions and the Barycentric coordinates are not explicitly
 formed.
 
  -- : IDX = tsearch (X, Y, T, XI, YI)
      Search for the enclosing Delaunay convex hull.
 
      For ‘T = delaunay (X, Y)’, finds the index in T containing the
      points ‘(XI, YI)’.  For points outside the convex hull, IDX is NaN.
 
DONTPRINTYET       See also: Seedelaunay XREFdelaunay, *notedelaunayn:
DONTPRINTYET       See also: Seedelaunay XREFdelaunay, Seedelaunayn

      XREFdelaunayn.
 
  -- : IDX = tsearchn (X, T, XI)
  -- : [IDX, P] = tsearchn (X, T, XI)
      Search for the enclosing Delaunay convex hull.
 
      For ‘T = delaunayn (X)’, finds the index in T containing the points
      XI.  For points outside the convex hull, IDX is NaN.
 
      If requested ‘tsearchn’ also returns the Barycentric coordinates P
      of the enclosing triangles.
 
DONTPRINTYET       See also: Seedelaunay XREFdelaunay, *notedelaunayn:
DONTPRINTYET       See also: Seedelaunay XREFdelaunay, Seedelaunayn

      XREFdelaunayn.
 
    An example of the use of ‘tsearch’ can be seen with the simple
 triangulation
 
      X = [-1; -1; 1; 1];
      Y = [-1; 1; -1; 1];
      TRI = [1, 2, 3; 2, 3, 4];
 
 consisting of two triangles defined by TRI.  We can then identify which
 triangle a point falls in like
 
      tsearch (X, Y, TRI, -0.5, -0.5)
      ⇒ 1
      tsearch (X, Y, TRI, 0.5, 0.5)
      ⇒ 2
 
 and we can confirm that a point doesn’t lie within one of the triangles
 like
 
      tsearch (X, Y, TRI, 2, 2)
      ⇒ NaN
 
    The ‘dsearch’ and ‘dsearchn’ find the closest point in a tessellation
 to the desired point.  The desired point does not necessarily have to be
 in the tessellation, and even if it the returned point of the
 tessellation does not have to be one of the vertexes of the N-simplex
 within which the desired point is found.
 
  -- : IDX = dsearch (X, Y, TRI, XI, YI)
  -- : IDX = dsearch (X, Y, TRI, XI, YI, S)
      Return the index IDX of the closest point in ‘X, Y’ to the elements
      ‘[XI(:), YI(:)]’.
 
      The variable S is accepted for compatibility but is ignored.
 
      See also: Seedsearchn XREFdsearchn, Seetsearch XREFtsearch.
 
  -- : IDX = dsearchn (X, TRI, XI)
  -- : IDX = dsearchn (X, TRI, XI, OUTVAL)
  -- : IDX = dsearchn (X, XI)
  -- : [IDX, D] = dsearchn (...)
      Return the index IDX of the closest point in X to the elements XI.
 
      If OUTVAL is supplied, then the values of XI that are not contained
      within one of the simplices TRI are set to OUTVAL.  Generally, TRI
      is returned from ‘delaunayn (X)’.
 
      See also: Seedsearch XREFdsearch, Seetsearch XREFtsearch.
 
    An example of the use of ‘dsearch’, using the above values of X, Y
 and TRI is
 
      dsearch (X, Y, TRI, -2, -2)
      ⇒ 1
 
    If you wish the points that are outside the tessellation to be
 flagged, then ‘dsearchn’ can be used as
 
      dsearchn ([X, Y], TRI, [-2, -2], NaN)
      ⇒ NaN
      dsearchn ([X, Y], TRI, [-0.5, -0.5], NaN)
      ⇒ 1
 
 where the point outside the tessellation are then flagged with ‘NaN’.