GNU Octave  4.4.1
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-2018 Andreas Stahel
4 
5 This file is part of Octave.
6 
7 Octave is free software: you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by
9 the Free Software Foundation, either version 3 of the License, or
10 (at your option) any later version.
11 
12 Octave is distributed in the hope that it will be useful, but
13 WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with Octave; see the file COPYING. If not, see
19 <https://www.gnu.org/licenses/>.
20 
21 */
22 
23 // Author: Andreas Stahel <Andreas.Stahel@hta-bi.bfh.ch>
24 
25 #if defined (HAVE_CONFIG_H)
26 # include "config.h"
27 #endif
28 
29 #include <cmath>
30 
31 #include "lo-ieee.h"
32 
33 #include "defun.h"
34 #include "error.h"
35 #include "ovl.h"
36 
37 inline double max (double a, double b, double c)
38 {
39  if (a < b)
40  return (b < c ? c : b);
41  else
42  return (a < c ? c : a);
43 }
44 
45 inline double min (double a, double b, double c)
46 {
47  if (a > b)
48  return (b > c ? c : b);
49  else
50  return (a > c ? c : a);
51 }
52 
53 #define REF(x,k,i) x(static_cast<octave_idx_type>(elem((k), (i))) - 1)
54 
55 // for large data set the algorithm is very slow
56 // one should presort (how?) either the elements of the points of evaluation
57 // to cut down the time needed to decide which triangle contains the
58 // given point
59 
60 // e.g., build up a neighbouring triangle structure and use a simplex-like
61 // method to traverse it
62 
63 DEFUN (tsearch, args, ,
64  doc: /* -*- texinfo -*-
65 @deftypefn {} {@var{idx} =} tsearch (@var{x}, @var{y}, @var{t}, @var{xi}, @var{yi})
66 Search for the enclosing Delaunay convex hull.
67 
68 For @code{@var{t} = delaunay (@var{x}, @var{y})}, finds the index in @var{t}
69 containing the points @code{(@var{xi}, @var{yi})}. For points outside the
70 convex hull, @var{idx} is NaN.
71 @seealso{delaunay, delaunayn}
72 @end deftypefn */)
73 {
74  if (args.length () != 5)
75  print_usage ();
76 
77  const double eps = 1.0e-12;
78 
79  const ColumnVector x (args(0).vector_value ());
80  const ColumnVector y (args(1).vector_value ());
81  const Matrix elem (args(2).matrix_value ());
82  const ColumnVector xi (args(3).vector_value ());
83  const ColumnVector yi (args(4).vector_value ());
84 
85  const octave_idx_type nelem = elem.rows ();
86 
87  ColumnVector minx (nelem);
88  ColumnVector maxx (nelem);
89  ColumnVector miny (nelem);
90  ColumnVector maxy (nelem);
91  for (octave_idx_type k = 0; k < nelem; k++)
92  {
93  minx(k) = min (REF (x, k, 0), REF (x, k, 1), REF (x, k, 2)) - eps;
94  maxx(k) = max (REF (x, k, 0), REF (x, k, 1), REF (x, k, 2)) + eps;
95  miny(k) = min (REF (y, k, 0), REF (y, k, 1), REF (y, k, 2)) - eps;
96  maxy(k) = max (REF (y, k, 0), REF (y, k, 1), REF (y, k, 2)) + eps;
97  }
98 
99  const octave_idx_type np = xi.numel ();
100  ColumnVector values (np);
101 
102  double x0, y0, a11, a12, a21, a22, det;
103  x0 = y0 = 0.0;
104  a11 = a12 = a21 = a22 = 0.0;
105  det = 0.0;
106 
107  octave_idx_type k = nelem; // k is a counter of elements
108  for (octave_idx_type kp = 0; kp < np; kp++)
109  {
110  const double xt = xi(kp);
111  const double yt = yi(kp);
112 
113  // check if last triangle contains the next point
114  if (k < nelem)
115  {
116  const double dx1 = xt - x0;
117  const double dx2 = yt - y0;
118  const double c1 = (a22 * dx1 - a21 * dx2) / det;
119  const double c2 = (-a12 * dx1 + a11 * dx2) / det;
120  if (c1 >= -eps && c2 >= -eps && (c1 + c2) <= (1 + eps))
121  {
122  values(kp) = double(k+1);
123  continue;
124  }
125  }
126 
127  // it doesn't, so go through all elements
128  for (k = 0; k < nelem; k++)
129  {
130  octave_quit ();
131 
132  if (xt >= minx(k) && xt <= maxx(k) && yt >= miny(k) && yt <= maxy(k))
133  {
134  // element inside the minimum rectangle: examine it closely
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  const double dx1 = xt - x0;
145  const double dx2 = yt - y0;
146  const double c1 = (a22 * dx1 - a21 * dx2) / det;
147  const double c2 = (-a12 * dx1 + a11 * dx2) / det;
148  if ((c1 >= -eps) && (c2 >= -eps) && ((c1 + c2) <= (1 + eps)))
149  {
150  values(kp) = double(k+1);
151  break;
152  }
153  } //endif # examine this element closely
154  } //endfor # each element
155 
156  if (k == nelem)
157  values(kp) = lo_ieee_nan_value ();
158 
159  } //endfor # kp
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 */
#define REF(x, k, i)
Definition: tsearch.cc:53
OCTINTERP_API void print_usage(void)
Definition: defun.cc:54
for large enough k
Definition: lu.cc:617
#define DEFUN(name, args_name, nargout_name, doc)
Macro to define a builtin function.
Definition: defun.h:53
double lo_ieee_nan_value(void)
Definition: lo-ieee.cc:91
double max(double a, double b, double c)
Definition: tsearch.cc:37
nd example oindent opens the file binary numeric values will be read assuming they are stored in IEEE format with the least significant bit and then converted to the native representation Opening a file that is already open simply opens it again and returns a separate file id It is not an error to open a file several though writing to the same file through several different file ids may produce unexpected results The possible values of text mode reading and writing automatically converts linefeeds to the appropriate line end character for the you may append a you must also open the file in binary mode The parameter conversions are currently only supported for and permissions will be set to and then everything is written in a single operation This is very efficient and improves performance c
Definition: file-io.cc:587
cell array If invoked with two or more scalar integer or a vector of integer values
Definition: ov-cell.cc:1241
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
Definition: cellfun.cc:400
double min(double a, double b, double c)
Definition: tsearch.cc:45
static int elem
Definition: __contourc__.cc:47
Definition: dMatrix.h:36
#define eps(C)
OCTAVE_EXPORT octave_value_list isa nd deftypefn *return ovl(args(0).isinteger())
the element is set to zero In other the statement xample y
Definition: data.cc:5264
b
Definition: cellfun.cc:400
static const double xi[33]
Definition: quadcc.cc:62
F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T const F77_DBLE * x