Mercurial > octave-nkf

comparison libinterp/corefcn/tsearch.cc @ 19912:12ecb7212b44

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
move some files without external dependencies from dldfcn to corefcn * __dsearchn__.cc, __ichol__.cc, __ilu__.cc, tsearch.cc: Move from dldfcn to corefcn directory. Use DEFUN instead of DEFUN_DLD. * libinterp/corefcn/module.mk, libinterp/dldfcn/module-files: Update.
author John W. Eaton <jwe@octave.org>
date Fri, 27 Feb 2015 19:44:28 -0500
parents libinterp/dldfcn/tsearch.cc@4197fc428c7d
children 4f45eaf83908
comparison
equal deleted inserted replaced
19911:f799bf70350f 19912:12ecb7212b44
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. For @code{@var{t} =\n\
66 delaunay (@var{x}, @var{y})}, finds the index in @var{t} containing the\n\
67 points @code{(@var{xi}, @var{yi})}. For points outside the convex hull,\n\
68 @var{idx} is NaN.\n\
69 @seealso{delaunay, delaunayn}\n\
70 @end deftypefn")
71 {
72 const double eps=1.0e-12;
73
74 octave_value_list retval;
75 const int nargin = args.length ();
76 if (nargin != 5)
77 {
78 print_usage ();
79 return retval;
80 }
81
82 const ColumnVector x (args(0).vector_value ());
83 const ColumnVector y (args(1).vector_value ());
84 const Matrix elem (args(2).matrix_value ());
85 const ColumnVector xi (args(3).vector_value ());
86 const ColumnVector yi (args(4).vector_value ());
87
88 if (error_state)
89 return retval;
90
91 const octave_idx_type nelem = elem.rows ();
92
93 ColumnVector minx (nelem);
94 ColumnVector maxx (nelem);
95 ColumnVector miny (nelem);
96 ColumnVector maxy (nelem);
97 for (octave_idx_type k = 0; k < nelem; k++)
98 {
99 minx(k) = min (REF (x, k, 0), REF (x, k, 1), REF (x, k, 2)) - eps;
100 maxx(k) = max (REF (x, k, 0), REF (x, k, 1), REF (x, k, 2)) + eps;
101 miny(k) = min (REF (y, k, 0), REF (y, k, 1), REF (y, k, 2)) - eps;
102 maxy(k) = max (REF (y, k, 0), REF (y, k, 1), REF (y, k, 2)) + eps;
103 }
104
105 const octave_idx_type np = xi.length ();
106 ColumnVector values (np);
107
108 double x0, y0, a11, a12, a21, a22, det;
109 x0 = y0 = 0.0;
110 a11 = a12 = a21 = a22 = 0.0;
111 det = 0.0;
112
113 octave_idx_type k = nelem; // k is a counter of elements
114 for (octave_idx_type kp = 0; kp < np; kp++)
115 {
116 const double xt = xi(kp);
117 const double yt = yi(kp);
118
119 // check if last triangle contains the next point
120 if (k < nelem)
121 {
122 const double dx1 = xt - x0;
123 const double dx2 = yt - y0;
124 const double c1 = (a22 * dx1 - a21 * dx2) / det;
125 const double c2 = (-a12 * dx1 + a11 * dx2) / det;
126 if (c1 >= -eps && c2 >= -eps && (c1 + c2) <= (1 + eps))
127 {
128 values(kp) = double(k+1);
129 continue;
130 }
131 }
132
133 // it doesn't, so go through all elements
134 for (k = 0; k < nelem; k++)
135 {
136 OCTAVE_QUIT;
137 if (xt >= minx(k) && xt <= maxx(k) && yt >= miny(k) && yt <= maxy(k))
138 {
139 // element inside the minimum rectangle: examine it closely
140 x0 = REF (x, k, 0);
141 y0 = REF (y, k, 0);
142 a11 = REF (x, k, 1) - x0;
143 a12 = REF (y, k, 1) - y0;
144 a21 = REF (x, k, 2) - x0;
145 a22 = REF (y, k, 2) - y0;
146 det = a11 * a22 - a21 * a12;
147
148 // solve the system
149 const double dx1 = xt - x0;
150 const double dx2 = yt - y0;
151 const double c1 = (a22 * dx1 - a21 * dx2) / det;
152 const double c2 = (-a12 * dx1 + a11 * dx2) / det;
153 if ((c1 >= -eps) && (c2 >= -eps) && ((c1 + c2) <= (1 + eps)))
154 {
155 values(kp) = double(k+1);
156 break;
157 }
158 } //endif # examine this element closely
159 } //endfor # each element
160
161 if (k == nelem)
162 values(kp) = lo_ieee_nan_value ();
163
164 } //endfor # kp
165
166 retval(0) = values;
167
168 return retval;
169 }
170
171 /*
172 %!shared x, y, tri
173 %! x = [-1;-1;1];
174 %! y = [-1;1;-1];
175 %! tri = [1, 2, 3];
176 %!assert (tsearch (x,y,tri,-1,-1), 1)
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/3, -1/3), 1)
180 %!assert (tsearch (x,y,tri, 1, 1), NaN)
181
182 %!error tsearch ()
183 */