Mercurial > octave
annotate src/DLD-FUNCTIONS/tsearch.cc @ 9245:16f53d29049f
update copyright notices
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Fri, 22 May 2009 10:46:00 -0400 |
| parents | 7c02ec148a3c |
| children | 2c279308f6ab |
| 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 | |
| 6823 | 25 #include <iostream> |
| 26 #include <fstream> | |
| 27 #include <string> | |
| 28 | |
| 29 #include "oct.h" | |
| 30 #include "parse.h" | |
| 31 #include "lo-ieee.h" | |
| 7231 | 32 #include "lo-math.h" |
| 6823 | 33 |
| 34 inline double max(double a, double b, double c) | |
| 35 { | |
| 36 if (a < b) | |
| 37 return ( b < c ? c : b ); | |
| 38 else | |
| 39 return ( a < c ? c : a ); | |
| 40 } | |
| 41 | |
| 42 inline double min(double a, double b, double c) | |
| 43 { | |
| 44 if (a > b) | |
| 45 return ( b > c ? c : b ); | |
| 46 else | |
| 47 return ( a > c ? c : a ); | |
| 48 } | |
| 49 | |
| 50 #define REF(x,k,i) x(static_cast<octave_idx_type>(elem((k), (i))) - 1) | |
| 51 | |
| 52 // for large data set the algorithm is very slow | |
| 53 // one should presort (how?) either the elements of the points of evaluation | |
| 54 // to cut down the time needed to decide which triangle contains the | |
| 55 // given point | |
| 56 | |
| 57 // e.g., build up a neighbouring triangle structure and use a simplex-like | |
| 58 // method to traverse it | |
| 59 | |
| 60 DEFUN_DLD (tsearch, args, , | |
| 61 "-*- texinfo -*-\n\ | |
| 62 @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
|
63 Searches for the enclosing Delaunay convex hull. For @code{@var{t} =\n\ |
| 6823 | 64 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
|
65 points @code{(@var{xi}, @var{yi})}. For points outside the convex hull,\n\ |
| 6823 | 66 @var{idx} is NaN.\n\ |
| 67 @seealso{delaunay, delaunayn}\n\ | |
| 68 @end deftypefn") | |
| 69 { | |
| 70 const double eps=1.0e-12; | |
| 71 | |
| 72 octave_value_list retval; | |
| 73 const int nargin = args.length (); | |
| 74 if ( nargin != 5 ) { | |
| 75 print_usage (); | |
| 76 return retval; | |
| 77 } | |
| 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 if (error_state) return retval; | |
| 86 | |
| 87 const octave_idx_type nelem = elem.rows(); | |
| 88 | |
| 89 ColumnVector minx(nelem); | |
| 90 ColumnVector maxx(nelem); | |
| 91 ColumnVector miny(nelem); | |
| 92 ColumnVector maxy(nelem); | |
| 93 for(octave_idx_type k = 0; k < nelem; k++) | |
| 94 { | |
| 95 minx(k) = min(REF(x,k,0), REF(x,k,1), REF(x,k,2)) - eps; | |
| 96 maxx(k) = max(REF(x,k,0), REF(x,k,1), REF(x,k,2)) + eps; | |
| 97 miny(k) = min(REF(y,k,0), REF(y,k,1), REF(y,k,2)) - eps; | |
| 98 maxy(k) = max(REF(y,k,0), REF(y,k,1), REF(y,k,2)) + eps; | |
| 99 } | |
| 100 | |
| 101 const octave_idx_type np = xi.length(); | |
| 102 ColumnVector values(np); | |
| 103 | |
| 104 double x0 = 0.0, y0 = 0.0; | |
| 105 double a11 = 0.0, a12 = 0.0, a21 = 0.0, a22 = 0.0, 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 if (xt >= minx(k) && xt <= maxx(k) && | |
| 132 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 retval(0) = values; | |
| 162 | |
| 163 return retval; | |
| 164 } | |
| 165 | |
| 166 /* | |
| 167 %!shared x, y, tri | |
| 168 %! x = [-1;-1;1]; | |
| 169 %! y = [-1;1;-1]; | |
| 170 %! tri = [1, 2, 3]; | |
| 171 %!error (tsearch()) | |
| 172 %!assert (tsearch (x,y,tri,-1,-1), 1) | |
| 173 %!assert (tsearch (x,y,tri, 1,-1), 1) | |
| 174 %!assert (tsearch (x,y,tri,-1, 1), 1) | |
| 175 %!assert (tsearch (x,y,tri,-1/3, -1/3), 1) | |
| 176 %!assert (tsearch (x,y,tri, 1, 1), NaN) | |
| 177 | |
| 178 */ |