# HG changeset patch
# User dbateman
# Date 1187944049 0
# Node ID 9fddcc5860653d78fb1cc5a82b9c6edb62137e93
# Parent  5c4702203cc4e6822e07b55d3d75b28f54b0060d
[project @ 2007-08-24 08:27:27 by dbateman]

diff -r 5c4702203cc4 -r 9fddcc586065 ChangeLog
--- a/ChangeLog	Thu Aug 23 19:43:08 2007 +0000
+++ b/ChangeLog	Fri Aug 24 08:27:29 2007 +0000
@@ -1,3 +1,10 @@
+2007-08-24  David Bateman  <dbateman@free.fr>
+
+        * configure.in: Probe for the use of Qhull.
+        * aclocal.m4 (AC_CHECK_QHULL_VERSION): Macro to check whether
+        Qhull needs a version number.
+        * Makeconf.in: Add QHULL_LIBS.
+
 2007-08-23  John W. Eaton  <jwe@octave.org>
 
 	* aclocal.m4 (OCTAVE_PROG_SED): Don't clobber value from environment.
diff -r 5c4702203cc4 -r 9fddcc586065 Makeconf.in
--- a/Makeconf.in	Thu Aug 23 19:43:08 2007 +0000
+++ b/Makeconf.in	Fri Aug 24 08:27:29 2007 +0000
@@ -207,6 +207,7 @@
 LIBREADLINE = @LIBREADLINE@
 TERMLIBS = @TERMLIBS@
 
+QHULL_LIBS = @QHULL_LIBS@
 REGEX_LIBS = @REGEX_LIBS@
 BLAS_LIBS = @BLAS_LIBS@
 FFTW_LIBS = @FFTW_LIBS@
diff -r 5c4702203cc4 -r 9fddcc586065 aclocal.m4
--- a/aclocal.m4	Thu Aug 23 19:43:08 2007 +0000
+++ b/aclocal.m4	Fri Aug 24 08:27:29 2007 +0000
@@ -987,3 +987,40 @@
   if test "$octave_cv_hdf5_dll" = yes; then
     AC_DEFINE(_HDF5USEDLL_, 1, [Define if using HDF5 dll (Win32)])
   fi])
+dnl
+dnl Check for the QHull version.
+dnl
+AC_DEFUN(AC_CHECK_QHULL_VERSION,
+[AC_MSG_CHECKING([for qh_qhull in -lqhull with qh_version])
+AC_CACHE_VAL(octave_cv_lib_qhull_version,  [
+cat > conftest.c <<EOF
+#include <stdio.h>
+char *qh_version = "version";
+char qh_qhull();
+int
+main(argc, argv)
+  int argc;
+  char **argv;
+{
+  qh_qhull();
+  return 0;
+}
+EOF
+
+octave_qhull_try="${CC-cc} $CFLAGS $CPPFLAGS $LDFLAGS conftest.c -o conftest -lqhull $LIBS"
+if AC_TRY_EVAL(octave_qhull_try) && test -s conftest ; then
+    octave_cv_lib_qhull_version=yes
+else
+    octave_cv_lib_qhull_version=no
+fi
+rm -f conftest.c conftest.o conftest
+])dnl
+if test "$octave_cv_lib_qhull_version" = "yes"; then
+  AC_MSG_RESULT(yes)
+  ifelse([$1], , , [$1])
+else
+  AC_MSG_RESULT(no)
+  ifelse([$2], , , [$2])
+fi
+])
+
diff -r 5c4702203cc4 -r 9fddcc586065 configure.in
--- a/configure.in	Thu Aug 23 19:43:08 2007 +0000
+++ b/configure.in	Fri Aug 24 08:27:29 2007 +0000
@@ -29,7 +29,7 @@
 EXTERN_CXXFLAGS="$CXXFLAGS"
 
 AC_INIT
-AC_REVISION($Revision: 1.570 $)
+AC_REVISION($Revision: 1.571 $)
 AC_PREREQ(2.57)
 AC_CONFIG_SRCDIR([src/octave.cc])
 AC_CONFIG_HEADER(config.h)
@@ -424,6 +424,24 @@
   ;;
 esac
 
+### Check for the QHull library
+AC_SUBST(QHULL_LIBS)
+AC_CHECK_HEADER(qhull/qhull_a.h, have_qhull=yes, have_qhull=no)
+if test "$have_qhull" = yes; then
+  AC_CHECK_LIB(qhull, qh_qhull, have_qhull=yes, have_qhull=no)
+  if test "$have_qhull" != yes; then
+    AC_CHECK_QHULL_VERSION(have_qhull=yes, have_qhull=no)
+    AC_DEFINE(NEED_QHULL_VERSION, 1, [Define if the QHull library needs a wh_version variable defined.])
+  fi
+fi
+if test "$have_qhull" = yes; then
+  AC_DEFINE(HAVE_QHULL, 1, [Define if the QHull library is used.])
+  QHULL_LIBS="-lqhull"
+else
+  warn_qhull="Qhull library not found --- This will result in loss of functionality of some geometry functions."
+  AC_MSG_WARN($warn_qhull)
+fi
+
 ### Check for pcre/regex library.
 AC_SUBST(REGEX_LIBS)
 WITH_PCRE_CONFIG=no
diff -r 5c4702203cc4 -r 9fddcc586065 doc/interpreter/geometry.txi
--- a/doc/interpreter/geometry.txi	Thu Aug 23 19:43:08 2007 +0000
+++ b/doc/interpreter/geometry.txi	Fri Aug 24 08:27:29 2007 +0000
@@ -1,8 +1,215 @@
 @c Copyright (C) 2007 John W. Eaton
+@c Copyright (C) 2007 David Bateman
 @c This is part of the Octave manual.
 @c For copying conditions, see the file gpl.texi.
 
 @node Geometry
 @chapter Geometry
 
+@menu
+* Delaunay Triangulation::
+* Voronoi Diagrams::
+* Convex Hull::
+* Plotting the Triangulation::
+* Interpolation on Scattered Data::
+@end menu
+
+@node Delaunay Triangulation
+@section Delaunay Triangulation
+
+@DOCSTRING(delaunay)
+
+The 3- and N-dimensional extension of the Delaunay triangulation are
+given by @code{delaunay3} and @code{delaunayn} respectively.  
+@code{delaunay3} returns a set of tetrahedra that satisfy the
+Delaunay circum-circle criteria.  Similarly, @code{delaunayn} returns the
+N-dimensional simplex satisfying the Delaunay circum-circle criteria.  
+The N-dimensional extension of a triangulation is called a tesselation.
+
+@DOCSTRING(delaunay3)
+
+@DOCSTRING(delaunayn)
+
+@menu
+* Identifying points in Triangulation::
+@end menu
+
+@node Identifying points in Triangulation
+@subsection Identifying points in Triangulation
+
+It is often necessary to identify whether a particular point in the
+N-dimensional space is within the Delaunay tesselation of a set of
+points in this N-dimensional space, and if so which N-Simplex contains
+the point and which point in the tesselation is closest to the desired
+point.  The functions @code{tsearch} and @code{dsearch} perform this
+function in a triangulation, and @code{tsearchn} and @code{dsearchn} in
+an N-dimensional tesselation.
+
+To identify whether a particular point represented by a vector @var{p}
+falls within one of the simplices of an N-Simplex, we can write the
+Cartesian coordinates of the point in a parametric form with respect to
+the N-Simplex.  This parametric form is called the Barycentric
+Coordinates of the point.  If the points defining the N-Simplex are given
+by @code{@var{N} + 1} vectors @var{t}(@var{i},:), then the Barycentric
+coordinates defining the point @var{p} is given by
+
+@example
+@var{p} = sum (@var{beta}(1:@var{N}+1) * @var{t}(1:@var{N}+1),:)
+@end example
+
+@noindent
+where there are @code{@var{N} + 1} values @code{@var{beta}(@var{i})}
+that together as a vector represent the Barycentric coordinates of the
+point @var{p}.  To ensure a unique solution for the values of
+@code{@var{beta}(@var{i})} an additional criteria of
+
+@example
+sum (@var{beta}(1:@var{N}+1)) == 1
+@end example
+
+@noindent
+is imposed, and we can therefore write the above as
+
+@example
+@var{p} - @var{t}(end, :) = @var{beta}(1:end-1) * (@var{t}(1:end-1, :)
+      - ones(@var{N}, 1) * @var{t}(end, :)
+@end example
+
+@noindent
+Solving for @var{beta} we can then write
+
+@example
+@var{beta}(1:end-1) = (@var{p} - @var{t}(end, :)) / (@var{t}(1:end-1, :)
+      - ones(@var{N}, 1) * @var{t}(end, :))
+@var{beta}(end) = sum(@var{beta}(1:end-1))
+@end example
+
+@noindent
+which gives the formula for the conversion of the Cartesian coordinates
+of the point @var{p} to the Barycentric coordinates @var{beta}.  An
+important property of the Barycentric coordinates is that for all points
+in the N-Simplex
+
+@example
+0 <= @var{beta}(@var{i}) <= 1
+@end example
+
+@noindent
+Therefore, the test in @code{tsearch} and @code{tsearchn} essentially
+only needs to express each point in terms of the Barycentric coordinates
+of each of the simplices of the N-Simplex and test the values of
+@var{beta}. This is exactly the implementation used in
+@code{tsearchn}. @code{tsearch} is optimized for 2-dimensions and the
+Barycentric coordinates are not explicitly formed.
+
+@DOCSTRING(tsearch)
+
+@DOCSTRING(tsearchn)
+
+An example of the use of @code{tsearch} can be seen with the simple
+triangulation
+
+@example
+@group
+@var{x} = [-1; -1; 1; 1];
+@var{y} = [-1; 1; -1; 1];
+@var{tri} = [1, 2, 3; 2, 3, 1];
+@end group
+@end example
+
+@noindent
+consisting of two triangles defined by @var{tri}. We can then identify
+which triangle a point falls in like
+
+@example
+@group
+tsearch (@var{x}, @var{y}, @var{tri}, -0.5, -0.5)
+@result{} 1
+tsearch (@var{x}, @var{y}, @var{tri}, 0.5, 0.5)
+@result{} 2
+@end group
+@end example
+
+@noindent
+and we can confirm that a point doesn't lie within one of the triangles like
+
+@example
+@group
+tsearch (@var{x}, @var{y}, @var{tri}, 2, 2)
+@result{} NaN
+@end group
+@end example
+
+The @code{dsearch} and @code{dsearchn} find the closest point in a
+tessellation to the desired point. The desired point does not
+necessarily have to be in the tessellation, and even if it the returned
+point of the tessellation does not have to be one of the vertices of the
+N-simplex within which the desired point is found.
+
+@DOCSTRING(dsearch)
+
+@DOCSTRING(dsearchn)
+
+An example of the use of @code{dsearch}, using the above values of
+@var{x}, @var{y} and @var{tri} is
+
+@example
+@group
+dsearch (@var{x}, @var{y}, @var{tri}, -2, -2)
+@result{} 1
+@end group
+@end example
+
+If you wish the points that are outside the tessellation to be flagged,
+then @code{dsearchn} can be used as
+
+@example
+@group
+dsearchn ([@var{x}, @var{y}], @var{tri}, [-2, -2], NaN)
+@result{} NaN
+dsearchn ([@var{x}, @var{y}], @var{tri}, [-0.5, -0.5], NaN)
+@result{} 1
+@end group
+@end example
+
+@noindent
+where the point outside the tessellation are then flagged with @code{NaN}.
+
+@node Voronoi Diagrams
+@section Voronoi Diagrams
+
+A Voronoi diagram or Voronoi tessellation of a set of points @var{s} in
+an N-dimensional space, is the tessellation of the N-dimensional space
+such that all points in @code{@var{v}(@var{p})}, a partitions of the
+tessellation where @var{p} is a member of @var{s}, are closer to @var{p}
+than any other point in @var{s}.  The Voronoi diagram is related to the
+Delaunay triangulation of a set of points.
+
+@DOCSTRING(voronoi)
+
+@DOCSTRING(voronoin)
+
+@node Convex Hull
+@section Convex Hull
+
+@DOCSTRING(convhull)
+
+@DOCSTRING(convhulln)
+
+@node Plotting the Triangulation
+@section Plotting the Triangulation
+
+@DOCSTRING(triplot)
+
+@DOCSTRING(trimesh)
+
 @DOCSTRING(polyarea)
+
+@node Interpolation on Scattered Data
+@section Interpolation on Scattered Data
+
+@DOCSTRING(griddata)
+
+@DOCSTRING(griddata3)
+
+@DOCSTRING(griddatan)
diff -r 5c4702203cc4 -r 9fddcc586065 doc/interpreter/octave.texi
--- a/doc/interpreter/octave.texi	Thu Aug 23 19:43:08 2007 +0000
+++ b/doc/interpreter/octave.texi	Fri Aug 24 08:27:29 2007 +0000
@@ -446,9 +446,18 @@
 * Set Operations::
 
 Interpolation
+
 * One-dimensional Interpolation::
 * Multi-dimensional Interpolation::
 
+Geometry
+
+* Delaunay Triangulation::
+* Voronoi Diagrams::
+* Convex Hull::
+* Plotting the Triangulation::
+* Interpolation on Scattered Data::
+
 Control Theory
 
 * sysstruct::                   
diff -r 5c4702203cc4 -r 9fddcc586065 liboctave/CSparse.cc
--- a/liboctave/CSparse.cc	Thu Aug 23 19:43:08 2007 +0000
+++ b/liboctave/CSparse.cc	Fri Aug 24 08:27:29 2007 +0000
@@ -536,6 +536,20 @@
   return *this;
 }
 
