Mercurial > octave-nkf
changeset 6823:9fddcc586065
[project @ 2007-08-24 08:27:27 by dbateman]
line wrap: on
line diff
--- 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.
--- 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@
--- 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 +]) +
--- 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
--- 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)
--- 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::
--- 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)
--- 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);
--- 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.
--- 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 (); }
--- 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) {
--- 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);
--- 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 {
--- 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);
--- 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
--- 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
--- 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 \
--- /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
--- /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]);
--- /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*');
--- /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]) +
--- /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
--- /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);
--- /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);
--- /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');
--- /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 +
--- /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)
--- /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);
--- /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);
--- /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])
--- /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
--- /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
--- 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.
--- /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; +}
--- /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; +}
--- /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; +}
--- /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; +}
--- /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) + +*/
--- 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)
--- 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 ()) {
--- 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")
--- 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)
--- 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)
--- 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 ();
--- 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
--- 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")