[SciPy-user] [Fwd: 3D interpolation over irregular data]

481 views
Skip to first unread message

mark starnes

unread,
Jul 30, 2008, 11:11:30 AM7/30/08
to scipy...@scipy.org
Hi everyone,

I've looked through the list here and in Numpy-users, and checked the
'net but can't find an answer to this problem (with luck, I've missed
something obvious!).

I've an array of velocities at 80,000 points, irregularly spaced (from a
CFD analysis). I'd like to generate the interpolated velocity at any
position in the domain, to map the data to an acoustics analysis on a
different mesh.

I tried a least squares approach but the errors are too large using
polynomials and trigonometric functions. My conclusion is that I need a
nearest-neighbour type interpolation routine.

Is there such a routine in Scipy or Numpy? From inspection of Matlab's
functions 'interp3' may be similar to what I would like, but I can't test.

If there's nothing available, I'll end up converting one of the
submitted scripts at the url,

http://www.mathworks.com/matlabcentral/fileexchange/loadCategory.do?objectId=14&objectType=Category


Thanks in advance,

Mark.

_______________________________________________
SciPy-user mailing list
SciPy...@scipy.org
http://projects.scipy.org/mailman/listinfo/scipy-user

José María García Pérez

unread,
Jul 30, 2008, 11:18:02 AM7/30/08
to SciPy Users List
Hi Mark,
Try Radial Basis Functions. But I'm not sure if it's implemented at SciPy/Numpy. I implemented myself long time ago, but better if there is something already done within the project.
Regards,
José M.

2008/7/30 mark starnes <m.sta...@imperial.ac.uk>

Zachary Pincus

unread,
Jul 30, 2008, 11:55:09 AM7/30/08
to SciPy Users List
Mark:

> My conclusion is that I need a
> nearest-neighbour type interpolation routine.

José:
> Try Radial Basis Functions.

This was my thought as well.

Scipy has a handy implementation:
http://www.scipy.org/Cookbook/RadialBasisFunctions

Zach

James Turner

unread,
Jul 30, 2008, 2:46:25 PM7/30/08
to SciPy Users List
Hi Mark,

> I've an array of velocities at 80,000 points, irregularly spaced (from a
> CFD analysis). I'd like to generate the interpolated velocity at any
> position in the domain, to map the data to an acoustics analysis on a
> different mesh.
>
> I tried a least squares approach but the errors are too large using
> polynomials and trigonometric functions. My conclusion is that I need a
> nearest-neighbour type interpolation routine.

If you don't have any joy with the other suggestions, there are some
useful routines by Renka in Netlib, at least one of which is 3D,
though I believe they are also piecewise polynomials (and in ForTran,
of course, so need wrapping).

http://www.netlib.org/toms/

I also have a program implementing the Hubble Telescope's "Drizzle"
algorithm in 3D, which is similar to nearest neighbour or linear,
but unfortunately it's written in an obscure language that would
make it difficult to use from Python unless you are prepared to
install STScI's PyRAF. I also suspect it will only work well for
data that are mildly irregular and it may be fiddly to set up for
your application, but anyway, just so you have the option...
The paper is at http://xxx.lanl.gov/abs/astro-ph/9708242.

Are the points really irregular in 3D, or just in 1-2 of the
dimensions? If any dimension is regularly spaced, you might be able
to simplify the problem by resampling in 1-2D at a time?

Stefan might also have suggestions ;-). But maybe the radial basis
functions will do the trick... I'll have to look at those.

Cheers,

James.

mark starnes

unread,
Jul 30, 2008, 4:15:32 PM7/30/08
to SciPy Users List
José, Zachary,

I seem to be having some success running the example from the link
posted! Thanks for your rapid and helpful responses.

Best regards,

Mark.

mark starnes

unread,
Jul 30, 2008, 4:32:31 PM7/30/08
to SciPy Users List
Hi James,

I'm having a go with the radial basis functions, as they seem to be
working with my Scipy (surprisingly!). Thanks for your suggestions;
I'll move onto them if I fail with this approach.

The points are partly on a structured grid, though not regular, and
partly an unstructured grid. I was hoping not to treat the regions
differently.

Thanks for your response, though. I'll be sure to write back if the
radial functions don't work.

Best regards,

Mark.

Stéfan van der Walt

unread,
Jul 30, 2008, 5:11:36 PM7/30/08
to SciPy Users List
Hi Mark

2008/7/30 mark starnes <m.sta...@imperial.ac.uk>:


> Hi everyone,
>
> I've looked through the list here and in Numpy-users, and checked the
> 'net but can't find an answer to this problem (with luck, I've missed
> something obvious!).
>
> I've an array of velocities at 80,000 points, irregularly spaced (from a
> CFD analysis). I'd like to generate the interpolated velocity at any
> position in the domain, to map the data to an acoustics analysis on a
> different mesh.
>
> I tried a least squares approach but the errors are too large using
> polynomials and trigonometric functions. My conclusion is that I need a
> nearest-neighbour type interpolation routine.

Also take a look at

http://www.scipy.org/Cookbook/Matplotlib/Gridding_irregularly_spaced_data

Robert Kern's delaunay package does natural neighbour interpolation, IIRC.

Regards
Stéfan

mark starnes

