[SciPy-User] [SciPy-user] 3D interpolation

[Burak1327]

Hi,

I have a 3D array of data with dimensions (40x40x40). The majority of data
consists of zeros
with a relative handful of negative numbers (all between -3115 to -3020).
I've been attempting to interpolate this data to a grid of 80x80x80
dimensions.
I want to interpolate the existing data for points in between the existing
points.
The newly interpolated results seem to be incorect. There are now large
positive numbers and
the new minimum value is much lower than the from the coarser grid dataset.
I probably have
misunderstood the purpose of the coordinates variable in map_coordinates().

I've listed a snippet of the code below:

# Read the PES data
pes40 = ReadPES("pes-0.5-noopt.xsf", 40)

# Interpolation
newx,newy,newz = mgrid[0:40:0.5, 0:40:0.5, 0:40:0.5]
coords = array([newx, newy, newz])
pes80 = ndimage.map_coordinates(pes40, coords)

Thanks
Burak
--
View this message in context: http://old.nabble.com/3D-interpolation-tp26919717p26919717.html
Sent from the Scipy-User mailing list archive at Nabble.com.

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

[denis]

Burak,
there's a short 2d example in
http://advice.mechanicalkern.com/question/17/getting-started-with-2d-interpolation-in-scipy
;
does that help ? Let me know please
cheers
-- denis

[Vincent Schut]

On 12/27/2009 12:29 AM, Burak1327 wrote:
>
> Hi,
>
> I have a 3D array of data with dimensions (40x40x40). The majority of data
> consists of zeros
> with a relative handful of negative numbers (all between -3115 to -3020).
> I've been attempting to interpolate this data to a grid of 80x80x80
> dimensions.
> I want to interpolate the existing data for points in between the existing
> points.
> The newly interpolated results seem to be incorect. There are now large
> positive numbers and
> the new minimum value is much lower than the from the coarser grid dataset.
> I probably have
> misunderstood the purpose of the coordinates variable in map_coordinates().

IIRC map_coordinates by default uses 3rd order splines. This may lead to
under/overshoots, which is probably what you're seeing. Try adding
'order=1' to your map_coordinates call, it will then use linear
interpolation and values should never be lower than pes40.min() or
higher than pes40.max().
Your results will be a bit less smooth, though, than with higher order
splines. What I usually do to eliminate this is something like the
following pseudocode:

spline3 = map_coordinates(data, coords, order=3) # 3rd order spline
spline1 = map_coordinates(data, coords, order=1) # linear
spline0 = map_coordinates(data, coords, order=0) # nearest-neighbour
undershoots = spline3 < spline0
overshoots = spline3 > spline0
result = numpy.where(undershoots | overshoots, spline1, spline3)

this will give you 3rd order spline where there are no under/overshoots,
and linear interpolated replacements where 3rd order gives under/overshoots.

Vincent.

>
> I've listed a snippet of the code below:
>
> # Read the PES data
> pes40 = ReadPES("pes-0.5-noopt.xsf", 40)
>
> # Interpolation
> newx,newy,newz = mgrid[0:40:0.5, 0:40:0.5, 0:40:0.5]
> coords = array([newx, newy, newz])
> pes80 = ndimage.map_coordinates(pes40, coords)
>
> Thanks
> Burak

[denis]

Vincent, you're right, splines can overshoot, though not by a lot:
0 1 1 0 overshoots by 20 % in the code snippet below.
Burak, does order=1, just averaging 8 neighbors, work ?

There are many kinds of cubic splines, including
-- global cubic, which map_coordinates seems to use
-- local B-spline: does not interpolate, go through the input points,
but can't overshoot, stays in the convex hull of the input points
-- local Catmull-Rom: interpolates, but can overshoot (over shoots and
leaves :)
("Local" means that e.g. in 1d, out[ 3 to 4 ] is calculated from
the 4 nearest input points at 2 3 4 5,
or equivalently that in[3] affects only out[ 1 to 5 ]: a desirable
property.)

Different corners of Scipy have different spline routines, some with
doc.
To see exactly what they do, you can only plot their impulse response
like so.

cheers
-- denis

""" plot map_coordinates( 1d spike )
"""

from __future__ import division
import sys
import numpy as np
from scipy.ndimage import map_coordinates
import pylab as pl

N = 10
exec( "\n".join( sys.argv[1:] )) # N= ...
np.set_printoptions( 2, threshold=100, suppress=True ) # .2f

a = np.r_[ 5*[0], 1., 5*[0], 1,1, 5*[0], 5*[1] ] # 0 1 1 0 -> 1.2
na = len(a) - 1
x = np.linspace( 0, na, na*N + 1 )
print "a:", a.astype(int)

for order in (3,):
z = map_coordinates( a, [x], order=order )
print "z[%d] halfint:" % order, z[N//2::N]
pl.plot( x, z )

pl.show()

[Burak1327]

Thank you everyone for your replies.

I tried Vincent's advice, and it seems to work. Though, the overshoot is
much more
than 20% in this case. Below I listed the minimum and maximum values
of the data set, after interpolation at different orders:

Raw Data
Minimium Value : -3118.532224
Maximum Value : 0.0 (There next value is about -3115)

Order = 0
Minimium Value : -3118.532224
Maximum Value : 0.0

Order = 1
Minimium Value : -3118.532224
Maximum Value : 0.0

Order = 2
Minimium Value : -3905.68081297
Maximum Value : 851.160340484

Order = 3
Minimium Value : -4077.10199215
Maximum Value : 1028.48609365

Thanks
Burak

--
View this message in context: http://old.nabble.com/3D-interpolation-tp26919717p26954588.html

Sent from the Scipy-User mailing list archive at Nabble.com.

_______________________________________________