Use iterative finite differences for diffusion, instead of gaussian blur

Allows variable diffusion coefficients

generate.py

@@ -1,30 +1,54 @@
 #!/usr/bin/env python3
 
 import numpy as np
-from noise import snoise2
+from noise import snoise2, snoise3
 import os
 import sys
 
 import terrainlib
 
-def noisemap(X, Y, scale=0.01, vscale=1.0, offset=0.0, log=False, **params):
-    # Determine noise offset randomly
-    xbase = np.random.randint(8192)-4096
-    ybase = np.random.randint(8192)-4096
+class noisemap:
+    def __init__(self, X, Y, scale=0.01, vscale=1.0, tscale=1.0, offset=0.0, log=None, xbase=None, ybase=None, **params):
+        # Determine noise offset randomly
+        if xbase is None:
+            xbase = np.random.randint(8192)-4096
+        if ybase is None:
+            ybase = np.random.randint(8192)-4096
+        self.xbase = xbase
+        self.ybase = ybase
+        self.X = X
+        self.Y = Y
+        self.scale = scale
+        if log:
+            vscale /= offset
+        self.vscale = vscale
+        self.tscale = tscale
+        self.offset = offset
+        self.log = log
+        self.params = params
 
-    if log:
-        vscale /= offset
+    def get2d(self):
+        n = np.zeros((self.X, self.Y))
+        for x in range(self.X):
+            for y in range(self.Y):
+                n[x,y] = snoise2(x/self.scale + self.xbase, y/self.scale + self.ybase, **self.params)
 
-    # Generate the noise
-    n = np.zeros((X, Y))
-    for x in range(X):
-        for y in range(Y):
-            n[x,y] = snoise2(x/scale + xbase, y/scale + ybase, **params)
+        if self.log:
+            return np.exp(n*self.vscale) * self.offset
+        else:
+            return n*self.vscale + self.offset
 
-    if log:
-        return np.exp(n*vscale) * offset
-    else:
-        return n*vscale + offset
+    def get3d(self, t=0):
+        t /= self.tscale
+        n = np.zeros((self.X, self.Y))
+        for x in range(self.X):
+            for y in range(self.Y):
+                n[x,y] = snoise3(x/self.scale + self.xbase, y/self.scale + self.ybase, t, **self.params)
+
+        if self.log:
+            return np.exp(n*self.vscale) * self.offset
+        else:
+            return n*self.vscale + self.offset
 
 ### PARSE COMMAND-LINE ARGUMENTS
 argc = len(sys.argv)
@@ -87,6 +111,7 @@ sea_level = float(get_setting('sea_level', 0.0))
 sea_level_variations = float(get_setting('sea_level_variations', 0.0))
 sea_level_variations_time = float(get_setting('sea_level_variations_time', 1.0))
 flex_radius = float(get_setting('flex_radius', 20.0))
+tectonics_time = float(get_setting('tectonics_time', 0.0))
 flow_method = get_setting('flow_method', 'semirandom')
 
 time = float(get_setting('time', 10.0))
@@ -111,17 +136,43 @@ params_sealevel = {
     "lacunarity" : 2,
 }
 
+params_K = {
+    "offset" : K,
+    "vscale" : K,
+    "scale" : 400,
+    "octaves" : 1,
+    "persistence" : 0.5,
+    "lacunarity" : 2,
+    "log" : True,
+}
+
+params_m = {
+    "offset" : m,
+    "vscale" : m*0.2,
+    "scale" : 400,
+    "octaves" : 1,
+    "persistence" : 0.5,
+    "lacunarity" : 2,
+    "log" : False,
+}
+
 if sea_level_variations != 0.0:
     sea_ybase = np.random.randint(8192)-4096
     sea_level_ref = snoise2(time * (1-1/niter) / sea_level_variations, sea_ybase, **params_sealevel) * sea_level_variations
     params['offset'] -= (sea_level_ref + sea_level)
 
