Landform Generation using perlin noise and systematic distortions.

landgen.py

import math, numpy, random

def north(matrix):
  return numpy.concatenate((matrix[1:,:], matrix[:1,:]))

def east(matrix):
	return numpy.concatenate((matrix[:,-1:], matrix[:,:-1]), 1)

def south(matrix):
	return numpy.concatenate((matrix[-1:,:], matrix[:-1,:]))

def west(matrix):
	return numpy.concatenate((matrix[:,1:], matrix[:,:1]), 1)

class TerrainChunk:

	def __init__(self, position, shape=(16,16)):
		self.shape = shape
		self.z = numpy.zeros(shape)
		self.sealevel = 0.5
		self.position = position

	def smooth(self):
		"""Smooth on the Moore neighborhood; apply this filter:
		1 2 1
		2 4 2
		1 2 1
		"""
		n = north(self.z)
		s = south(self.z)
		self.z = (self.z * 4 + (n + s + east(self.z) + west(self.z)) * 2 +
							(west(n) + east(n) + west(s) + east(s))) / 16
		return self

	def meadowize(self):
		"""Fill in local minima so that the entire above-sea-level portion of
		the map flows downhill to the sea. Not sure if this is worth anything."""
		W,H = self.z.shape
		def vnn(x,y):
			VNN = (-1,0),(1,0),(0,-1),(0,1)	# Von Neumann neighborhood vectors
			return [((x+dx)%W,(y+dy)%H) for dx,dy in VNN]
		blessing = numpy.ones(self.z.shape) * self.z.max() + 1.0
		blessed = [(x,y) for y in range(H) for x in range(W)
							 if self.z[x,y] <= self.sealevel]
		for x,y in blessed:
			blessing[x,y] = self.z[x,y] # self.sealevel
		todo = set()
		for x,y in blessed:
			for x1,y1 in vnn(x,y):
				if blessing[x1,y1] > self.sealevel:
					todo.add((x1,y1))
		while todo:
			x,y = todo.pop()
			b = max(self.z[x,y], min([blessing[xn,yn] for xn,yn in vnn(x,y)]))
			if b < blessing[x,y]:
				blessing[x,y] = b
				todo |= set([(x1,y1) for (x1,y1) in vnn(x,y)
										 if max(b,self.z[x1,y1]) < blessing[x1,y1]])
		self.z = blessing
		#print blessing
		return self

	def hydro_erode(self, iters, rainfall=0.02, solubility=0.1, evap=0.3):
		def drip(rain, speed=0.5):
			tz = self.z + rain
			NSEW = north,south,east,west
			drop = numpy.array([(tz - x(tz)).clip(0.0, 999.0) for x in NSEW])
			flowing = numpy.minimum(drop.max(0), rain) * speed
			flowto = drop.argmax(0)
			nada = numpy.zeros(self.z.shape)
			result = 0-flowing
			
			for i in range(4):
				dflow = numpy.where(flowto-i, nada, flowing)
				result += (south,north,west,east)[i](dflow)
			return result

		z = self.z + 0
		rain = z * 0
		#rain = z * 0 + rainfall
		for i in range(iters):
			# rainfall (which dissolves when is falls)
			rain += rainfall
			#z -= rainfall * solubility
			# flow
			d = drip(z, rain)
			z += d * solubility
			rain += d
			# evaporate
			#z += evap * solubility * rain
			#rain -= evap * rain
			# final evap
			#z += rain * solubility
		self.z = z
		return self

	def fillPerlinNoise(self, Nrange=None, B=256):
		shape = self.shape

		def perlin_precomps(dimensions, B=256):
			"Build random permutation and gradient vectors"

			def permutation(n):
				result = range(n)
				for i in range(n):
					j = random.randrange(n)
					t = result[j]
					result[j] = result[i]
					result[i] = t
				return result

			def normalize_vector(v):
				s = math.sqrt(sum([a*a for a in v]))
				return [a/s for a in v]

			p = permutation(B)
			g = [normalize_vector([(random.random() * 2.0) - 1.0 for j in range(dimensions)])
					 for i in range(B)]
			# avoid need for modulo op
			p += p + p[:2]	
			g += g + g[:2]
			return p,g

		def perlin_noise2d(point, B, p, g):
			def setup(v):
				b0 = int(v % B)
				b1 = (b0+1) % B
				r0 = v % 1.0
				r1 = r0 - 1.0
				return b0,b1,r0,r1
			((bx0,bx1,rx0,rx1), (by0,by1,ry0,ry1)) = [setup(coord) for coord in point]

			i = p[bx0]
			j = p[bx1]

			b00 = p[i + by0]; b10 = p[j + by0]
			b01 = p[i + by1]; b11 = p[j + by1]

			sx = s_curve(rx0)
			sy = s_curve(ry0)

			u = numpy.dot([rx0, ry0], g[b00])
			v = numpy.dot([rx1, ry0], g[b10])
			a = u + sx * (v - u)	# linear interpolation in x

			u = numpy.dot([rx0, ry1], g[b01])
			v = numpy.dot([rx1, ry1], g[b11])
			b = u + sx * (v - u)	# linear interpolation in x

			return a + sy * (b - a)	# linear interpolation in y


		p,g = perlin_precomps(2, B)

		for n in Nrange or range(1, 6+1):
			R = 2.0 ** n	# could speed up the 0 case by losing interpolation!
			for j in range(self.z.shape[1]):
				for i in range(self.z.shape[0]):
					x = i + self.position[0]
					y = j + self.position[1]
					# put abs() around here for fire effect
					self.z[i,j] += R * perlin_noise2d((x / R, y / R), B, p, g)
		return self

	def reset(self):
		self.z = numpy.zeros(self.z.shape)
		return self

	def normalize(self):
		#if self.z.max() > self.z.min():
		self.z = (self.z - self.z.min()) / (self.z.max() - self.z.min())
		#else:
		#	self.reset()
		return self

	def thermally_erode(self, talus=0.01, erosion=0.5, fill=False):
		nsew = [(self.z - x(self.z)).clip(fill and -999.0 or 0.0, 999.0)
						for x in north,south,east,west]
		nsew = [(x - x.clip(-talus, talus)) for x in nsew]
		self.z = self.z - sum(nsew) * (erosion/4.0)
		return self

	def makepretty(self):
		self.normalize()
		self.hydro_erode(random.randrange(1,10)).normalize()
		self.meadowize().normalize()
		self.smooth().smooth().normalize()
		talus_stuff = (self.z-north(self.z)).clip(0,10)
		mean = talus_stuff.mean()
		std = talus_stuff.mean()
		talus = random.uniform(mean,mean+std)
		self.thermally_erode(talus=talus).normalize()
		return self

	def export(self, size=16, max=(30,10)):
		max_height, min_height = max

		ret = []
		for x in xrange(size):
			for y in xrange(size):
				z = self.z[x,y]
				oz = (z + 1) / 2 * max_height + min_height
				ret.append((x,y,int(oz)))

		return ret


if __name__ == "__main__":	
	terrain = TerrainChunk((0,0))
	terrain.fillPerlinNoise()
	terrain.makepretty()
	terrain = terrain.export()
	f = open('f.txt','w')
	w = []
	for x in range(16):
		for y in range(16):
			w.append(str(int(terrain[x + y][2])))
		f.write(' '.join(w) + "\n")