+SparseComplexMatrix&
+SparseComplexMatrix::insert (const SparseMatrix& a, const Array<octave_idx_type>& indx)
+{
+  SparseComplexMatrix tmp (a);
+  return insert (tmp /*a*/, indx);
+}
+
+SparseComplexMatrix&
+SparseComplexMatrix::insert (const SparseComplexMatrix& a, const Array<octave_idx_type>& indx)
+{
+  MSparse<Complex>::insert (a, indx);
+  return *this;
+}
+
 SparseComplexMatrix
 SparseComplexMatrix::concat (const SparseComplexMatrix& rb, 
 			     const Array<octave_idx_type>& ra_idx)
diff -r 5c4702203cc4 -r 9fddcc586065 liboctave/CSparse.h
--- a/liboctave/CSparse.h	Thu Aug 23 19:43:08 2007 +0000
+++ b/liboctave/CSparse.h	Fri Aug 24 08:27:29 2007 +0000
@@ -51,6 +51,9 @@
 
   SparseComplexMatrix (octave_idx_type r, octave_idx_type c) : MSparse<Complex> (r, c) { }
 
+  SparseComplexMatrix (const dim_vector& dv, octave_idx_type nz = 0) : 
+    MSparse<Complex> (dv, nz) { }
+
   explicit SparseComplexMatrix (octave_idx_type r, octave_idx_type c, Complex val) 
     : MSparse<Complex> (r, c, val) { }
 
@@ -107,6 +110,8 @@
 
   SparseComplexMatrix& insert (const SparseComplexMatrix& a, octave_idx_type r, octave_idx_type c);
   SparseComplexMatrix& insert (const SparseMatrix& a, octave_idx_type r, octave_idx_type c);
+  SparseComplexMatrix& insert (const SparseComplexMatrix& a, const Array<octave_idx_type>& indx);
+  SparseComplexMatrix& insert (const SparseMatrix& a, const Array<octave_idx_type>& indx);
 
   SparseComplexMatrix concat (const SparseComplexMatrix& rb,
 			      const Array<octave_idx_type>& ra_idx);
diff -r 5c4702203cc4 -r 9fddcc586065 liboctave/ChangeLog
--- a/liboctave/ChangeLog	Thu Aug 23 19:43:08 2007 +0000
+++ b/liboctave/ChangeLog	Fri Aug 24 08:27:29 2007 +0000
@@ -1,3 +1,27 @@
+2007-08-24  David Bateman  <dbateman@free.fr>
+
+        * MSparse.h (MSparse<T>& insert (const Sparse<T>&, 
+        const Array<octave_idx_type>&)): New method.
+        (MSparse (const dim_vector&, octave_idx_type)): Ditto.
+        * dSparse.h (SparseMatrix& SparseMatrix::insert (const
+        SparseMatrix&, const Array<octave_idx_type>&)): ditto.
+        (SparseMatrix (const dim_vector&, octave_idx_type)): ditto.
+        * dSparse.cc (SparseMatrix& SparseMatrix::insert (const
+        SparseMatrix&, const Array<octave_idx_type>&)): ditto.
+        * boolSparse.h (SparseBoolMatrix& SparseBoolMatrix::insert (const
+        SparseBoolMatrix&, const Array<octave_idx_type>&)): ditto.
+        * boolSparse.cc (SparseBoolMatrix& SparseBoolMatrix::insert (const
+        SparseBoolMatrix&, const Array<octave_idx_type>&)): ditto.
+        * CSparse.h (SparseComplexMatrix& SparseComplexMatrix::insert (const
+        SparseMatrix&, const Array<octave_idx_type>&),
+        SparseComplexMatrix& SparseComplexMatrix::insert (const
+        SparseComplexMatrix&, const Array<octave_idx_type>&)): ditto.
+        (SparseComplexMatrix (const dim_vector&, octave_idx_type)): ditto.
+        * CSparse.cc (SparseComplexMatrix& SparseComplexMatrix::insert (const
+        SparseMatrix&, const Array<octave_idx_type>&),
+        SparseComplexMatrix& SparseComplexMatrix::insert (const
+        SparseComplexMatrix&, const Array<octave_idx_type>&)): ditto.
+
 2007-08-19  David Bateman  <dbateman@free.fr>
 
 	* Sparse.cc (Sparse<T>::permute): Avoid shadowing warning.
diff -r 5c4702203cc4 -r 9fddcc586065 liboctave/MSparse.h
--- a/liboctave/MSparse.h	Thu Aug 23 19:43:08 2007 +0000
+++ b/liboctave/MSparse.h	Fri Aug 24 08:27:29 2007 +0000
@@ -45,6 +45,9 @@
 
   MSparse (octave_idx_type n, octave_idx_type m) : Sparse<T> (n, m) { }
 
+  MSparse (const dim_vector& dv, octave_idx_type nz = 0) : 
+    Sparse<T> (dv, nz) { }
+
   MSparse (const MSparse<T>& a) : Sparse<T> (a) { }
 
   MSparse (const MSparse<T>& a, const dim_vector& dv) : Sparse<T> (a, dv) { }
@@ -79,6 +82,12 @@
     return *this;
   }
 
+  MSparse<T>& insert (const Sparse<T>& a, const Array<octave_idx_type>& indx)
+  {
+    Sparse<T>::insert (a, indx);
+    return *this;
+  }
+
   MSparse<T> transpose (void) const { return Sparse<T>::transpose (); }
 
   MSparse<T> squeeze (void) const { return Sparse<T>::squeeze (); }
diff -r 5c4702203cc4 -r 9fddcc586065 liboctave/boolSparse.cc
--- a/liboctave/boolSparse.cc	Thu Aug 23 19:43:08 2007 +0000
+++ b/liboctave/boolSparse.cc	Fri Aug 24 08:27:29 2007 +0000
@@ -73,6 +73,13 @@
   return *this;
 }
 
+SparseBoolMatrix&
+SparseBoolMatrix::insert (const SparseBoolMatrix& a, const Array<octave_idx_type>& indx)
+{
+  Sparse<bool>::insert (a, indx);
+  return *this;
+}
+
 SparseBoolMatrix
 SparseBoolMatrix::concat (const SparseBoolMatrix& rb, const Array<octave_idx_type>& ra_idx)
 {
diff -r 5c4702203cc4 -r 9fddcc586065 liboctave/boolSparse.h
--- a/liboctave/boolSparse.h	Thu Aug 23 19:43:08 2007 +0000
+++ b/liboctave/boolSparse.h	Fri Aug 24 08:27:29 2007 +0000
@@ -40,6 +40,9 @@
   explicit SparseBoolMatrix (octave_idx_type r, octave_idx_type c, bool val) 
     : Sparse<bool> (r, c, val) { }
 
+  SparseBoolMatrix (const dim_vector& dv, octave_idx_type nz = 0) : 
+    Sparse<bool> (dv, nz) { }
+
   SparseBoolMatrix (const Sparse<bool>& a) : Sparse<bool> (a) { }
 
   SparseBoolMatrix (const SparseBoolMatrix& a) : Sparse<bool> (a) { }
@@ -79,6 +82,8 @@
 
   SparseBoolMatrix& insert (const SparseBoolMatrix& a, octave_idx_type r, octave_idx_type c);
 
+  SparseBoolMatrix& insert (const SparseBoolMatrix& a, const Array<octave_idx_type>& indx);
+
   SparseBoolMatrix concat (const SparseBoolMatrix& rb, 
 			   const Array<octave_idx_type>& ra_idx);
 
diff -r 5c4702203cc4 -r 9fddcc586065 liboctave/dSparse.cc
--- a/liboctave/dSparse.cc	Thu Aug 23 19:43:08 2007 +0000
+++ b/liboctave/dSparse.cc	Fri Aug 24 08:27:29 2007 +0000
@@ -216,6 +216,13 @@
   return *this;
 }
 
+SparseMatrix&
+SparseMatrix::insert (const SparseMatrix& a, const Array<octave_idx_type>& indx)
+{
+  MSparse<double>::insert (a, indx);
+  return *this;
+}
+
 SparseMatrix
 SparseMatrix::max (int dim) const
 {
diff -r 5c4702203cc4 -r 9fddcc586065 liboctave/dSparse.h
--- a/liboctave/dSparse.h	Thu Aug 23 19:43:08 2007 +0000
+++ b/liboctave/dSparse.h	Fri Aug 24 08:27:29 2007 +0000
@@ -50,6 +50,9 @@
 
   SparseMatrix (octave_idx_type r, octave_idx_type c) : MSparse<double> (r, c) { }
 
+  SparseMatrix (const dim_vector& dv, octave_idx_type nz = 0) : 
+    MSparse<double> (dv, nz) { }
+
   explicit SparseMatrix (octave_idx_type r, octave_idx_type c, double val) 
     : MSparse<double> (r, c, val) { }
 
@@ -98,6 +101,8 @@
 
   SparseMatrix& insert (const SparseMatrix& a, octave_idx_type r, octave_idx_type c);
 
+  SparseMatrix& insert (const SparseMatrix& a, const Array<octave_idx_type>& indx);
+
   SparseMatrix concat (const SparseMatrix& rb, const Array<octave_idx_type>& ra_idx);
   SparseComplexMatrix concat (const SparseComplexMatrix& rb,
 			      const Array<octave_idx_type>& ra_idx);
diff -r 5c4702203cc4 -r 9fddcc586065 scripts/ChangeLog
--- a/scripts/ChangeLog	Thu Aug 23 19:43:08 2007 +0000
+++ b/scripts/ChangeLog	Fri Aug 24 08:27:29 2007 +0000
@@ -1,3 +1,16 @@
+2007-08-24  David Bateman  <dbateman@free.fr>
+
+        * geometry/convhull.m, geometry/delaunay.m, geometry/delaunay3.m, 
+        geometry/griddata.m, geometry/voronoi.m, geometry/voronoin.m: New 
+        functions ported from octave-forge.
+        * geometry/delaunayn.m, geometry/dsearch.m, geometry/dsearchn.m,
+        geometry/griddata3.m, geometry/griddatan.m, geometry/trimesh.m, 
+        geometry/triplot.m, geometry/tsearchn.m:
+        New functions.
+        * geometry/voronoi.m: Remove duplicate edges from Voronoi diagram.
+        * geometry/Makefile.in (SOURCES): Add functions above.
+        * configure.in (AC_CONFIG_FILES): Add new file geometry/Makefile.
+        
 2007-08-23  John W. Eaton  <jwe@octave.org>
 
 	* pkg/pkg.m: Avoid using installed_packages for both function and
diff -r 5c4702203cc4 -r 9fddcc586065 scripts/Makefile.in
--- a/scripts/Makefile.in	Thu Aug 23 19:43:08 2007 +0000
+++ b/scripts/Makefile.in	Fri Aug 24 08:27:29 2007 +0000
@@ -29,7 +29,7 @@
 	configure.in configure mkinstalldirs mkdoc mkpkgadd gethelp.cc \
 	skip-autoheader move-if-change) DOCSTRINGS
 
-SUBDIRS = audio control deprecated elfun finance general image io \
+SUBDIRS = audio control deprecated elfun finance general geometry image io \
 	linear-algebra miscellaneous optimization path pkg plot polynomial \
 	quaternion set signal sparse specfun special-matrix startup \
 	statistics strings testfun time
diff -r 5c4702203cc4 -r 9fddcc586065 scripts/configure.in
--- a/scripts/configure.in	Thu Aug 23 19:43:08 2007 +0000
+++ b/scripts/configure.in	Fri Aug 24 08:27:29 2007 +0000
@@ -32,8 +32,8 @@
 	  control/base/Makefile control/hinf/Makefile \
 	  control/obsolete/Makefile control/system/Makefile \
 	  control/util/Makefile deprecated/Makefile elfun/Makefile \
-	  finance/Makefile general/Makefile image/Makefile io/Makefile \
-	  linear-algebra/Makefile miscellaneous/Makefile \
+	  finance/Makefile general/Makefile geometry/Makefile image/Makefile \
+	  io/Makefile linear-algebra/Makefile miscellaneous/Makefile \
 	  optimization/Makefile path/Makefile pkg/Makefile plot/Makefile \
 	  polynomial/Makefile quaternion/Makefile set/Makefile \
 	  signal/Makefile sparse/Makefile specfun/Makefile \
