GNU Octave  4.0.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
tsearch.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2002-2015 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 the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11 
12 Octave is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 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 <http://www.gnu.org/licenses/>.
20 
21 */
22 
23 // Author: Andreas Stahel <Andreas.Stahel@hta-bi.bfh.ch>
24 
25 #ifdef HAVE_CONFIG_H
26 #include <config.h>
27 #endif
28 
29 #include "lo-ieee.h"
30 #include "lo-math.h"
31 
32 #include "defun.h"
33 #include "error.h"
34 #include "oct-obj.h"
35 
36 inline double max (double a, double b, double c)
37 {
38  if (a < b)
39  return (b < c ? c : b);
40  else
41  return (a < c ? c : a);
42 }
43 
44 inline double min (double a, double b, double c)
45 {
46  if (a > b)
47  return (b > c ? c : b);
48  else
49  return (a > c ? c : a);
50 }
51 
52 #define REF(x,k,i) x(static_cast<octave_idx_type>(elem((k), (i))) - 1)
53 
54 // for large data set the algorithm is very slow
55 // one should presort (how?) either the elements of the points of evaluation
56 // to cut down the time needed to decide which triangle contains the
57 // given point
58 
59 // e.g., build up a neighbouring triangle structure and use a simplex-like
60 // method to traverse it
61 
62 DEFUN (tsearch, args, ,
63  "-*- texinfo -*-\n\
64 @deftypefn {Built-in Function} {@var{idx} =} tsearch (@var{x}, @var{y}, @var{t}, @var{xi}, @var{yi})\n\
65 Search for the enclosing Delaunay convex hull.\n\
66 \n\
67 For @code{@var{t} = delaunay (@var{x}, @var{y})}, finds the index in @var{t}\n\
68 containing the points @code{(@var{xi}, @var{yi})}. For points outside the\n\
69 convex hull, @var{idx} is NaN.\n\
70 @seealso{delaunay, delaunayn}\n\
71 @end deftypefn")
72 {
73  const double eps=1.0e-12;
74 
75  octave_value_list retval;
76  const int nargin = args.length ();
77  if (nargin != 5)
78  {
79  print_usage ();
80  return retval;
81  }
82 
83  const ColumnVector x (args(0).vector_value ());
84  const ColumnVector y (args(1).vector_value ());
85  const Matrix elem (args(2).matrix_value ());
86  const ColumnVector xi (args(3).vector_value ());
87  const ColumnVector yi (args(4).vector_value ());
88 
89  if (error_state)
90  return retval;
91 
92  const octave_idx_type nelem = elem.rows ();
93 
94  ColumnVector minx (nelem);
95  ColumnVector maxx (nelem);
96  ColumnVector miny (nelem);
97  ColumnVector maxy (nelem);
98  for (octave_idx_type k = 0; k < nelem; k++)
99  {
100  minx(k) = min (REF (x, k, 0), REF (x, k, 1), REF (x, k, 2)) - eps;
101  maxx(k) = max (REF (x, k, 0), REF (x, k, 1), REF (x, k, 2)) + eps;
102  miny(k) = min (REF (y, k, 0), REF (y, k, 1), REF (y, k, 2)) - eps;
103  maxy(k) = max (REF (y, k, 0), REF (y, k, 1), REF (y, k, 2)) + eps;
104  }
105 
106  const octave_idx_type np = xi.length ();
107  ColumnVector values (np);
108 
109  double x0, y0, a11, a12, a21, a22, det;
110  x0 = y0 = 0.0;
111  a11 = a12 = a21 = a22 = 0.0;
112  det = 0.0;
113 
114  octave_idx_type k = nelem; // k is a counter of elements
115  for (octave_idx_type kp = 0; kp < np; kp++)
116  {
117  const double xt = xi(kp);
118  const double yt = yi(kp);
119 
120  // check if last triangle contains the next point
121  if (k < nelem)
122  {
123  const double dx1 = xt - x0;
124  const double dx2 = yt - y0;
125  const double c1 = (a22 * dx1 - a21 * dx2) / det;
126  const double c2 = (-a12 * dx1 + a11 * dx2) / det;
127  if (c1 >= -eps && c2 >= -eps && (c1 + c2) <= (1 + eps))
128  {
129  values(kp) = double(k+1);
130  continue;
131  }
132  }
133 
134  // it doesn't, so go through all elements
135  for (k = 0; k < nelem; k++)
136  {
137  OCTAVE_QUIT;
138  if (xt >= minx(k) && xt <= maxx(k) && yt >= miny(k) && yt <= maxy(k))
139  {
140  // element inside the minimum rectangle: examine it closely
141  x0 = REF (x, k, 0);
142  y0 = REF (y, k, 0);
143  a11 = REF (x, k, 1) - x0;
144  a12 = REF (y, k, 1) - y0;
145  a21 = REF (x, k, 2) - x0;
146  a22 = REF (y, k, 2) - y0;
147  det = a11 * a22 - a21 * a12;
148 
149  // solve the system
150  const double dx1 = xt - x0;
151  const double dx2 = yt - y0;
152  const double c1 = (a22 * dx1 - a21 * dx2) / det;
153  const double c2 = (-a12 * dx1 + a11 * dx2) / det;
154  if ((c1 >= -eps) && (c2 >= -eps) && ((c1 + c2) <= (1 + eps)))
155  {
156  values(kp) = double(k+1);
157  break;
158  }
159  } //endif # examine this element closely
160  } //endfor # each element
161 
162  if (k == nelem)
163  values(kp) = lo_ieee_nan_value ();
164 
165  } //endfor # kp
166 
167  retval(0) = values;
168 
169  return retval;
170 }
171 
172 /*
173 %!shared x, y, tri
174 %! x = [-1;-1;1];
175 %! y = [-1;1;-1];
176 %! tri = [1, 2, 3];
177 %!assert (tsearch (x,y,tri,-1,-1), 1)
178 %!assert (tsearch (x,y,tri, 1,-1), 1)
179 %!assert (tsearch (x,y,tri,-1, 1), 1)
180 %!assert (tsearch (x,y,tri,-1/3, -1/3), 1)
181 %!assert (tsearch (x,y,tri, 1, 1), NaN)
182 
183 %!error tsearch ()
184 */
#define REF(x, k, i)
Definition: tsearch.cc:52
OCTINTERP_API void print_usage(void)
Definition: defun.cc:51
octave_idx_type length(void) const
Definition: oct-obj.h:89
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:44
double lo_ieee_nan_value(void)
Definition: lo-ieee.cc:126
double max(double a, double b, double c)
Definition: tsearch.cc:36
octave_idx_type rows(void) const
Definition: Array.h:313
double min(double a, double b, double c)
Definition: tsearch.cc:44
int error_state
Definition: error.cc:101
static int elem
Definition: __contourc__.cc:49
Definition: dMatrix.h:35
#define eps(C)
octave_idx_type length(void) const
Number of elements in the array.
Definition: Array.h:267
static MArray< double > const octave_idx_type const octave_idx_type octave_idx_type octave_idx_type octave_idx_type c1
#define OCTAVE_QUIT
Definition: quit.h:130
static const double xi[33]
Definition: quadcc.cc:56
F77_RET_T const double * x