Geometry
********

Much of the geometry code in Octave is based on the Qhull library (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, http://www.qhull.org). Some of the documentation for Qhull, particularly for the options that can be passed to delaunay, voronoi and convhull, etc., is relevant to Octave users.

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 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.

 -- TRI =: delaunay (X, Y)
 -- TETR =: delaunay (X, Y, Z)
 -- TRI =: delaunay (X)
 -- TRI =: delaunay (..., OPTIONS)
     Compute the Delaunay triangulation for a 2-D or 3-D set of points.

     For 2-D sets, the return value TRI is a set of triangles which satisfies the Delaunay circum-circle criterion, i.e., no data point from [X, Y] is within the circum-circle of the defining triangle. The set of triangles TRI is a matrix of size [n, 3].  Each row defines a triangle and the three columns are the three vertices of the triangle. The value of TRI(i,j) is an index into X and Y for the location of the j-th vertex of the i-th triangle.

     For 3-D sets, the return value TETR is a set of tetrahedrons which satisfies the Delaunay circum-circle criterion, i.e., no data point from [X, Y, Z] is within the circum-circle of the defining tetrahedron.  The set of tetrahedrons is a matrix of size [n, 4].  Each row defines a tetrahedron and the four columns are the four vertices of the tetrahedron.  The value of TETR(i,j) is an index into X, Y, Z for the location of the j-th vertex of the i-th tetrahedron.

     The input X may also be a matrix with two or three columns where the first column contains x-data, the second y-data, and the optional third column contains z-data.

     An optional final argument, which must be a string or cell array of strings, contains options passed to the underlying qhull command. See the documentation for the Qhull library for details http://www.qhull.org/html/qh-quick.htm#options. The default options are {"Qt", "Qbb", "Qc"}. If Qhull fails for 2-D input the triangulation is attempted again with the options {"Qt", "Qbb", "Qc", "Qz"} which may result in reduced accuracy.

     If OPTIONS is not present or [] then the default arguments are used.  Otherwise, OPTIONS replaces the default argument list. To append user options to the defaults it is necessary to repeat the default arguments in OPTIONS.  Use a null string to pass no arguments.

          x = rand (1, 10);
          y = rand (1, 10);
          tri = delaunay (x, y);
          triplot (tri, x, y);
          hold on;
          plot (x, y, "r*");
          axis ([0,1,0,1]);

     delaunayn, convhull, voronoi, triplot, trimesh, tetramesh, trisurf

For 3-D inputs delaunay returns a set of tetrahedra that satisfy the Delaunay circum-circle criteria.  Similarly, delaunayn returns the N-dimensional simplex satisfying the Delaunay circum-circle criteria. The N-dimensional extension of a triangulation is called a tessellation.

 -- T =: delaunayn (PTS)
 -- T =: delaunayn (PTS, OPTIONS)
     Compute the Delaunay triangulation for an N-dimensional set of points.

     The Delaunay triangulation is a tessellation of the convex hull of a set of points such that no N-sphere defined by the N-triangles contains any other points from the set.

     The input matrix PTS of size [n, dim] contains n points in a space of dimension dim.  The return matrix T has size [m, dim+1].  Each row of T contains a set of indices back into the original set of points PTS which describes a simplex of dimension dim.  For example, a 2-D simplex is a triangle and 3-D simplex is a tetrahedron.

     An optional second argument, which must be a string or cell array of strings, contains options passed to the underlying qhull command.  See the documentation for the Qhull library for details http://www.qhull.org/html/qh-quick.htm#options. The default options depend on the dimension of the input:

       * 2-D and 3-D: OPTIONS = {"Qt", "Qbb", "Qc"}

       * 4-D and higher: OPTIONS = {"Qt", "Qbb", "Qc", "Qx"}

     If Qhull fails for 2-D input the triangulation is attempted again with the options {"Qt", "Qbb", "Qc", "Qz"} which may result in reduced accuracy.

     If OPTIONS is not present or [] then the default arguments are used.  Otherwise, OPTIONS replaces the default argument list. To append user options to the defaults it is necessary to repeat the default arguments in OPTIONS.  Use a null string to pass no arguments.

     delaunay, convhulln, voronoin, trimesh, tetramesh

An example of a Delaunay triangulation of a set of points is

     rand ("state", 1);
     x = rand (1, 10);
     y = rand (1, 10);
     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*");

Plotting the Triangulation
--------------------------

Octave has the functions triplot, trimesh, and trisurf to plot the Delaunay triangulation of a 2-dimensional set of points. tetramesh will plot the triangulation of a 3-dimensional set of points.

 -- triplot: (TRI, X, Y)
 -- triplot: (TRI, X, Y, LINESPEC)
 -- H =: triplot (...)
     Plot a 2-D triangular mesh.

     TRI is typically the output of a Delaunay triangulation over the grid of X, Y.  Every row of TRI represents one triangle and contains three indices into [X, Y] which are the vertices of the triangles in the x-y plane.

     The linestyle to use for the plot can be defined with the argument LINESPEC of the same format as the plot command.

     The optional return value H is a graphics handle to the created patch object. plot, trimesh, trisurf, delaunay

 -- trimesh: (TRI, X, Y, Z, C)
 -- trimesh: (TRI, X, Y, Z)
 -- trimesh: (TRI, X, Y)
 -- trimesh: (..., PROP, VAL, ...)
 -- H =: trimesh (...)
     Plot a 3-D triangular wireframe mesh.

     In contrast to mesh, which plots a mesh using rectangles, trimesh plots the mesh using triangles.

     TRI is typically the output of a Delaunay triangulation over the grid of X, Y.  Every row of TRI represents one triangle and contains three indices into [X, Y] which are the vertices of the triangles in the x-y plane.  Z determines the height above the plane of each vertex.  If no Z input is given then the triangles are plotted as a 2-D figure.

     The color of the trimesh is computed by linearly scaling the Z values to fit the range of the current colormap.  Use caxis and/or change the colormap to control the appearance.

     Optionally, the color of the mesh can be specified independently of Z by supplying C, which is a vector for colormap data, or a matrix with three columns for RGB data.  The number of colors specified in C must either equal the number of vertices in Z or the number of triangles in TRI.

     Any property/value pairs are passed directly to the underlying patch object. The full list of properties is documented at Patch Properties.

     The optional return value H is a graphics handle to the created patch object. mesh, tetramesh, triplot, trisurf, delaunay, patch, hidden

 -- trisurf: (TRI, X, Y, Z, C)
 -- trisurf: (TRI, X, Y, Z)
 -- trisurf: (..., PROP, VAL, ...)
 -- H =: trisurf (...)
     Plot a 3-D triangular surface.

     In contrast to surf, which plots a surface mesh using rectangles, trisurf plots the mesh using triangles.

     TRI is typically the output of a Delaunay triangulation over the grid of X, Y.  Every row of TRI represents one triangle and contains three indices into [X, Y] which are the vertices of the triangles in the x-y plane.  Z determines the height above the plane of each vertex.

     The color of the trisurf is computed by linearly scaling the Z values to fit the range of the current colormap.  Use caxis and/or change the colormap to control the appearance.

     Optionally, the color of the mesh can be specified independently of Z by supplying C, which is a vector for colormap data, or a matrix with three columns for RGB data.  The number of colors specified in C must either equal the number of vertices in Z or the number of triangles in TRI.  When specifying the color at each vertex the triangle will be colored according to the color of the first vertex only (see patch documentation and the "FaceColor" property when set to "flat").

     Any property/value pairs are passed directly to the underlying patch object. The full list of properties is documented at Patch Properties.

     The optional return value H is a graphics handle to the created patch object. surf, triplot, trimesh, delaunay, patch, shading

 -- tetramesh: (T, X)
 -- tetramesh: (T, X, C)
 -- tetramesh: (..., PROPERTY, VAL, ...)
 -- H =: tetramesh (...)
     Display the tetrahedrons defined in the m-by-4 matrix T as 3-D patches.

     T is typically the output of a Delaunay triangulation of a 3-D set of points.  Every row of T contains four indices into the n-by-3 matrix X of the vertices of a tetrahedron.  Every row in X represents one point in 3-D space.

     The vector C specifies the color of each tetrahedron as an index into the current colormap.  The default value is 1:m where m is the number of tetrahedrons; the indices are scaled to map to the full range of the colormap.  If there are more tetrahedrons than colors in the colormap then the values in C are cyclically repeated.

     Calling tetramesh (..., "property", "value", ...) passes all property/value pairs directly to the patch function as additional arguments. The full list of properties is documented at Patch Properties.

     The optional return value H is a vector of patch handles where each handle represents one tetrahedron in the order given by T. A typical use case for H is to turn the respective patch "visible" property "on" or "off".

     Type demo tetramesh to see examples on using tetramesh. trimesh, delaunay, delaunayn, patch

The difference between triplot, and trimesh or trisurf, is that the former only plots the 2-dimensional triangulation itself, whereas the second two plot the value of a function f (X, Y).  An example of the use of the triplot function is

     rand ("state", 2)
     x = rand (20, 1);
     y = rand (20, 1);
     tri = delaunay (x, y);
     triplot (tri, x, y);

which plots the Delaunay triangulation of a set of random points in 2-dimensions.

Identifying Points in Triangulation
-----------------------------------

It is often necessary to identify whether a particular point in the N-dimensional space is within the Delaunay tessellation of a set of points in this N-dimensional space, and if so which N-simplex contains the point and which point in the tessellation is closest to the desired point.  The functions tsearch and dsearch perform this function in a triangulation, and tsearchn and dsearchn in an N-dimensional tessellation.

To identify whether a particular point represented by a vector 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 N + 1 vectors T(I,:), then the Barycentric coordinates defining the point P are given by

     P = BETA * T

where BETA contains N + 1 values that together as a vector represent the Barycentric coordinates of the point P.  To ensure a unique solution for the values of BETA an additional criteria of

     sum (BETA) == 1

is imposed, and we can therefore write the above as

     P - T(end, :) = BETA(1:end-1) * (T(1:end-1, :)
                     - ones (N, 1) * T(end, :)

Solving for BETA we can then write

     BETA(1:end-1) = (P - T(end, :)) /
                     (T(1:end-1, :) - ones (N, 1) * T(end, :))
     BETA(end) = sum (BETA(1:end-1))

which gives the formula for the conversion of the Cartesian coordinates of the point P to the Barycentric coordinates BETA.  An important property of the Barycentric coordinates is that for all points in the N-simplex

     0 <= BETA(I) <= 1

Therefore, the test in tsearch and 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 BETA.  This is exactly the implementation used in tsearchn.  tsearch is optimized for 2-dimensions and the Barycentric coordinates are not explicitly formed.

 -- IDX =: tsearch (X, Y, T, XI, YI)
     Search for the enclosing Delaunay convex hull.

     For T = delaunay (X, Y), finds the index in T containing the points (XI, YI).  For points outside the convex hull, IDX is NaN.

     Programming Note: The algorithm is O(M*N) for locating M points in N triangles.  Performance is typically much faster if the points to be located are in a single continuous path; a point is first checked against the region its predecessor was found in, speeding up lookups for points along a continuous path.

     delaunay, delaunayn

 -- IDX =: tsearchn (X, T, XI)
 -- [IDX, P] =: tsearchn (X, T, XI)
     Find the simplexes enclosing the given points.

     tsearchn is typically used with delaunayn: T = delaunayn (X) returns a set of simplexes t, then tsearchn returns the row index of T containing each point of XI.  For points outside the convex hull, IDX is NaN.

     If requested, tsearchn also returns the barycentric coordinates P of the enclosing simplexes.

     delaunay, delaunayn, tsearch

An example of the use of tsearch can be seen with the simple triangulation

     X = [-1; -1; 1; 1];
     Y = [-1; 1; -1; 1];
     TRI = [1, 2, 3; 2, 3, 4];

consisting of two triangles defined by TRI.  We can then identify which triangle a point falls in like

     tsearch (X, Y, TRI, -0.5, -0.5)
     => 1
     tsearch (X, Y, TRI, 0.5, 0.5)
     => 2

and we can confirm that a point doesn't lie within one of the triangles like

     tsearch (X, Y, TRI, 2, 2)
     => NaN

The dsearch and 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.

 -- IDX =: dsearch (X, Y, TRI, XI, YI)
 -- IDX =: dsearch (X, Y, TRI, XI, YI, S)
     Return the index IDX of the closest point in X, Y to the elements [XI(:), YI(:)].

     The variable S is accepted for compatibility but is ignored. dsearchn, tsearch

 -- IDX =: dsearchn (X, TRI, XI)
 -- IDX =: dsearchn (X, TRI, XI, OUTVAL)
 -- IDX =: dsearchn (X, XI)
 -- [IDX, D] =: dsearchn (...)
     Return the index IDX of the closest point in X to the elements XI.

     If OUTVAL is supplied, then the values of XI that are not contained within one of the simplices TRI are set to OUTVAL. Generally, TRI is returned from delaunayn (X).

     The optional output D contains a column vector of distances between the query points XI and the nearest simplex points X.

     dsearch, tsearch

An example of the use of dsearch, using the above values of X, Y and TRI is

     dsearch (X, Y, TRI, -2, -2)
     => 1

If you wish the points that are outside the tessellation to be flagged, then dsearchn can be used as

     dsearchn ([X, Y], TRI, [-2, -2], NaN)
     => NaN
     dsearchn ([X, Y], TRI, [-0.5, -0.5], NaN)
     => 1

where the point outside the tessellation are then flagged with NaN.

Voronoi Diagrams
================

A Voronoi diagram or Voronoi tessellation of a set of points S in an N-dimensional space, is the tessellation of the N-dimensional space such that all points in V(P), a partitions of the tessellation where P is a member of S, are closer to P than any other point in S.  The Voronoi diagram is related to the Delaunay triangulation of a set of points, in that the vertices of the Voronoi tessellation are the centers of the circum-circles of the simplices of the Delaunay tessellation.

 -- voronoi: (X, Y)
 -- voronoi: (X, Y, OPTIONS)
 -- voronoi: (..., "linespec")
 -- voronoi: (HAX, ...)
 -- H =: voronoi (...)
 -- [VX, VY] =: voronoi (...)
     Plot the Voronoi diagram of points (X, Y).

     The Voronoi facets with points at infinity are not drawn.

     The OPTIONS argument, which must be a string or cell array of strings, contains options passed to the underlying qhull command. See the documentation for the Qhull library for details http://www.qhull.org/html/qh-quick.htm#options.

     If "linespec" is given it is used to set the color and line style of the plot.

     If an axes graphics handle HAX is supplied then the Voronoi diagram is drawn on the specified axes rather than in a new figure.

     If a single output argument is requested then the Voronoi diagram will be plotted and a graphics handle H to the plot is returned.

     [VX, VY] = voronoi (...) returns the Voronoi vertices instead of plotting the diagram.

          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");

     voronoin, delaunay, convhull

 -- [C, F] =: voronoin (PTS)
 -- [C, F] =: voronoin (PTS, OPTIONS)
     Compute N-dimensional Voronoi facets.

     The input matrix PTS of size [n, dim] contains n points in a space of dimension dim.

     C contains the points of the Voronoi facets.  The list F contains, for each facet, the indices of the Voronoi points.

     An optional second argument, which must be a string or cell array of strings, contains options passed to the underlying qhull command.  See the documentation for the Qhull library for details http://www.qhull.org/html/qh-quick.htm#options.

     The default options depend on the dimension of the input:

       * 2-D and 3-D: OPTIONS = {"Qbb"}

       * 4-D and higher: OPTIONS = {"Qbb", "Qx"}

     If OPTIONS is not present or [] then the default arguments are used.  Otherwise, OPTIONS replaces the default argument list. To append user options to the defaults it is necessary to repeat the default arguments in OPTIONS.  Use a null string to pass no arguments.

     voronoi, convhulln, delaunayn

An example of the use of voronoi is

     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");

Additional information about the size of the facets of a Voronoi diagram, and which points of a set of points is in a polygon can be had with the polyarea and inpolygon functions respectively.

 -- A =: polyarea (X, Y)
 -- A =: polyarea (X, Y, DIM)

     Determine area of a polygon by triangle method.

     The variables X and Y define the vertex pairs, and must therefore have the same shape.  They can be either vectors or arrays.  If they are arrays then the columns of X and Y are treated separately and an area returned for each.

     If the optional DIM argument is given, then polyarea works along this dimension of the arrays X and Y.

An example of the use of polyarea might be

     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

Facets of the Voronoi diagram with a vertex at infinity have infinity area.  A simplified version of polyarea for rectangles is available with rectint

 -- AREA =: rectint (A, B)
     Compute area or volume of intersection of rectangles or N-D boxes.

     Compute the area of intersection of rectangles in A and rectangles in B.  N-dimensional boxes are supported in which case the volume, or hypervolume is computed according to the number of dimensions.

     2-dimensional rectangles are defined as [xpos ypos width height] where xpos and ypos are the position of the bottom left corner.  Higher dimensions are supported where the coordinates for the minimum value of each dimension follow the length of the box in that dimension, e.g., [xpos ypos zpos kpos ... width height depth k_length ...].

     Each row of A and B define a rectangle, and if both define multiple rectangles, then the output, AREA, is a matrix where the i-th row corresponds to the i-th row of a and the j-th column corresponds to the j-th row of b.

     polyarea

 -- IN =: inpolygon (X, Y, XV, YV)
 -- [IN, ON] =: inpolygon (X, Y, XV, YV)

     For a polygon defined by vertex points (XV, YV), return true if the points (X, Y) are inside (or on the boundary) of the polygon; Otherwise, return false.

     The input variables X and Y, must have the same dimension.

     The optional output ON returns true if the points are exactly on the polygon edge, and false otherwise. delaunay

An example of the use of inpolygon might be

     randn ("state", 2);
     x = randn (100, 1);
     y = randn (100, 1);
     vx = cos (pi * [-1 : 0.1: 1]);
     vy = sin (pi * [-1 : 0.1 : 1]);
     in = inpolygon (x, y, vx, vy);
     plot (vx, vy, x(in), y(in), "r+", x(!in), y(!in), "bo");
     axis ([-2, 2, -2, 2]);

Convex Hull
===========

The convex hull of a set of points is the minimum convex envelope containing all of the points.  Octave has the functions convhull and convhulln to calculate the convex hull of 2-dimensional and N-dimensional sets of points.

 -- H =: convhull (X, Y)
 -- H =: convhull (X, Y, Z)
 -- H =: convhull (X)
 -- H =: convhull (..., OPTIONS)
 -- [H, V] =: convhull (...)
     Compute the convex hull of a 2-D or 3-D set of points.

     The hull H is a linear index vector into the original set of points that specifies which points form the enclosing hull.  For 2-D inputs only, the output is ordered in a counterclockwise manner around the hull.

     The input X may also be a matrix with two or three columns where the first column contains x-data, the second y-data, and the optional third column contains z-data.

     An optional final argument, which must be a string or cell array of strings, contains options passed to the underlying qhull command. See the documentation for the Qhull library for details http://www.qhull.org/html/qh-quick.htm#options. The default option is {"Qt"}.

     If OPTIONS is not present or [] then the default arguments are used.  Otherwise, OPTIONS replaces the default argument list. To append user options to the defaults it is necessary to repeat the default arguments in OPTIONS.  Use a null string to pass no arguments.

     If the second output V is requested the volume of the enclosing convex hull is calculated.

     convhulln, delaunay, voronoi

 -- H =: convhulln (PTS)
 -- H =: convhulln (PTS, OPTIONS)
 -- [H, V] =: convhulln (...)
     Compute the convex hull of the set of points PTS.

     PTS is a matrix of size [n, dim] containing n points in a space of dimension dim.

     The hull H is an index vector into the set of points and specifies which points form the enclosing hull.

     An optional second argument, which must be a string or cell array of strings, contains options passed to the underlying qhull command.  See the documentation for the Qhull library for details http://www.qhull.org/html/qh-quick.htm#options. The default options depend on the dimension of the input:

       * 2D, 3D, 4D: OPTIONS = {"Qt"}

       * 5D and higher: OPTIONS = {"Qt", "Qx"}

     If OPTIONS is not present or [] then the default arguments are used.  Otherwise, OPTIONS replaces the default argument list. To append user options to the defaults it is necessary to repeat the default arguments in OPTIONS.  Use a null string to pass no arguments.

     If the second output V is requested the volume of the enclosing convex hull is calculated. convhull, delaunayn, voronoin

An example of the use of convhull is

     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]);

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 delaunay or delaunayn.  Then the simplices in to which the desired points are found are identified.  Finally the vertices of the simplices are used to interpolate to the desired points.  The functions that perform this interpolation are griddata, griddata3 and griddatan.

 -- ZI =: griddata (X, Y, Z, XI, YI)
 -- ZI =: griddata (X, Y, Z, XI, YI, METHOD)
 -- [XI, YI, ZI] =: griddata (...)
 -- VI =: griddata (X, Y, Z, V, XI, YI, ZI)
 -- VI =: griddata (X, Y, Z, V, XI, YI, ZI, METHOD)
 -- VI =: griddata (X, Y, Z, V, XI, YI, ZI, METHOD, OPTIONS)

     Interpolate irregular 2-D and 3-D source data at specified points.

     For 2-D interpolation, the inputs X and Y define the points where the function Z = f (X, Y) is evaluated. The inputs X, Y, Z are either vectors of the same length, or the unequal vectors X, Y are expanded to a 2-D grid with meshgrid and Z is a 2-D matrix matching the resulting size of the X-Y grid.

     The interpolation points are (XI, YI).  If, and only if, XI is a row vector and YI is a column vector, then meshgrid will be used to create a mesh of interpolation points.

     For 3-D interpolation, the inputs X, Y, and Z define the points where the function V = f (X, Y, Z) is evaluated.  The inputs X, Y, Z are either vectors of the same length, or if they are of unequal length, then they are expanded to a 3-D grid with meshgrid.  The size of the input V must match the size of the original data, either as a vector or a matrix.

     The optional input interpolation METHOD can be "nearest", "linear", or for 2-D data "v4".  When the method is "nearest", the output VI will be the closest point in the original data (X, Y, Z) to the query point (XI, YI, ZI).  When the method is "linear", the output VI will be a linear interpolation between the two closest points in the original source data in each dimension.  For 2-D cases only, the "v4" method is also available which implements a biharmonic spline interpolation.  If METHOD is omitted or empty, it defaults to "linear".

     For 3-D interpolation, the optional argument OPTIONS is passed directly to Qhull when computing the Delaunay triangulation used for interpolation.  For more information on the defaults and how to pass different values, see delaunayn.

     Programming Notes: If the input is complex the real and imaginary parts are interpolated separately.  Interpolation is normally based on a Delaunay triangulation.  Any query values outside the convex hull of the input points will return NaN.  However, the "v4" method does not use the triangulation and will return values outside the original data (extrapolation). griddata3, griddatan, delaunay

 -- VI =: griddata3 (X, Y, Z, V, XI, YI, ZI)
 -- VI =: griddata3 (X, Y, Z, V, XI, YI, ZI, METHOD)
 -- VI =: griddata3 (X, Y, Z, V, XI, YI, ZI, METHOD, OPTIONS)

     Interpolate irregular 3-D source data at specified points.

     The inputs X, Y, and Z define the points where the function V = f (X, Y, Z) is evaluated.  The inputs X, Y, Z are either vectors of the same length, or if they are of unequal length, then they are expanded to a 3-D grid with meshgrid.  The size of the input V must match the size of the original data, either as a vector or a matrix.

     The interpolation points are specified by XI, YI, ZI.

     The optional input interpolation METHOD can be "nearest" or "linear".  When the method is "nearest", the output VI will be the closest point in the original data (X, Y, Z) to the query point (XI, YI, ZI).  When the method is "linear", the output VI will be a linear interpolation between the two closest points in the original source data in each dimension. If METHOD is omitted or empty, it defaults to "linear".

     The optional argument OPTIONS is passed directly to Qhull when computing the Delaunay triangulation used for interpolation.  See delaunayn for information on the defaults and how to pass different values.

     Programming Notes: If the input is complex the real and imaginary parts are interpolated separately.  Interpolation is based on a Delaunay triangulation and any query values outside the convex hull of the input points will return NaN. griddata, griddatan, delaunayn

 -- YI =: griddatan (X, Y, XI)
 -- YI =: griddatan (X, Y, XI, METHOD)
 -- YI =: griddatan (X, Y, XI, METHOD, OPTIONS)

     Interpolate irregular source data X, Y at points specified by XI.

     The input X is an MxN matrix representing M points in an N-dimensional space.  The input Y is a single-valued column vector (Mx1) representing a function evaluated at the points X, i.e., Y = fcn (X).  The input XI is a list of points for which the function output YI should be approximated through interpolation.  XI must have the same number of columns (N) as X so that the dimensionality matches.

     The optional input interpolation METHOD can be "nearest" or "linear".  When the method is "nearest", the output YI will be the closest point in the original data X to the query point XI.  When the method is "linear", the output YI will be a linear interpolation between the two closest points in the original source data.  If METHOD is omitted or empty, it defaults to "linear".

     The optional argument OPTIONS is passed directly to Qhull when computing the Delaunay triangulation used for interpolation.  See delaunayn for information on the defaults and how to pass different values.

     Example

          ## Evaluate sombrero() function at irregular data points
          x = 16*gallery ("uniformdata", [200,1], 1) - 8;
          y = 16*gallery ("uniformdata", [200,1], 11) - 8;
          z = sin (sqrt (x.^2 + y.^2)) ./ sqrt (x.^2 + y.^2);
          ## Create a regular grid and interpolate data
          [xi, yi] = ndgrid (linspace (-8, 8, 50));
          zi = griddatan ([x, y], z, [xi(:), yi(:)]);
          zi = reshape (zi, size (xi));
          ## Plot results
          clf ();
          plot3 (x, y, z, "or");
          hold on
          surf (xi, yi, zi);
          legend ("Original Data", "Interpolated Data");

     Programming Notes: If the input is complex the real and imaginary parts are interpolated separately.  Interpolation is based on a Delaunay triangulation and any query values outside the convex hull of the input points will return NaN.  For 2-D and 3-D data additional interpolation methods are available by using the griddata function. griddata, griddata3, delaunayn

