GNU Octave  8.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
tsearch.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 2002-2023 The Octave Project Developers
4 //
5 // See the file COPYRIGHT.md in the top-level directory of this
6 // distribution or <https://octave.org/copyright/>.
7 //
8 // This file is part of Octave.
9 //
10 // Octave is free software: you can redistribute it and/or modify it
11 // under the terms of the GNU General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // Octave is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU General Public License for more details.
19 //
20 // You should have received a copy of the GNU General Public License
21 // along with Octave; see the file COPYING. If not, see
22 // <https://www.gnu.org/licenses/>.
23 //
24 ////////////////////////////////////////////////////////////////////////
25 
26 #if defined (HAVE_CONFIG_H)
27 # include "config.h"
28 #endif
29 
30 #include <cmath>
31 
32 #include "lo-ieee.h"
33 
34 #include "defun.h"
35 #include "error.h"
36 #include "ovl.h"
37 
39 
40 inline double max (double a, double b, double c)
41 {
42  return (a > b) ? (a > c ? a : c) : (b > c ? b : c);
43 }
44 
45 inline double min (double a, double b, double c)
46 {
47  return (a < b) ? (a < c ? a : c) : (b < c ? b : c);
48 }
49 
50 #define REF(x,k,i) x(static_cast<octave_idx_type> (elem((k), (i))) - 1)
51 
52 // for large data set the algorithm is very slow
53 // one should presort (how?) either the elements of the points of evaluation
54 // to cut down the time needed to decide which triangle contains the
55 // given point
56 
57 // e.g., build up a neighbouring triangle structure and use a simplex-like
58 // method to traverse it
59 
60 DEFUN (tsearch, args, ,
61  doc: /* -*- texinfo -*-
62 @deftypefn {} {@var{idx} =} tsearch (@var{x}, @var{y}, @var{t}, @var{xi}, @var{yi})
63 Search for the enclosing Delaunay convex hull.
64 
65 For @code{@var{t} = delaunay (@var{x}, @var{y})}, finds the index in @var{t}
66 containing the points @code{(@var{xi}, @var{yi})}. For points outside the
67 convex hull, @var{idx} is NaN.
68 @seealso{delaunay, delaunayn}
69 @end deftypefn */)
70 {
71  if (args.length () != 5)
72  print_usage ();
73 
74  const double eps = 1.0e-12;
75 
76  const ColumnVector x (args(0).vector_value ());
77  const ColumnVector y (args(1).vector_value ());
78  const Matrix elem (args(2).matrix_value ());
79  const ColumnVector xi (args(3).vector_value ());
80  const ColumnVector yi (args(4).vector_value ());
81 
82  const octave_idx_type nelem = elem.rows ();
83 
84  ColumnVector minx (nelem);
85  ColumnVector maxx (nelem);
86  ColumnVector miny (nelem);
87  ColumnVector maxy (nelem);
88  for (octave_idx_type k = 0; k < nelem; k++)
89  {
90  minx(k) = min (REF (x, k, 0), REF (x, k, 1), REF (x, k, 2)) - eps;
91  maxx(k) = max (REF (x, k, 0), REF (x, k, 1), REF (x, k, 2)) + eps;
92  miny(k) = min (REF (y, k, 0), REF (y, k, 1), REF (y, k, 2)) - eps;
93  maxy(k) = max (REF (y, k, 0), REF (y, k, 1), REF (y, k, 2)) + eps;
94  }
95 
96  const octave_idx_type np = xi.numel ();
97  ColumnVector values (np);
98 
99  double x0 = 0.0, y0 = 0.0;
100  double a11 = 0.0, a12 = 0.0, a21 = 0.0, a22 = 0.0, det = 0.0;
101  double xt = 0.0, yt = 0.0;
102  double dx1 = 0.0, dx2 = 0.0, c1 = 0.0, c2 = 0.0;
103 
104  octave_idx_type k = nelem; // k is more than just an index variable.
105 
106  for (octave_idx_type kp = 0; kp < np; kp++) // for each point
107  {
108  xt = xi (kp);
109  yt = yi (kp);
110 
111  // Check if point (xt,yt) is in the triangle that was last examined.
112  // This is for inputs where points are in contiguous order,
113  // like when the points are sampled from a continuous path.
114  if (k < nelem) // This check will be false for the very first point.
115  {
116  // If we are here, then x0, y0, det all exist from before.
117  dx1 = xt - x0;
118  dx2 = yt - y0;
119  c1 = (a22 * dx1 - a21 * dx2) / det;
120  c2 = (-a12 * dx1 + a11 * dx2) / det;
121  if (c1 >= -eps && c2 >= -eps && (c1 + c2) <= 1 + eps)
122  {
123  values (kp) = k+1;
124  continue;
125  }
126  }
127 
128  // The point is not in the same triangle, so go through all triangles.
129  for (k = 0; k < nelem; k++)
130  {
131  if (xt >= minx(k) && xt <= maxx(k) && yt >= miny(k) && yt <= maxy(k))
132  {
133  // Point is inside the triangle's bounding rectangle:
134  // See if it's inside the triangle itself.
135  x0 = REF (x, k, 0);
136  y0 = REF (y, k, 0);
137  a11 = REF (x, k, 1) - x0;
138  a12 = REF (y, k, 1) - y0;
139  a21 = REF (x, k, 2) - x0;
140  a22 = REF (y, k, 2) - y0;
141  det = a11 * a22 - a21 * a12;
142 
143  // solve the system
144  dx1 = xt - x0;
145  dx2 = yt - y0;
146  c1 = (a22 * dx1 - a21 * dx2) / det;
147  c2 = (-a12 * dx1 + a11 * dx2) / det;
148  if (c1 >= -eps && c2 >= -eps && (c1 + c2) <= 1 + eps)
149  {
150  values (kp) = k+1;
151  break;
152  }
153  } //end see if it's inside the triangle itself
154  } //end for each triangle
155 
156  if (k == nelem)
157  values (kp) = lo_ieee_nan_value ();
158 
159  } //end for each point
160 
161  return ovl (values);
162 }
163 
164 /*
165 %!shared x, y, tri
166 %! x = [-1;-1;1];
167 %! y = [-1;1;-1];
168 %! tri = [1, 2, 3];
169 %!assert (tsearch (x,y,tri,-1,-1), 1)
170 %!assert (tsearch (x,y,tri, 1,-1), 1)
171 %!assert (tsearch (x,y,tri,-1, 1), 1)
172 %!assert (tsearch (x,y,tri,-1/3, -1/3), 1)
173 %!assert (tsearch (x,y,tri, 1, 1), NaN)
174 
175 %!error tsearch ()
176 */
177 
OCTAVE_END_NAMESPACE(octave)
static int elem
Definition: __contourc__.cc:54
Definition: dMatrix.h:42
OCTAVE_BEGIN_NAMESPACE(octave) static octave_value daspk_fcn
T eps(const T &x)
Definition: data.cc:4942
OCTINTERP_API void print_usage(void)
Definition: defun-int.h:72
#define DEFUN(name, args_name, nargout_name, doc)
Macro to define a builtin function.
Definition: defun.h:56
double lo_ieee_nan_value(void)
Definition: lo-ieee.cc:84
F77_RET_T const F77_DBLE * x
octave_value_list ovl(const OV_Args &... args)
Construct an octave_value_list with less typing.
Definition: ovl.h:211
static const double xi[33]
Definition: quadcc.cc:68
double min(double a, double b, double c)
Definition: tsearch.cc:45
double max(double a, double b, double c)
Definition: tsearch.cc:40
#define REF(x, k, i)
Definition: tsearch.cc:50