octave: Multi-dimensional Interpolation

 
 29.2 Multi-dimensional Interpolation
 ====================================
 
 There are three multi-dimensional interpolation functions in Octave,
 with similar capabilities.  Methods using Delaunay tessellation are
 described in SeeInterpolation on Scattered Data.
 
  -- : ZI = interp2 (X, Y, Z, XI, YI)
  -- : ZI = interp2 (Z, XI, YI)
  -- : ZI = interp2 (Z, N)
  -- : ZI = interp2 (Z)
  -- : ZI = interp2 (..., METHOD)
  -- : ZI = interp2 (..., METHOD, EXTRAP)
 
      Two-dimensional interpolation.
 
      Interpolate reference data X, Y, Z to determine ZI at the
      coordinates XI, YI.  The reference data X, Y can be matrices, as
      returned by ‘meshgrid’, in which case the sizes of X, Y, and Z must
      be equal.  If X, Y are vectors describing a grid then ‘length (X)
      == columns (Z)’ and ‘length (Y) == rows (Z)’.  In either case the
      input data must be strictly monotonic.
 
      If called without X, Y, and just a single reference data matrix Z,
      the 2-D region ‘X = 1:columns (Z), Y = 1:rows (Z)’ is assumed.
      This saves memory if the grid is regular and the distance between
      points is not important.
 
      If called with a single reference data matrix Z and a refinement
      value N, then perform interpolation over a grid where each original
      interval has been recursively subdivided N times.  This results in
      ‘2^N-1’ additional points for every interval in the original grid.
      If N is omitted a value of 1 is used.  As an example, the interval
      [0,1] with ‘N==2’ results in a refined interval with points at [0,
      1/4, 1/2, 3/4, 1].
 
      The interpolation METHOD is one of:
 
      "nearest"
           Return the nearest neighbor.
 
      "linear" (default)
           Linear interpolation from nearest neighbors.
 
      "pchip"
           Piecewise cubic Hermite interpolating
           polynomial—shape-preserving interpolation with smooth first
           derivative.
 
      "cubic"
           Cubic interpolation (same as "pchip").
 
      "spline"
           Cubic spline interpolation—smooth first and second derivatives
           throughout the curve.
 
      EXTRAP is a scalar number.  It replaces values beyond the endpoints
      with EXTRAP.  Note that if EXTRAP is used, METHOD must be specified
      as well.  If EXTRAP is omitted and the METHOD is "spline", then the
      extrapolated values of the "spline" are used.  Otherwise the
      default EXTRAP value for any other METHOD is "NA".
 
      See also: Seeinterp1 XREFinterp1, Seeinterp3 XREFinterp3,
      Seeinterpn XREFinterpn, Seemeshgrid XREFmeshgrid.
 
  -- : VI = interp3 (X, Y, Z, V, XI, YI, ZI)
  -- : VI = interp3 (V, XI, YI, ZI)
  -- : VI = interp3 (V, N)
  -- : VI = interp3 (V)
  -- : VI = interp3 (..., METHOD)
  -- : VI = interp3 (..., METHOD, EXTRAPVAL)
 
      Three-dimensional interpolation.
 
      Interpolate reference data X, Y, Z, V to determine VI at the
      coordinates XI, YI, ZI.  The reference data X, Y, Z can be
      matrices, as returned by ‘meshgrid’, in which case the sizes of X,
      Y, Z, and V must be equal.  If X, Y, Z are vectors describing a
      cubic grid then ‘length (X) == columns (V)’, ‘length (Y) == rows
      (V)’, and ‘length (Z) == size (V, 3)’.  In either case the input
      data must be strictly monotonic.
 
      If called without X, Y, Z, and just a single reference data matrix
      V, the 3-D region ‘X = 1:columns (V), Y = 1:rows (V), Z = 1:size
      (V, 3)’ is assumed.  This saves memory if the grid is regular and
      the distance between points is not important.
 
      If called with a single reference data matrix V and a refinement
      value N, then perform interpolation over a 3-D grid where each
      original interval has been recursively subdivided N times.  This
      results in ‘2^N-1’ additional points for every interval in the
      original grid.  If N is omitted a value of 1 is used.  As an
      example, the interval [0,1] with ‘N==2’ results in a refined
      interval with points at [0, 1/4, 1/2, 3/4, 1].
 
      The interpolation METHOD is one of:
 
      "nearest"
           Return the nearest neighbor.
 
      "linear" (default)
           Linear interpolation from nearest neighbors.
 
      "cubic"
           Piecewise cubic Hermite interpolating
           polynomial—shape-preserving interpolation with smooth first
           derivative (not implemented yet).
 
      "spline"
           Cubic spline interpolation—smooth first and second derivatives
           throughout the curve.
 
      EXTRAPVAL is a scalar number.  It replaces values beyond the
      endpoints with EXTRAPVAL.  Note that if EXTRAPVAL is used, METHOD
      must be specified as well.  If EXTRAPVAL is omitted and the METHOD
      is "spline", then the extrapolated values of the "spline" are used.
      Otherwise the default EXTRAPVAL value for any other METHOD is "NA".
 
      See also: Seeinterp1 XREFinterp1, Seeinterp2 XREFinterp2,
      Seeinterpn XREFinterpn, Seemeshgrid XREFmeshgrid.
 
  -- : VI = interpn (X1, X2, ..., V, Y1, Y2, ...)
  -- : VI = interpn (V, Y1, Y2, ...)
  -- : VI = interpn (V, M)
  -- : VI = interpn (V)
  -- : VI = interpn (..., METHOD)
  -- : VI = interpn (..., METHOD, EXTRAPVAL)
 
      Perform N-dimensional interpolation, where N is at least two.
 
      Each element of the N-dimensional array V represents a value at a
      location given by the parameters X1, X2, ..., XN.  The parameters
      X1, X2, ..., XN are either N-dimensional arrays of the same size as
      the array V in the "ndgrid" format or vectors.  The parameters Y1,
      etc.  respect a similar format to X1, etc., and they represent the
      points at which the array VI is interpolated.
 
      If X1, ..., XN are omitted, they are assumed to be ‘x1 = 1 : size
      (V, 1)’, etc.  If M is specified, then the interpolation adds a
      point half way between each of the interpolation points.  This
      process is performed M times.  If only V is specified, then M is
      assumed to be ‘1’.
 
      The interpolation METHOD is one of:
 
      "nearest"
           Return the nearest neighbor.
 
      "linear" (default)
           Linear interpolation from nearest neighbors.
 
      "pchip"
           Piecewise cubic Hermite interpolating
           polynomial—shape-preserving interpolation with smooth first
           derivative (not implemented yet).
 
      "cubic"
           Cubic interpolation (same as "pchip" [not implemented yet]).
 
      "spline"
           Cubic spline interpolation—smooth first and second derivatives
           throughout the curve.
 
      The default method is "linear".
 
      EXTRAPVAL is a scalar number.  It replaces values beyond the
      endpoints with EXTRAPVAL.  Note that if EXTRAPVAL is used, METHOD
      must be specified as well.  If EXTRAPVAL is omitted and the METHOD
      is "spline", then the extrapolated values of the "spline" are used.
      Otherwise the default EXTRAPVAL value for any other METHOD is "NA".
 
      See also: Seeinterp1 XREFinterp1, Seeinterp2 XREFinterp2,
