--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/DLD-FUNCTIONS/tsearch.cc	Fri Aug 24 08:27:29 2007 +0000
@@ -0,0 +1,177 @@
+/*
+
+Copyright (C) 2002 Andreas Stahel <Andreas.Stahel@hta-bi.bfh.ch>
+
+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 2, 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, write to the Free
+Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+02110-1301, USA.
+
+*/
+
+#include <iostream>
+#include <fstream>
+#include <cmath>
+#include <string>
+
+#include "oct.h"
+#include "parse.h"
+#include "lo-ieee.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_DLD (tsearch, args, ,
+	"-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {@var{idx} =} tsearch (@var{x}, @var{y}, @var{t}, @var{xi}, @var{yi})\n\
+Searches for the enclosing Delaunay convex hull. For @code{@var{t} =\n\
+delaunay (@var{x}, @var{y})}, finds the index in @var{t} containing the\n\
+points @code{(@var{xi}, @var{yi})}. For points outside the convex hull,\n\
+@var{idx} is NaN.\n\
+@seealso{delaunay, delaunayn}\n\
+@end deftypefn")
+{
+  const double eps=1.0e-12;
+
+  octave_value_list retval;
+  const int nargin = args.length ();
+  if ( nargin != 5 ) {
+    print_usage ();
+    return retval;
+  }
+  
+  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());
+
+  if (error_state) return retval;
+
+  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.length();
+  ColumnVector values(np);
+
+  double x0 = 0.0, y0 = 0.0;
+  double a11 = 0.0, a12 = 0.0, a21 = 0.0, 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
+  
+  retval(0) = values;
+  
+  return retval;
+}
+
+/*
+%!shared x, y, tri
+%! x = [-1;-1;1];
+%! y = [-1;1;-1];
+%! tri = [1, 2, 3];
+%!error (tsearch())
+%!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)
+
+*/