-n = noisemap(mapsize+1, mapsize+1, **params)
+if tectonics_time == 0.0:
+    n = noisemap(mapsize+1, mapsize+1, **params).get2d()
+else:
+    terrain_noisemap = noisemap(mapsize+1, mapsize+1, tscale=tectonics_time, **params)
+    n = terrain_noisemap.get3d()
+m_map = noisemap(mapsize+1, mapsize+1, **params_m).get2d()
+K_map = noisemap(mapsize+1, mapsize+1, **params_K).get2d()
 
 ### COMPUTE LANDSCAPE EVOLUTION
 # Initialize landscape evolution model
 print('Initializing model')
-model = terrainlib.EvolutionModel(n, K=K, m=m, d=d, sea_level=sea_level, flex_radius=flex_radius, flow_method=flow_method)
+model = terrainlib.EvolutionModel(n, K=K_map, m=m_map, d=K_map, sea_level=sea_level, flex_radius=flex_radius, flow_method=flow_method)
 terrainlib.update(model.dem, model.lakes, t=5, sea_level=model.sea_level, title='Initializing...')
 
 dt = time/niter
@@ -141,6 +192,9 @@ for i in range(niter):
     terrainlib.update(model.dem, model.lakes, sea_level=model.sea_level, title=disp_niter)
     print('Advection')
     model.advection(dt)
+    if tectonics_time != 0.0:
+        print('Isostasy reference redefinition')
+        model.define_isostasy(terrain_noisemap.get3d((i+1)*dt))
     print('Isostatic equilibration')
     model.adjust_isostasy()
 

terrainlib/erosion.py

@@ -1,5 +1,6 @@
 import numpy as np
 import scipy.ndimage as im
+import scipy.signal as si
 from .rivermapper import flow
 
 def advection(dem, dirs, rivers, time, K=1, m=0.5, sea_level=0):
@@ -47,11 +48,31 @@ def advection(dem, dirs, rivers, time, K=1, m=0.5, sea_level=0):
 
     return dem_new
 
-def diffusion(dem, time, d=1):
-    radius = d * time**.5
-    if radius == 0:
-        return dem
-    return im.gaussian_filter(dem, radius, mode='reflect') # Diffusive erosion is a simple Gaussian blur
+second_derivative_matrix = np.array([
+    [0., 0.25, 0.],
+    [0.25,-1., 0.25],
+    [0., 0.25, 0.],
+])
+
+diff_max = 1.0
+
+def diffusion(dem, time, d=1.0):
+    if isinstance(d, np.ndarray):
+        dmax = d.max()
+    else:
+        dmax = d
+    diff = time * dmax
+    print(diff)
+    niter = int(diff//diff_max) + 1
+    ddiff = d * (time / niter)
+
+    #print('{:d} iterations'.format(niter))
+    for i in range(niter):
+        dem[1:-1,1:-1] += si.convolve2d(dem, second_derivative_matrix, mode='valid') * ddiff
+        #print('iteration {:d}'.format(i+1))
+
+    return dem
+    #return im.gaussian_filter(dem, radius, mode='reflect') # Diffusive erosion is a simple Gaussian blur
 
 class EvolutionModel:
     def __init__(self, dem, K=1, m=0.5, d=1, sea_level=0, flow=False, flex_radius=100, flow_method='semirandom'):
@@ -59,7 +80,10 @@ class EvolutionModel:
         #self.bedrock = dem
         self.K = K
         self.m = m
-        self.d = d
+        if isinstance(d, np.ndarray):
+            self.d = d[1:-1,1:-1]
+        else:
+            self.d = d
         self.sea_level = sea_level
         self.flex_radius = flex_radius
         self.define_isostasy()
@@ -86,8 +110,10 @@ class EvolutionModel:
         self.dem = diffusion(self.dem, time, d=self.d)
         self.flow_uptodate = False
 
-    def define_isostasy(self):
-        self.ref_isostasy = im.gaussian_filter(self.dem, self.flex_radius, mode='reflect') # Define a blurred version of the DEM that will be considered as the reference isostatic elevation.
+    def define_isostasy(self, dem=None):
+        if dem is None:
+            dem = self.dem
+        self.ref_isostasy = im.gaussian_filter(dem, self.flex_radius, mode='reflect') # Define a blurred version of the DEM that will be considered as the reference isostatic elevation.
 
     def adjust_isostasy(self, rate=1):
         isostasy = im.gaussian_filter(self.dem, self.flex_radius, mode='reflect') # Calculate blurred DEM