diff -r 5c4702203cc4 -r 9fddcc586065 scripts/geometry/Makefile.in
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/geometry/Makefile.in	Fri Aug 24 08:27:29 2007 +0000
@@ -0,0 +1,69 @@
+#
+# Makefile for octave's scripts/geometry directory
+#
+# John W. Eaton
+# jwe@bevo.che.wisc.edu
+# University of Wisconsin-Madison
+# Department of Chemical Engineering
+
+TOPDIR = ../..
+
+script_sub_dir = geometry
+
+srcdir = @srcdir@
+top_srcdir = @top_srcdir@
+VPATH = @srcdir@
+
+include $(TOPDIR)/Makeconf
+
+INSTALL = @INSTALL@
+INSTALL_PROGRAM = @INSTALL_PROGRAM@
+INSTALL_DATA = @INSTALL_DATA@
+
+SOURCES = convhull.m delaunay3.m delaunayn.m delaunay.m dsearch.m dsearchn.m \
+	griddata.m griddata3.m griddatan.m trimesh.m triplot.m tsearchn.m \
+	voronoi.m voronoin.m
+
+DISTFILES = $(addprefix $(srcdir)/,Makefile.in $(SOURCES))
+
+FCN_FILES = $(addprefix $(srcdir)/, $(SOURCES))
+FCN_FILES_NO_DIR = $(notdir $(FCN_FILES))
+
+all: PKG_ADD
+.PHONY: all
+
+install install-strip:
+	$(do-script-install)
+.PHONY: install install-strip
+
+uninstall:
+	$(do-script-uninstall)
+.PHONY: uninstall
+
+clean:
+.PHONY: clean
+
+PKG_ADD: $(FCN_FILES)
+	@echo "making PKG_ADD"
+	@$(do-mkpkgadd)
+
+tags: $(SOURCES)
+	ctags $(SOURCES)
+
+TAGS: $(SOURCES)
+	etags $(SOURCES)
+
+mostlyclean: clean
+.PHONY: mostlyclean
+
+distclean: clean
+	rm -f Makefile PKG_ADD
+.PHONY: distclean
+
+maintainer-clean: distclean
+	rm -f tags TAGS
+.PHONY: maintainer-clean
+
+dist:
+	ln $(DISTFILES) ../../`cat ../../.fname`/scripts/geometry
+.PHONY: dist
diff -r 5c4702203cc4 -r 9fddcc586065 scripts/geometry/convhull.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/geometry/convhull.m	Fri Aug 24 08:27:29 2007 +0000
@@ -0,0 +1,87 @@
+## Copyright (C) 2000  Kai Habel
+##
+## 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.
+
+## -*- texinfo -*-
+## @deftypefn {Loadable Function} {@var{H} =} convhull (@var{x}, @var{y})
+## @deftypefnx {Loadable Function} {@var{H} =} convhull (@var{x}, @var{y}, @var{opt})
+## Returns the index vector to the points of the enclosing convex hull.  The
+## data points are defined by the x and y vectors.
+##
+## A third optional argument, which must be a string, contains extra options
+## passed to the underlying qhull command.  See the documentation for the 
+## Qhull library for details.
+##
+## @seealso{delaunay, convhulln}
+## @end deftypefn
+
+## Author:	Kai Habel <kai.habel@gmx.de>
+
+function H = convhull (x,y,opt)
+
+  if ((nargin != 2) && (nargin != 3))
+    print_usage ();
+  endif
+
+  if (isvector(x) && isvector(y) && (length(x) == length(y)))
+    if (nargin == 2)
+      i = convhulln([x(:), y(:)]);
+    elseif (ischar(opt) || iscell (opt))
+      i = convhulln([x(:), y(:)], opt);
+    else
+      error("third argument must be a string or cell array of strings");
+    endif
+  else
+    error("first two input arguments must be vectors of same size");
+  endif
+
+  n = rows(i);
+  i = i'(:);
+  H = zeros(n + 1,1);
+
+  H(1) = i(1);
+  next_i = i(2);
+  i(2) = 0;
+  for k = 2:n
+    next_idx = find (i == next_i);
+
+    if (rem (next_idx, 2) == 0)
+      H(k) = i(next_idx);
+      next_i = i(next_idx - 1);
+      i(next_idx - 1) = 0;
+    else
+      H(k) = i(next_idx);
+      next_i = i(next_idx + 1);
+      i(next_idx + 1) = 0;
+    endif
+  endfor
+
+  H(n + 1) = H(1);
+endfunction
+
+%!test
+%! x = -3:0.5:3;
+%! y = abs (sin (x));
+%! assert (convhull (x, y, {"s","Qci","Tcv","Pp"}), [1;7;13;12;11;10;4;3;2;1])
+
+%!demo
+%! x = -3:0.05:3;
+%! y = abs (sin (x));
+%! k = convhull (x, y);
+%! plot (x(k),y(k),'r-',x,y,'b+');
+%! axis ([-3.05, 3.05, -0.05, 1.05]);
diff -r 5c4702203cc4 -r 9fddcc586065 scripts/geometry/delaunay.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/geometry/delaunay.m	Fri Aug 24 08:27:29 2007 +0000
@@ -0,0 +1,89 @@
+## Copyright (C) 1999,2000  Kai Habel
+##
+## 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.
+
+## -*- texinfo -*-
+## @deftypefn {Loadable Function} {@var{tri}=} delaunay (@var{x}, @var{y})
+## @deftypefnx {Loadable Function} {@var{tri}=} delaunay (@var{x}, @var{y}, @var{opt})
+## The return matrix of size [n, 3] contains a set triangles which are
+## described by the indices to the data point x and y vector.
+## The triangulation satisfies the Delaunay circumcircle criterion.
+## No other data point is in the circumcircle of the defining triangle.
+##
+## A third optional argument, which must be a string, contains extra options
+## passed to the underlying qhull command.  See the documentation for the 
+## Qhull library for details.
+##
+## @example
+## @group
+## x = rand(1,10);
+## y = rand(size(x));
+## T = delaunay(x,y);
+## X = [ x(T(:,1)); x(T(:,2)); x(T(:,3)); x(T(:,1)) ];
+## Y = [ y(T(:,1)); y(T(:,2)); y(T(:,3)); y(T(:,1)) ];
+## axis([0,1,0,1]);
+## plot(X,Y,'b',x,y,'r*');
+## @end group
+## @end example
+## @seealso{voronoi, delaunay3, delaunayn}
+## @end deftypefn
+
+## Author:	Kai Habel <kai.habel@gmx.de>
+
+function ret = delaunay (x,y,opt)
+
+  if ((nargin != 2) && (nargin != 3))
+    print_usage ();
+  endif
+  
+  if (isvector(x) && isvector(y) && (length(x) == length(y)))
+    if (nargin == 2)
+      tri = delaunayn([x(:), y(:)]);
+    elseif ischar(opt)
+      tri = delaunayn([x(:), y(:)], opt);
+    else
+      error("third argument must be a string");
+    endif
+  else
+    error("first two input arguments must be vectors of same size");
+  endif
+
+  if nargout == 0
+    x = x(:).'; y = y(:).';
+    X = [ x(tri(:,1)); x(tri(:,2)); x(tri(:,3)); x(tri(:,1)) ];
+    Y = [ y(tri(:,1)); y(tri(:,2)); y(tri(:,3)); y(tri(:,1)) ];
+    plot(X, Y, 'b', x, y, 'r*');
+  else
+    ret = tri;
+  endif
+endfunction
+
+%!test
+%! x = [-1, 0, 1, 0, 0];
+%! y = [0, 1, 0, -1, 0];
+%! assert (sortrows (sort (delaunay (x, y), 2)), [1,2,5;1,4,5;2,3,5;3,4,5])
+
+%!demo
+%! rand ('state', 1);
+%! x = rand(1,10);
+%! y = rand(size(x));
+%! T = delaunay(x,y);
+%! X = [ x(T(:,1)); x(T(:,2)); x(T(:,3)); x(T(:,1)) ];
+%! Y = [ y(T(:,1)); y(T(:,2)); y(T(:,3)); y(T(:,1)) ];
+%! axis([0,1,0,1]);
+%! plot(X,Y,'b',x,y,'r*');
diff -r 5c4702203cc4 -r 9fddcc586065 scripts/geometry/delaunay3.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/geometry/delaunay3.m	Fri Aug 24 08:27:29 2007 +0000
@@ -0,0 +1,59 @@
+## Copyright (C) 1999,2000  Kai Habel
+##
+## 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.
+
+## -*- texinfo -*-
+## @deftypefn {Loadable Function} {@var{T} =} delaunay3 (@var{x}, @var{y}, @var{z})
+## @deftypefnx {Loadable Function} {@var{T} =} delaunay3 (@var{x}, @var{y}, @var{z}, @var{opt})
+## A matrix of size [n, 4] is returned. Each row contains a 
+## set of tetrahedron which are
+## described by the indices to the data point vectors (x,y,z).
+##
+## A fourth optional argument, which must be a string or cell array of strings,
+## contains extra options passed to the underlying qhull command.  See the 
+## documentation for the Qhull library for details.
+## @seealso{delaunay,delaunayn}
+## @end deftypefn
+
+## Author:	Kai Habel <kai.habel@gmx.de>
+
+function tetr = delaunay3 (x,y,z,opt)
+
+  if ((nargin != 3) && (nargin != 4))
+    print_usage ();
+  endif
+
+  if (isvector(x) && isvector(y) &&isvector(z) && ...
+      (length(x) == length(y)) && (length(x) == length(z)))
+    if (nargin == 3)
+      tetr = delaunayn([x(:),y(:),z(:)]);
+    elseif (ischar(opt) || iscell (opt))
+      tetr = delaunayn([x(:),y(:),z(:)], opt);
+    else
+      error("delaunay3: fourth argument must be a string or cell array of strings");
+    endif
+  else
+    error("delaunay3: first three input arguments must be vectors of same size");
+  endif
+
+endfunction
+
+%!test
+%! x = [-1, -1, 1, 0, -1]; y = [-1, 1, 1, 0, -1]; z = [0, 0, 0, 1, 1];
+%! assert (sortrows (sort (delaunay3 (x, y, z), 2)), [1,2,3,4;1,2,4,5])
+
diff -r 5c4702203cc4 -r 9fddcc586065 scripts/geometry/delaunayn.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/geometry/delaunayn.m	Fri Aug 24 08:27:29 2007 +0000
@@ -0,0 +1,76 @@
+## Copyright (C) 2007  David Bateman
+##
+## 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.
+
+## -*- texinfo -*-
+## @deftypefn {Loadable Function} {@var{T} =} delaunayn (@var{P})
+## @deftypefnx {Loadable Function} {@var{T} =} delaunayn (@var{P}, @var{opt})
+## Form the Delaunay triangulation for a set of points.
+## The Delaunay triangulation is a tessellation of the convex hull of the
+## points such that no n-sphere defined by the n-triangles contains
+## any other points from the set.\n
+## The input matrix @var{P} of size [n, dim] contains n points in a space
+## of dimension dim. The return matrix @var{T} has the size [m, dim+1]. It
+## contains for each row a set of indices to the points, which describes a
+## simplex of dimension dim.  For example, a 2d simplex is a triangle and 3d
+## simplex is a tetrahedron.
+## 
+## Extra options for the underlying Qhull command can be specified by the
+## second argument. This argument is a cell array of strings. The default
+## options depend on the dimension of the input: 
+## 
+## @itemize 
+## @item  2D and 3D: @var{opt} = @{'Qt','Qbb','Qc'@}  
+## @item  4D and higher: @var{opt} = @{'Qt','Qbb','Qc','Qz'@} 
+## @end itemize 
+## 
+## If @var{opt} is [], then the default arguments are used. If @var{opt}
+## is @{'@w{}'@}, then none of the default arguments are used by Qhull. 
+## See the Qhull documentation for the available options. 
+## 
+## All options can also be specified as single string, for example
+## 'Qt Qbb Qc Qz'. 
+## 
+## @end deftypefn
+
+function t = delaunayn (x, varargin)
+  if (nargin < 1)
+    print_usage()
+  endif
+
+  t = __delaunayn__ (x, varargin {:});
+
+  ## Try to remove the zero volume simplices. The volume of the i-th simplex is
+  ## given by abs(det(x(t(i,1:end-1),:)-x(t(i,2:end),:)))/prod(1:n) 
+  ## (reference http://en.wikipedia.org/wiki/Simplex). Any simplex with a 
+  ## relative volume less than some arbitrary criteria is rejected. The 
+  ## criteria we use is the volume of the simplex corresponding to an 
+  ## orthogonal simplex is equal edge length all equal to the edge length of 
+  ## the original simplex. If the relative volume is 1e3*eps then the simplex
+  ## is rejected. Note division of the two volumes means that the factor 
+  ## prod(1:n) is dropped.
+  idx = [];
+  [nt, n] = size(t);
+  for i = 1 : nt
+    X = x(t(i,1:end-1),:) - x(t(i,2:end),:);
+    if (abs (det (X)) /  sqrt (sum (X .^ 2, 2)) < 1e3 * eps)
+     idx = [idx, i];
+    endif
+  endfor
+  t(idx,:) = [];
+endfunction
diff -r 5c4702203cc4 -r 9fddcc586065 scripts/geometry/dsearch.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/geometry/dsearch.m	Fri Aug 24 08:27:29 2007 +0000
@@ -0,0 +1,41 @@
+## Copyright (C) 2007  David Bateman
+##
+## 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.
+
+## -*- texinfo -*-
+## @deftypefn {Loadable Function} {@var{idx} =} dsearch (@var{x}, @var{y}, @var{tri}, @var{xi}, @var{yi})
+## @deftypefnx {Loadable Function} {@var{idx} =} dsearch (@var{x}, @var{y}, @var{tri}, @var{xi}, @var{yi}, @var{s})
+## Returns the index @var{idx} or the closest point in @code{@var{x}, @var{y})}
+## to the elements @code{[@var{xi}(:), @var{yi}(:)]}. The variable @var{s} is
+## accepted but ignored for compatibility.
+## @seealso{dsearchn, tsearch}
+## @end deftypefn
+
+function idx = dsearch (x, y, t, xi, yi, s)
+  if (nargin < 5 || nargin > 6)
+    print_usage();
+  endif
+  idx = __dsearchn__ ([x(:),y(:)], [xi(:), yi(:)]);
+endfunction
+
+%!shared x, y, tri
+%! x = [-1;-1;1];
+%! y = [-1;1;-1];
+%! tri = [1,2,3]; 
+%!assert (dsearch(x,y,tri,1,1/3), 3);
+%!assert (dsearch(x,y,tri,1/3,1), 2);
diff -r 5c4702203cc4 -r 9fddcc586065 scripts/geometry/dsearchn.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/geometry/dsearchn.m	Fri Aug 24 08:27:29 2007 +0000
@@ -0,0 +1,58 @@
+## Copyright (C) 2007  David Bateman
+##
+## 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.
+
+## -*- texinfo -*-
+## @deftypefn {Loadable Function} {@var{idx} =} dsearchn (@var{x}, @var{tri}, @var{xi})
+## @deftypefnx {Loadable Function} {@var{idx} =} dsearchn (@var{x}, @var{tri}, @var{xi}, @var{outval})
+## @deftypefnx {Loadable Function} {@var{idx} =} dsearchn (@var{x}, @var{xi})
+## @deftypefnx {Loadable Function} {[@var{idx}, @var{d}] =} dsearchn (@dots{})
+## Returns the index @var{idx} or the closest point in @var{x} to the elements
+## @var{xi}. If @var{outval} is supplied, then the values of @var{xi} that are
+## not contained within one of the simplicies @var{tri} are set to 
+## @var{outval}. Generally, @var{tri} is returned from @code{delaunayn 
+## (@var{x})}.
+## @seealso{dsearch, tsearch}
+## @end deftypefn
+
+function [idx, d] = dsearchn (x, t, xi, outval)
+  if (nargin < 2 || nargin > 4)
+    print_usage ();
+  endif
+
+  if (nargin == 2)
+    [idx, d] = __dsearchn__ (x, t);
+  else
+    [idx, d] = __dsearchn__ (x, xi);
+    if (nargin == 4)
+      idx2 = isnan (tsearchn (x, t, xi));
+      idx(idx2) = outval;
+      d(idx2) = outval;
+    endif
+  endif
+endfunction
+
+%!shared x, tri
+%! x = [-1,-1;-1,1;1,-1]; 
+%! tri = [1,2,3]; 
+%!assert (dsearchn(x,tri,[1,1/3]), 3);
+%!assert (dsearchn(x,tri,[1,1/3],NaN), NaN);
+%!assert (dsearchn(x,tri,[1,1/3],NA), NA);
+%!assert (dsearchn(x,tri,[1/3,1]), 2);
+%!assert (dsearchn(x,tri,[1/3,1],NaN), NaN);
+%!assert (dsearchn(x,tri,[1/3,1],NA), NA);
diff -r 5c4702203cc4 -r 9fddcc586065 scripts/geometry/griddata.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/geometry/griddata.m	Fri Aug 24 08:27:29 2007 +0000
@@ -0,0 +1,155 @@
+## Copyright (C) 1999,2000  Kai Habel
+##
+## 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.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{zi} =} griddata (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi}, @var{method})
+## @deftypefnx {Function File} {[@var{xi}, @var{yi}, @var{zi}] =} griddata (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi}, @var{method})
+## 
+## Generate a regular mesh from irregular data using interpolation.
+## The function is defined by @code{@var{z} = f (@var{x}, @var{y})}.
+## The interpolation points are all @code{(@var{xi}, @var{yi})}.  If
+## @var{xi}, @var{yi} are vectors then they are made into a 2D mesh.
+##
+## The interpolation method can be 'nearest', 'cubic' or 'linear'.
+## If method is omitted it defaults to 'linear'.
+## @seealso{delaunay}
+## @end deftypefn
+
+## Author:	Kai Habel <kai.habel@gmx.de>
+## Adapted-by:  Alexander Barth <barth.alexander@gmail.com>
+##              xi and yi are not "meshgridded" if both are vectors 
+##              of the same size (for compatibility)
+
+function [rx, ry, rz] = griddata (x,y,z,xi,yi,method)
+	
+  if (nargin == 5)
+    method="linear";
+  endif
+  if (nargin < 5 || nargin > 7) 
+    print_usage();
+  endif
+
+  if (ischar(method))
+    method=tolower(method);
+  endif
+  if (!all( (size(x)==size(y)) & (size(x)==size(z))))
+    error("griddata: x,y,z must be vectors of same length");
+  endif
+  
+  ## meshgrid xi and yi if they are vectors unless they
+  ## are vectors of the same length 
+  if (isvector(xi) && isvector(yi) && numel(xi) ~= numel(yi)) 
+    [xi,yi]=meshgrid(xi,yi);
+  endif
+
+  if (any(size(xi) != size(yi)))
+    error("griddata: xi and yi must be vectors or matrices of same size");
+  endif
+
+  [nr,nc] = size(xi);
+  
+  ## triangulate data
+  tri = delaunay(x,y);
+  zi = nan(size(xi));
+  
+  if strcmp(method,"cubic")
+    error("griddata(...,'cubic') cubic interpolation not yet implemented\n")
+
+  elseif strcmp(method,'nearest')
+    ## search index of nearest point
+    idx = dsearch(x,y,tri,xi,yi);
+    valid = !isnan(idx);
+    zi(valid) = z(idx(valid));
+
+  elseif strcmp(method,'linear')
+    ## search for every point the enclosing triangle
+    tri_list= tsearch(x,y,tri,xi(:),yi(:));
+
+    ## only keep the points within triangles.
+    valid = !isnan(reshape(tri_list,size(xi)));
+    tri_list = tri_list(!isnan(tri_list));
+    nr_t = rows(tri_list);
+
+    ## assign x,y,z for each point of triangle
+    x1 = x(tri(tri_list,1));y1=y(tri(tri_list,1));z1=z(tri(tri_list,1));
+    x2 = x(tri(tri_list,2));y2=y(tri(tri_list,2));z2=z(tri(tri_list,2));
+    x3 = x(tri(tri_list,3));y3=y(tri(tri_list,3));z3=z(tri(tri_list,3));
+
+    ## calculate norm vector
+    N = cross([x2-x1, y2-y1, z2-z1],[x3-x1, y3-y1, z3-z1]);
+    N_norm = sqrt(sumsq(N,2));
+    N = N ./ N_norm(:,[1,1,1]);
+    
+    ## calculate D of plane equation
+    ## Ax+By+Cz+D=0;
+    D = -(N(:,1) .* x1 + N(:,2) .* y1 + N(:,3) .* z1);
+    
+    ## calculate zi by solving plane equation for xi,yi
+    zi(valid) = -(N(:,1).*xi(valid) + N(:,2).*yi(valid) + D ) ./ N(:,3);
+    
+  else
+    error("griddata: unknown interpolation method");
+  endif
+
+  if nargout == 3
+    rx = xi;
+    ry = yi;
+    rz = zi;
+  elseif nargout == 1
+    rx = zi;
+  elseif nargout == 0
+    mesh(xi, yi, zi);
+  endif
+endfunction
+
+%!test
+%! [xx,yy]=meshgrid(linspace(-1,1,32));
+%! x = xx(:);
+%! x = x + 10 * (2 * round(rand(size(x))) - 1) * eps;
+%! y = yy(:);
+%! y = y + 10 * (2 * round(rand(size(y))) - 1) * eps;
+%! z = sin(2*(x.^2+y.^2));
+%! zz = griddata(x,y,z,xx,yy,'linear');
+%! zz2 = sin(2*(xx.^2+yy.^2));
+%! zz2(isnan(zz)) = NaN;
+%! assert (zz, zz2, 100 * eps)
+
+%!demo
+%! x=2*rand(100,1)-1;
+%! y=2*rand(size(x))-1;
+%! z=sin(2*(x.^2+y.^2));
+%! [xx,yy]=meshgrid(linspace(-1,1,32));
+%! griddata(x,y,z,xx,yy);
+%! title('nonuniform grid sampled at 100 points');
+
+%!demo
+%! x=2*rand(1000,1)-1;
+%! y=2*rand(size(x))-1;
+%! z=sin(2*(x.^2+y.^2));
+%! [xx,yy]=meshgrid(linspace(-1,1,32));
+%! griddata(x,y,z,xx,yy);
+%! title('nonuniform grid sampled at 1000 points');
+
+%!demo
+%! x=2*rand(1000,1)-1;
+%! y=2*rand(size(x))-1;
+%! z=sin(2*(x.^2+y.^2));
+%! [xx,yy]=meshgrid(linspace(-1,1,32));
+%! griddata(x,y,z,xx,yy,'nearest');
+%! title('nonuniform grid sampled at 1000 points with nearest neighbor');
diff -r 5c4702203cc4 -r 9fddcc586065 scripts/geometry/griddata3.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/geometry/griddata3.m	Fri Aug 24 08:27:29 2007 +0000
@@ -0,0 +1,57 @@
+## Copyright (C) 2007  David Bateman  <dbateman@free.fr>
+##
+## 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.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{vi} =} griddata3 (@var{x}, @var{y}, @var{z}, @var{v} @var{xi}, @var{yi}, @var{zi}, @var{method}, @var{options})
+## 
+## Generate a regular mesh from irregular data using interpolation.
+## The function is defined by @code{@var{y} = f (@var{x},@var{y},@var{z})}.
+## The interpolation points are all @var{xi}.  
+##
+## The interpolation method can be 'nearest' or 'linear'. If method is 
+## omitted it defaults to 'linear'.
+## @seealso{griddata, delaunayn}
+## @end deftypefn
+
+function [yi] = griddata3 (x,y,z,v,xi,yi,zi,method,varargin)
+	
+  if (nargin < 7) 
+    print_usage();
+  endif
+
+  if (!all( (size(x)==size(y)) & (size(x)==size(z)) & ...
+	   (size(x)==size(v))))
+    error("griddata3: x,y,z,v must be vectors of same length");
+  endif
+
+  ## meshgrid xi, yi and zi if they are vectors unless they
+  ## are vectors of the same length 
+  if (isvector(xi) && isvector(yi) && isvector(zi) && ...
+      (numel(xi) != numel(yi) || numel(xi) != numel(zi))) 
+    [xi,yi,zi] = meshgrid(xi,yi,zi);
+  endif
+
+  if (any(size(xi) != size(yi)) || any(size(xi) != size(zi)))
+    error("griddata: xi, yi and zi must be vectors or matrices of same size");
+  endif
+
+  vi = gridata ([x(:),y(:),z(:)], v(:), [xi(:),yi(:),zi(:)], varargin{:});
+  vi = reshape (vi, size(xi));
+endfunction
+
diff -r 5c4702203cc4 -r 9fddcc586065 scripts/geometry/griddatan.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/geometry/griddatan.m	Fri Aug 24 08:27:29 2007 +0000
@@ -0,0 +1,106 @@
+## Copyright (C) 2007  David Bateman  <dbateman@free.fr>
+##
+## 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.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{yi} =} griddatan (@var{x}, @var{y}, @var{xi}, @var{method}, @var{options})
+## 
+## Generate a regular mesh from irregular data using interpolation.
+## The function is defined by @code{@var{y} = f (@var{x})}.
+## The interpolation points are all @var{xi}.  
+##
+## The interpolation method can be 'nearest' or 'linear'. If method is 
+## omitted it defaults to 'linear'.
+## @seealso{griddata, delaunayn}
+## @end deftypefn
+
+function [yi] = griddatan (x,y,xi,method,varargin)
+	
+  if (nargin == 3)
+    method="linear";
+  endif
+  if (nargin < 3) 
+    print_usage();
+  endif
+
+  if (ischar(method))
+    method=tolower(method);
+  endif
+
+  [m, n] = size(x);
+  [mi, ni] = size(xi);
+  
+  if (n != ni || size(y,1) != m || size(y,2) != 1)
+    error ("griddatan: dimensional mismatch");
+  endif
+
+  ## triangulate data
+  ## tri = delaunayn(x, varargin{:});
+  tri = delaunayn(x);
+
+  yi = nan(mi, 1);
+  
+  if strcmp(method,'nearest')
+    ## search index of nearest point
+    idx = dsearchn(x,tri,xi);
+    valid = !isnan(idx);
+    yi(valid) = y(idx(valid));
+
+  elseif strcmp(method,'linear')
+    ## search for every point the enclosing triangle
+    [tri_list, bary_list] = tsearchn(x,tri,xi);
+
+    ## only keep the points within triangles.
+    valid = !isnan(tri_list);
+    tri_list = tri_list(!isnan(tri_list));
+    bary_list = bary_list(!isnan(tri_list), :);
+    nr_t = rows(tri_list);
+
+    ## assign x,y for each point of simplex
+    xt =  reshape (x(tri(tri_list,:),:),[nr_t, n+1, n]);
+    yt = y(tri(tri_list,:));
+
+    ## Use barycentric coordinate of point to calculate yi
+    yi(valid) = sum (y(tri(tri_list,:)) .* bary_list, 2);
+
+  else
+    method
+    error("griddatan: unknown interpolation method");
+  endif
+
+endfunction
+
+%!test
+%! [xx,yy] = meshgrid(linspace(-1,1,32));
+%! xi = [xx(:), yy(:)];
+%! x = (2 * rand(100,2) - 1);
+%! x = [x;1,1;1,-1;-1,-1;-1,1];
+%! y = sin(2*(sum(x.^2,2)));
+%! zz = griddatan(x,y,xi,'linear');
+%! zz2 = griddata(x(:,1),x(:,2),y,xi(:,1),xi(:,2),'linear');
+%! assert (zz, zz2, 1e-10)
+
+%!test
+%! [xx,yy] = meshgrid(linspace(-1,1,32));
+%! xi = [xx(:), yy(:)];
+%! x = (2 * rand(100,2) - 1);
+%! x = [x;1,1;1,-1;-1,-1;-1,1];
+%! y = sin(2*(sum(x.^2,2)));
+%! zz = griddatan(x,y,xi,'nearest');
+%! zz2 = griddata(x(:,1),x(:,2),y,xi(:,1),xi(:,2),'nearest');
+%! assert (zz, zz2, 1e-10)
diff -r 5c4702203cc4 -r 9fddcc586065 scripts/geometry/trimesh.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/geometry/trimesh.m	Fri Aug 24 08:27:29 2007 +0000
@@ -0,0 +1,64 @@
+## Copyright (C) 2007  David Bateman
+##
+## 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.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} trimesh (@var{tri}, @var{x}, @var{y}, @var{z})
+## @deftypefnx {Function File} {@var{h} = } trimesh (@dots{})
+## Plot a triangular mesh in 3D. The variable @var{tri} is the triangular
+## meshing of the points @code{(@var{x}, @var{y}, @var{z})} which is returned 
+## from @code{delaunay3}. The output argument @var{h} is the graphic handle
+## to the plot.
+## @seealso{triplot, delaunay3}
+## @end deftypefn
+
+function h = trimesh (tri, x, y, z, varargin)
+
+  if (nargin < 3)
+    print_usage();
+  endif
+
+  if (nargin == 3)
+    triplot (tri, x, y);
+  elseif (ischar (z))
+    triplot (tri, x, y, z, varargin{:});
+  else
+    idx = tri(:, [1,2,3,1]).';
+    nt = size(tri, 1);
+    ## FIXME We should really use patch instead of plot3, but we don't
+    ## have a patch function, and probably won't in 3D that works with
+    ## gnuplot
+    if (nargout > 0)
+      h = plot3([x(idx); NaN*ones(1, nt)](:), ...
+		[y(idx); NaN*ones(1, nt)](:), ...
+		[z(idx); NaN*ones(1, nt)](:), varargin{:});
+    else
+      plot3([x(idx); NaN*ones(1, nt)](:), ...
+	    [y(idx); NaN*ones(1, nt)](:), ...
+	    [z(idx); NaN*ones(1, nt)](:), varargin{:});
+    endif
+  endif
+endfunction
+
+%!demo
+%! rand ('state', 10)
+%! x = 3 - 6 * rand (1, 50);
+%! y = 3 - 6 * rand (1, 50);
+%! z = peaks (x, y);
+%! tri = delaunay (x, y);
+%! trimesh (tri, x, y, z);
diff -r 5c4702203cc4 -r 9fddcc586065 scripts/geometry/triplot.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/geometry/triplot.m	Fri Aug 24 08:27:29 2007 +0000
@@ -0,0 +1,54 @@
+## Copyright (C) 2007  David Bateman
+##
+## 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.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} triplot (@var{tri}, @var{x}, @var{y})
+## @deftypefnx {Function File} {} triplot (@var{tri}, @var{x}, @var{y}, @var{linespec})
+## @deftypefnx {Function File} {@var{h} = } triplot (@dots{})
+## Plot a triangular mesh in 2D. The variable @var{tri} is the triangular
+## meshing of the points @code{(@var{x}, @var{y})} which is returned from
+## @code{delaunay}. If given, the @var{linespec} determines the properties
+## to use for the lines. The output argument @var{h} is the graphic handle
+## to the plot.
+## @seealso{plot, trimesh, delaunay}
+## @end deftypefn
+
+function h = triplot (tri, x, y, varargin)
+
+  if (nargin < 3)
+    print_usage();
+  endif
+
+  idx = tri (:, [1, 2, 3, 1]).';
+  nt = size(tri, 1);
+  if (nargout > 0)
+    h = plot([x(idx); NaN*ones(1, nt)](:), ...
+	     [y(idx); NaN*ones(1, nt)](:), varargin{:});
+  else
+    plot([x(idx); NaN*ones(1, nt)](:), ...
+	 [y(idx); NaN*ones(1, nt)](:), varargin{:});
+  endif
+endfunction
+
+%!demo
+%! rand ('state', 2)
+%! x = rand (20, 1);
+%! y = rand (20, 1);
+%! tri = delaunay (x, y);
+%! triplot (tri, x, y);
diff -r 5c4702203cc4 -r 9fddcc586065 scripts/geometry/tsearchn.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/geometry/tsearchn.m	Fri Aug 24 08:27:29 2007 +0000
@@ -0,0 +1,108 @@
+## Copyright (C) 2007  David Bateman
+##
+## 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.
+
+## -*- texinfo -*-
+## @deftypefn {Loadable Function} {[@var{idx}, @var{p}] =} tsearchn (@var{x}, @var{t}, @var{xi})
+## Searches for the enclosing Delaunay convex hull. For @code{@var{t} =
+## delaunayn (@var{x})}, finds the index in @var{t} containing the
+## points @var{xi}. For points outside the convex hull, @var{idx} is NaN.
+## If requested @code{tsearchn} also returns the barycentric coorinates @var{p}
+## of the enclosing triangles.
+## @seealso{delaunay, delaunayn}
+## @end deftypefn
+
+function [idx, p] = tsearchn (x, t, xi)
+  if (nargin != 3)
+    print_usage ();
+  endif
+
+  nt = size (t, 1);
+  [m, n] = size (x);
+  mi = size (xi, 1);
+  idx = nan (mi, 1);
+  p = nan (mi, n + 1);
+
+  ni = [1:mi].';
+  for i = 1 : nt
+    ## Only calculate the Barycentric coordinates for points that have not
+    ## already been found in a triangle.
+    b = cart2bary (x (t (i, :), :), xi(ni,:));
+
+    ## Our points xi are in the current triangle if
+    ## (all(b >= 0) && all (b <= 1)). However as we impose that 
+    ## sum(b,2) == 1 we only need to test all(b>=0). Note need to add
+    ## a small margin for rounding errors
+    intri = all (b >= -1e-12, 2);
+    idx(ni(intri)) = i;
+    p(ni(intri),:) = b(intri, :);
+    ni(intri) = [];
+  endfor
+endfunction
+
+function Beta = cart2bary (T, P)
+  ## Conversion of Cartesian to Barycentric coordinates.
+  ## Given a reference simplex in N dimensions represented by a
+  ## (N+1)-by-(N) matrix, and arbitrary point P in cartesion coordinates,
+  ## represented by a N-by-1 row vector can be written as
+  ##
+  ## P = Beta * T
+  ##
+  ## Where Beta is a N+1 vector of the barycentric coordinates. A criteria
+  ## on Beta is that
+  ##
+  ## sum (Beta) == 1
+  ##
+  ## and therefore we can write the above as
+  ##
+  ## P - T(end, :) = Beta(1:end-1) * (T(1:end-1,:) - ones(N,1) * T(end,:))
+  ##
+  ## and then we can solve for Beta as
+  ##
+  ## Beta(1:end-1) = (P - T(end,:)) / (T(1:end-1,:) - ones(N,1) * T(end,:))
+  ## Beta(end) = sum(Beta)
+  ##
+  ## Note below is generalize for multiple values of P, one per row.
+  [M, N] = size (P);
+  Beta = (P - ones (M,1) * T(end,:)) / (T(1:end-1,:) - ones(N,1) * T(end,:));
+  Beta (:,end+1) = 1 - sum(Beta, 2);
+endfunction
+
+%!shared x, tri
+%! x = [-1,-1;-1,1;1,-1];
+%! tri = [1, 2, 3];
+%!test
+%! [idx, p] = tsearchn (x,tri,[-1,-1]);
+%! assert (idx, 1)
+%! assert (p, [1,0,0], 1e-12)
+%!test
+%! [idx, p] = tsearchn (x,tri,[-1,1]);
+%! assert (idx, 1)
+%! assert (p, [0,1,0], 1e-12)
+%!test
+%! [idx, p] = tsearchn (x,tri,[1,-1]);
+%! assert (idx, 1)
+%! assert (p, [0,0,1], 1e-12)
+%!test
+%! [idx, p] = tsearchn (x,tri,[-1/3,-1/3]);
+%! assert (idx, 1)
+%! assert (p, [1/3,1/3,1/3], 1e-12)
+%!test
+%! [idx, p] = tsearchn (x,tri,[1,1]);
+%! assert (idx, NaN)
+%! assert (p, [NaN, NaN, NaN])
diff -r 5c4702203cc4 -r 9fddcc586065 scripts/geometry/voronoi.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/geometry/voronoi.m	Fri Aug 24 08:27:29 2007 +0000
@@ -0,0 +1,156 @@
+## Copyright (C) 2000  Kai Habel
+##
+## 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.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} voronoi (@var{x}, @var{y})
+## @deftypefnx {Function File} {} voronoi (@var{x}, @var{y}, "plotstyle")
+## @deftypefnx {Function File} {} voronoi (@var{x}, @var{y}, "plotstyle", @var{options})
+## @deftypefnx {Function File} {[@var{vx}, @var{vy}] =} voronoi (@dots{})
+## plots voronoi diagram of points @code{(@var{x}, @var{y})}.
+## The voronoi facets with points at infinity are not drawn.
+## [@var{vx}, @var{vy}] = voronoi(...) returns the vertices instead plotting the
+## diagram. plot (@var{vx}, @var{vy}) shows the voronoi diagram.
+##
+## A fourth optional argument, which must be a string, contains extra options
+## passed to the underlying qhull command.  See the documentation for the
+## Qhull library for details.
+##
+## @example
+## @group
+##   x = rand(10,1); y = rand(size(x));
+##   h = convhull(x,y);
+##   [vx,vy] = voronoi(x,y);
+##   plot(vx, vy, "-b", x, y, "o", x(h), y(h), "-g")
+##   legend("", "points", "hull");
+## @end group
+## @end example
+##
+## @seealso{voronoin, delaunay, convhull}
+## @end deftypefn
+
+## Author: Kai Habel <kai.habel@gmx.de>
+## First Release: 20/08/2000
+
+## 2002-01-04 Paul Kienzle <pkienzle@users.sf.net>
+## * limit the default graph to the input points rather than the whole diagram
+## * provide example
+## * use unique(x,"rows") rather than __unique_rows__
+
+## 2003-12-14 Rafael Laboissiere <rafael@laboissiere.net>
+## Added optional fourth argument to pass options to the underlying
+## qhull command
+
+function [vvx, vvy] = voronoi (varargin)
+
+  if (nargin < 1)
+    print_usage ();
+  endif
+
+  narg = 1;
+  if (isscalar (varargin{1}) && ishandle (varargin{1}))
+    handl = varargin{1};
+    narg++;
+    obj = get (handl);
+    if (! strcmp (obj.type, "axes"))
+      error ("voronoi: expecting first argument to be an axes object");
+    endif
+  else
+    if (nargout < 2)    
+      handl = gca();
+    endif
+  endif
+
+  if (nargin < 1 + narg || nargin > 3 + narg)
+    print_usage ();
+  endif
+
+  x = varargin{narg++};
+  y = varargin{narg++};
+  
+  opts = {};
+  if (narg <= nargin) 
+    if (iscell (varargin{narg}))
+      opts = varargin(narg++);
+    elseif (ismatrix (varargin{narg}))
+      ## Accept but ignore the triangulation
+      narg++;
+    endif
+  endif
+
+  linespec = {"b"};
+  if (narg <= nargin) 
+    if (ischar (varargin{narg}))
+      linespec = varargin(narg);
+    endif
+  endif
+
+  lx = length (x);
+  ly = length (y);
+
+  if (lx != ly)
+    error ("voronoi: arguments must be vectors of same length");
+  endif
+
+  [p, c, infi] = __voronoi__ ([x(:),y(:)], opts{:});
+
+  idx = find (!infi);
+
+  ll = length (idx);
+  k = 0;r = 1;
+
+  for i = 1:ll
+    k += length (c{idx(i)});
+  endfor
+
+  edges = zeros (2, k);
+
+  for i=1:ll
+    fac = c {idx(i)};
+    lf = length(fac);
+    fac = [fac, fac(1)];
+    fst = fac(1:length(fac) - 1);
+    sec = fac(2:length(fac));
+    edges(:,r:r+lf-1) = [fst; sec];
+    r += lf;
+  endfor
+
+  ## Identify the unique edges of the Voronoi diagram
+  edges = sortrows (sort (edges).').';
+  edges = edges (:,[(edges(1,1:end-1) != edges(1,2:end) | ...
+		     edges(2,1:end-1) != edges(2,2:end)), true]);
+
+  ## Get points of the diagram
+  vx = reshape(p(edges, 1), size(edges));
+  vy = reshape(p(edges, 2), size(edges));
+
+  if (nargout < 2)    
+    lim = [min(x(:)), max(x(:)), min(y(:)), max(y(:))];
+    h = plot (handl, vx, vy, linespec{:}, x, y, '+');
+    axis(lim+0.1*[[-1,1]*(lim(2)-lim(1)),[-1,1]*(lim(4)-lim(3))]);
+    if (nargout == 1)
+      vxx = h;
+    endif
+  elseif (nargout == 2)
+    vvx = vx;
+    vvy = vy;
+  else
+    error ("voronoi: only two or zero output arguments supported")
+  endif
+
+endfunction
diff -r 5c4702203cc4 -r 9fddcc586065 scripts/geometry/voronoin.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/geometry/voronoin.m	Fri Aug 24 08:27:29 2007 +0000
@@ -0,0 +1,60 @@
+## Copyright (C) 2000  Kai Habel
+##
+## 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.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{C}, @var{F}] =} voronoin (@var{pts})
+## @deftypefnx {Function File} {[@var{C}, @var{F}] =} voronoin (@var{pts}, @var{options})
+## computes n- dimensinal voronoi facets.  The input matrix @var{pts}
+## of size [n, dim] contains n points of dimension dim.
+## @var{C} contains the points of the voronoi facets. The list @var{F}
+## contains for each facet the indices of the voronoi points.
+##
+## A second optional argument, which must be a string, contains extra options
+## passed to the underlying qhull command.  See the documentation for the
+## Qhull library for details.
+## @seealso{voronoin, delaunay, convhull}
+## @end deftypefn
+
+## Author: Kai Habel <kai.habel@gmx.de>
+## First Release: 20/08/2000
+
+## 2003-12-14 Rafael Laboissiere <rafael@laboissiere.net>
+## Added optional second argument to pass options to the underlying
+## qhull command
+
+function [C, F] = voronoin (pts, opt)
+
+  if ((nargin != 1) && (nargin != 2))
+    print_usage ();
+  endif
+
+  [np, dims] = size (pts);
+  if (np > dims)
+    if (nargin == 1)
+      [C, F, infi] = __voronoi__ (pts);
+    elseif ischar(opt)
+      [C, F, infi] = __voronoi__ (pts, opt);
+    else
+      error("second argument must be a string");
+    endif
+
+  else
+    error ("voronoin: number of points must be greater than their dimension")
+  endif
+endfunction
diff -r 5c4702203cc4 -r 9fddcc586065 src/ChangeLog
--- a/src/ChangeLog	Thu Aug 23 19:43:08 2007 +0000
+++ b/src/ChangeLog	Fri Aug 24 08:27:29 2007 +0000
@@ -1,3 +1,36 @@
+2007-08-24  David Bateman  <dbateman@free.fr>
+
+        * interpreter/geometry.txi: Document new functions.
+        * interpreter/octave.texi: Update indexing of geometry items.
+
+2007-08-24  David Bateman  <dbateman@free.fr>
+
+        * ov-bool-sparse.cc (DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA): Class
+        is now logical.
+        * ov-re-sparse.cc, ov-cx-sparse.cc 
+        (DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA): Class is now double.
+        * ov-mapper.cc (octave_mapper::apply): Use is_sparse_type method
+        rather than comparing class name.
+        * ls-mat5.cc (save_mat5_element_length, save_mat5_binary_element):
+        ditto.
+        * pt-mat.cc (tree_matrix:rvalue): sparse matrices are now of class
+        "double" or "logical", create new single type concat clause for
+        them.
+        * mex.cc (get_class_id): No longer need to special case sparse
+        matrices.
+
+        * DLD-FUNCTIONS/__delaunayn__.cc, DLD-FUNCTIONS/convhulln.cc, 
+        DLD-FUNCTIONS/tsearch.cc, DLD-FUNCTIONS/__voronoi__.cc: New
+        functions ported from octave-forge.
+        * DLD-FUCTIONS/__dsearchn__.cc: New file.
+        * DLD-FUNCTIONS/__voronoi__.cc: Return point at infinity and
+        include it in the polygrons of the Voronoi diagram for
+        compatibility.
+        * Makefile.in: Add specific build targets for __delanayn__.cc,
+        convhulln.cc and __voronoi__.cc to link to Qhull.
+        (DLD_SRC): Add new functions.
+        (OCTAVE_LIBS): Add QHULL_LIBS
+
 2007-08-22  David Bateman  <dbateman@free.fr>
 
 	* variables.cc (Fmunlock): Call munlock and not mlock.
diff -r 5c4702203cc4 -r 9fddcc586065 src/DLD-FUNCTIONS/__delaunayn__.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/DLD-FUNCTIONS/__delaunayn__.cc	Fri Aug 24 08:27:29 2007 +0000
@@ -0,0 +1,228 @@
+/*
+
+Copyright (C) 2000  Kai Habel
+
+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.
+
+*/
+
+/*
+  16. July 2000 - Kai Habel: first release
+
+  25. September 2002 - Changes by Rafael Laboissiere <rafael@laboissiere.net>
+
+  * Added Qbb option to normalize the input and avoid crashes in Octave.
+  * delaunayn accepts now a second (optional) argument that must be a string
+  containing extra options to the qhull command.
+  * Fixed doc string.  The dimension of the result matrix is [m, dim+1], and
+  not [n, dim-1].
+
+  6. June 2006: Changes by Alexander Barth <abarth@marine.usf.edu>
+
+  * triangulate non-simplicial facets 
+  * allow options to be specified as cell array of strings
+  * change the default options (for compatibility with matlab)
+*/
+
+#include <iostream>
+#include <string>
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+#include "oct.h"
+#include "ov-cell.h"
+
+#ifdef HAVE_QHULL
+extern "C" {
+#include <qhull/qhull_a.h>
+}
+
+#ifdef NEED_QHULL_VERSION
+char qh_version[] = "__delaunayn__.oct 2007-08-21";
+#endif
+FILE *outfile = stdout;
+FILE *errfile = stderr;
+char flags[250];
+#endif
+
+DEFUN_DLD (__delaunayn__, args, ,
+	   "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {@var{T} =} __delaunayn__ (@var{P}[, @var{opt}])\n\
+Internal function for delaunayn.\n\
+@end deftypefn")
+
+{
+  octave_value_list retval;
+#ifdef HAVE_QHULL
+  retval(0) = 0.0;
+  std::string options = "";
+
+  int nargin = args.length();
+  if (nargin < 1 || nargin > 2)
+    {
+      print_usage ();
+      return retval;
+    }
+
+  Matrix p(args(0).matrix_value());
+  const octave_idx_type dim = p.columns();
+  const octave_idx_type n = p.rows();
+
+  // default options
+  if (dim <= 3)
+    options = "Qt Qbb Qc";
+  else
+    options = "Qt Qbb Qc Qx";
+
+
+  if (nargin == 2)
+    {
+    if (args(1).is_empty())
+      {
+	// keep default options
+      }
+    else if ( args(1).is_string ()) 
+      {
+	// option string is directly provided
+	options = args(1).string_value();
+      }
+    else if ( args(1).is_cell ()) 
+      {
+	options = "";
+
+	Cell c =  args(1).cell_value();
+	for (octave_idx_type i = 0; i < c.numel(); i++)
+	  {
+
+	    if (!c.elem(i).is_string()) 
+	      {
+		error ("__delaunayn__: all options must be strings");
+		return retval;
+	      }
+
+	    options = options + c.elem(i).string_value() + " ";
+	  }
+      }
+    else 
+      {
+	error ("__delaunayn__: second argument must be a string, cell of stringsor empty");
+	return retval;
+      }
+    } 
+
+  //octave_stdout << "options " << options << std::endl;
+
+  if (n > dim) 
+    {
+      p = p.transpose();
+      double *pt_array = p.fortran_vec();
+      boolT ismalloc = False;
+
+      sprintf(flags,"qhull d %s",options.c_str());
+
+      // If you want some debugging information replace the NULL
+      // pointer with outfile.
+
+      int exitcode =  qh_new_qhull (dim, n, pt_array, ismalloc,flags, 
+				    NULL, errfile);
+
+      if (exitcode)
+	{
+	  error("__delaunayn__: qhull failed.");
+	  return retval;
+	}
+
+      // triangulate non-simplicial facets
+      qh_triangulate(); 
+
+      facetT *facet;
+      vertexT *vertex, **vertexp;
+      octave_idx_type nf = 0, i = 0;
+
+      FORALLfacets
+	{
+	  if (!facet->upperdelaunay)
+	    nf++;
+
+	  // Double check
+	  if (!facet->simplicial) 
+	    {
+	      error("__delaunayn__: Qhull returned non-simplicial facets.\n",
+		    "Try delaunayn with different options.");
+	      break;
+	    }
+	}
+
+    if (!error_state) 
+      {
+	Matrix simpl(nf, dim+1);
+	FORALLfacets
+	  {
+	    if (!facet->upperdelaunay) 
+	      {
+		octave_idx_type j = 0;
+		FOREACHvertex_ (facet->vertices)
+		  {
+		    // if delaunayn crashes, enable this check
+#if 0
+		    if (j > dim)
+		      {
+			error("__delaunayn__: internal error. Qhull returned non-simplicial facets.");
+			return retval;
+		      }
+#endif
+
+		    simpl(i, j++) = 1 + qh_pointid(vertex->point);
+		  }
+		i++;
+	      }
+	  }
+	retval(0) = simpl;
+      }
+    
+
+    qh_freeqhull(!qh_ALL);
+    //free long memory
+
+    int curlong, totlong;
+    qh_memfreeshort (&curlong, &totlong);
+    //free short memory and memory allocator
+
+    if (curlong || totlong)
+      {
+	warning("__delaunay__: did not free %d bytes of long memory (%d pieces)",
+		totlong, curlong);
+      }
+    } 
+  else if (n == dim + 1) 
+    {
+      // one should check if nx points span a simplex
+      // I will look at this later.
+      RowVector vec(n);
+      for (octave_idx_type i = 0; i < n; i++) 
+	{
+	  vec(i) = i + 1.0;
+	}
+      retval(0) = vec;
+    }
+#else
+  error ("__delaunayn__: not available in this version of Octave");
+#endif
+  return retval;
+}
diff -r 5c4702203cc4 -r 9fddcc586065 src/DLD-FUNCTIONS/__dsearchn__.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/DLD-FUNCTIONS/__dsearchn__.cc	Fri Aug 24 08:27:29 2007 +0000
@@ -0,0 +1,108 @@
+/*
+
+Copyright (C) 2007 David Bateman <dbateman@free.fr>
+
+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>
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "oct.h"
+
+DEFUN_DLD (__dsearchn__, args, ,
+	"-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {[@var{idx}, @var{d}] =} dsearch (@var{x}, @var{xi})\n\
+Internal function for dsearchn.\n\
+@end deftypefn")
+{
+  int nargin = args.length();
+  octave_value_list retval;
+
+  if (nargin != 2)
+    {
+      print_usage ();
+      return retval;
+    }
+
+  Matrix x = args(0).matrix_value().transpose ();
+  Matrix xi = args(1).matrix_value().transpose ();
+
+  if (! error_state)
+    {
+      if (x.rows() != xi.rows() || x.columns() < 1)
+	error ("__dsearch__: dimensional mismatch.");
+      else
+	{
+	  octave_idx_type n = x.rows();
+	  octave_idx_type nx = x.columns();
+	  octave_idx_type nxi = xi.columns();
+
+	  ColumnVector idx (nxi);
+	  double *pidx = idx.fortran_vec ();
+	  ColumnVector dist (nxi);
+	  double *pdist = dist.fortran_vec ();
+
+#define DIST(dd, y, yi, m) \
+  dd = 0.; \
+  for (octave_idx_type k = 0; k < m; k++) \
+   { \
+     double yd = y[k] - yi[k]; \
+     dd += yd * yd; \
+   } \
+  dd = sqrt (dd);
+
+	  const double *pxi = xi.fortran_vec ();
+	  for (octave_idx_type i = 0; i < nxi; i++)
+	    {
+	      double d0;
+	      const double *px = x.fortran_vec ();
+	      DIST(d0, px, pxi, n);
+	      *pidx = 1.;
+	      for (octave_idx_type j = 1; j < nx; j++)
+		{
+		  px += n;
+		  double d;
+		  DIST (d, px, pxi, n);
+		  if (d < d0)
+		    {
+		      d0 = d;
+		      *pidx = static_cast<double>(j + 1);
+		    }
+		  OCTAVE_QUIT;
+ 		}
+
+	      *pdist++ = d0; 
+	      pidx++;
+	      pxi += n;
+	    }
+
+	  retval(1) = dist;
+	  retval(0) = idx;
+	}
+    }
+
+  return retval;
+}
diff -r 5c4702203cc4 -r 9fddcc586065 src/DLD-FUNCTIONS/__voronoi__.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/DLD-FUNCTIONS/__voronoi__.cc	Fri Aug 24 08:27:29 2007 +0000
@@ -0,0 +1,228 @@
+/*
+
+Copyright (C) 2000  Kai Habel
+
+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.
+
+*/
+
+/*
+20. Augiust 2000 - Kai Habel: first release
+*/
+
+/*
+2003-12-14 Rafael Laboissiere <rafael@laboissiere.net>
+Added optional second argument to pass options to the underlying
+qhull command
+*/
+
+#include <iostream>
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+#include "oct.h"
+#include "Cell.h"
+
+#ifdef HAVE_QHULL
+extern "C" {
+#include <qhull/qhull_a.h>
+}
+
+#ifdef NEED_QHULL_VERSION
+char qh_version[] = "__voronoi__.oct 2007-07-24";
+#endif
+
+FILE *outfile = stdout;
+FILE *errfile = stderr;
+char flags[250];
+const char *options;
+#endif
+
+DEFUN_DLD (__voronoi__, args, ,
+        "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {@var{tri} =} __voronoi__ (@var{point})\n\
+@deftypefnx {Loadable Function} {@var{tri} =} __voronoi__ (@var{point}, @var{options})\n\
+Internal function for voronoi.\n\
+@end deftypefn")
+{
+  octave_value_list retval;
+#ifdef HAVE_QHULL
+  retval(0) = 0.0;
+
+  int nargin = args.length();
+  if (nargin < 1 || nargin > 2) 
+    {
+      print_usage ();
+      return retval;
+    }
+
+  if (nargin == 2) 
+    {
+      if (!args (1).is_string ()) 
+	{
+	  error ("__voronoi__: second argument must be a string");
+	  return retval;
+	}
+      options = args (1).string_value().c_str();
+    }
+  else
+    options = "";
+
+  Matrix p(args(0).matrix_value());
+
+  const octave_idx_type dim = p.columns();
+  const octave_idx_type np = p.rows();
+  p = p.transpose();
+
+  double *pt_array = p.fortran_vec();
+  //double  pt_array[dim * np];
+  //for (int i = 0; i < np; i++) 
+  //  {
+  //    for (int j = 0; j < dim; j++)
+  //	  {
+  //	    pt_array[j+i*dim] = p(i,j);
+  //	  }
+  //  }
+
+  boolT ismalloc = false;
+
+  // hmm  lot's of options for qhull here
+  sprintf(flags,"qhull v Fv T0 %s",options);
+  if (!qh_new_qhull (dim, np, pt_array, ismalloc, flags, NULL, errfile)) 
+    {
+
+      // If you want some debugging information replace the NULL
+      // pointer with outfile.
+      facetT *facet;
+      vertexT *vertex;
+
+      octave_idx_type i = 0, n = 1, k = 0, m = 0, fidx = 0, j = 0, r = 0;
+      OCTAVE_LOCAL_BUFFER(octave_idx_type, ni, np);
+      
+      for (i = 0; i < np; i++) 
+	ni[i] = 0;
+      qh_setvoronoi_all();
+      bool infinity_seen = false;
+      facetT *neighbor,**neighborp;
+      coordT *voronoi_vertex;
+      FORALLfacets 
+	{
+	  facet->seen = false;
+	}
+      FORALLvertices 
+	{
+	  if (qh hull_dim == 3)
+	    qh_order_vertexneighbors(vertex);
+	  infinity_seen = false;
+	  FOREACHneighbor_(vertex)
+	    {
+	      if (!neighbor->upperdelaunay)
+		{
+		  if (!neighbor->seen)
+		    {
+		      n++;
+		      neighbor->seen=True;
+		    }
+		  ni[k]++;
+		}
+	      else if (!infinity_seen)
+		{
+		  infinity_seen = true;
+		  ni[k]++;
+		}
+	    }
+	  k++;
+	}
+
+      Matrix v(n, dim);
+      for (octave_idx_type d = 0; d < dim; d++)
+	v(0,d) = octave_Inf;
+
+      boolMatrix AtInf(np, 1);
+      for (i = 0; i < np; i++) 
+	AtInf(i) = false;
+      octave_value_list F;
+      k = 0;
+      i = 0;
+      FORALLfacets 
+	{
+	  facet->seen = false;
+	}
+      FORALLvertices
+	{
+	  if (qh hull_dim == 3)
+	    qh_order_vertexneighbors(vertex);
+	  infinity_seen = false;
+	  RowVector facet_list(ni[k++]);
+	  m = 0;
+	  FOREACHneighbor_(vertex)
+	    {
+	      if (neighbor->upperdelaunay)
+		{
+		  if (!infinity_seen)
+		    {
+		      infinity_seen = true;
+		      facet_list(m++) = 1;
+		      AtInf(j) = true;
+		    }
+		} 
+	      else 
+		{
+		  if (!neighbor->seen)
+		    {
+		      voronoi_vertex = neighbor->center;
+		      fidx = neighbor->id;
+		      i++;
+		      for (octave_idx_type d = 0; d < dim; d++)
+			{
+			  v(i,d) = *(voronoi_vertex++);
+			}
+		      neighbor->seen = true;
+		      neighbor->visitid = i;
+		    }
+		  facet_list(m++)=neighbor->visitid + 1;
+		}
+	    }
+	  F(r++) = facet_list;
+	  j++;
+	}
+
+      Cell C(r, 1);
+      for (i = 0; i < r; i++)
+	C.elem (i) = F(i);
+
+      retval(0) = v;
+      retval(1) = C;
+      retval(2) = AtInf;
+
+      //free long memory
+      qh_freeqhull(!qh_ALL);
+
+      //free short memory and memory allocator
+      int curlong, totlong;
+      qh_memfreeshort (&curlong, &totlong);
+
+      if (curlong || totlong) {
+	warning("__voronoi__: did not free %d bytes of long memory (%d pieces)", totlong, curlong);
+      }
+    }
+#else
+  error ("__voronoi__: not available in this version of Octave");
+#endif
+  return retval;
+}
diff -r 5c4702203cc4 -r 9fddcc586065 src/DLD-FUNCTIONS/convhulln.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/DLD-FUNCTIONS/convhulln.cc	Fri Aug 24 08:27:29 2007 +0000
@@ -0,0 +1,185 @@
+/*
+
+Copyright (C) 2000  Kai Habel
+
+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.
+
+*/
+
+/*
+29. July 2000 - Kai Habel: first release
+2002-04-22 Paul Kienzle
+* Use warning(...) function rather than writing to cerr
+2006-05-01 Tom Holroyd
+* add support for consistent winding in all dimensions; output is
+* guaranteed to be simplicial.
+*/
+
+#include <sstream>
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+#include "oct.h"
+#include "Cell.h"
+
+#ifdef HAVE_QHULL
+#if defined(HAVE__SNPRINTF) && !defined(HAVE_SNPRINTF)
+#define snprintf _snprintf
+#endif
+
+extern "C" {
+#include <qhull/qhull_a.h>
+}
+
+#ifdef NEED_QHULL_VERSION
+char qh_version[] = "convhulln.oct 2007-07-24";
+#endif
+char flags[250];
+#endif
+
+DEFUN_DLD (convhulln, args, ,
+	"-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {@var{H} =} convhulln (@var{p})\n\
+@deftypefnx {Loadable Function} {@var{H} =} convhulln (@var{p}, @var{opt})\n\
+Returns an index vector to the points of the enclosing convex hull.\n\
+The input matrix of size [n, dim] contains n points of dimension dim.\n\n\
+If a second optional argument is given, it must be a string or cell array\n\
+of strings containing options for the underlying qhull command.  (See\n\
+the Qhull documentation for the available options.)  The default options\n\
+are \"s Qci Tcv\".\n\n\
+@seealso{convhull, delaunayn}\n\
+@end deftypefn")
+{
+  octave_value_list retval;
+#ifdef HAVE_QHULL
+  std::string options;
+
+  int nargin = args.length();
+  if (nargin < 1 || nargin > 2) 
+    {
+      print_usage ();
+      return retval;
+    }
+
+  if (nargin == 2) 
+    {
+      if (args (1).is_string () ) 
+	options = args (1).string_value();
+      else if (args (1).is_cell ())
+	{
+	  Cell c = args (1). cell_value ();
+	  options = "";
+	  for (octave_idx_type i = 0; i < c.numel(); i++)
+	    {
+	      if (!c.elem(i).is_string()) 
+		{
+		  error ("convhulln: second argument must be a string or cell array of strings");
+		  return retval;
+	      }
+
+	      options = options + c.elem(i).string_value() + " ";
+	    }
+	}
+      else
+	{
+	  error ("convhulln: second argument must be a string or cell array of strings");
+	  return retval;
+	}
+    }
+  else
+    // turn on some consistency checks
+    options = "s Qci Tcv";
+
+  Matrix p(args(0).matrix_value());
+
+  const octave_idx_type dim = p.columns();
+  const octave_idx_type n = p.rows();
+  p = p.transpose();
+
+  double *pt_array = p.fortran_vec();
+
+  boolT ismalloc = False;
+
+  // hmm, lots of options for qhull here
+  // QJ guarantees that the output will be triangles
+  snprintf(flags, sizeof(flags), "qhull QJ %s", options.c_str());
+
+  if (!qh_new_qhull (dim, n, pt_array, ismalloc, flags, NULL, stderr)) 
+    {
+      // If you want some debugging information replace the NULL
+      // pointer with stdout
+
+      vertexT *vertex,**vertexp;
+      facetT *facet;
+      setT *vertices;
+      unsigned int nf = qh num_facets;
+
+      Matrix idx(nf, dim);
+
+      octave_idx_type j, i = 0;
+      FORALLfacets 
+	{
+	  j = 0;
+	  if (!facet->simplicial)
+	    // should never happen with QJ
+	    error("convhulln: non-simplicial facet");
+
+	  if (dim == 3) 
+	    {
+	      vertices = qh_facet3vertex (facet);
+	      FOREACHvertex_(vertices)
+		idx(i, j++) = 1 + qh_pointid(vertex->point);
+	      qh_settempfree(&vertices);
+	    } 
+	  else 
+	    {
+	      if (facet->toporient ^ qh_ORIENTclock) 
+		{
+		  FOREACHvertex_(facet->vertices)
+		    idx(i, j++) = 1 + qh_pointid(vertex->point);
+		} 
+	      else 
+		{
+		  FOREACHvertexreverse12_(facet->vertices)
+		    idx(i, j++) = 1 + qh_pointid(vertex->point);
+		}
+	    }
+	  if (j < dim)
+	    // likewise but less fatal
+	    warning("facet %d only has %d vertices",i,j);
+	  i++;
+	}
+      retval(0) = octave_value(idx);
+    }
+  qh_freeqhull (!qh_ALL);
+  //free long memory
+  int curlong, totlong;
+  qh_memfreeshort (&curlong, &totlong);
+  //free short memory and memory allocator
+
+  if (curlong || totlong) 
+    {
+      warning("convhulln: did not free %d bytes of long memory (%d pieces)",
+	      totlong, curlong);
+    }
+#else
+  error ("convhulln: not available in this version of Octave");
+#endif
+  return retval;
+}
diff -r 5c4702203cc4 -r 9fddcc586065 src/DLD-FUNCTIONS/tsearch.cc
--- /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)
+
+*/
diff -r 5c4702203cc4 -r 9fddcc586065 src/Makefile.in
--- a/src/Makefile.in	Thu Aug 23 19:43:08 2007 +0000
+++ b/src/Makefile.in	Fri Aug 24 08:27:29 2007 +0000
@@ -48,17 +48,18 @@
 	LSODE-opts.cc NLEqn-opts.cc Quad-opts.cc
 
 DLD_XSRC := balance.cc besselj.cc betainc.cc cellfun.cc chol.cc \
