#!/usr/bin/env python3

from numpy import *
from PIL import Image, ImageDraw, ImageFont
import random

random.seed("Perlin")
fnt = ImageFont.truetype('Courier Prime.ttf', 100)

XDIM, YDIM = 100, 100     # gradient mesh size (must be big enough for highest octave)
XDP, YDP = 1024, 1024     # Perlin noise array size

phi = zeros((YDIM, XDIM))   # gradient phase angles
perlin = zeros((YDP, XDP))  # Perlin noise array

# choose random phase angles for gradients
for y in range(YDIM):
    for x in range(XDIM):
        phi[y,x] = 2 * pi * random.random()

def grad(x0, y0, dx, dy):
    "Compute the dot product between the gradient and position vector"
    gx, gy = cos(phi[y0,x0]), sin(phi[y0,x0])
    return gx * dx + gy * dy

def inter(a, b, w):
    "Interpolate between a and b using weight w"
    return (b - a) * ((w * (w * 6 - 15) + 10) * w * w * w) + a  # smootherstep

# compute Perlin noise in several octaves and add it together
for nscale, SCALE in enumerate((512, 256, 128, 64, 32, 16)):
    for yi in range(YDP):
        for xi in range(XDP):
            x, y = xi / SCALE, yi / SCALE
            x0, y0 = int(x), int(y)
            dx, dy = x - x0, y - y0

            g1 = grad(x0, y0, dx, dy)
            g2 = grad(x0 + 1, y0, -1 + dx, dy)
            g3 = grad(x0, y0 + 1, dx, -1 + dy)
            g4 = grad(x0 + 1, y0 + 1, -1 + dx, -1 + dy)

            perlin[yi,xi] += (.5**nscale) * inter(inter(g1, g2, dx), inter(g3, g4, dx), dy)

    p2 = (perlin * 128 + 128).astype(uint8)
    im = Image.fromarray(p2, mode = "L")
    d = ImageDraw.Draw(im)
    d.text((0, 0), "%u" % (nscale + 1), font = fnt, fill = 255, align = "left")
    im.save("out%u.png" % nscale)
    print("saved out%u.png" % nscale)

# Create GIF animation with ImageMagick:
# convert -delay 150 out*.png -loop 0 perlin_animation_6_octaves.gif