Sorry this is such a specialized question, but I went through posts in
the numerical analysis/OR groups and haven't found a solution yet:

What I am looking for is multidimensional cubic spline that:
1) Works for an arbitrary number of dimensions.  Though I will know it
at compile time, so template arguments are reasonable.
2) Modern C++ code if possible, but I could wrap solid C code if
necessary.
3) Can return the gradient.  Hessian if possible, but I guess i could
go finite differences for this if absolutely necessary.  The algorithm
for derivatives at nodes doesn't matter much to me.
4) Can take const pointers to data (or equivalent C++ concepts) and
set of extents.  And it doesn't copy all of the grid and node values
into its own data structure.  Too inefficient for me since my grid may
big and the interpolator needs to be generated fairly frequently.
5) Passing in a regular grid is fine, I don't need anything sparse.
Though a solution to interpolation of irregular grids would be really
nice.
6) Ideally, the types for the grid axis could be different along each
axis as long as the type supports some kind of linear space concept.
But I can live without this for now.
7) If there is similar code that fits the bill for chebyshev or some
other polynomial interpolation, I am willing to take it.  I am not all
that picky on the interpolation method.
8) LGPL or more liberal licensing.

As a note, I found 
http://quantlib.org/reference/class_quant_lib_1_1_multi_cubic_spline
which seems to work.  The problems are:
1) It uses its own data structures and hence would require me to copy
the entire grid, grid values.  In the longterm, this is a deal killer
for me.  My current data structure is boost::multi_array, but I would
love it to be generic enough to work with any "recursive" iterator
concept where you can keep getting new iterators until you run out of
dimensions.
2) It has no tricks to deal with grids of only a few points in one or
more dimensions (reverting to linear interpolation, etc).  This is a
real pain for my purposes, and I needed to contrive some points at the
beginning to work.  But this necessitates even more expensive data
copying to pad data.
3) It doesn't have the hessian, which I will need for my optimizer.
4) It can't evaluate at the boundaries, which seems a pretty nasty bug
since I will frequently have corner solutions.
5) I am not smart enough to change its elegant code to fix these
problems.

Thanks,
Jesse

-- 
      [ See http://www.gotw.ca/resources/clcm.htm for info about ]
      [ comp.lang.c++.moderated.    First time posters: Do this! ]

On 24 Jan., 20:12, Jesse Perla <jessepe...@gmail.com> wrote:
> What I am looking for is multidimensional cubic spline

I can't suggest a library that fulfills all your requirements.  But I
can recommend a very good paper on it.  It even includes C code if I
remember correctly:

http://bigwww.epfl.ch/publications/thevenaz9901.pdf

It is a detailed coverage of various interpolation and approximation
approaches for 2D (medical) images (including spectral properties of
reconstruction filters and speed/quality graphs).  Extending their
cubic B-Spline code to higher dimensions is straight forward.

BTW: You may want to skip the weights computation part.  It's just a
filter that amplifies high (spatial) frequencies.  In my applications
I skipped that step and got better results even though the curve
doesn't interpolate anymore.  It depends on what you're trying to do.

[snip]

> 5) Passing in a regular grid is fine, I don't need anything sparse.
> Though a solution to interpolation of irregular grids would be really
> nice.

Oh!  I don't think there's even a method thats widely agreed upon for
reconstructing irregularly sampled data.  For regularly sampled data
all methods somehow try to approximate sinc interpolation.

[snip]


Cheers!
SG

-- 
      [ See http://www.gotw.ca/resources/clcm.htm for info about ]
      [ comp.lang.c++.moderated.    First time posters: Do this! ]