-	ccolamd.cc colamd.cc colloc.cc conv2.cc daspk.cc dasrt.cc \
-	dassl.cc det.cc dispatch.cc eig.cc expm.cc fft.cc fft2.cc \
-	fftn.cc fftw.cc filter.cc find.cc fsolve.cc \
+	ccolamd.cc colamd.cc colloc.cc conv2.cc convhulln.cc daspk.cc \
+	dasrt.cc dassl.cc det.cc dispatch.cc eig.cc expm.cc \
+	fft.cc fft2.cc fftn.cc fftw.cc filter.cc find.cc fsolve.cc \
 	gammainc.cc gcd.cc getgrent.cc getpwent.cc getrusage.cc \
 	givens.cc hess.cc inv.cc kron.cc lpsolve.cc lsode.cc \
 	lu.cc luinc.cc matrix_type.cc md5sum.cc minmax.cc pinv.cc qr.cc \
 	quad.cc qz.cc rand.cc regexp.cc schur.cc sort.cc sparse.cc \
 	spchol.cc spdet.cc spfind.cc spkron.cc splu.cc spparms.cc spqr.cc \
-	sqrtm.cc svd.cc syl.cc symrcm.cc time.cc urlwrite.cc __contourc__.cc \
-	__gnuplot_raw__.l __glpk__.cc __lin_interpn__.cc __pchip_deriv__.cc \
-	__qp__.cc
+	sqrtm.cc svd.cc syl.cc symrcm.cc time.cc tsearch.cc urlwrite.cc \
+	__contourc__.cc __delaunayn__.cc __dsearchn__.cc __gnuplot_raw__.l \
+	__glpk__.cc __lin_interpn__.cc __pchip_deriv__.cc __qp__.cc \
+	__voronoi__.cc
 
 DLD_SRC := $(addprefix DLD-FUNCTIONS/, $(DLD_XSRC))
 
