Re: Latest port of Qhull/goemetry functions from octave-forge

[David Bateman]

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