#!/usr/bin/env python

################################################################################
# A landscape generator which employs some basic techniques for generating     #
# noise and interpolating images.                                              #
# Author: David J. Oftedal.                                                    #
#                                                                              #
# Depends on mlab: http://github.enthought.com/mayavi/mayavi/                  #
# and NumPy:       http://numpy.org/                                           #
#                                                                              #
# Special thanks to the authors of the following:                              #           
# http://www.decarpentier.nl/scape-procedural-basics                           #
# http://freespace.virgin.net/hugo.elias/models/m_perlin.htm                   #
################################################################################

from math import *
from numpy import *
from mayavi.mlab import *

# Linear interpolation between two points.
def linear(p1, p2, x):
	return p1*(1-x) + p2*x

# Cosine interpolation between two points.
def cosine(p1, p2, x):
	x = (1-cos(x * pi))/2.0
	return p1*(1-x) + p2*x

# Interpolate a matrix in two dimensions.
# The edges of the new matrix include the first and last element exactly as they
# were in the original matrix; this is especially simple when indexing from 0.
def interpolate2d(old, newx, newy, interpolate):

	new = zeros((newy, len(old[0])))

	# Interpolate the rows first.
	oldmaxy = float(len(old)) -1.0
	newmaxy = float(newy) -1.0

	for y in range(newy):
		try:
			j = y / (newmaxy/oldmaxy)
		except:
			# Debug code to investigate division by 0.
			# print str(y) + "\t" + str(newmaxy) + "/" + str(oldmaxy)
			# When oldmaxy is 0, j should be 0.
			j = 0
		k = max(int(floor(j)),0) # Just to be on the safe side
		l = min(int(ceil(j)),len(old)-1)
		for x in range(len(old[k])):
			val = interpolate(old[k][x], old[l][x], (j-k))
			#	Debug code to test for off-by-one/floating-point rounding errors:
			#	print str(j) + "->" + str(i) + ": old[" + str(k) + "][" + str(m)+ "], old[" + str(l)+"][" +str(m)+ "]. max = old[" + str(len(old)-1)+"][" +str(len(old[0])-1)+ "]."
			new[y][x] = val

	old = new
	new = zeros((newy, newx))

	# Then interpolate the columns.
	oldmaxx = float(len(old[0])) -1.0
	newmaxx = float(newx) -1.0

	for y in range(newy):
		for x in range(newx):
			try:
				j = x / (newmaxx/oldmaxx)
			except:
				# print str(x) + "\t" + str(newmaxx) + "/" + str(oldmaxx)
				j = 0
			k = max(int(floor(j)),0) # Just to be on the safe side
			l = min(int(ceil(j)),len(old)-1)
			val = interpolate(old[y][k], old[y][l], (j-k))
			new[y][x] = val	

	return new

# Generate a neat sequence of the numbers from 1-25.
def testinterpolation():
	return interpolate2d([[1, 5],[21,25]], 5, 5, linear)

# Generate some random noise.
def noise(x, y, scaledx, scaledy, minval, maxval):
	tmp = random.random_integers(minval, maxval, (y, x))
	return interpolate2d(tmp, scaledx, scaledy, cosine)

# Generate landscape-like noise by adding together noise at different scales.
def genlandscape(x, y, height, divscale, divheight):
	tmpx = 1
	tmpy = 1

	landscape = zeros((y, x))
	
	while tmpx <= x and tmpy <= y and height > 0:
		landscape += noise(int(tmpx), int(tmpy), x, y, -height, height)
		height /= divheight
		tmpx *= divscale
		tmpy *= divscale

	return landscape

# Visualize a random landscape with mlab.
a = array(genlandscape(800,600,400,1.97,2))
s = surf(a)
show()