@@ -242,7 +243,7 @@
     $(SPECIAL_MATH_LIB) $(LIBCRUFT) \
     $(LIBPLPLOT) $(LIBGLOB) $(LIBDLFCN)
 else
-  OCTAVE_LIBS = $(LIBOCTINTERP) $(LIBOCTAVE) \
+  OCTAVE_LIBS = $(LIBOCTINTERP) $(LIBOCTAVE) $(QHULL_LIBS) \
     $(GLPK_LIBS) $(REGEX_LIBS) $(SPECIAL_MATH_LIB) $(LIBCRUFT) \
     $(LIBPLPLOT) $(LIBGLOB) $(LIBDLFCN)
 endif
@@ -601,6 +602,12 @@
 
 ifeq ($(ENABLE_DYNAMIC_LINKING), true)
   ifdef CXXPICFLAG
+    convhulln.oct : pic/convhulln.o octave$(EXEEXT)
+	  $(DL_LD) $(DL_LDFLAGS) -o $@ $< $(OCT_LINK_DEPS) $(QHULL_LIBS)
+    __delaunayn__.oct : pic/__delaunayn__.o octave$(EXEEXT)
+	  $(DL_LD) $(DL_LDFLAGS) -o $@ $< $(OCT_LINK_DEPS) $(QHULL_LIBS)
+    __voronoi__.oct : pic/__voronoi__.o octave$(EXEEXT)
+	  $(DL_LD) $(DL_LDFLAGS) -o $@ $< $(OCT_LINK_DEPS) $(QHULL_LIBS)
     regexp.oct : pic/regexp.o octave$(EXEEXT)
 	  $(DL_LD) $(DL_LDFLAGS) -o $@ $< $(OCT_LINK_DEPS) $(REGEX_LIBS)
     urlwrite.oct : pic/urlwrite.o octave$(EXEEXT)
