#!/usr/bin/env python

from __future__ import division
import numpy as np
from scipy import interpolate
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


def noise(nnodes, amplitude):
    x = np.linspace(0, 1, nnodes)
    y = x
    z = amplitude * np.random.uniform(-1, +1, (nnodes, nnodes))
    return interpolate.RectBivariateSpline(x, y, z)


class Noise(object):
    def __init__(self, size, details=1):
        self._f = [noise(2**(size+n), 1/2**n) for n in range(details)]

    def __call__(self, x, y):
        return sum([f(x, y, grid=True) for f in self._f])


class ColorMap(object):
    def __init__(self, colors):
        self._r = interpolate.interp1d(colors[:, 0], colors[:, 1])
        self._g = interpolate.interp1d(colors[:, 0], colors[:, 2])
        self._b = interpolate.interp1d(colors[:, 0], colors[:, 3])

    def __call__(self, z):
        m, n = z.shape
        result = np.empty((m, n, 3))
        result[:, :, 0] = self._r(z)
        result[:, :, 1] = self._g(z)
        result[:, :, 2] = self._b(z)
        return result


f = Noise(3, 4)
x = np.linspace(0, 1, 1024)
y = x
z = f(x, y)
z /= np.max(abs(z))
x, y = np.meshgrid(x, y)

with plt.style.context('ggplot'):
    cmap = np.array([[-1.00, 0.13, 0.25, 0.51],
                     [-0.10, 0.30, 0.51, 0.99],
                     [-0.08, 0.99, 0.87, 0.54],
                     [ 0.00, 0.68, 0.45, 0.30],
                     [ 0.10, 0.54, 0.75, 0.24],
                     [ 0.45, 0.13, 0.51, 0.13],
                     [ 0.65, 0.13, 0.51, 0.13],
                     [ 0.80, 1.00, 1.00, 1.00],
                     [ 1.00, 1.00, 1.00, 1.00]])
    cmap = ColorMap(cmap)

    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    ax.set_zlim(0, +2)
    ax.plot_surface(x, y, np.clip(z, 0, 1), facecolors=cmap(z))
    plt.savefig('src05.png')