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