unread,
Jul 31, 2008, 4:13:26 AM7/31/08
to SciPy Users List
Hi Stéfan, thanks for the tip. I tried this approach before, but
couldn't get it to work, due to my installation (I broke it, got it
working again, then wasn't brave enough to tinker). However, I got the
radial basis functions working by modifying my PYTHONPATH, rather than
with a rebuild, so may try the Delauney / natgrid approaches too. The
radial basis function seems to have trouble with more than 5000 points
so I may need other options.

Thanks for taking the time to help!

Best regards,

Mark.

mark starnes

unread,
Jul 31, 2008, 4:58:05 AM7/31/08
to SciPy Users List
Hi again, Stéfan.

Is Robert Kern's package limited to two-dimensional data? I've had a
look and can't see any three-dimensional options.

Best regards,

Mark.

Stéfan van der Walt wrote:

José María García Pérez

unread,
Jul 31, 2008, 5:18:50 AM7/31/08
to SciPy Users List
Mark,
From my experience working with RBF, they work pretty well even when you use few points for the interpolation. They track non linear behavior very well.
I have worked with RBF with very big FEM models (200000-500000 grid points) and with more than 3D (in other disciplines), but I don't take all the points at the same. What I would do is to use for example the 100 nearest points to the geometric point where you want to interpolate (probably with 10 would be enough). That's something you can try: test using 10, 20, 50, 100 points, and you will see that the difference is small pretty soon, and for sure smaller that the error you may expect from a CFD simulation.
Hope this tip helps!
José M.


2008/7/31 mark starnes <m.sta...@imperial.ac.uk>

mark starnes

unread,
Jul 31, 2008, 7:12:10 AM7/31/08
to SciPy Users List
Hi José,

Thanks for your comments. Knowing that other people split the
interpolation into regions is useful. I'll be attempting an approach
similar to that you suggested; splitting the domain, not around each
desired point but into regions containing sets of points.

I didn't expect the effect you mention about accuracy being related to
the number of points. Looks like I need to read up on the theory....

Thanks again,

Mark.


José María García Pérez wrote:
> Mark,
> From my experience working with RBF, they work pretty well even when you
> use few points for the interpolation. They track non linear behavior
> very well.
> I have worked with RBF with very big FEM models (200000-500000 grid
> points) and with more than 3D (in other disciplines), but I don't take
> all the points at the same. What I would do is to use for example the
> 100 nearest points to the geometric point where you want to interpolate
> (probably with 10 would be enough). That's something you can try: test
> using 10, 20, 50, 100 points, and you will see that the difference is
> small pretty soon, and for sure smaller that the error you may expect
> from a CFD simulation.
> Hope this tip helps!
> José M.
>
>
> 2008/7/31 mark starnes <m.sta...@imperial.ac.uk

> <mailto:m.sta...@imperial.ac.uk>>


>
> Hi again, Stéfan.
>
> Is Robert Kern's package limited to two-dimensional data? I've had a
> look and can't see any three-dimensional options.
>
> Best regards,
>
> Mark.
>
> Stéfan van der Walt wrote:
> > Hi Mark
> >
> > 2008/7/30 mark starnes <m.sta...@imperial.ac.uk

> <mailto:m.sta...@imperial.ac.uk>>:


> >> Hi everyone,
> >>
> >> I've looked through the list here and in Numpy-users, and checked the
> >> 'net but can't find an answer to this problem (with luck, I've missed
> >> something obvious!).
> >>
> >> I've an array of velocities at 80,000 points, irregularly spaced
> (from a
> >> CFD analysis). I'd like to generate the interpolated velocity at any
> >> position in the domain, to map the data to an acoustics analysis on a
> >> different mesh.
> >>
> >> I tried a least squares approach but the errors are too large using
> >> polynomials and trigonometric functions. My conclusion is that I
> need a
> >> nearest-neighbour type interpolation routine.
> >
> > Also take a look at
> >
> >
> http://www.scipy.org/Cookbook/Matplotlib/Gridding_irregularly_spaced_data
> >
> > Robert Kern's delaunay package does natural neighbour
> interpolation, IIRC.
> >
> > Regards
> > Stéfan
> > _______________________________________________
> > SciPy-user mailing list

> > SciPy...@scipy.org <mailto:SciPy...@scipy.org>


> > http://projects.scipy.org/mailman/listinfo/scipy-user
> >
> _______________________________________________
> SciPy-user mailing list

> SciPy...@scipy.org <mailto:SciPy...@scipy.org>
> http://projects.scipy.org/mailman/listinfo/scipy-user
>
>
>
> ------------------------------------------------------------------------

Stéfan van der Walt

unread,
Jul 31, 2008, 8:32:38 AM7/31/08
to SciPy Users List
2008/7/31 mark starnes <m.sta...@imperial.ac.uk>:

> Is Robert Kern's package limited to two-dimensional data? I've had a
> look and can't see any three-dimensional options.

Sorry, Mark, my mistake, it doesn't.

The only Python wrapper for 3-dimensional Delaunay I'm aware of is

from enthought.tvtk.api import tvtk
d = tvtk.Delaunay3D()

Unfortunately, I don't have much more detail on how to use it.

Reply all
Reply to author
Forward
0 new messages