[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Latest port of Qhull/goemetry functions from octave-forge
From: |
David Bateman |
Subject: |
Re: Latest port of Qhull/goemetry functions from octave-forge |
Date: |
Sat, 25 Aug 2007 00:37:34 +0200 |
User-agent: |
Thunderbird 1.5.0.7 (X11/20060921) |
John W. Eaton wrote:
> On 24-Aug-2007, David Bateman wrote:
>
> | John W. Eaton wrote:
> | > I agree. Please check in these changes.
> | >
> | >
> | Ok, though I went a bit further than the previous patch in that I
> | started documenting and adding examples to the manual, modified
> | __voronoi__.cc to return the point at infinity for compatibility, and
> | made voronoi.m only return the unique edges of the Voronoi diagram. I
> | still have to make voronoi.m return the lines to infinity for
> | compatibility, though I now see how it could be done. I'll also finish
> | the documentation. This might take a while though :-)
>
> Thanks. I also made some additional style changes to the .m files.
>
> jwe
>
Here are the rest of the documentation that I wanted to add. The only
thing still missing is the edges of the Voronoi diagram with one point
at infinity.
D.
*** ./doc/interpreter/geometry.txi.orig10 2007-08-24 22:16:38.703857718
+0200
--- ./doc/interpreter/geometry.txi 2007-08-25 00:36:55.652179055 +0200
***************
*** 6,22 ****
@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
--- 6,39 ----
@node Geometry
@chapter Geometry
+ Much of geometry code in Octave is based on the QHull @footnote{Barber,
+ C.B., Dobkin, D.P., and Huhdanpaa, H.T., "The Quickhull algorithm for
+ convex hulls," ACM Trans. on Mathematical Software, 22(4):469-483, Dec
+ 1996, @url{http://www.qhull.org}}. Some of the documentation for Qhull,
+ particularly for the options that can be passed to @code{delaunay},
+ @code{voronoi} and @code{convhull}, etc, is relevant to Octave users.
+
@menu
* Delaunay Triangulation::
* Voronoi Diagrams::
* Convex Hull::
* Interpolation on Scattered Data::
@end menu
@node Delaunay Triangulation
@section Delaunay Triangulation
+ The Delaunay triangulation is constructed from a set of
+ circum-circles. These circum-circles are chosen so that there are at
+ least three of the points in the set to triangulation on the
+ circumference of the circum-circle. None of the points in the set of
+ points falls within any of the circum-circles.
+
+ In general there are only three points on the circumference of any
+ circum-circle. However, in the some cases, and in particular for the
+ case of a regular grid, 4 or more points can be on a single
+ circum-circle. In this case the Delaunay triangulation is not unique.
+
@DOCSTRING(delaunay)
The 3- and N-dimensional extension of the Delaunay triangulation are
***************
*** 30,39 ****
--- 47,118 ----
@DOCSTRING(delaunayn)
+ An example of a Delaunay triangulation of a set of points is
+
+ @example
+ @group
+ rand ("state", 2);
+ x = rand (10, 1);
+ y = rand (10, 1);
+ 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
+
+ @ifnotinfo
+ @noindent
+ The result of which can be seen in @ref{fig:delaunay}.
+
+ @float Figure,fig:delaunay
+ @image{delaunay,8cm}
+ @caption{Delaunay triangulation of a random set of points}
+ @end float
+ @end ifnotinfo
+
@menu
+ * Plotting the Triangulation::
* Identifying points in Triangulation::
@end menu
+ @node Plotting the Triangulation
+ @subsection Plotting the Triangulation
+
+ Octave has the functions @code{triplot} and @code{trimesh} to plot the
+ Delaunay triangulation of a 2-dimensional set of points.
+
+ @DOCSTRING(triplot)
+
+ @DOCSTRING(trimesh)
+
+ The difference between @code{triplot} and @code{trimesh} is that the
+ former only plots the 2-dimensional triangulation itself, whereas the
+ second plots the value of some function @code{f (@var{x}, @var{y})}.
+ An example of the use of the @code{triplot} function is
+
+ @example
+ @group
+ rand ("state", 2)
+ x = rand (20, 1);
+ y = rand (20, 1);
+ tri = delaunay (x, y);
+ triplot (tri, x, y);
+ @end group
+ @end example
+
+ that plot the Delaunay triangulation of a set of random points in
+ 2-dimensions.
+ @ifnotinfo
+ The output of the above can be seen in @ref{fig:triplot}.
+
+ @float Figure,fig:triplot
+ @image{triplot,8cm}
+ @caption{Delaunay triangulation of a random set of points}
+ @end float
+ @end ifnotinfo
+
@node Identifying points in Triangulation
@subsection Identifying points in Triangulation
***************
*** 143,149 ****
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)
--- 222,228 ----
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 vertexes of the
N-simplex within which the desired point is found.
@DOCSTRING(dsearch)
***************
*** 183,215 ****
such that all points in @address@hidden(@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)
--- 262,399 ----
such that all points in @address@hidden(@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, in that the vertexes of the
! Voronoi tessellation are the center's of the circum-circles of the
! simplicies of the Delaunay tessellation.
@DOCSTRING(voronoi)
@DOCSTRING(voronoin)
+ An example of the use of @code{voronoi} is
+
+ @example
+ @group
+ rand("state",9);
+ x = rand(10,1);
+ y = rand(10,1);
+ tri = delaunay (x, y);
+ [vx, vy] = voronoi (x, y, tri);
+ triplot (tri, x, y, "b");
+ hold on;
+ plot (vx, vy, "r");
+ @end group
+ @end example
+
+ @ifnotinfo
+ @noindent
+ The result of which can be seen in @ref{fig:voronoi}. Note that the
+ circum-circle of one of the triangles has been added to this figure, to
+ make the relationship between the Delaunay tessellation and the Voronoi
+ diagram clearer.
+
+ @float Figure,fig:voronoi
+ @image{voronoi,8cm}
+ @caption{Delaunay triangulation and Voronoi diagram of a random set of points}
+ @end float
+ @end ifnotinfo
+
+ Additional information about the size of the facets of a Voronoi diagram
+ can be had with the @code{polyarea} function.
+
+ @DOCSTRING(polyarea)
+
+ An example of the use of @code{polyarea} might be
+
+ @example
+ @group
+ rand ("state", 2);
+ x = rand (10, 1);
+ y = rand (10, 1);
+ [c, f] = voronoin ([x, y]);
+ af = zeros (size(f));
+ for i = 1 : length (f)
+ af(i) = polyarea (c (f @{i, :@}, 1), c (f @{i, :@}, 2));
+ endfor
+ @end group
+ @end example
+
+ Facets of the Voronoi diagram with a vertex at infinity have infinity area.
+
@node Convex Hull
@section Convex Hull
+ The convex hull of a set of points, is the minimum convex envelope
+ containing all of the points. Octave has the functions @code{convhull}
+ and @code{convhulln} to calculate the convec hull of 2-dimensional and
+ N-dimensional sets of points.
+
@DOCSTRING(convhull)
@DOCSTRING(convhulln)
! An example of the use of @code{convhull} is
! @example
! @group
! 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]);
! @end group
! @end example
! @ifnotinfo
! @noindent
! The output of the above can be seen in @ref{fig:convhull}.
! @float Figure,fig:convhull
! @image{convhull,8cm}
! @caption{The convex hull of a simple set of points}
! @end float
! @end ifnotinfo
@node Interpolation on Scattered Data
@section Interpolation on Scattered Data
+ An important use of the Delaunay tessellation is that it can be used to
+ interpolate from scattered data to an arbitrary set of points. To do
+ this the N-simplex of the known set of points is calculated with
+ @code{delaunay}, @code{delaunay3} or @code{delaunayn}. Then the
+ simplicies in to which the desired points are found are
+ identified. Finally the vertices of the simplicies are used to
+ interpolate to the desired points. The functions that perform this
+ interpolation are @code{griddata}, @code{griddata3} and
+ @code{griddatan}.
+
@DOCSTRING(griddata)
@DOCSTRING(griddata3)
@DOCSTRING(griddatan)
+
+ An example of the use of the @code{griddata} function is
+
+ @example
+ @group
+ rand("state",1);
+ 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);
+ @end group
+ @end example
+
+ @noindent
+ that interpolates from a random scattering of points, to a uniform
+ grid.
+ @ifnotinfo
+ The output of the above can be seen in @ref{fig:griddata}.
+
+ @float Figure,fig:griddata
+ @image{griddata,8cm}
+ @caption{Interpolation from a scattered data to a regular grid}
+ @end float
+ @end ifnotinfo
*** ./doc/interpreter/octave.texi.orig10 2007-08-24 22:16:38.885848535
+0200
--- ./doc/interpreter/octave.texi 2007-08-24 22:30:03.226265393 +0200
***************
*** 455,461 ****
* Delaunay Triangulation::
* Voronoi Diagrams::
* Convex Hull::
- * Plotting the Triangulation::
* Interpolation on Scattered Data::
Control Theory
--- 455,460 ----
*** ./doc/interpreter/geometryimages.m.orig10 2007-08-24 18:09:41.974281530
+0200
--- ./doc/interpreter/geometryimages.m 2007-08-25 00:30:07.641765295 +0200
***************
*** 32,37 ****
--- 32,112 ----
[xx,yy]=meshgrid(linspace(-1,1,32));
griddata(x,y,z,xx,yy);
print (strcat (nm, ".", typ), strcat ("-d", typ))
+ elseif (strcmp (nm, "convhull"))
+ 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]);
+ print (strcat (nm, ".", typ), strcat ("-d", typ))
+ elseif (strcmp (nm, "delaunay"))
+ rand ("state", 1);
+ x = rand (10, 1);
+ y = rand (10, 1);
+ 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*");
+ print (strcat (nm, ".", typ), strcat ("-d", typ))
+ else
+ error ("unrecognized plot requested");
+ endif
+ bury_output ();
+ endfunction
+
+ function [r, c] = tri2circ (tri, xx, yy)
+ x = xx(tri);
+ y = yy(tri);
+ m = (y(1:end-1) - y(2:end)) ./ (x(1:end-1) - x(2:end));
+ xc = (prod(m) .* (y(1) - y(end)) + m(end)*(x(1)+x(2)) - m(1)*(x(2)+x(3)))
...
+ ./ (2 * (m(end) - m(1)));
+ yc = - (xc - (x(2) + x(3))./2) ./ m(end) + (y(2) + y(3)) / 2;
+ c = [xc, yc];
+ r = sqrt ((xc - x(1)).^2 + (yc - y(1)).^2);
+ endfunction
+
+ ## Use this function before plotting commands and after every call to
+ ## print since print() resets output to stdout (unfortunately, gnpulot
+ ## can't pop output as it can the terminal type).
+ function bury_output ()
+ f = figure (1);
+ set (f, "visible", "off");
+ endfunction
+ function geometryimages (nm, typ)
+ bury_output ();
+ if (strcmp (nm, "voronoi"))
+ rand("state",9);
+ x = rand(10,1);
+ y = rand(10,1);
+ tri = delaunay (x, y);
+ [vx, vy] = voronoi (x, y, tri);
+ triplot (tri, x, y, "b");
+ hold on;
+ plot (vx, vy, "r");
+ [r, c] = tri2circ (tri(end,:), x, y);
+ pc = [-1:0.01:1];
+ xc = r * sin(pi*pc) + c(1);
+ yc = r * cos(pi*pc) + c(2);
+ plot (xc, yc, "g-", "LineWidth", 3);
+ axis([0, 1, 0, 1]);
+ legend ("Delaunay Triangulation", "Voronoi Diagram");
+ print (strcat (nm, ".", typ), strcat ("-d", typ))
+ elseif (strcmp (nm, "triplot"))
+ rand ("state", 2)
+ x = rand (20, 1);
+ y = rand (20, 1);
+ tri = delaunay (x, y);
+ triplot (tri, x, y);
+ print (strcat (nm, ".", typ), strcat ("-d", typ))
+ elseif (strcmp (nm, "griddata"))
+ rand("state",1);
+ 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);
+ print (strcat (nm, ".", typ), strcat ("-d", typ))
else
error ("unrecognized plot requested");
endif
*** ./doc/interpreter/Makefile.in.orig10 2007-08-24 22:29:52.271818101
+0200
--- ./doc/interpreter/Makefile.in 2007-08-25 00:20:04.696187024 +0200
***************
*** 18,24 ****
INSTALL_PROGRAM = @INSTALL_PROGRAM@
INSTALL_DATA = @INSTALL_DATA@
! SCRIPT_SOURCES = sparseimages.m interpimages.m
EXAMPLE_FILES_NODIR = \
addtwomatrices.cc \
--- 18,24 ----
INSTALL_PROGRAM = @INSTALL_PROGRAM@
INSTALL_DATA = @INSTALL_DATA@
! SCRIPT_SOURCES = sparseimages.m interpimages.m geometryimages.m
EXAMPLE_FILES_NODIR = \
addtwomatrices.cc \
***************
*** 43,48 ****
--- 43,53 ----
EXAMPLE_FILES = $(addprefix $(top_srcdir)/examples/, $(EXAMPLE_FILES_NODIR))
+ GEOMETRYIMAGES = voronoi triplot griddata convhull delaunay
+ GEOMETRYIMAGES_EPS = $(addsuffix .eps, $(GEOMETRYIMAGES))
+ GEOMETRYIMAGES_PDF = $(addsuffix .pdf, $(GEOMETRYIMAGES))
+ GEOMETRYIMAGES_PNG = $(addsuffix .png, $(GEOMETRYIMAGES))
+
INTERPIMAGES = interpft interpn interpderiv1 interpderiv2
INTERPIMAGES_EPS = $(addsuffix .eps, $(INTERPIMAGES))
INTERPIMAGES_PDF = $(addsuffix .pdf, $(INTERPIMAGES))
***************
*** 54,62 ****
SPARSEIMAGES_PNG = $(addsuffix .png, $(SPARSEIMAGES_1))
SPARSEIMAGES_TXT = $(addsuffix .txt, $(SPARSEIMAGES_1))
! IMAGES_EPS = $(SPARSEIMAGES_EPS) $(INTERPIMAGES_EPS)
! IMAGES_PDF = $(SPARSEIMAGES_PDF) $(INTERPIMAGES_PDF)
! IMAGES_PNG = $(SPARSEIMAGES_PNG) $(INTERPIMAGES_PNG)
IMAGES_TXT = $(SPARSEIMAGES_TXT)
HTML_IMAGES_PNG = $(addprefix HTML/, $(IMAGES_PNG))
--- 59,70 ----
SPARSEIMAGES_PNG = $(addsuffix .png, $(SPARSEIMAGES_1))
SPARSEIMAGES_TXT = $(addsuffix .txt, $(SPARSEIMAGES_1))
! IMAGES_EPS = $(SPARSEIMAGES_EPS) $(INTERPIMAGES_EPS) \
! $(GEOMETRYIMAGES_EPS)
! IMAGES_PDF = $(SPARSEIMAGES_PDF) $(INTERPIMAGES_PDF) \
! $(GEOMETRYIMAGES_PDF)
! IMAGES_PNG = $(SPARSEIMAGES_PNG) $(INTERPIMAGES_PNG) \
! $(GEOMETRYIMAGES_PNG)
IMAGES_TXT = $(SPARSEIMAGES_TXT)
HTML_IMAGES_PNG = $(addprefix HTML/, $(IMAGES_PNG))
***************
*** 210,216 ****
--eval "$(notdir $(basename $<)) ('$(notdir $(basename $@))', '$(patsubst
.%,%, $(suffix $@))'); sleep (1);"
endef
! $(INTERPIMAGES_EPS) $(INTERPIMAGES_PNG) $(INTERPIMAGES_TXT): interpimages.m
$(run-octave)
$(SPARSEIMAGES_EPS) $(SPARSEIMAGES_PNG) $(SPARSEIMAGES_TXT): sparseimages.m
--- 218,227 ----
--eval "$(notdir $(basename $<)) ('$(notdir $(basename $@))', '$(patsubst
.%,%, $(suffix $@))'); sleep (1);"
endef
! $(GEOMETRYIMAGES_EPS) $(GEOMETRYIMAGES_PNG) $(GEOMETRYIMAGES_TXT):
geometryimages.m
! $(run-octave)
!
! $(INTERPIMAGES_EPS) $(INTEgRPIMAGES_PNG) $(INTERPIMAGES_TXT): interpimages.m
$(run-octave)
$(SPARSEIMAGES_EPS) $(SPARSEIMAGES_PNG) $(SPARSEIMAGES_TXT): sparseimages.m
*** ./scripts/geometry/trimesh.m.orig10 2007-08-24 18:02:07.000000000 +0200
--- ./scripts/geometry/trimesh.m 2007-08-24 22:30:03.228265292 +0200
***************
*** 21,29 ****
## @deftypefn {Function File} {} trimesh (@var{tri}, @var{x}, @var{y},
@var{z})
## @deftypefnx {Function File} address@hidden = } 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
--- 21,30 ----
## @deftypefn {Function File} {} trimesh (@var{tri}, @var{x}, @var{y},
@var{z})
## @deftypefnx {Function File} address@hidden = } trimesh (@dots{})
## Plot a triangular mesh in 3D. The variable @var{tri} is the triangular
! ## meshing of the points @code{(@var{x}, @var{y})} which is returned
! ## from @code{delaunay}. The variable @var{z} is value at the point
! ## @code{(@var{x}, @var{y})}. The output argument @var{h} is the graphic
! ## handle to the plot.
## @seealso{triplot, delaunay3}
## @end deftypefn