@@ -608,6 +615,12 @@
     __glpk__.oct : pic/__glpk__.o octave$(EXEEXT)
 	  $(DL_LD) $(DL_LDFLAGS) -o $@ $< $(OCT_LINK_DEPS) $(GLPK_LIBS)
   else
+    convhulln.oct : convhulln.o octave$(EXEEXT)
+	  $(DL_LD) $(DL_LDFLAGS) -o $@ $< $(OCT_LINK_DEPS) $(QHULL_LIBS)
+    __delaunayn__.oct : __delaunayn__.o octave$(EXEEXT)
+	  $(DL_LD) $(DL_LDFLAGS) -o $@ $< $(OCT_LINK_DEPS) $(QHULL_LIBS)
+    __voronoi__.oct : __voronoi__.o octave$(EXEEXT)
+	  $(DL_LD) $(DL_LDFLAGS) -o $@ $< $(OCT_LINK_DEPS) $(QHULL_LIBS)
     regexp.oct : regexp.o octave$(EXEEXT)
 	  $(DL_LD) $(DL_LDFLAGS) -o $@ $< $(OCT_LINK_DEPS) $(REGEX_LIBS)
     urlwrite.oct : urlwrite.o octave$(EXEEXT)
diff -r 5c4702203cc4 -r 9fddcc586065 src/ls-mat5.cc
--- a/src/ls-mat5.cc	Thu Aug 23 19:43:08 2007 +0000
+++ b/src/ls-mat5.cc	Fri Aug 24 08:27:29 2007 +0000
@@ -1725,7 +1725,7 @@
       if (chm.nelem () > 2)
 	ret += PAD (2 * chm.nelem ());
     }
