Mercurial > octave-nkf

comparison src/DLD-FUNCTIONS/tsearch.cc @ 10154:40dfc0c99116

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
DLD-FUNCTIONS/*.cc: untabify
author John W. Eaton <jwe@octave.org>
date Wed, 20 Jan 2010 17:33:41 -0500
parents 2c279308f6ab
children c48b7048e720
comparison
equal deleted inserted replaced
10153:2c28f9d0360f 10154:40dfc0c99116
63 63
64 // e.g., build up a neighbouring triangle structure and use a simplex-like 64 // e.g., build up a neighbouring triangle structure and use a simplex-like
65 // method to traverse it 65 // method to traverse it
66 66
67 DEFUN_DLD (tsearch, args, , 67 DEFUN_DLD (tsearch, args, ,
68 "-*- texinfo -*-\n\ 68 "-*- texinfo -*-\n\
69 @deftypefn {Loadable Function} {@var{idx} =} tsearch (@var{x}, @var{y}, @var{t}, @var{xi}, @var{yi})\n\ 69 @deftypefn {Loadable Function} {@var{idx} =} tsearch (@var{x}, @var{y}, @var{t}, @var{xi}, @var{yi})\n\
70 Searches for the enclosing Delaunay convex hull. For @code{@var{t} =\n\ 70 Searches for the enclosing Delaunay convex hull. For @code{@var{t} =\n\
71 delaunay (@var{x}, @var{y})}, finds the index in @var{t} containing the\n\ 71 delaunay (@var{x}, @var{y})}, finds the index in @var{t} containing the\n\
72 points @code{(@var{xi}, @var{yi})}. For points outside the convex hull,\n\ 72 points @code{(@var{xi}, @var{yi})}. For points outside the convex hull,\n\
73 @var{idx} is NaN.\n\ 73 @var{idx} is NaN.\n\
117 const double xt = xi(kp); 117 const double xt = xi(kp);
118 const double yt = yi(kp); 118 const double yt = yi(kp);
119 119
120 // check if last triangle contains the next point 120 // check if last triangle contains the next point
121 if (k < nelem) 121 if (k < nelem)
122 { 122 {
123 const double dx1 = xt - x0; 123 const double dx1 = xt - x0;
124 const double dx2 = yt - y0; 124 const double dx2 = yt - y0;
125 const double c1 = (a22 * dx1 - a21 * dx2) / det; 125 const double c1 = (a22 * dx1 - a21 * dx2) / det;
126 const double c2 = (-a12 * dx1 + a11 * dx2) / det; 126 const double c2 = (-a12 * dx1 + a11 * dx2) / det;
127 if ( c1 >= -eps && c2 >= -eps && (c1 + c2) <= (1 + eps)) 127 if ( c1 >= -eps && c2 >= -eps && (c1 + c2) <= (1 + eps))
128 { 128 {
129 values(kp) = double(k+1); 129 values(kp) = double(k+1);
130 continue; 130 continue;
131 } 131 }
132 } 132 }
133 133
134 // it doesn't, so go through all elements 134 // it doesn't, so go through all elements
135 for (k = 0; k < nelem; k++) 135 for (k = 0; k < nelem; k++)
136 { 136 {
137 OCTAVE_QUIT; 137 OCTAVE_QUIT;
138 if (xt >= minx(k) && xt <= maxx(k) && 138 if (xt >= minx(k) && xt <= maxx(k) &&
139 yt >= miny(k) && yt <= maxy(k) ) 139 yt >= miny(k) && yt <= maxy(k) )
140 { 140 {
141 // element inside the minimum rectangle: examine it closely 141 // element inside the minimum rectangle: examine it closely
142 x0 = REF(x,k,0); 142 x0 = REF(x,k,0);
143 y0 = REF(y,k,0); 143 y0 = REF(y,k,0);
144 a11 = REF(x,k,1)-x0; 144 a11 = REF(x,k,1)-x0;
145 a12 = REF(y,k,1)-y0; 145 a12 = REF(y,k,1)-y0;
146 a21 = REF(x,k,2)-x0; 146 a21 = REF(x,k,2)-x0;
147 a22 = REF(y,k,2)-y0; 147 a22 = REF(y,k,2)-y0;
148 det = a11 * a22 - a21 * a12; 148 det = a11 * a22 - a21 * a12;
149 149
150 // solve the system 150 // solve the system
151 const double dx1 = xt - x0; 151 const double dx1 = xt - x0;
152 const double dx2 = yt - y0; 152 const double dx2 = yt - y0;
153 const double c1 = (a22 * dx1 - a21 * dx2) / det; 153 const double c1 = (a22 * dx1 - a21 * dx2) / det;
154 const double c2 = (-a12 * dx1 + a11 * dx2) / det; 154 const double c2 = (-a12 * dx1 + a11 * dx2) / det;
155 if ((c1 >= -eps) && (c2 >= -eps) && ((c1 + c2) <= (1 + eps))) 155 if ((c1 >= -eps) && (c2 >= -eps) && ((c1 + c2) <= (1 + eps)))
156 { 156 {
157 values(kp) = double(k+1); 157 values(kp) = double(k+1);
158 break; 158 break;
159 } 159 }
160 } //endif # examine this element closely 160 } //endif # examine this element closely
161 } //endfor # each element 161 } //endfor # each element
162 162
163 if (k == nelem) 163 if (k == nelem)
164 values(kp) = lo_ieee_nan_value (); 164 values(kp) = lo_ieee_nan_value ();
165 165
166 } //endfor # kp 166 } //endfor # kp
167 167
168 retval(0) = values; 168 retval(0) = values;
169 169