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

[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