# HG changeset patch # User jwe # Date 1171525948 0 # Node ID 08a261edd738f80be864f2dae498f8ab8b6f5223 # Parent 1ef6ac75795108f0c304ab3f1b8f6bbdf13b77e5 [project @ 2007-02-15 07:52:28 by jwe] diff -r 1ef6ac757951 -r 08a261edd738 src/DLD-FUNCTIONS/interpn.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/DLD-FUNCTIONS/interpn.cc Thu Feb 15 07:52:28 2007 +0000 @@ -0,0 +1,261 @@ +/* + +Copyright (C) 2007 Alexander Barth + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, write to the Free +Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +02110-1301, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include "dMatrix.cc" + +#include "defun-dld.h" +#include "error.h" +#include "oct-obj.h" + +// equivalent to isvector.m + +bool +isvector (const NDArray& array) +{ + const dim_vector dv = array.dims (); + return dv.length () == 2 && (dv(0) == 1 || dv(1) == 1); +} + +// lookup a value in a sorted table (lookup.m) +octave_idx_type +lookup (const double *x, octave_idx_type n, double y) +{ + octave_idx_type j, j0, j1; + + if (y > x[n-1] || y < x[0]) + return -1; + +#ifdef EXHAUSTIF + for (j = 0; j < n - 1; j++) + { + if (x[j] <= y && y <= x[j+1]) + return j; + } +#else + j0 = 0; + j1 = n - 1; + + while (true) + { + j = (j0+j1)/2; + + if (y <= x[j+1]) + { + if (x[j] <= y) + return j; + + j1 = j; + } + + if (x[j] <= y) + j0 = j; + } + +#endif +} + +// n-dimensional linear interpolation + +void +lin_interpn (int n, const octave_idx_type *size, const octave_idx_type *scale, + octave_idx_type Ni, double extrapval, const double **x, + const double *v, const double **y, double *vi) +{ + bool out = false; + int bit; + + OCTAVE_LOCAL_BUFFER (double, coef, 2*n); + OCTAVE_LOCAL_BUFFER (octave_idx_type, index, n); + + // loop over all points + for (octave_idx_type m = 0; m < Ni; m++) + { + // loop over all dimensions + for (int i = 0; i < n; i++) + { + index[i] = lookup (x[i], size[i], y[i][m]); + out = index[i] == -1; + + if (out) + break; + else + { + octave_idx_type j = index[i]; + coef[2*i+1] = (y[i][m] - x[i][j])/(x[i][j+1] - x[i][j]); + coef[2*i] = 1 - coef[2*i+1]; + } + } + + + if (out) + vi[m] = extrapval; + else + { + vi[m] = 0; + + // loop over all corners of hypercube (1<> j & 1; + l += scale[j] * (index[j] + bit); + c *= coef[2*j+bit]; + } + + vi[m] += c * v[l]; + } + } + } +} + +DEFUN_DLD (interpn, args, , + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{vi} =} interpn (@var{x1}, @var{x2}, @dots{}, @var{xn}, @var{v}, @var{y1}, @var{y2}, @dots{}, @var{yn})\n\ +Perform @var{n}-dimensional interpolation. Each element of then\n\ +@var{n}-dimensional array @var{v} represents a value at a location\n\ +given by the parameters @var{x1}, @var{x2},...,@var{xn}. The parameters\n\ +@var{x1}, @var{x2}, @dots{}, @var{xn} are either @var{n}-dimensional\n\ +arrays of the same size as the array @var{v} in the \"ndgrid\" format\n\ +or vectors. The parameters @var{y1}, @var{y2}, @dots{}, @var{yn} are\n\ +all @var{n}-dimensional arrays of the same size and represent the\n\ +points at which the array @var{vi} is interpolated.\n\ +\n\ +This function only performs linear interpolation.\n\ +@seealso{interp1, interp2, ndgrid}\n\ +@end deftypefn") +{ + octave_value retval; + + int nargin = args.length (); + + if (nargin < 2 || nargin % 2 == 0) + { + print_usage (); + return retval; + } + + // dimension of the problem + int n = (nargin-1)/2; + + OCTAVE_LOCAL_BUFFER (NDArray, X, n); + OCTAVE_LOCAL_BUFFER (NDArray, Y, n); + + OCTAVE_LOCAL_BUFFER (const double *, x, n); + OCTAVE_LOCAL_BUFFER (const double *, y, n); + OCTAVE_LOCAL_BUFFER (octave_idx_type, scale, n); + OCTAVE_LOCAL_BUFFER (octave_idx_type, size, n); + + const NDArray V = args(n).array_value (); + NDArray Vi = NDArray (args(n+1).dims ()); + + if (error_state) + { + print_usage (); + return retval; + } + + const double *v = V.data (); + double *vi = Vi.fortran_vec (); + octave_idx_type Ni = Vi.numel (); + + double extrapval = octave_NaN; + + for (int i = 0; i < n; i++) + { + X[i] = args(i).array_value (); + Y[i] = args(n+i+1).array_value (); + + if (error_state) + { + print_usage (); + return retval; + } + + y[i] = Y[i].data (); + size[i] = V.dims()(i); + + if (Y[0].dims () != Y[i].dims ()) + { + error ("interpn: incompatible size of argument number %d", n+i+2); + return retval; + } + } + + // offset in memory of each dimension + + scale[0] = 1; + + for (int i = 1; i < n; i++) + scale[i] = scale[i-1] * size[i-1]; + + // tests if X[0] is a vector, if yes, assume that all elements of X are + // in the ndgrid format. + + if (! isvector (X[0])) + { + for (int i = 0; i < n; i++) + { + if (X[i].dims () != V.dims ()) + { + error ("interpn: incompatible size of argument number %d", i+1); + return retval; + } + else + { + NDArray tmp = NDArray (dim_vector (size[i], 1)); + + for (octave_idx_type j = 0; j < size[i]; j++) + tmp(j) = X[i](scale[i]*j); + + X[i] = tmp; + } + } + } + + for (int i = 0; i < n; i++) + { + if (! isvector (X[i]) && X[i].numel () != size[i]) + { + error ("interpn: incompatible size of argument number %d", i+1); + return retval; + } + else + x[i] = X[i].data (); + } + + lin_interpn (n, size, scale, Ni, extrapval, x, v, y, vi); + + retval = Vi; + + return retval; +}