A quick comment about array shapes would be handy here or in the docstring. If I've got it right:
 - interp_points = (n_samples, n_sample_dims)
 - data = always N-dimensional, where N = cube.ndim
 - result = (n_samples, ...non-sampled dimensions...) ?