DONTPRINTYET       Seeinterp3 XREFinterp3, Seespline XREFspline, *notendgrid:
DONTPRINTYET       Seeinterp3 XREFinterp3, Seespline XREFspline, Seendgrid

      XREFndgrid.
 
    A significant difference between ‘interpn’ and the other two
 multi-dimensional interpolation functions is the fashion in which the
 dimensions are treated.  For ‘interp2’ and ‘interp3’, the y-axis is
 considered to be the columns of the matrix, whereas the x-axis
 corresponds to the rows of the array.  As Octave indexes arrays in
 column major order, the first dimension of any array is the columns, and
 so ‘interpn’ effectively reverses the ’x’ and ’y’ dimensions.  Consider
 the example,
 
      x = y = z = -1:1;
      f = @(x,y,z) x.^2 - y - z.^2;
      [xx, yy, zz] = meshgrid (x, y, z);
      v = f (xx,yy,zz);
      xi = yi = zi = -1:0.1:1;
      [xxi, yyi, zzi] = meshgrid (xi, yi, zi);
      vi = interp3 (x, y, z, v, xxi, yyi, zzi, "spline");
      [xxi, yyi, zzi] = ndgrid (xi, yi, zi);
      vi2 = interpn (x, y, z, v, xxi, yyi, zzi, "spline");
      mesh (zi, yi, squeeze (vi2(1,:,:)));
 
 where ‘vi’ and ‘vi2’ are identical.  The reversal of the dimensions is
 treated in the ‘meshgrid’ and ‘ndgrid’ functions respectively.