An example of the use of the griddata function is

     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));
     zz = griddata (x, y, z, xx, yy);
     mesh (xx, yy, zz);

that interpolates from a random scattering of points, to a uniform grid.

Vector Rotation Matrices
========================

Also included in Octave's geometry functions are primitive functions to enable vector rotations in 3-dimensional space.  Separate functions are provided for rotation about each of the principle axes, X, Y, and Z. According to Euler's rotation theorem, any arbitrary rotation, R, of any vector, P, can be expressed as a product of the three principle rotations:

     p' = Rp = Rz*Ry*Rx*p

 -- T =: rotx (ANGLE)

     rotx returns the 3x3 transformation matrix corresponding to an active rotation of a vector about the x-axis by the specified ANGLE, given in degrees, where a positive angle corresponds to a counterclockwise rotation when viewing the y-z plane from the positive x side.

     The form of the transformation matrix is:

               | 1      0           0      |
           T = | 0  cos(ANGLE) -sin(ANGLE) |
               | 0  sin(ANGLE)  cos(ANGLE) |

     This rotation matrix is intended to be used as a left-multiplying matrix when acting on a column vector, using the notation V = T*U. For example, a vector, U, pointing along the positive y-axis, rotated 90-degrees about the x-axis, will result in a vector pointing along the positive z-axis:

          >> u = [0 1 0]'
          u =
             0
             1
             0

          >> T = rotx (90)
          T =
             1.00000   0.00000   0.00000
             0.00000   0.00000  -1.00000
             0.00000   1.00000   0.00000

          >> v = T*u
          v =
             0.00000
             0.00000
             1.00000

     roty, rotz

 -- T =: roty (ANGLE)

     roty returns the 3x3 transformation matrix corresponding to an active rotation of a vector about the y-axis by the specified ANGLE, given in degrees, where a positive angle corresponds to a counterclockwise rotation when viewing the z-x plane from the positive y side.

     The form of the transformation matrix is:

               |  cos(ANGLE)  0  sin(ANGLE) |
           T = |      0       1      0      |
               | -sin(ANGLE)  0  cos(ANGLE) |

     This rotation matrix is intended to be used as a left-multiplying matrix when acting on a column vector, using the notation V = T*U. For example, a vector, U, pointing along the positive z-axis, rotated 90-degrees about the y-axis, will result in a vector pointing along the positive x-axis:

            >> u = [0 0 1]'
             u =
                0
                0
                1

             >> T = roty (90)
             T =
                0.00000   0.00000   1.00000
                0.00000   1.00000   0.00000
               -1.00000   0.00000   0.00000

             >> v = T*u
             v =
                1.00000
                0.00000
                0.00000

     rotx, rotz

 -- T =: rotz (ANGLE)

     rotz returns the 3x3 transformation matrix corresponding to an active rotation of a vector about the z-axis by the specified ANGLE, given in degrees, where a positive angle corresponds to a counterclockwise rotation when viewing the x-y plane from the positive z side.

     The form of the transformation matrix is:

               | cos(ANGLE) -sin(ANGLE) 0 |
           T = | sin(ANGLE)  cos(ANGLE) 0 |
               |     0           0      1 |

     This rotation matrix is intended to be used as a left-multiplying matrix when acting on a column vector, using the notation V = T*U. For example, a vector, U, pointing along the positive x-axis, rotated 90-degrees about the z-axis, will result in a vector pointing along the positive y-axis:

            >> u = [1 0 0]'
             u =
                1
                0
                0

             >> T = rotz (90)
             T =
                0.00000  -1.00000   0.00000
                1.00000   0.00000   0.00000
                0.00000   0.00000   1.00000

             >> v = T*u
             v =
                0.00000
                1.00000
                0.00000

     rotx, roty