-  else if (cname == "sparse")
+  else if (tc.is_sparse_type ())
     {
       if (tc.is_complex_type ())
 	{
@@ -1916,7 +1916,7 @@
     flags |= MAT_FILE_UINT32_CLASS;
   else if (cname == "uint64")
     flags |= MAT_FILE_UINT64_CLASS;
-  else if (cname == "sparse")
+  else if (tc.is_sparse_type ())
     {
       flags |= MAT_FILE_SPARSE_CLASS;
       if (tc.is_complex_type ())
@@ -2013,7 +2013,7 @@
 	  os.write (padbuf, paddedlength - len);
 	}
     }
-  else if (cname == "sparse")
+  else if (tc.is_sparse_type ())
     {
       if (tc.is_complex_type ())
 	{
diff -r 5c4702203cc4 -r 9fddcc586065 src/mex.cc
--- a/src/mex.cc	Thu Aug 23 19:43:08 2007 +0000
+++ b/src/mex.cc	Fri Aug 24 08:27:29 2007 +0000
@@ -390,13 +390,6 @@
       id = mxCHAR_CLASS;
     else if (cn == "double")
       id = mxDOUBLE_CLASS;
-    else if (cn == "sparse")
-      {
-	if (val.is_bool_type ())
-	  id = mxLOGICAL_CLASS;
-	else
-	  id = mxDOUBLE_CLASS;
-      }
     else if (cn == "single")
       id = mxSINGLE_CLASS;
     else if (cn == "int8")
diff -r 5c4702203cc4 -r 9fddcc586065 src/ov-bool-sparse.cc
--- a/src/ov-bool-sparse.cc	Thu Aug 23 19:43:08 2007 +0000
+++ b/src/ov-bool-sparse.cc	Fri Aug 24 08:27:29 2007 +0000
@@ -47,7 +47,7 @@
 
 DEFINE_OCTAVE_ALLOCATOR (octave_sparse_bool_matrix);
 
-DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA (octave_sparse_bool_matrix, "sparse bool matrix", "sparse");
+DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA (octave_sparse_bool_matrix, "sparse bool matrix", "logical");
 
 static octave_base_value *
 default_numeric_conversion_function (const octave_base_value& a)
diff -r 5c4702203cc4 -r 9fddcc586065 src/ov-cx-sparse.cc
--- a/src/ov-cx-sparse.cc	Thu Aug 23 19:43:08 2007 +0000
+++ b/src/ov-cx-sparse.cc	Fri Aug 24 08:27:29 2007 +0000
@@ -46,7 +46,7 @@
 
 DEFINE_OCTAVE_ALLOCATOR (octave_sparse_complex_matrix);
 
-DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA (octave_sparse_complex_matrix, "sparse complex matrix", "sparse");
+DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA (octave_sparse_complex_matrix, "sparse complex matrix", "double");
 
 octave_base_value *
 octave_sparse_complex_matrix::try_narrowing_conversion (void)
diff -r 5c4702203cc4 -r 9fddcc586065 src/ov-mapper.cc
--- a/src/ov-mapper.cc	Thu Aug 23 19:43:08 2007 +0000
+++ b/src/ov-mapper.cc	Fri Aug 24 08:27:29 2007 +0000
@@ -241,7 +241,7 @@
 	    error ("%s: unable to handle real arguments",
 		   name().c_str ());
 	}
-      else if (arg.class_name () == "sparse")
+      else if (arg.is_sparse_type ())
 	{
 	  const SparseMatrix m = arg.sparse_matrix_value ();
 
@@ -309,7 +309,7 @@
 	    error ("%s: unable to handle complex arguments",
 		   name().c_str ());
 	}
-      else if (arg.class_name () == "sparse")
+      else if (arg.is_sparse_type ())
 	{
 	  SparseComplexMatrix cm = arg.sparse_complex_matrix_value ();
 
diff -r 5c4702203cc4 -r 9fddcc586065 src/ov-re-sparse.cc
--- a/src/ov-re-sparse.cc	Thu Aug 23 19:43:08 2007 +0000
+++ b/src/ov-re-sparse.cc	Fri Aug 24 08:27:29 2007 +0000
@@ -46,7 +46,7 @@
 
 DEFINE_OCTAVE_ALLOCATOR (octave_sparse_matrix);
 
-DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA (octave_sparse_matrix, "sparse matrix", "sparse");
+DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA (octave_sparse_matrix, "sparse matrix", "double");
 
 idx_vector
 octave_sparse_matrix::index_vector (void) const
diff -r 5c4702203cc4 -r 9fddcc586065 src/pt-mat.cc
--- a/src/pt-mat.cc	Thu Aug 23 19:43:08 2007 +0000
+++ b/src/pt-mat.cc	Fri Aug 24 08:27:29 2007 +0000
@@ -791,10 +791,20 @@
 
       if (result_type == "double")
 	{
-	  if (all_real_p)
-	    DO_SINGLE_TYPE_CONCAT (NDArray, array_value);
+	  if (any_sparse_p)
+	    {	    
+	      if (all_real_p)
+		DO_SINGLE_TYPE_CONCAT (SparseMatrix, sparse_matrix_value);
+	      else
+		DO_SINGLE_TYPE_CONCAT (SparseComplexMatrix, sparse_complex_matrix_value);
+	    }
 	  else
-	    DO_SINGLE_TYPE_CONCAT (ComplexNDArray, complex_array_value);
+	    {
+	      if (all_real_p)
+		DO_SINGLE_TYPE_CONCAT (NDArray, array_value);
+	      else
+		DO_SINGLE_TYPE_CONCAT (ComplexNDArray, complex_array_value);
+	    }
 	}
 #if 0
       else if (result_type == "single")
@@ -812,7 +822,12 @@
 	  retval = octave_value (result, true, type);
 	}
       else if (result_type == "logical")
-	DO_SINGLE_TYPE_CONCAT (boolNDArray, bool_array_value);
+	{
+	  if (any_sparse_p)
+	    DO_SINGLE_TYPE_CONCAT (SparseBoolMatrix, sparse_bool_matrix_value);
+	  else
+	    DO_SINGLE_TYPE_CONCAT (boolNDArray, bool_array_value);
+	}
       else if (result_type == "int8")
 	DO_SINGLE_TYPE_CONCAT (int8NDArray, int8_array_value);
       else if (result_type == "int16")