Any differences in 3d interpolation between MATLAB and Numpy/Scipy?

Question

I'm a MATLAB user and I'm trying to translate some code in Python as an assignment. Since I noticed some differences between the two languages in 3d interpolation results from my original code, I am trying to address the issue by analysing a simple example.

I set a 2x2x2 matrix (named blocc below) with some values, and its coordinates in three vectors (X,Y,Z). Given a query point, I use 3D-linear interpolation to find the intepolated value. Again,I get different results in MATLAB and Python (code below).

Python

import numpy as np
import scipy.interpolate as si

X,Y,Z =(np.array([1, 2]),np.array([1, 2]),np.array([1, 2]))

a = np.ones((2,2,1))
b = np.ones((2,2,1))*2

blocc = np.concatenate((a,b),axis=2) # Matrix with values

blocc[1,0,0]=7
blocc[0,1,1]=7

qp = np.array([2,1.5,1.5]) #My query point

value=si.interpn((X,Y,Z),blocc,qp,'linear')

print(value)

Here I get value=3

MATLAB

blocc = zeros(2,2,2);
blocc(:,:,1) = ones(2,2);
blocc(:,:,2) = ones(2,2)*2;

blocc(2,1,1)=7;
blocc(1,2,2)=7;

X=[1,2];
Y=[1,2];
Z=[1,2];

qp = [2 1.5 1.5];

value=interp3(X,Y,Z,blocc,qp(1),qp(2),qp(3),'linear')

And here value=2.75

I can't understand why: I think there is something I don't get about how does interpolation and/or matrix indexing work in Python. Can you please make it clear for me? Thanks!

Answer 1

Apparently, for MATLAB when X, Y and Z are vectors, then it considers that the order of the dimensions in the values array is (Y, X, Z). From the documentation:

V — Sample values
array

Sample values, specified as a real or complex array. The size requirements for V depend on the size of X, Y, and Z:

• If X, Y, and Z are arrays representing a full grid (in meshgrid format), then the size of V matches the size of X, Y, or Z .
• If X, Y, and Z are grid vectors, then size(V) = [length(Y) length(X) length(Z)].

If V contains complex numbers, then interp3 interpolates the real and imaginary parts separately.

Example: rand(10,10,10)

Data Types: single | double
Complex Number Support: Yes

This means that, to get the same result in Python, you just need to swap the first and second values in the query:

qp = np.array([1.5, 2, 1.5])
f = si.interpn((X, Y, Z), blocc, qp, 'linear')
print(f)
# [2.75]