Mercurial > hg > octave-avbm
changeset 6832:3c500bc71e14
[project @ 2007-08-25 00:35:33 by dbateman]
author | dbateman |
---|---|
date | Sat, 25 Aug 2007 00:35:43 +0000 |
parents | af63f57a19ae |
children | e8a18d380097 |
files | doc/ChangeLog doc/interpreter/Makefile.in doc/interpreter/geometry.txi doc/interpreter/geometryimages.m doc/interpreter/octave.texi scripts/geometry/trimesh.m |
diffstat | 6 files changed, 355 insertions(+), 17 deletions(-) [+] |
line wrap: on
line diff
--- 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.
--- 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
--- 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
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
--- 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
--- 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