import numpy
import scipy.ndimage.interpolation

exposures = [0, 1, 10, 20, 25]
num_samples = 10

values2d = numpy.arange((20), dtype = numpy.float32).reshape(5, 2, 2)
input_values = numpy.array([[10, 11], [12, 13]])

vmin2d, vmax2d = values2d.min(axis = 0), values2d.max(axis = 0)
uniform_values2d = numpy.zeros((num_samples, 2, 2))
uniform_exposures2d = numpy.zeros((num_samples, 2, 2))
for i in xrange(2):
    for j in xrange(2):
        uniform_values2d[:, i, j] = numpy.linspace(vmin2d[i, j],
                vmax2d[i, j], num_samples)
        uniform_exposures2d[:, i, j] = numpy.interp(uniform_values2d[:, i, j],
                values2d[:, i, j], exposures)
indices2d = (num_samples - 1) * (input_values - vmin2d) / (vmax2d - vmin2d)

map_input = numpy.zeros((3, 2, 2))
map_input[:1] = indices2d
map_input[1:] = numpy.indices((2, 2))
exposure_estimate2d = scipy.ndimage.interpolation.map_coordinates(
        uniform_exposures2d, map_input, order = 1, cval = -1)

print "Initial values:\n",values2d
print "\nUniformly-sampled:\n",uniform_values2d
print "\nMapped to exposure-time indices:\n",indices2d
print "\nEstimated exposure times:\n",exposure_estimate2d