# HG changeset patch
# User dbateman
# Date 1188002143 0
# Node ID 3c500bc71e148da9fe7f1a5880ce74e9617ea817
# Parent  af63f57a19aeb86b44f72e09559c773da7ea5fc8
[project @ 2007-08-25 00:35:33 by dbateman]

diff --git a/doc/ChangeLog b/doc/ChangeLog
--- a/doc/ChangeLog
+++ b/doc/ChangeLog
@@ -1,3 +1,13 @@
+2007-08-25  David Bateman  <dbateman@free.fr>
+
+        * interpreter/geometry.txi: Add examples and explanatory text.
+        * interpreter/octave.texi: Update indexing of geometry functions.
+	* interpreter/geometryimage.m: New script to create images for
+	geometry chapter.
+	* interpreter/Makefile.in (SCRIPT_SORCES): add geometryimages.m
+	(GEOMETRYIMAGES*): New variables.
+	(IMAGES_EPS, IMAGES_PDF, IMAGES_PNG): Add the GEOMETRYIMAGES*.
+
 2007-08-24  David Bateman  <dbateman@free.fr>
 
         * interpreter/geometry.txi: Document new functions.
diff --git a/doc/interpreter/Makefile.in b/doc/interpreter/Makefile.in
--- a/doc/interpreter/Makefile.in
+++ b/doc/interpreter/Makefile.in
@@ -18,7 +18,7 @@
 INSTALL_PROGRAM = @INSTALL_PROGRAM@
 INSTALL_DATA = @INSTALL_DATA@
 
-SCRIPT_SOURCES = sparseimages.m interpimages.m
+SCRIPT_SOURCES = sparseimages.m interpimages.m geometryimages.m
 
 EXAMPLE_FILES_NODIR = \
   addtwomatrices.cc \
@@ -43,6 +43,11 @@
 
 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,9 +59,12 @@
 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_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,7 +218,10 @@
     --eval "$(notdir $(basename $<)) ('$(notdir $(basename $@))', '$(patsubst .%,%, $(suffix $@))'); sleep (1);"
 endef
 
-$(INTERPIMAGES_EPS) $(INTERPIMAGES_PNG) $(INTERPIMAGES_TXT): interpimages.m
+$(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
diff --git a/doc/interpreter/geometry.txi b/doc/interpreter/geometry.txi
--- a/doc/interpreter/geometry.txi
+++ b/doc/interpreter/geometry.txi
@@ -6,17 +6,34 @@
 @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::
-* Plotting the Triangulation::
 * 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,10 +47,72 @@
 
 @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,7 +222,7 @@
 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
+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,33 +262,138 @@
 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.
+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)
 
-@node Plotting the Triangulation
-@section Plotting the Triangulation
+An example of the use of @code{convhull} is
 
-@DOCSTRING(triplot)
+@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
 
-@DOCSTRING(trimesh)
+@ifnotinfo
+@noindent
+The output of the above can be seen in @ref{fig:convhull}.
 
-@DOCSTRING(polyarea)
+@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
diff --git a/doc/interpreter/geometryimages.m b/doc/interpreter/geometryimages.m
new file mode 100644
--- /dev/null
+++ b/doc/interpreter/geometryimages.m
@@ -0,0 +1,133 @@
+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))    
+  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
+  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
diff --git a/doc/interpreter/octave.texi b/doc/interpreter/octave.texi
--- a/doc/interpreter/octave.texi
+++ b/doc/interpreter/octave.texi
@@ -455,7 +455,6 @@
 * Delaunay Triangulation::
 * Voronoi Diagrams::
 * Convex Hull::
-* Plotting the Triangulation::
 * Interpolation on Scattered Data::
 
 Control Theory
diff --git a/scripts/geometry/trimesh.m b/scripts/geometry/trimesh.m
--- a/scripts/geometry/trimesh.m
+++ b/scripts/geometry/trimesh.m
@@ -21,9 +21,10 @@
 ## @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.
+## 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