Given a colormap sampled at various values, compute splines that
    interpolate in R, G, and B (separately) for fast evaluation of the
    colormap for arbitrary float values. We have primitive support for some
    transformations, though these are generally best done upstream of the
    color mapping code.

    samples   - Shape (4, n). samples[0,:] are the normalized values at
                which the map is sampled, hopefully ranging uniformly
                between 0 and 1. samples[1:4,:] are the RGB values of
                the colormap. (They don't need to actually be RGB, but
                there need to be three of them.)
    transform - One of 'none', 'reverse', or 'sqrt'.
    fitfactor - Sets the tightness of the spline interpolation.

    Returns: a function `map` following `map(n) -> rgb`, where if `n` has
    shape S the result has shape shape (S + (3,)), following a spline
    interpolation from the sampled values.

def approx_colormap (samples, transform='none', fitfactor=1.):
    """Given a colormap sampled at various values, compute splines that
    interpolate in R, G, and B (separately) for fast evaluation of the
    colormap for arbitrary float values. We have primitive support for some
    transformations, though these are generally best done upstream of the
    color mapping code.

    samples   - Shape (4, n). samples[0,:] are the normalized values at
                which the map is sampled, hopefully ranging uniformly
                between 0 and 1. samples[1:4,:] are the RGB values of
                the colormap. (They don't need to actually be RGB, but
                there need to be three of them.)
    transform - One of 'none', 'reverse', or 'sqrt'.
    fitfactor - Sets the tightness of the spline interpolation.

    Returns: a function `map` following `map(n) -> rgb`, where if `n` has
    shape S the result has shape shape (S + (3,)), following a spline
    interpolation from the sampled values.

    """
    import scipy.interpolate as SI

    values = samples[0]
    if transform == 'none':
        pass
    elif transform == 'reverse':
        samples = samples[:,::-1]
    elif transform == 'sqrt':
        values = np.sqrt (values)
    else:
        raise ValueError ('unknown transformation: ' + str (transform))

    nsamp = samples.shape[1]
    rspline = SI.splrep (values, samples[R+1], s=fitfactor/nsamp)
    gspline = SI.splrep (values, samples[G+1], s=fitfactor/nsamp)
    bspline = SI.splrep (values, samples[B+1], s=fitfactor/nsamp)

    def colormap (values):
        values = np.asarray (values)
        mapped = np.empty (values.shape + (3,))

        flatvalues = values.flatten ()
        flatmapped = mapped.reshape (flatvalues.shape + (3,))

        flatmapped[:,R] = SI.splev (flatvalues, rspline)
        flatmapped[:,G] = SI.splev (flatvalues, gspline)
        flatmapped[:,B] = SI.splev (flatvalues, bspline)

        return mapped

    return colormap