# Various common functions. from PIL import Image import collections import csv from matplotlib.colors import LightSource, LinearSegmentedColormap import matplotlib.pyplot as plt import numpy as np import scipy as sp import scipy.spatial # Open CSV file as a dict. def read_csv(csv_path): with open(csv_path, 'r') as csv_file: return list(csv.DictReader(csv_file)) # Renormalizes the values of `x` to `bounds` def normalize(x, bounds=(0, 1)): return np.interp(x, (x.min(), x.max()), bounds) # Fourier-based power law noise with frequency bounds. def fbm(shape, p, lower=-np.inf, upper=np.inf): freqs = tuple(np.fft.fftfreq(n, d=1.0 / n) for n in shape) freq_radial = np.hypot(*np.meshgrid(*freqs)) envelope = (np.power(freq_radial, p, where=freq_radial!=0) * (freq_radial > lower) * (freq_radial < upper)) envelope[0][0] = 0.0 phase_noise = np.exp(2j * np.pi * np.random.rand(*shape)) return normalize(np.real(np.fft.ifft2(np.fft.fft2(phase_noise) * envelope))) # Returns each value of `a` with coordinates offset by `offset` (via complex # values). The values at the new coordiantes are the linear interpolation of # neighboring values in `a`. def sample(a, offset): shape = np.array(a.shape) delta = np.array((offset.real, offset.imag)) coords = np.array(np.meshgrid(*map(range, shape))) - delta lower_coords = np.floor(coords).astype(int) upper_coords = lower_coords + 1 coord_offsets = coords - lower_coords lower_coords %= shape[:, np.newaxis, np.newaxis] upper_coords %= shape[:, np.newaxis, np.newaxis] result = lerp(lerp(a[lower_coords[1], lower_coords[0]], a[lower_coords[1], upper_coords[0]], coord_offsets[0]), lerp(a[upper_coords[1], lower_coords[0]], a[upper_coords[1], upper_coords[0]], coord_offsets[0]), coord_offsets[1]) return result # Takes each value of `a` and offsets them by `delta`. Treats each grid point # like a unit square. def displace(a, delta): fns = { -1: lambda x: -x, 0: lambda x: 1 - np.abs(x), 1: lambda x: x, } result = np.zeros_like(a) for dx in range(-1, 2): wx = np.maximum(fns[dx](delta.real), 0.0) for dy in range(-1, 2): wy = np.maximum(fns[dy](delta.imag), 0.0) result += np.roll(np.roll(wx * wy * a, dy, axis=0), dx, axis=1) return result # Returns the gradient of the gaussian blur of `a` encoded as a complex number. def gaussian_gradient(a, sigma=1.0): [fy, fx] = np.meshgrid(*(np.fft.fftfreq(n, 1.0 / n) for n in a.shape)) sigma2 = sigma**2 g = lambda x: ((2 * np.pi * sigma2) ** -0.5) * np.exp(-0.5 * (x / sigma)**2) dg = lambda x: g(x) * (x / sigma2) fa = np.fft.fft2(a) dy = np.fft.ifft2(np.fft.fft2(dg(fy) * g(fx)) * fa).real dx = np.fft.ifft2(np.fft.fft2(g(fy) * dg(fx)) * fa).real return 1j * dx + dy # Simple gradient by taking the diff of each cell's horizontal and vertical # neighbors. def simple_gradient(a): dx = 0.5 * (np.roll(a, 1, axis=0) - np.roll(a, -1, axis=0)) dy = 0.5 * (np.roll(a, 1, axis=1) - np.roll(a, -1, axis=1)) return 1j * dx + dy # Loads the terrain height array (and optionally the land mask from the given # file. def load_from_file(path): result = np.load(path) if type(result) == np.lib.npyio.NpzFile: return (result['height'], result['land_mask']) else: return (result, None) # Saves the array as a PNG image. Assumes all input values are [0, 1] def save_as_png(a, path): image = Image.fromarray(np.round(a * 255).astype('uint8')) image.save(path) # Creates a hillshaded RGB array of heightmap `a`. _TERRAIN_CMAP = LinearSegmentedColormap.from_list('my_terrain', [ (0.00, (0.15, 0.3, 0.15)), (0.25, (0.3, 0.45, 0.3)), (0.50, (0.5, 0.5, 0.35)), (0.80, (0.4, 0.36, 0.33)), (1.00, (1.0, 1.0, 1.0)), ]) def hillshaded(a, land_mask=None, angle=270): if land_mask is None: land_mask = np.ones_like(a) ls = LightSource(azdeg=angle, altdeg=30) land = ls.shade(a, cmap=_TERRAIN_CMAP, vert_exag=10.0, blend_mode='overlay')[:, :, :3] water = np.tile((0.25, 0.35, 0.55), a.shape + (1,)) return lerp(water, land, land_mask[:, :, np.newaxis]) # Linear interpolation of `x` to `y` with respect to `a` def lerp(x, y, a): return (1.0 - a) * x + a * y # Returns a list of grid coordinates for every (x, y) position bounded by # `shape` def make_grid_points(shape): [Y, X] = np.meshgrid(np.arange(shape[0]), np.arange(shape[1])) grid_points = np.column_stack([X.flatten(), Y.flatten()]) return grid_points # Returns a list of points sampled within the bounds of `shape` and with a # minimum spacing of `radius`. # NOTE: This function is fairly slow, given that it is implemented with almost # no array operations. def poisson_disc_sampling(shape, radius, retries=16): grid = {} points = [] # The bounds of `shape` are divided into a grid of cells, each of which can # contain a maximum of one point. cell_size = radius / np.sqrt(2) cells = np.ceil(np.divide(shape, cell_size)).astype(int) offsets = [(0, 0), (0, -1), (0, 1), (-1, 0), (1, 0), (-1, -1), (-1, 1), (1, -1), (1, 1), (-2, 0), (2, 0), (0, -2), (0, 2)] to_cell = lambda p: (p / cell_size).astype('int') # Returns true if there is a point within `radius` of `p`. def has_neighbors_in_radius(p): cell = to_cell(p) for offset in offsets: cell_neighbor = (cell[0] + offset[0], cell[1] + offset[1]) if cell_neighbor in grid: p2 = grid[cell_neighbor] diff = np.subtract(p2, p) if np.dot(diff, diff) <= radius * radius: return True return False # Adds point `p` to the cell grid. def add_point(p): grid[tuple(to_cell(p))] = p q.append(p) points.append(p) q = collections.deque() first = shape * np.random.rand(2) add_point(first) while len(q) > 0: point = q.pop() # Make `retries` attemps to find a point within [radius, 2 * radius] from # `point`. for _ in range(retries): diff = 2 * radius * (2 * np.random.rand(2) - 1) r2 = np.dot(diff, diff) new_point = diff + point if (new_point[0] >= 0 and new_point[0] < shape[0] and new_point[1] >= 0 and new_point[1] < shape[1] and not has_neighbors_in_radius(new_point) and r2 > radius * radius and r2 < 4 * radius * radius): add_point(new_point) num_points = len(points) # Return points list as a numpy array. return np.concatenate(points).reshape((num_points, 2)) # Returns an array in which all True values of `mask` contain the distance to # the nearest False value. def dist_to_mask(mask): border_mask = (np.maximum.reduce([ np.roll(mask, 1, axis=0), np.roll(mask, -1, axis=0), np.roll(mask, -1, axis=1), np.roll(mask, 1, axis=1)]) * (1 - mask)) border_points = np.column_stack(np.where(border_mask > 0)) kdtree = sp.spatial.cKDTree(border_points) grid_points = make_grid_points(mask.shape) return kdtree.query(grid_points)[0].reshape(mask.shape) # Generates worley noise with points separated by `spacing`. def worley(shape, spacing): points = poisson_disc_sampling(shape, spacing) coords = np.floor(points).astype(int) mask = np.zeros(shape, dtype=bool) mask[coords[:, 0], coords[:, 1]] = True return normalize(dist_to_mask(mask)) # Peforms a gaussian blur of `a`. def gaussian_blur(a, sigma=1.0): freqs = tuple(np.fft.fftfreq(n, d=1.0 / n) for n in a.shape) freq_radial = np.hypot(*np.meshgrid(*freqs)) sigma2 = sigma**2 g = lambda x: ((2 * np.pi * sigma2) ** -0.5) * np.exp(-0.5 * (x / sigma)**2) kernel = g(freq_radial) kernel /= kernel.sum() return np.fft.ifft2(np.fft.fft2(a) * np.fft.fft2(kernel)).real