Mercurial > octave
view libinterp/corefcn/tsearch.cc @ 20939:b17fda023ca6
maint: Use new C++ archetype in more files.
Place input validation first in files.
Move declaration of retval down in function to be closer to point of usage.
Eliminate else clause after if () error.
Use "return ovl()" where it makes sense.
* find.cc, gammainc.cc, gcd.cc, getgrent.cc, getpwent.cc, givens.cc,
graphics.cc, help.cc, hess.cc, hex2num.cc, input.cc, kron.cc, load-path.cc,
load-save.cc, lookup.cc, mappers.cc, matrix_type.cc, mgorth.cc, nproc.cc,
ordschur.cc, pager.cc, pinv.cc, pr-output.cc, profiler.cc, psi.cc, quad.cc,
rcond.cc, regexp.cc, schur.cc, sighandlers.cc, sparse.cc, str2double.cc,
strfind.cc, strfns.cc, sub2ind.cc, svd.cc, sylvester.cc, symtab.cc,
syscalls.cc, sysdep.cc, time.cc, toplev.cc, tril.cc, tsearch.cc, typecast.cc,
urlwrite.cc, utils.cc, variables.cc, __delaunayn__.cc, __eigs__.cc,
__glpk__.cc, __magick_read__.cc, __osmesa_print__.cc, __voronoi__.cc, amd.cc,
audiodevinfo.cc, audioread.cc, chol.cc, colamd.cc, dmperm.cc, fftw.cc, qr.cc,
symbfact.cc, symrcm.cc, ov-bool-mat.cc, ov-cell.cc, ov-class.cc,
ov-classdef.cc, ov-fcn-handle.cc, ov-fcn-inline.cc, ov-flt-re-mat.cc,
ov-java.cc, ov-null-mat.cc, ov-oncleanup.cc, ov-re-mat.cc, ov-struct.cc,
ov-typeinfo.cc, ov-usr-fcn.cc, ov.cc, octave.cc:
Use new C++ archetype in more files.
author | Rik <rik@octave.org> |
---|---|
date | Fri, 18 Dec 2015 15:37:22 -0800 |
parents | 1142cf6abc0d |
children | 48b2ad5ee801 |
line wrap: on
line source
/* Copyright (C) 2002-2015 Andreas Stahel This file is part of Octave. Octave is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version. Octave is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with Octave; see the file COPYING. If not, see <http://www.gnu.org/licenses/>. */ // Author: Andreas Stahel <Andreas.Stahel@hta-bi.bfh.ch> #ifdef HAVE_CONFIG_H #include <config.h> #endif #include "lo-ieee.h" #include "lo-math.h" #include "defun.h" #include "error.h" #include "oct-obj.h" inline double max (double a, double b, double c) { if (a < b) return (b < c ? c : b); else return (a < c ? c : a); } inline double min (double a, double b, double c) { if (a > b) return (b > c ? c : b); else return (a > c ? c : a); } #define REF(x,k,i) x(static_cast<octave_idx_type>(elem((k), (i))) - 1) // for large data set the algorithm is very slow // one should presort (how?) either the elements of the points of evaluation // to cut down the time needed to decide which triangle contains the // given point // e.g., build up a neighbouring triangle structure and use a simplex-like // method to traverse it DEFUN (tsearch, args, , "-*- texinfo -*-\n\ @deftypefn {} {@var{idx} =} tsearch (@var{x}, @var{y}, @var{t}, @var{xi}, @var{yi})\n\ Search for the enclosing Delaunay convex hull.\n\ \n\ For @code{@var{t} = delaunay (@var{x}, @var{y})}, finds the index in @var{t}\n\ containing the points @code{(@var{xi}, @var{yi})}. For points outside the\n\ convex hull, @var{idx} is NaN.\n\ @seealso{delaunay, delaunayn}\n\ @end deftypefn") { if (args.length () != 5) print_usage (); const double eps = 1.0e-12; const ColumnVector x (args(0).vector_value ()); const ColumnVector y (args(1).vector_value ()); const Matrix elem (args(2).matrix_value ()); const ColumnVector xi (args(3).vector_value ()); const ColumnVector yi (args(4).vector_value ()); const octave_idx_type nelem = elem.rows (); ColumnVector minx (nelem); ColumnVector maxx (nelem); ColumnVector miny (nelem); ColumnVector maxy (nelem); for (octave_idx_type k = 0; k < nelem; k++) { minx(k) = min (REF (x, k, 0), REF (x, k, 1), REF (x, k, 2)) - eps; maxx(k) = max (REF (x, k, 0), REF (x, k, 1), REF (x, k, 2)) + eps; miny(k) = min (REF (y, k, 0), REF (y, k, 1), REF (y, k, 2)) - eps; maxy(k) = max (REF (y, k, 0), REF (y, k, 1), REF (y, k, 2)) + eps; } const octave_idx_type np = xi.numel (); ColumnVector values (np); double x0, y0, a11, a12, a21, a22, det; x0 = y0 = 0.0; a11 = a12 = a21 = a22 = 0.0; det = 0.0; octave_idx_type k = nelem; // k is a counter of elements for (octave_idx_type kp = 0; kp < np; kp++) { const double xt = xi(kp); const double yt = yi(kp); // check if last triangle contains the next point if (k < nelem) { const double dx1 = xt - x0; const double dx2 = yt - y0; const double c1 = (a22 * dx1 - a21 * dx2) / det; const double c2 = (-a12 * dx1 + a11 * dx2) / det; if (c1 >= -eps && c2 >= -eps && (c1 + c2) <= (1 + eps)) { values(kp) = double(k+1); continue; } } // it doesn't, so go through all elements for (k = 0; k < nelem; k++) { OCTAVE_QUIT; if (xt >= minx(k) && xt <= maxx(k) && yt >= miny(k) && yt <= maxy(k)) { // element inside the minimum rectangle: examine it closely x0 = REF (x, k, 0); y0 = REF (y, k, 0); a11 = REF (x, k, 1) - x0; a12 = REF (y, k, 1) - y0; a21 = REF (x, k, 2) - x0; a22 = REF (y, k, 2) - y0; det = a11 * a22 - a21 * a12; // solve the system const double dx1 = xt - x0; const double dx2 = yt - y0; const double c1 = (a22 * dx1 - a21 * dx2) / det; const double c2 = (-a12 * dx1 + a11 * dx2) / det; if ((c1 >= -eps) && (c2 >= -eps) && ((c1 + c2) <= (1 + eps))) { values(kp) = double(k+1); break; } } //endif # examine this element closely } //endfor # each element if (k == nelem) values(kp) = lo_ieee_nan_value (); } //endfor # kp return ovl (values); } /* %!shared x, y, tri %! x = [-1;-1;1]; %! y = [-1;1;-1]; %! tri = [1, 2, 3]; %!assert (tsearch (x,y,tri,-1,-1), 1) %!assert (tsearch (x,y,tri, 1,-1), 1) %!assert (tsearch (x,y,tri,-1, 1), 1) %!assert (tsearch (x,y,tri,-1/3, -1/3), 1) %!assert (tsearch (x,y,tri, 1, 1), NaN) %!error tsearch () */