1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 2002-2021 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
max(double a,double b,double c)38 inline double max (double a, double b, double c)
39 {
40 if (a < b)
41 return (b < c ? c : b);
42 else
43 return (a < c ? c : a);
44 }
45
min(double a,double b,double c)46 inline double min (double a, double b, double c)
47 {
48 if (a > b)
49 return (b > c ? c : b);
50 else
51 return (a > c ? c : a);
52 }
53
54 #define REF(x,k,i) x(static_cast<octave_idx_type> (elem((k), (i))) - 1)
55
56 // for large data set the algorithm is very slow
57 // one should presort (how?) either the elements of the points of evaluation
58 // to cut down the time needed to decide which triangle contains the
59 // given point
60
61 // e.g., build up a neighbouring triangle structure and use a simplex-like
62 // method to traverse it
63
64 DEFUN (tsearch, args, ,
65 doc: /* -*- texinfo -*-
66 @deftypefn {} {@var{idx} =} tsearch (@var{x}, @var{y}, @var{t}, @var{xi}, @var{yi})
67 Search for the enclosing Delaunay convex hull.
68
69 For @code{@var{t} = delaunay (@var{x}, @var{y})}, finds the index in @var{t}
70 containing the points @code{(@var{xi}, @var{yi})}. For points outside the
71 convex hull, @var{idx} is NaN.
72 @seealso{delaunay, delaunayn}
73 @end deftypefn */)
74 {
75 if (args.length () != 5)
76 print_usage ();
77
78 const double eps = 1.0e-12;
79
80 const ColumnVector x (args(0).vector_value ());
81 const ColumnVector y (args(1).vector_value ());
82 const Matrix elem (args(2).matrix_value ());
83 const ColumnVector xi (args(3).vector_value ());
84 const ColumnVector yi (args(4).vector_value ());
85
86 const octave_idx_type nelem = elem.rows ();
87
88 ColumnVector minx (nelem);
89 ColumnVector maxx (nelem);
90 ColumnVector miny (nelem);
91 ColumnVector maxy (nelem);
92 for (octave_idx_type k = 0; k < nelem; k++)
93 {
94 minx(k) = min (REF (x, k, 0), REF (x, k, 1), REF (x, k, 2)) - eps;
95 maxx(k) = max (REF (x, k, 0), REF (x, k, 1), REF (x, k, 2)) + eps;
96 miny(k) = min (REF (y, k, 0), REF (y, k, 1), REF (y, k, 2)) - eps;
97 maxy(k) = max (REF (y, k, 0), REF (y, k, 1), REF (y, k, 2)) + eps;
98 }
99
100 const octave_idx_type np = xi.numel ();
101 ColumnVector values (np);
102
103 double x0, y0, a11, a12, a21, a22, det;
104 x0 = y0 = 0.0;
105 a11 = a12 = a21 = a22 = 0.0;
106 det = 0.0;
107
108 octave_idx_type k = nelem; // k is a counter of elements
109 for (octave_idx_type kp = 0; kp < np; kp++)
110 {
111 const double xt = xi(kp);
112 const double yt = yi(kp);
113
114 // check if last triangle contains the next point
115 if (k < nelem)
116 {
117 const double dx1 = xt - x0;
118 const double dx2 = yt - y0;
119 const double c1 = (a22 * dx1 - a21 * dx2) / det;
120 const double c2 = (-a12 * dx1 + a11 * dx2) / det;
121 if (c1 >= -eps && c2 >= -eps && (c1 + c2) <= (1 + eps))
122 {
123 values(kp) = double(k+1);
124 continue;
125 }
126 }
127
128 // it doesn't, so go through all elements
129 for (k = 0; k < nelem; k++)
130 {
131 octave_quit ();
132
133 if (xt >= minx(k) && xt <= maxx(k) && yt >= miny(k) && yt <= maxy(k))
134 {
135 // element inside the minimum rectangle: examine it closely
136 x0 = REF (x, k, 0);
137 y0 = REF (y, k, 0);
138 a11 = REF (x, k, 1) - x0;
139 a12 = REF (y, k, 1) - y0;
140 a21 = REF (x, k, 2) - x0;
141 a22 = REF (y, k, 2) - y0;
142 det = a11 * a22 - a21 * a12;
143
144 // solve the system
145 const double dx1 = xt - x0;
146 const double dx2 = yt - y0;
147 const double c1 = (a22 * dx1 - a21 * dx2) / det;
148 const double c2 = (-a12 * dx1 + a11 * dx2) / det;
149 if ((c1 >= -eps) && (c2 >= -eps) && ((c1 + c2) <= (1 + eps)))
150 {
151 values(kp) = double(k+1);
152 break;
153 }
154 } //endif # examine this element closely
155 } //endfor # each element
156
157 if (k == nelem)
158 values(kp) = lo_ieee_nan_value ();
159
160 } //endfor # kp
161
162 return ovl (values);
163 }
164
165 /*
166 %!shared x, y, tri
167 %! x = [-1;-1;1];
168 %! y = [-1;1;-1];
169 %! tri = [1, 2, 3];
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, 1), 1)
173 %!assert (tsearch (x,y,tri,-1/3, -1/3), 1)
174 %!assert (tsearch (x,y,tri, 1, 1), NaN)
175
176 %!error tsearch ()
177 */
178