Mercurial > octave
annotate src/DLD-FUNCTIONS/tsearch.cc @ 10154:40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Wed, 20 Jan 2010 17:33:41 -0500 |
parents | 2c279308f6ab |
children | c48b7048e720 |
rev | line source |
---|---|
6823 | 1 /* |
2 | |
9245 | 3 Copyright (C) 2002, 2007, 2009 Andreas Stahel |
6823 | 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 | |
7016 | 9 Free Software Foundation; either version 3 of the License, or (at your |
10 option) any later version. | |
6823 | 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 | |
7016 | 18 along with Octave; see the file COPYING. If not, see |
19 <http://www.gnu.org/licenses/>. | |
6823 | 20 |
21 */ | |
22 | |
7016 | 23 // Author: Andreas Stahel <Andreas.Stahel@hta-bi.bfh.ch> |
24 | |
9786
2c279308f6ab
fix includes in some src/DLD-FUNCTIONS files
John W. Eaton <jwe@octave.org>
parents:
9245
diff
changeset
|
25 #ifdef HAVE_CONFIG_H |
2c279308f6ab
fix includes in some src/DLD-FUNCTIONS files
John W. Eaton <jwe@octave.org>
parents:
9245
diff
changeset
|
26 #include <config.h> |
2c279308f6ab
fix includes in some src/DLD-FUNCTIONS files
John W. Eaton <jwe@octave.org>
parents:
9245
diff
changeset
|
27 #endif |
2c279308f6ab
fix includes in some src/DLD-FUNCTIONS files
John W. Eaton <jwe@octave.org>
parents:
9245
diff
changeset
|
28 |
6823 | 29 #include <iostream> |
30 #include <fstream> | |
31 #include <string> | |
32 | |
33 #include "lo-ieee.h" | |
7231 | 34 #include "lo-math.h" |
6823 | 35 |
9786
2c279308f6ab
fix includes in some src/DLD-FUNCTIONS files
John W. Eaton <jwe@octave.org>
parents:
9245
diff
changeset
|
36 #include "defun-dld.h" |
2c279308f6ab
fix includes in some src/DLD-FUNCTIONS files
John W. Eaton <jwe@octave.org>
parents:
9245
diff
changeset
|
37 #include "error.h" |
2c279308f6ab
fix includes in some src/DLD-FUNCTIONS files
John W. Eaton <jwe@octave.org>
parents:
9245
diff
changeset
|
38 #include "oct-obj.h" |
2c279308f6ab
fix includes in some src/DLD-FUNCTIONS files
John W. Eaton <jwe@octave.org>
parents:
9245
diff
changeset
|
39 #include "parse.h" |
2c279308f6ab
fix includes in some src/DLD-FUNCTIONS files
John W. Eaton <jwe@octave.org>
parents:
9245
diff
changeset
|
40 |
6823 | 41 inline double max(double a, double b, double c) |
42 { | |
43 if (a < b) | |
44 return ( b < c ? c : b ); | |
45 else | |
46 return ( a < c ? c : a ); | |
47 } | |
48 | |
49 inline double min(double a, double b, double c) | |
50 { | |
51 if (a > b) | |
52 return ( b > c ? c : b ); | |
53 else | |
54 return ( a > c ? c : a ); | |
55 } | |
56 | |
57 #define REF(x,k,i) x(static_cast<octave_idx_type>(elem((k), (i))) - 1) | |
58 | |
59 // for large data set the algorithm is very slow | |
60 // one should presort (how?) either the elements of the points of evaluation | |
61 // to cut down the time needed to decide which triangle contains the | |
62 // given point | |
63 | |
64 // e.g., build up a neighbouring triangle structure and use a simplex-like | |
65 // method to traverse it | |
66 | |
67 DEFUN_DLD (tsearch, args, , | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
68 "-*- texinfo -*-\n\ |
6823 | 69 @deftypefn {Loadable Function} {@var{idx} =} tsearch (@var{x}, @var{y}, @var{t}, @var{xi}, @var{yi})\n\ |
9064
7c02ec148a3c
Check grammar on all .cc files
Rik <rdrider0-list@yahoo.com>
parents:
7231
diff
changeset
|
70 Searches for the enclosing Delaunay convex hull. For @code{@var{t} =\n\ |
6823 | 71 delaunay (@var{x}, @var{y})}, finds the index in @var{t} containing the\n\ |
9064
7c02ec148a3c
Check grammar on all .cc files
Rik <rdrider0-list@yahoo.com>
parents:
7231
diff
changeset
|
72 points @code{(@var{xi}, @var{yi})}. For points outside the convex hull,\n\ |
6823 | 73 @var{idx} is NaN.\n\ |
74 @seealso{delaunay, delaunayn}\n\ | |
75 @end deftypefn") | |
76 { | |
77 const double eps=1.0e-12; | |
78 | |
79 octave_value_list retval; | |
80 const int nargin = args.length (); | |
81 if ( nargin != 5 ) { | |
82 print_usage (); | |
83 return retval; | |
84 } | |
85 | |
86 const ColumnVector x(args(0).vector_value()); | |
87 const ColumnVector y(args(1).vector_value()); | |
88 const Matrix elem(args(2).matrix_value()); | |
89 const ColumnVector xi(args(3).vector_value()); | |
90 const ColumnVector yi(args(4).vector_value()); | |
91 | |
92 if (error_state) return retval; | |
93 | |
94 const octave_idx_type nelem = elem.rows(); | |
95 | |
96 ColumnVector minx(nelem); | |
97 ColumnVector maxx(nelem); | |
98 ColumnVector miny(nelem); | |
99 ColumnVector maxy(nelem); | |
100 for(octave_idx_type k = 0; k < nelem; k++) | |
101 { | |
102 minx(k) = min(REF(x,k,0), REF(x,k,1), REF(x,k,2)) - eps; | |
103 maxx(k) = max(REF(x,k,0), REF(x,k,1), REF(x,k,2)) + eps; | |
104 miny(k) = min(REF(y,k,0), REF(y,k,1), REF(y,k,2)) - eps; | |
105 maxy(k) = max(REF(y,k,0), REF(y,k,1), REF(y,k,2)) + eps; | |
106 } | |
107 | |
108 const octave_idx_type np = xi.length(); | |
109 ColumnVector values(np); | |
110 | |
111 double x0 = 0.0, y0 = 0.0; | |
112 double a11 = 0.0, a12 = 0.0, a21 = 0.0, a22 = 0.0, 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) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
122 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
123 const double dx1 = xt - x0; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
124 const double dx2 = yt - y0; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
125 const double c1 = (a22 * dx1 - a21 * dx2) / det; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
126 const double c2 = (-a12 * dx1 + a11 * dx2) / det; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
127 if ( c1 >= -eps && c2 >= -eps && (c1 + c2) <= (1 + eps)) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
128 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
129 values(kp) = double(k+1); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
130 continue; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
131 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
132 } |
6823 | 133 |
134 // it doesn't, so go through all elements | |
135 for (k = 0; k < nelem; k++) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
136 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
137 OCTAVE_QUIT; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
138 if (xt >= minx(k) && xt <= maxx(k) && |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
139 yt >= miny(k) && yt <= maxy(k) ) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
140 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
141 // element inside the minimum rectangle: examine it closely |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
142 x0 = REF(x,k,0); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
143 y0 = REF(y,k,0); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
144 a11 = REF(x,k,1)-x0; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
145 a12 = REF(y,k,1)-y0; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
146 a21 = REF(x,k,2)-x0; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
147 a22 = REF(y,k,2)-y0; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
148 det = a11 * a22 - a21 * a12; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
149 |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
150 // solve the system |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
151 const double dx1 = xt - x0; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
152 const double dx2 = yt - y0; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
153 const double c1 = (a22 * dx1 - a21 * dx2) / det; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
154 const double c2 = (-a12 * dx1 + a11 * dx2) / det; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
155 if ((c1 >= -eps) && (c2 >= -eps) && ((c1 + c2) <= (1 + eps))) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
156 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
157 values(kp) = double(k+1); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
158 break; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
159 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
160 } //endif # examine this element closely |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
161 } //endfor # each element |
6823 | 162 |
163 if (k == nelem) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9786
diff
changeset
|
164 values(kp) = lo_ieee_nan_value (); |
6823 | 165 |
166 } //endfor # kp | |
167 | |
168 retval(0) = values; | |
169 | |
170 return retval; | |
171 } | |
172 | |
173 /* | |
174 %!shared x, y, tri | |
175 %! x = [-1;-1;1]; | |
176 %! y = [-1;1;-1]; | |
177 %! tri = [1, 2, 3]; | |
178 %!error (tsearch()) | |
179 %!assert (tsearch (x,y,tri,-1,-1), 1) | |
180 %!assert (tsearch (x,y,tri, 1,-1), 1) | |
181 %!assert (tsearch (x,y,tri,-1, 1), 1) | |
182 %!assert (tsearch (x,y,tri,-1/3, -1/3), 1) | |
183 %!assert (tsearch (x,y,tri, 1, 1), NaN) | |
184 | |
185 */ |