ant landscape: update/various fixes

10 files changed, 475 insertions, 539 deletions

diff --git a/ant_landscape/ErosionR.txt b/ant_landscape/ErosionR.txt
index 1c461945..1dae7aa0 100644
--- a/ant_landscape/ErosionR.txt
+++ b/ant_landscape/ErosionR.txt
@@ -1,3 +1,2 @@
 http://blog.michelanders.nl/search/label/erosion
-https://github.com/nerk987/ErosionR
-https://blenderartists.org/forum/showthread.php?323808-Simulating-erosion-in-Blender
\ No newline at end of file
+https://github.com/nerk987/ErosionR
\ No newline at end of file
diff --git a/ant_landscape/__init__.py b/ant_landscape/__init__.py
index d0f1528d..de916966 100644
--- a/ant_landscape/__init__.py
+++ b/ant_landscape/__init__.py
@@ -17,7 +17,7 @@
 # ##### END GPL LICENSE BLOCK #####
 
 # Another Noise Tool - Suite (W.I.P.)
-# Jim Hazevoet 5/2017
+# Jimmy Hazevoet 5/2017
 
 bl_info = {
     "name": "A.N.T.Landscape",
@@ -67,7 +67,9 @@ from .ant_functions import (
 # Menu's and panels
 
 def menu_func_eroder(self, context):
-    self.layout.operator('mesh.eroder', text="Landscape Eroder", icon='SMOOTHCURVE')
+    ob = bpy.context.active_object
+    if ob and (ob.ant_landscape.keys() and not ob.ant_landscape['sphere_mesh']):
+        self.layout.operator('mesh.eroder', text="Landscape Eroder", icon='SMOOTHCURVE')
 
 
 def menu_func_landscape(self, context):
@@ -107,8 +109,9 @@ class AntLandscapeToolsPanel(bpy.types.Panel):
         ob = context.active_object
         col = layout.column()
         col.operator('mesh.ant_displace', text="Mesh Displace", icon="RNDCURVE")
-        col.operator('mesh.eroder', text="Landscape Eroder", icon='SMOOTHCURVE')
         col.operator('mesh.ant_slope_map', icon='GROUP_VERTEX')
+        if ob.ant_landscape.keys() and not ob.ant_landscape['sphere_mesh']:
+            col.operator('mesh.eroder', text="Landscape Eroder", icon='SMOOTHCURVE')
 
 
 # Landscape Settings / Properties:
@@ -299,6 +302,12 @@ class AntNoiseSettingsPanel(bpy.types.Panel):
             col.prop(ant, "gain")
             col.separator()
             col.prop(ant, "vl_basis_type")
+        elif ant.noise_type == "rocks_noise":
+            col.prop(ant, "noise_depth")
+            col.prop(ant, "distortion")
+            col.separator()
+            row = col.row(align=True)
+            row.prop(ant, "hard_noise", expand=True)
         elif ant.noise_type == "slick_rock":
             col.prop(ant, "noise_depth")
             col.prop(ant, "dimension")
@@ -514,9 +523,10 @@ class AntLandscapePropertiesGroup(bpy.types.PropertyGroup):
                 ('vl_hTerrain', "vlNoise hTerrain", "A.N.T: vlNoise hTerrain", 12),
                 ('distorted_heteroTerrain', "Distorted hTerrain", "A.N.T distorted hTerrain", 13),
                 ('double_multiFractal', "Double MultiFractal", "A.N.T: double multiFractal", 14),
-                ('slick_rock', "Slick Rock", "A.N.T: slick rock", 15),
-                ('planet_noise', "Planet Noise", "Planet Noise by: Farsthary", 16),
-                ('blender_texture', "Blender Texture - Texture Nodes", "Blender texture data block", 17)]
+                ('rocks_noise', "Noise Rocks", "A.N.T: turbulence variation", 15),
+                ('slick_rock', "Slick Rock", "A.N.T: slick rock", 16),
+                ('planet_noise', "Planet Noise", "Planet Noise by: Farsthary", 17),
+                ('blender_texture', "Blender Texture - Texture Nodes", "Blender texture data block", 18)]
             )
     basis_type = EnumProperty(
             name="Noise Basis",
@@ -735,8 +745,8 @@ class AntLandscapePropertiesGroup(bpy.types.PropertyGroup):
                 ("1", "Smooth", "Smooth transitions", 1),
                 ("2", "Sharp Sub", "Sharp substract transitions", 2),
                 ("3", "Sharp Add", "Sharp add transitions", 3),
-                ("4", "Posterize", "Posterize", 4),
-                ("5", "Posterize Mix", "Posterize mixed", 5)]
+                ("4", "Quantize", "Quantize", 4),
+                ("5", "Quantize Mix", "Quantize mixed", 5)]
             )
     water_plane = BoolProperty(
             name="Water Plane",
diff --git a/ant_landscape/add_mesh_ant_landscape.py b/ant_landscape/add_mesh_ant_landscape.py
index c7de90f3..7fc4a268 100644
--- a/ant_landscape/add_mesh_ant_landscape.py
+++ b/ant_landscape/add_mesh_ant_landscape.py
@@ -17,7 +17,7 @@
 # ##### END GPL LICENSE BLOCK #####
 
 # Another Noise Tool - Landscape
-# Jim Hazevoet
+# Jimmy Hazevoet
 
 # import modules
 import bpy
@@ -47,7 +47,7 @@ from .ant_functions import (
 class AntAddLandscape(bpy.types.Operator):
     bl_idname = "mesh.landscape_add"
     bl_label = "Another Noise Tool - Landscape"
-    bl_description = "A.N.T. Add landscape mesh"
+    bl_description = "Add landscape mesh"
     bl_options = {'REGISTER', 'UNDO', 'PRESET'}
 
     ant_terrain_name = StringProperty(
@@ -191,9 +191,10 @@ class AntAddLandscape(bpy.types.Operator):
                 ('vl_hTerrain', "vlNoise hTerrain", "A.N.T: vlNoise hTerrain", 12),
                 ('distorted_heteroTerrain', "Distorted hTerrain", "A.N.T distorted hTerrain", 13),
                 ('double_multiFractal', "Double MultiFractal", "A.N.T: double multiFractal", 14),
-                ('slick_rock', "Slick Rock", "A.N.T: slick rock", 15),
-                ('planet_noise', "Planet Noise", "Planet Noise by: Farsthary", 16),
-                ('blender_texture', "Blender Texture - Texture Nodes", "Blender texture data block", 17)]
+                ('rocks_noise', "Noise Rocks", "A.N.T: turbulence variation", 15),
+                ('slick_rock', "Slick Rock", "A.N.T: slick rock", 16),
+                ('planet_noise', "Planet Noise", "Planet Noise by: Farsthary", 17),
+                ('blender_texture', "Blender Texture - Texture Nodes", "Blender texture data block", 18)]
             )
     basis_type = EnumProperty(
             name="Noise Basis",
@@ -412,8 +413,8 @@ class AntAddLandscape(bpy.types.Operator):
                 ("1", "Smooth", "Smooth transitions", 1),
                 ("2", "Sharp Sub", "Sharp substract transitions", 2),
                 ("3", "Sharp Add", "Sharp add transitions", 3),
-                ("4", "Posterize", "Posterize", 4),
-                ("5", "Posterize Mix", "Posterize mixed", 5)]
+                ("4", "Quantize", "Quantize", 4),
+                ("5", "Quantize Mix", "Quantize mixed", 5)]
             )
     water_plane = BoolProperty(
             name="Water Plane",
diff --git a/ant_landscape/ant_functions.py b/ant_landscape/ant_functions.py
index 5ccbeade..ce7751af 100644
--- a/ant_landscape/ant_functions.py
+++ b/ant_landscape/ant_functions.py
@@ -17,7 +17,7 @@
 # ##### END GPL LICENSE BLOCK #####
 
 # Another Noise Tool - Functions
-# Jim Hazevoet
+# Jimmy Hazevoet
 
 # ErosionR:
 # Michel Anders (varkenvarken), Ian Huish (nerk)
@@ -60,17 +60,11 @@ from bpy_extras import object_utils
 def create_mesh_object(context, verts, edges, faces, name):
     # Create new mesh
     mesh = bpy.data.meshes.new(name)
-
     # Make a mesh from a list of verts/edges/faces.
     mesh.from_pydata(verts, [], faces)
-
     # Update mesh geometry after adding stuff.
     mesh.update()
-
-    #new_ob = bpy.data.objects.new(name, mesh)
-    #context.scene.objects.link(new_ob)
     return object_utils.object_data_add(context, mesh, operator=None)
-    #return new_ob
 
 
 # Generate XY Grid
@@ -172,12 +166,7 @@ class AntVgSlopeMap(bpy.types.Operator):
             default=0.0,
             min=0.0,
             max=1.0,
-            description="Increase to select more vertices"
-            )
-    weight_mode = BoolProperty(
-            name="Enter WeightPaint Mode:",
-            default=True,
-            description="Enter weightpaint mode when done"
+            description="Increase to select more vertices on slopes"
             )
 
     @classmethod
@@ -192,8 +181,8 @@ class AntVgSlopeMap(bpy.types.Operator):
 
 
     def execute(self, context):
-        message = "Popup Values: %d, %f, %s, %s, %s" % \
-            (self.select_flat, self.select_range, self.group_name, self.z_method, self.weight_mode)
+        message = "Popup Values: %d, %f, %s, %s" % \
+            (self.select_flat, self.select_range, self.group_name, self.z_method)
         self.report({'INFO'}, message)
 
         bpy.ops.object.mode_set(mode='OBJECT')
@@ -223,8 +212,7 @@ class AntVgSlopeMap(bpy.types.Operator):
 
         vg_normal.name = self.group_name
 
-        if self.weight_mode:
-            bpy.ops.paint.weight_paint_toggle()
+        bpy.ops.paint.weight_paint_toggle()
         return {'FINISHED'}
 
 
@@ -395,9 +383,19 @@ def vl_hTerrain(coords, H, lacunarity, octaves, offset, basis, vlbasis, distort)
 # another turbulence
 def ant_turbulence(coords, depth, hardnoise, nbasis, amp, freq, distortion):
     x, y, z = coords
-    tv = turbulence_vector((x + 1, y + 2, z + 3), depth, hardnoise, nbasis, amp, freq)
-    d = (distortion * tv[0]) * 0.25
-    return (d + ((tv[0] - tv[1]) * (tv[2])**2))
+    t = turbulence_vector((x/2, y/2, z/2), depth, 0, nbasis, amp, freq) * 0.5 * distortion
+    return turbulence((t[0], t[1], t[2]), 2, hardnoise, 3) * 0.5 + 0.5
+
+
+# rocks noise
+def rocks_noise(coords, depth, hardnoise, nbasis, distortion):
+    x,y,z = coords
+    p = turbulence((x, y, z), 4, 0, 0) * 0.125 * distortion
+    xx, yy, zz = x, y, z
+    a = turbulence((xx + p, yy + p, zz), 2, 0, 7)
+    pa = a * 0.1875 * distortion
+    b = turbulence((x, y, z + pa), depth, hardnoise, nbasis)
+    return ((a + 0.5 * (b - a)) * 0.5 + 0.5)
 
 
 # shattered_hterrain:
@@ -570,13 +568,16 @@ def noise_gen(coords, props):
     elif ntype in [14, 'double_multiFractal']:
         value = double_multiFractal(ncoords, dimension, lacunarity, depth, offset, gain, nbasis, vlbasis)
 
-    elif ntype in [15, 'slick_rock']:
+    elif ntype in [15, 'rocks_noise']:
+        value = rocks_noise(ncoords, depth, hardnoise, nbasis, distortion)
+
+    elif ntype in [16, 'slick_rock']:
         value = slick_rock(ncoords,dimension, lacunarity, depth, offset, gain, distortion, nbasis, vlbasis)
 
-    elif ntype in [16, 'planet_noise']:
+    elif ntype in [17, 'planet_noise']:
         value = planet_noise(ncoords, depth, hardnoise, nbasis)[2] * 0.5 + 0.5
 
-    elif ntype in [17, 'blender_texture']:
+    elif ntype in [18, 'blender_texture']:
         if texture_name != "" and texture_name in bpy.data.textures:
             value = bpy.data.textures[texture_name].evaluate(ncoords)[3]
         else:
@@ -819,6 +820,12 @@ def draw_ant_noise(self, context, generate=True):
             col.prop(self, "gain")
             col.separator()
             col.prop(self, "vl_basis_type")
+        elif self.noise_type == "rocks_noise":
+            col.prop(self, "noise_depth")
+            col.prop(self, "distortion")
+            col.separator()
+            row = col.row(align=True)
+            row.prop(self, "hard_noise", expand=True)
         elif self.noise_type == "slick_rock":
             col.prop(self, "noise_depth")
             col.prop(self, "dimension")
@@ -840,9 +847,10 @@ def draw_ant_displace(self, context, generate=True):
     box = layout.box()
     box.prop(self, "show_displace_settings", toggle=True)
     if self.show_displace_settings:
-        col = box.column(align=False)
         if not generate:
+            col = box.column(align=False)
             col.prop(self, "direction", toggle=True)
+
         col = box.column(align=True)
         row = col.row(align=True).split(0.92, align=True)
         row.prop(self, "height")
@@ -961,7 +969,7 @@ def availableVertexGroupsOrNone(self, context):
 class Eroder(bpy.types.Operator):
     bl_idname = "mesh.eroder"
     bl_label = "ErosionR"
-    bl_description = "Apply various kinds of erosion to a landscape mesh. Also available in Weight Paint mode > Weights menu"
+    bl_description = "Apply various kinds of erosion to a square ANT-Landscape grid. Also available in Weight Paint mode > Weights menu"
     bl_options = {'REGISTER', 'UNDO', 'PRESET'}
 
     Iterations = IntProperty(
@@ -1022,7 +1030,8 @@ class Eroder(bpy.types.Operator):
             description="Total Rain amount",
             default=.01,
             min=0,
-            soft_max=1
+            soft_max=1,
+            precision=3
             )
     Kv = FloatProperty(
             name="Rain variance",
@@ -1131,9 +1140,7 @@ class Eroder(bpy.types.Operator):
     def execute(self, context):
 
         ob = context.active_object
-        #obwater = bpy.data.objects["water"]
         me = ob.data
-        #mewater = obwater.data
         self.stats.reset()
         try:
             vgActive = ob.vertex_groups.active.name
@@ -1180,14 +1187,14 @@ class Eroder(bpy.types.Operator):
             vgcapacity=ob.vertex_groups["capacity"]
         except:
             vgcapacity=ob.vertex_groups.new("capacity")
+
         g = Grid.fromBlenderMesh(me, vg, self.Ef)
 
         me = bpy.data.meshes.new(me.name)
-        #mewater = bpy.data.meshes.new(mewater.name)
 
-        self.counts['diffuse']=0
-        self.counts['avalanche']=0
-        self.counts['water']=0
+        self.counts['diffuse'] = 0
+        self.counts['avalanche'] = 0
+        self.counts['water'] = 0
         for i in range(self.Iterations):
             if self.IterRiver > 0:
                 for i in range(self.IterRiver):
@@ -1197,7 +1204,7 @@ class Eroder(bpy.types.Operator):
                 for k in range(self.IterDiffuse):
                     g.diffuse(self.Kd / 5, self.IterDiffuse, self.numexpr)
                     self.counts['diffuse']+=1
-            #if self.Kt < radians(90) and rand() < self.Pa:
+
             if self.Kt < radians(90) and self.Pa > 0:
                 for k in range(self.IterAva):
                     # since dx and dy are scaled to 1, tan(Kt) is the height for a given angle
@@ -1209,8 +1216,7 @@ class Eroder(bpy.types.Operator):
 
         g.toBlenderMesh(me)
         ob.data = me
-        #g.toWaterMesh(mewater)
-        #obwater.data = mewater
+
         if vg:
             for row in range(g.rainmap.shape[0]):
                 for col in range(g.rainmap.shape[1]):
@@ -1231,12 +1237,10 @@ class Eroder(bpy.types.Operator):
                 for col in range(g.rainmap.shape[1]):
                     i = row * g.rainmap.shape[1] + col
                     vgw.add([i],g.water[row,col]/g.watermax,'ADD')
-                    # vgw.add([i],g.water[row,col],'ADD')
         if vgscour:
             for row in range(g.rainmap.shape[0]):
                 for col in range(g.rainmap.shape[1]):
                     i = row * g.rainmap.shape[1] + col
-                    # vgscour.add([i],(g.scour[row,col]-g.scourmin)/(g.scourmax-g.scourmin),'ADD')
                     vgscour.add([i],g.scour[row,col]/max(g.scourmax, -g.scourmin),'ADD')
         if vgdeposit:
             for row in range(g.rainmap.shape[0]):
@@ -1247,13 +1251,11 @@ class Eroder(bpy.types.Operator):
             for row in range(g.rainmap.shape[0]):
                 for col in range(g.rainmap.shape[1]):
                     i = row * g.rainmap.shape[1] + col
-                    # vgflowrate.add([i],g.flowrate[row,col]/g.flowratemax,'ADD')
                     vgflowrate.add([i],g.flowrate[row,col],'ADD')
         if vgsediment:
             for row in range(g.rainmap.shape[0]):
                 for col in range(g.rainmap.shape[1]):
                     i = row * g.rainmap.shape[1] + col
-                    # vgsediment.add([i],g.sediment[row,col]/g.sedmax,'ADD')
                     vgsediment.add([i],g.sediment[row,col],'ADD')
         if vgsedimentpct:
             for row in range(g.rainmap.shape[0]):
@@ -1280,6 +1282,7 @@ class Eroder(bpy.types.Operator):
 
         return {'FINISHED'}
 
+
     def draw(self,context):
         layout = self.layout
 
@@ -1311,45 +1314,7 @@ class Eroder(bpy.types.Operator):
         col.prop(self, 'Kr')
         col.prop(self, 'Kv')
         col.prop(self, 'Kev')
-        #box2 = box.box()
-        #box2.prop(self, 'userainmap')
-        #box2.enabled = context.active_object.vertex_groups.active is not None
-        #box.prop(self, 'Ka')
-        col.prop(self, 'Ef')
 
-        #box = layout.box()
-        #box.label("Probabilities")
-        #box.prop(self, 'Pa')
-        #box.prop(self, 'Pw')
+        col.prop(self, 'Ef')
 
         layout.prop(self,'smooth')
-
-        #if numexpr_available:
-        #  layout.prop(self, 'numexpr')
-        #else:
-        #  box = layout.box()
-        #  box.alert=True
-        #  box.label("Numexpr not available. Will slow down large meshes")
-
-        #box = layout.box()
-        #box.prop(self,'showiterstats')
-        #if self.showiterstats:
-        #    row = box.row()
-        #    col1 = row.column()
-        #    col2 = row.column()
-        #    col1.label("Time"); col2.label("%.1f s"%self.stats.elapsedtime)
-        #    if self.stats.memstats_available:
-        #        col1.label("Memory"); col2.label("%.1f Mb"%(self.stats.maxmem/(1024.0*1024.0)))
-        #    col1.label("Diffusions"); col2.label("%d"% self.counts['diffuse'])
-        #    col1.label("Avalanches"); col2.label("%d"% self.counts['avalanche'])
-        #    col1.label("Water movements"); col2.label("%d"% self.counts['water'])
-        #box = layout.box()
-        #box.prop(self,'showmeshstats')
-        #if self.showmeshstats:
-        #    row = box.row()
-        #    col1 = row.column()
-        #    col2 = row.column()
-        #    for line in self.stats.meshstats.split('\n'):
-        #        label, value = line.split(':')
-        #        col1.label(label)
-        #        col2.label(value)
diff --git a/ant_landscape/ant_landscape_refresh.py b/ant_landscape/ant_landscape_refresh.py
index 010bbd1f..e3ff5366 100644
--- a/ant_landscape/ant_landscape_refresh.py
+++ b/ant_landscape/ant_landscape_refresh.py
@@ -17,7 +17,7 @@
 # ##### END GPL LICENSE BLOCK #####
 
 # Another Noise Tool - Landscape  Redraw - Regenerate
-# Jim Hazevoet
+# Jimmy Hazevoet
 
 
 # ------------------------------------------------------------
@@ -74,9 +74,12 @@ class AntLandscapeRefresh(bpy.types.Operator):
             if ob['use_vgroup']:
                 vertex_group = obj.vertex_groups.active
                 if vertex_group:
+                    gi = vertex_group.index
                     for v in mesh.vertices:
-                        v.co[2] = 0
-                        v.co[2] = vertex_group.weight(v.index) * noise_gen(v.co, prop)
+                        for g in v.groups:
+                            if g.group == gi:
+                                v.co[2] = 0.0
+                                v.co[2] = vertex_group.weight(v.index) * noise_gen(v.co, prop)
             else:
                 for v in mesh.vertices:
                     v.co[2] = 0
diff --git a/ant_landscape/eroder.py b/ant_landscape/eroder.py
index d36aef67..04e04bd7 100644
--- a/ant_landscape/eroder.py
+++ b/ant_landscape/eroder.py
@@ -2,7 +2,7 @@
 #
 #  erode.py  -- a script to simulate erosion of height fields
 #  (c) 2014 Michel J. Anders (varkenvarken)
-#  now with some modifications by Ian Huish (nerk)
+#  with some modifications by Ian Huish (nerk)
 #
 #  This program is free software; you can redistribute it and/or
 #  modify it under the terms of the GNU General Public License
@@ -20,37 +20,29 @@
 #
 # ##### END GPL LICENSE BLOCK #####
 
+
 from time import time
 import unittest
 import sys
 import os
-# import resource     # so much for platform independence, this only works on unix :-(
 from random import random as rand, shuffle
 import numpy as np
-#from .perlin import pnoise
 
 numexpr_available = False
-# Sorry, nerk can't handle numexpr at this time!
-#try:
-#    import numexpr as ne
-#    numexpr_available = True
-#except ImportError:
-#    pass
+
 
 def getmemsize():
-  return 0.0
-  #return resource.getrusage(resource.RUSAGE_SELF).ru_maxrss*resource.getpagesize()/(1024.0*1024.0)
-  
+    return 0.0
+
+
 def getptime():
-  #r = resource.getrusage(resource.RUSAGE_SELF)
-  #return r.ru_utime + r.ru_stime
-  return time()
-      
+    return time()
+
+
 class Grid:
 
     def __init__(self, size=10, dtype=np.single):
-        self.center = np.zeros([size,size], dtype)
-        #print("Centre\n", np.array_str(self.center,precision=3), file=sys.stderr)
+        self.center = np.zeros([size, size], dtype)
         self.water = None
         self.sediment = None
         self.scour = None
@@ -59,19 +51,20 @@ class Grid:
         self.sedimentpct = None
         self.capacity = None
         self.avalanced = None
-        self.minx=None
-        self.miny=None
-        self.maxx=None
-        self.maxy=None
-        self.zscale=1
-        self.maxrss=0.0
-        self.sequence=[0,1,2,3]
+        self.minx = None
+        self.miny = None
+        self.maxx = None
+        self.maxy = None
+        self.zscale = 1
+        self.maxrss = 0.0
+        self.sequence = [0, 1, 2, 3]
         self.watermax = 1.0
         self.flowratemax = 1.0
         self.scourmax = 1.0
         self.sedmax = 1.0
         self.scourmin = 1.0
- 
+
+
     def init_water_and_sediment(self):
         if self.water is None:
             self.water = np.zeros(self.center.shape, dtype=np.single)
@@ -88,9 +81,11 @@ class Grid:
         if self.avalanced is None:
             self.avalanced = np.zeros(self.center.shape, dtype=np.single)
 
+
     def __str__(self):
         return ''.join(self.__str_iter__(fmt="%.3f"))
-        
+
+
     def __str_iter__(self, fmt):
         for row in self.center[::]:
             values=[]
@@ -98,20 +93,24 @@ class Grid:
                 values.append(fmt%v)
             yield  ' '.join(values) + '\n'
 
+
     @staticmethod
     def fromFile(filename):
-        if filename == '-' : filename = sys.stdin
+        if filename == '-':
+            filename = sys.stdin
         g=Grid()
         g.center=np.loadtxt(filename,np.single)
         return g
 
+
     def toFile(self, filename, fmt="%.3f"):
         if filename == '-' : 
             filename = sys.stdout.fileno()
         with open(filename,"w") as f:
             for line in self.__str_iter__(fmt):
                 f.write(line)
-    
+
+
     def raw(self,format="%.3f"):
         fstr=format+" "+ format+" "+ format+" "
         a=self.center / self.zscale
@@ -134,6 +133,7 @@ class Grid:
                        fstr%(row1 ,col0 ,a[row+1][col  ])+
                        fstr%(row1 ,col1 ,a[row+1][col+1])+"\n")
 
+
     def toRaw(self, filename, infomap=None):
         with open(filename if type(filename) == str else sys.stdout.fileno() , "w") as f:
             f.writelines(self.raw())
@@ -141,85 +141,81 @@ class Grid:
             with open(os.path.splitext(filename)[0]+".inf" if type(filename) == str else sys.stdout.fileno() , "w") as f:
                 f.writelines("\n".join("%-15s: %s"%t for t in sorted(infomap.items())))
 
+
     @staticmethod
     def fromRaw(filename):
         """initialize a grid from a Blender .raw file.
         currenly suports just rectangular grids of all triangles
         """
-        g=Grid.fromFile(filename)
-        #  we assume tris and an axis aligned grid
-        g.center=np.reshape(g.center,(-1,3))
+        g = Grid.fromFile(filename)
+        # we assume tris and an axis aligned grid
+        g.center = np.reshape(g.center,(-1,3))
         g._sort()
         return g
 
+
     def _sort(self, expfact):
-        #  keep unique vertices only by creating a set and sort first on x then on y coordinate
-        #  using rather slow python sort but couldn;t wrap my head around np.lexsort
+        # keep unique vertices only by creating a set and sort first on x then on y coordinate
+        # using rather slow python sort but couldn;t wrap my head around np.lexsort
         verts = sorted(list({ tuple(t) for t in self.center[::] }))
-        x=set(c[0] for c in verts)
-        y=set(c[1] for c in verts)
-        nx=len(x)
-        ny=len(y)
-        self.minx=min(x)
-        self.maxx=max(x)
-        self.miny=min(y)
-        self.maxy=max(y)
-        xscale=(self.maxx-self.minx)/(nx-1)
-        yscale=(self.maxy-self.miny)/(ny-1)
+        x = set(c[0] for c in verts)
+        y = set(c[1] for c in verts)
+        nx = len(x)
+        ny = len(y)
+        self.minx = min(x)
+        self.maxx = max(x)
+        self.miny = min(y)
+        self.maxy = max(y)
+        xscale = (self.maxx-self.minx)/(nx-1)
+        yscale = (self.maxy-self.miny)/(ny-1)
         # note: a purely flat plane cannot be scaled 
-        if (yscale != 0.0) and (abs(xscale/yscale) - 1.0 > 1e-3) : raise ValueError("Mesh spacing not square %d x %d  %.4f x %4.f"%(nx,ny,xscale,yscale))
-        self.zscale=1.0
+        if (yscale != 0.0) and (abs(xscale/yscale) - 1.0 > 1e-3):
+            raise ValueError("Mesh spacing not square %d x %d  %.4f x %4.f"%(nx,ny,xscale,yscale))
+        self.zscale = 1.0
         if abs(yscale) > 1e-6 :
-            self.zscale=1.0/yscale
+            self.zscale = 1.0/yscale
 
-        #  keep just the z-values and null any ofsset
-        #  we might catch a reshape error that will occur if nx*ny != # of vertices (if we are not dealing with a heightfield but with a mesh with duplicate x,y coords, like an axis aligned cube
-        self.center=np.array([c[2] for c in verts],dtype=np.single).reshape(nx,ny)
-        self.center=(self.center-np.amin(self.center))*self.zscale
+        # keep just the z-values and null any ofsset
+        # we might catch a reshape error that will occur if nx*ny != # of vertices (if we are not dealing with a heightfield but with a mesh with duplicate x,y coords, like an axis aligned cube
+        self.center = np.array([c[2] for c in verts],dtype=np.single).reshape(nx,ny)
+        self.center = (self.center-np.amin(self.center))*self.zscale
         if self.rainmap is not None:
-            #rainmap = sorted(list({ tuple(t) for t in self.rainmap[::] }))
-            #self.rainmap=np.array([c[2] for c in rainmap],dtype=np.single).reshape(nx,ny)
             rmscale = np.max(self.center)
-            #self.rainmap = (self.center/rmscale) * np.exp(expfact*((self.center/rmscale)-1))
             self.rainmap = expfact + (1-expfact)*(self.center/rmscale)
 
+
     @staticmethod
     def fromBlenderMesh(me, vg, expfact):
-        g=Grid()
-        g.center=np.asarray(list(tuple(v.co) for v in me.vertices), dtype=np.single )
-        g.rainmap=None
-        print("VertexGroup\n",vg, file=sys.stderr)
+        g = Grid()
+        g.center = np.asarray(list(tuple(v.co) for v in me.vertices), dtype=np.single )
+        g.rainmap = None
         if vg is not None:
             for v in me.vertices:
                 vg.add([v.index],0.0,'ADD')
             g.rainmap=np.asarray(list( (v.co[0], v.co[1], vg.weight(v.index)) for v in me.vertices), dtype=np.single )
         g._sort(expfact)
-        #print("CentreMesh\n", np.array_str(g.center,precision=3), file=sys.stderr)
-        #print('rainmap',np.max(g.rainmap),np.min(g.rainmap))
         return g
 
-#    def rainmapcolor(me, vg):
-#        if vg is not None:
-#            for v in me.vertices:
-
 
     def setrainmap(self, rainmap):
         self.rainmap = rainmap
 
+
     def _verts(self, surface):
-        a=surface / self.zscale
-        minx=0.0 if self.minx is None else self.minx
-        miny=0.0 if self.miny is None else self.miny
-        maxx=1.0 if self.maxx is None else self.maxx
-        maxy=1.0 if self.maxy is None else self.maxy
-        dx=(maxx-minx)/(a.shape[0]-1)
-        dy=(maxy-miny)/(a.shape[1]-1)
+        a = surface / self.zscale
+        minx = 0.0 if self.minx is None else self.minx
+        miny = 0.0 if self.miny is None else self.miny
+        maxx = 1.0 if self.maxx is None else self.maxx
+        maxy = 1.0 if self.maxy is None else self.maxy
+        dx = (maxx - minx) / (a.shape[0] - 1)
+        dy = (maxy - miny) / (a.shape[1] - 1)
         for row in range(a.shape[0]):
-            row0=miny+row*dy
+            row0 = miny + row * dy
             for col in range(a.shape[1]):
-                col0=minx+col*dx
+                col0 = minx + col * dx
                 yield (row0 ,col0 ,a[row  ][col  ])
 
+
     def _faces(self):
         nrow, ncol = self.center.shape
         for row in range(nrow-1):
@@ -228,61 +224,70 @@ class Grid:
               yield (vi, vi+ncol, vi+1)
               yield (vi+1, vi+ncol, vi+ncol+1)      
 
-    def toBlenderMesh(self, me): # pass me as argument so that we don't need to import bpy and create a dependency
+
+    def toBlenderMesh(self, me):
+        # pass me as argument so that we don't need to import bpy and create a dependency
         # the docs state that from_pydata takes iterators as arguments but it will fail with generators because it does len(arg)
         me.from_pydata(list(self._verts(self.center)),[],list(self._faces()))
 
-    def toWaterMesh(self, me): # pass me as argument so that we don't need to import bpy and create a dependency
+
+    def toWaterMesh(self, me):
+        # pass me as argument so that we don't need to import bpy and create a dependency
         # the docs state that from_pydata takes iterators as arguments but it will fail with generators because it does len(arg)
         me.from_pydata(list(self._verts(self.water)),[],list(self._faces()))
 
+
     def peak(self, value=1):
         nx,ny = self.center.shape
         self.center[int(nx/2),int(ny/2)] += value
 
+
     def shelf(self, value=1):
         nx,ny = self.center.shape
         self.center[:nx/2] += value
 
+
     def mesa(self, value=1):
         nx,ny = self.center.shape
         self.center[nx/4:3*nx/4,ny/4:3*ny/4] += value
 
+
     def random(self, value=1):
         self.center += np.random.random_sample(self.center.shape)*value
 
+
     def neighborgrid(self):
-        self.up=np.roll(self.center,-1,0)
-        self.down=np.roll(self.center,1,0)
-        self.left=np.roll(self.center,-1,1)
-        self.right=np.roll(self.center,1,1)
+        self.up = np.roll(self.center,-1,0)
+        self.down = np.roll(self.center,1,0)
+        self.left = np.roll(self.center,-1,1)
+        self.right = np.roll(self.center,1,1)
+
 
     def zeroedge(self, quantity=None):
         c = self.center if quantity is None else quantity
-        c[0,:]=0
-        c[-1,:]=0
-        c[:,0]=0
-        c[:,-1]=0
+        c[0,:] = 0
+        c[-1,:] = 0
+        c[:,0] = 0
+        c[:,-1] = 0
+
 
     def diffuse(self, Kd, IterDiffuse, numexpr):
         self.zeroedge()
-        c     = self.center[1:-1,1:-1]
-        up    = self.center[ :-2,1:-1]
-        down  = self.center[2:  ,1:-1]
-        left  = self.center[1:-1, :-2]
+        c = self.center[1:-1,1:-1]
+        up = self.center[ :-2,1:-1]
+        down = self.center[2:  ,1:-1]
+        left = self.center[1:-1, :-2]
         right = self.center[1:-1,2:  ]
         if(numexpr and numexpr_available):
             self.center[1:-1,1:-1] = ne.evaluate('c + Kd * (up + down + left + right - 4.0 * c)')
         else:
             self.center[1:-1,1:-1] = c + (Kd/IterDiffuse) * (up + down + left + right - 4.0 * c)
-            print("diffuse: ", Kd)
         self.maxrss = max(getmemsize(), self.maxrss)
         return self.center
 
+
     def avalanche(self, delta, iterava, prob, numexpr):
         self.zeroedge()
-        #print(self.center)
-
         c     = self.center[1:-1,1:-1]
         up    = self.center[ :-2,1:-1]
         down  = self.center[2:  ,1:-1]
@@ -301,12 +306,12 @@ class Grid:
                  + where((right-c) < -delta,(right-c +delta)/2, 0)')
         else:
             sa = (
-            # incoming
+                # incoming
                    where((up   -c) > delta ,(up   -c -delta)/2, 0) 
                  + where((down -c) > delta ,(down -c -delta)/2, 0)
                  + where((left -c) > delta ,(left -c -delta)/2, 0)
                  + where((right-c) > delta ,(right-c -delta)/2, 0)
-            # outgoing
+                # outgoing
                  + where((up   -c) < -delta,(up   -c +delta)/2, 0) 
                  + where((down -c) < -delta,(down -c +delta)/2, 0)
                  + where((left -c) < -delta,(left -c +delta)/2, 0)
@@ -317,14 +322,16 @@ class Grid:
             self.avalanced[1:-1,1:-1] = self.avalanced[1:-1,1:-1] + sa/iterava 
             self.center[1:-1,1:-1] = c + sa/iterava
 
-        #print(self.center)
         self.maxrss = max(getmemsize(), self.maxrss)
         return self.center
 
+
     def rain(self, amount=1, variance=0, userainmap=False):
         self.water += (1.0 - np.random.random(self.water.shape) * variance) * (amount if ((self.rainmap is None) or (not userainmap)) else self.rainmap * amount) 
 
-    def spring(self, amount, px, py, radius): # px, py and radius are all fractions
+
+    def spring(self, amount, px, py, radius):
+        # px, py and radius are all fractions
         nx, ny = self.center.shape
         rx = max(int(nx*radius),1)
         ry = max(int(ny*radius),1)
@@ -332,320 +339,257 @@ class Grid:
         py = int(ny*py)
         self.water[px-rx:px+rx+1,py-ry:py+ry+1] += amount
 
+
     def river(self, Kc, Ks, Kdep, Ka, Kev, numexpr):
+        zeros = np.zeros
+        where = np.where
+        min = np.minimum
+        max = np.maximum
+        abs = np.absolute
+        arctan = np.arctan
+        sin = np.sin
+
+        center = (slice(   1,   -1,None),slice(   1,  -1,None))
+        up     = (slice(None,   -2,None),slice(   1,  -1,None))
+        down   = (slice(   2, None,None),slice(   1,  -1,None))
+        left   = (slice(   1,   -1,None),slice(None,  -2,None))
+        right  = (slice(   1,   -1,None),slice(   2,None,None))
+
+        water = self.water
+        rock = self.center
+        sediment = self.sediment
+        height = rock + water
+
+        # !! this gives a runtime warning for division by zero
+        verysmallnumber = 0.0000000001
+        water += verysmallnumber
+        sc = where(water > verysmallnumber, sediment / water, 0)
+
+        sdw = zeros(water[center].shape)
+        svdw = zeros(water[center].shape)
+        sds = zeros(water[center].shape)
+        angle = zeros(water[center].shape)
+        for d in (up,down,left,right):
+            if(numexpr and numexpr_available):
+                hdd = height[d]
+                hcc = height[center]
+                dw = ne.evaluate('hdd-hcc')
+                inflow = ne.evaluate('dw > 0')
+                wdd = water[d]
+                wcc = water[center]
+                dw = ne.evaluate('where(inflow, where(wdd<dw, wdd, dw), where(-wcc>dw, -wcc, dw))/4.0') # nested where() represent min() and max()
+                sdw  = ne.evaluate('sdw + dw')
+                scd  = sc[d]
+                scc  = sc[center]
+                rockd= rock[d]
+                rockc= rock[center]
+                sds  = ne.evaluate('sds + dw * where(inflow, scd, scc)')
+                svdw = ne.evaluate('svdw + abs(dw)')
+                angle= ne.evaluate('angle + arctan(abs(rockd-rockc))')
+            else:
+                dw = (height[d]-height[center])
+                inflow = dw > 0
+                dw = where(inflow, min(water[d], dw), max(-water[center], dw))/4.0
+                sdw  = sdw + dw
+                sds  = sds + dw * where(inflow, sc[d], sc[center])
+                svdw = svdw + abs(dw)
+                angle= angle + np.arctan(abs(rock[d]-rock[center]))
 
-      zeros = np.zeros
-      where = np.where
-      min = np.minimum
-      max = np.maximum
-      abs = np.absolute
-      arctan = np.arctan
-      sin = np.sin
-
-      center = (slice(   1,   -1,None),slice(   1,  -1,None))
-      #print("CentreSlice\n", np.array_str(center,precision=3), file=sys.stderr) 
-      up     = (slice(None,   -2,None),slice(   1,  -1,None))
-      down   = (slice(   2, None,None),slice(   1,  -1,None))
-      left   = (slice(   1,   -1,None),slice(None,  -2,None))
-      right  = (slice(   1,   -1,None),slice(   2,None,None))
-
-      water = self.water
-      rock = self.center
-      sediment = self.sediment
-      height = rock + water
-      sc = where(water>0, sediment/water, 0)  ##!! this gives a runtime warning for division by zero
-      sdw = zeros(water[center].shape)
-      svdw = zeros(water[center].shape)
-      sds = zeros(water[center].shape)
-      angle = zeros(water[center].shape)
-      #print(water[center])
-      for d in (up,down,left,right):
-          if(numexpr and numexpr_available):
-            hdd = height[d]
-            hcc = height[center]
-            dw = ne.evaluate('hdd-hcc')
-            inflow = ne.evaluate('dw > 0')
-            wdd = water[d]
+        if(numexpr and numexpr_available):
+            wcc = water[center]
+            scc = sediment[center]
+            rcc = rock[center]
+            water[center] = ne.evaluate('wcc + sdw')
+            sediment[center] = ne.evaluate('scc + sds')
+            sc = ne.evaluate('where(wcc>0, scc/wcc, 2000*Kc)')
+            fKc = ne.evaluate('Kc*sin(Ka*angle)*svdw')
+            ds = ne.evaluate('where(sc > fKc, -Kd * scc, Ks * svdw)')
+            rock[center] = ne.evaluate('rcc - ds')
+            # there isn't really a bottom to the rock but negative values look ugly
+            rock[center] = ne.evaluate('where(rcc<0,0,rcc)')
+            sediment[center] = ne.evaluate('scc + ds')
+        else:
             wcc = water[center]
-            dw = ne.evaluate('where(inflow, where(wdd<dw, wdd, dw), where(-wcc>dw, -wcc, dw))/4.0') # nested where() represent min() and max()
-            sdw  = ne.evaluate('sdw + dw')
-            scd  = sc[d]
-            scc  = sc[center]
-            rockd= rock[d]
-            rockc= rock[center]
-            sds  = ne.evaluate('sds + dw * where(inflow, scd, scc)')
-            svdw = ne.evaluate('svdw + abs(dw)')
-            angle= ne.evaluate('angle + arctan(abs(rockd-rockc))')
-          else:
-            dw = (height[d]-height[center])
-            inflow = dw > 0
-            dw = where(inflow, min(water[d], dw), max(-water[center], dw))/4.0
-            sdw  = sdw + dw
-            sds  = sds + dw * where(inflow, sc[d], sc[center])
-            svdw = svdw + abs(dw)
-            angle= angle + np.arctan(abs(rock[d]-rock[center]))
-
-      if(numexpr and numexpr_available):
-          wcc = water[center]
-          scc = sediment[center]
-          rcc = rock[center]
-          water[center] = ne.evaluate('wcc + sdw')
-          sediment[center] = ne.evaluate('scc + sds')
-          sc = ne.evaluate('where(wcc>0, scc/wcc, 2000*Kc)')
-          fKc = ne.evaluate('Kc*sin(Ka*angle)*svdw')
-          ds = ne.evaluate('where(sc > fKc, -Kd * scc, Ks * svdw)')
-          rock[center] = ne.evaluate('rcc - ds')
-          rock[center] = ne.evaluate('where(rcc<0,0,rcc)') # there isn't really a bottom to the rock but negative values look ugly
-          sediment[center] = ne.evaluate('scc + ds')
-      else:
-          wcc = water[center]
-          scc = sediment[center]
-          rcc = rock[center]
-          water[center] = wcc * (1-Kev) + sdw
-          sediment[center] = scc + sds
-          sc = where(wcc>0, scc/wcc, 2*Kc)
-          fKc = Kc*svdw
-          #fKc = Kc*np.sin(Ka*angle)*svdw*wcc
-          #ds = where(sc > fKc, -Kd * scc, Ks * svdw)
-          ds = where(fKc>sc,(fKc-sc)*Ks,(fKc-sc)*Kdep)*wcc
-          self.flowrate[center] = svdw
-          self.scour[center] = ds
-          self.sedimentpct[center] = sc
-          self.capacity[center] = fKc
-          #rock[center] = rcc - ds
-          #rock[center] = where(rcc<0,0,rcc) # there isn't really a bottom to the rock but negative values look ugly
-          sediment[center] = scc + ds + sds
-          #print("sdw", sdw[10,15])
+            scc = sediment[center]
+            rcc = rock[center]
+            water[center] = wcc * (1-Kev) + sdw
+            sediment[center] = scc + sds
+            sc = where(wcc > 0, scc / wcc, 2 * Kc)
+            fKc = Kc*svdw
+            ds = where(fKc > sc, (fKc - sc) * Ks, (fKc - sc) * Kdep) * wcc
+            self.flowrate[center] = svdw
+            self.scour[center] = ds
+            self.sedimentpct[center] = sc
+            self.capacity[center] = fKc
+            sediment[center] = scc + ds + sds
+
 
     def flow(self, Kc, Ks, Kz, Ka, numexpr):
+        zeros = np.zeros
+        where = np.where
+        min = np.minimum
+        max = np.maximum
+        abs = np.absolute
+        arctan = np.arctan
+        sin = np.sin
+
+        center = (slice(   1,   -1,None),slice(   1,  -1,None))
+        rock = self.center
+        ds = self.scour[center]    
+        rcc = rock[center]
+        rock[center] = rcc - ds * Kz
+        # there isn't really a bottom to the rock but negative values look ugly
+        rock[center] = where(rcc<0,0,rcc)
 
-      zeros = np.zeros
-      where = np.where
-      min = np.minimum
-      max = np.maximum
-      abs = np.absolute
-      arctan = np.arctan
-      sin = np.sin
-
-      center = (slice(   1,   -1,None),slice(   1,  -1,None))
-      #print("CentreSlice\n", np.array_str(center,precision=3), file=sys.stderr) 
-      #up     = (slice(None,   -2,None),slice(   1,  -1,None))
-      #down   = (slice(   2, None,None),slice(   1,  -1,None))
-      #left   = (slice(   1,   -1,None),slice(None,  -2,None))
-      #right  = (slice(   1,   -1,None),slice(   2,None,None))
-
-      #water = self.water
-      rock = self.center
-      #sediment = self.sediment
-      #height = rock + water
-      #sc = where(water>0, sediment/water, 0)  ##!! this gives a runtime warning for division by zero
-      #sdw = zeros(water[center].shape)
-      #svdw = zeros(water[center].shape)
-      #sds = zeros(water[center].shape)
-      #angle = zeros(water[center].shape)
-      #print(height[center])
-      #print(water[center])
-      #for d in (up,down,left,right):
-          #if(numexpr and numexpr_available):
-            #hdd = height[d]
-            #hcc = height[center]
-            #dw = ne.evaluate('hdd-hcc')
-            #inflow = ne.evaluate('dw > 0')
-            #wdd = water[d]
-            #wcc = water[center]
-            #dw = ne.evaluate('where(inflow, where(wdd<dw, wdd, dw), where(-wcc>dw, -wcc, dw))/4.0') # nested where() represent min() and max()
-            #sdw  = ne.evaluate('sdw + dw')
-            #scd  = sc[d]
-            #scc  = sc[center]
-            #rockd= rock[d]
-            #rockc= rock[center]
-            #sds  = ne.evaluate('sds + dw * where(inflow, scd, scc)')
-            #svdw = ne.evaluate('svdw + abs(dw)')
-            #angle= ne.evaluate('angle + arctan(abs(rockd-rockc))')
-          #else:
-            #dw = (height[d]-height[center])
-            #inflow = dw > 0
-            #dw = where(inflow, min(water[d], dw), max(-water[center], dw))/4.0
-            #sdw  = sdw + dw
-            #sds  = sds + dw * where(inflow, sc[d], sc[center])
-            #svdw = svdw + abs(dw)
-            #angle= angle + np.arctan(abs(rock[d]-rock[center]))
-
-      #if(numexpr and numexpr_available):
-          #wcc = water[center]
-          #scc = sediment[center]
-          #rcc = rock[center]
-          #water[center] = ne.evaluate('wcc + sdw')
-          #sediment[center] = ne.evaluate('scc + sds')
-          #sc = ne.evaluate('where(wcc>0, scc/wcc, 2000*Kc)')
-          #fKc = ne.evaluate('Kc*sin(Ka*angle)*svdw')
-          #ds = ne.evaluate('where(sc > fKc, -Kd * scc, Ks * svdw)')
-          #rock[center] = ne.evaluate('rcc - ds')
-          #rock[center] = ne.evaluate('where(rcc<0,0,rcc)') # there isn't really a bottom to the rock but negative values look ugly
-          #sediment[center] = ne.evaluate('scc + ds')
-      #else:
-          #wcc = water[center]
-          #scc = sediment[center]
-      ds = self.scour[center]    
-      rcc = rock[center]
-          #water[center] = wcc + sdw
-          #sediment[center] = scc + sds
-          #sc = where(wcc>0, scc/wcc, 2*Kc)
-          #fKc = Kc*np.sin(Ka*angle)*svdw
-          #ds = where(sc > fKc, -Kd * scc, Ks * svdw)
-      rock[center] = rcc - ds * Kz
-      rock[center] = where(rcc<0,0,rcc) # there isn't really a bottom to the rock but negative values look ugly
-          #sediment[center] = scc + ds
 
     def rivergeneration(self, rainamount, rainvariance, userainmap, Kc, Ks, Kdep, Ka, Kev, Kspring, Kspringx, Kspringy, Kspringr, numexpr):
         self.init_water_and_sediment()
         self.rain(rainamount, rainvariance, userainmap)
         self.zeroedge(self.water)
         self.zeroedge(self.sediment)
-        #self.spring(Kspring, Kspringx, Kspringy, Kspringr)
         self.river(Kc, Ks, Kdep, Ka, Kev, numexpr)
         self.watermax = np.max(self.water)
 
+
     def fluvial_erosion(self, rainamount, rainvariance, userainmap, Kc, Ks, Kdep, Ka, Kspring, Kspringx, Kspringy, Kspringr, numexpr):
-        #self.init_water_and_sediment()
-        #self.rain(rainamount, rainvariance, userainmap)
-        #self.zeroedge(self.water)
-        #self.zeroedge(self.sediment)
-        #self.spring(Kspring, Kspringx, Kspringy, Kspringr)
         self.flow(Kc, Ks, Kdep, Ka, numexpr)
         self.flowratemax = np.max(self.flowrate)
         self.scourmax = np.max(self.scour)
         self.scourmin = np.min(self.scour)
         self.sedmax = np.max(self.sediment)
-        print("DSMinMax", np.min(self.scour), np.max(self.scour))
+
 
     def analyze(self):
-      self.neighborgrid()
-      # just looking at up and left to avoid needless doubel calculations
-      slopes=np.concatenate((np.abs(self.left - self.center),np.abs(self.up - self.center)))
-      return '\n'.join(["%-15s: %.3f"%t for t in [
-        ('height average', np.average(self.center)),
-        ('height median', np.median(self.center)),
-        ('height max', np.max(self.center)),
-        ('height min', np.min(self.center)),
-        ('height std', np.std(self.center)),
-        ('slope average', np.average(slopes)),
-        ('slope median', np.median(slopes)),
-        ('slope max', np.max(slopes)),
-        ('slope min', np.min(slopes)),
-        ('slope std', np.std(slopes))
-        ]]
-      )
+        self.neighborgrid()
+        # just looking at up and left to avoid needless doubel calculations
+        slopes=np.concatenate((np.abs(self.left - self.center),np.abs(self.up - self.center)))
+        return '\n'.join(["%-15s: %.3f"%t for t in [
+                ('height average', np.average(self.center)),
+                ('height median', np.median(self.center)),
+                ('height max', np.max(self.center)),
+                ('height min', np.min(self.center)),
+                ('height std', np.std(self.center)),
+                ('slope average', np.average(slopes)),
+                ('slope median', np.median(slopes)),
+                ('slope max', np.max(slopes)),
+                ('slope min', np.min(slopes)),
+                ('slope std', np.std(slopes))
+                ]]
+            )
+
 
 class TestGrid(unittest.TestCase):
 
-  def test_diffuse(self):
-    g=Grid(5)
-    g.peak(1)
-    self.assertEqual(g.center[2,2],1.0)
-    g.diffuse(0.1, numexpr=False)
-    for n in [(2,1),(2,3),(1,2),(3,2)]:
-      self.assertAlmostEqual(g.center[n],0.1)
-    self.assertAlmostEqual(g.center[2,2],0.6)
-
-  def test_diffuse_numexpr(self):
-    g=Grid(5)
-    g.peak(1)
-    g.diffuse(0.1, numexpr=False)
-    h=Grid(5)
-    h.peak(1)
-    h.diffuse(0.1, numexpr=True)
-    self.assertEqual(list(g.center.flat),list(h.center.flat))
-
-  def test_avalanche_numexpr(self):
-    g=Grid(5)
-    g.peak(1)
-    g.avalanche(0.1, numexpr=False)
-    h=Grid(5)
-    h.peak(1)
-    h.avalanche(0.1, numexpr=True)
-    print(g)
-    print(h)
-    np.testing.assert_almost_equal(g.center,h.center)
+    def test_diffuse(self):
+        g = Grid(5)
+        g.peak(1)
+        self.assertEqual(g.center[2,2],1.0)
+        g.diffuse(0.1, numexpr=False)
+        for n in [(2,1),(2,3),(1,2),(3,2)]:
+            self.assertAlmostEqual(g.center[n],0.1)
+        self.assertAlmostEqual(g.center[2,2],0.6)
+
+
+    def test_diffuse_numexpr(self):
+        g = Grid(5)
+        g.peak(1)
+        g.diffuse(0.1, numexpr=False)
+        h = Grid(5)
+        h.peak(1)
+        h.diffuse(0.1, numexpr=True)
+        self.assertEqual(list(g.center.flat),list(h.center.flat))
+
+
+    def test_avalanche_numexpr(self):
+        g = Grid(5)
+        g.peak(1)
+        g.avalanche(0.1, numexpr=False)
+        h = Grid(5)
+        h.peak(1)
+        h.avalanche(0.1, numexpr=True)
+        print(g)
+        print(h)
+        np.testing.assert_almost_equal(g.center,h.center)
+
 
 if __name__ == "__main__":
 
-  import argparse
-
-  parser = argparse.ArgumentParser(description='Erode a terrain while assuming zero boundary conditions.')
-  parser.add_argument('-I', dest='iterations', type=int, default=1, help='the number of iterations')
-  parser.add_argument('-Kd', dest='Kd', type=float, default=0.01, help='Diffusion constant')
-  parser.add_argument('-Kh', dest='Kh', type=float, default=6, help='Maximum stable cliff height')
-  parser.add_argument('-Kp', dest='Kp', type=float, default=0.1, help='Avalanche probability for unstable cliffs')
-  parser.add_argument('-Kr', dest='Kr', type=float, default=0.1, help='Average amount of rain per iteration')
-  parser.add_argument('-Kspring', dest='Kspring', type=float, default=0.0, help='Average amount of wellwater per iteration')
-  parser.add_argument('-Kspringx', dest='Kspringx', type=float, default=0.5, help='relative x position of spring')
-  parser.add_argument('-Kspringy', dest='Kspringy', type=float, default=0.5, help='relative y position of spring')
-  parser.add_argument('-Kspringr', dest='Kspringr', type=float, default=0.02, help='radius of spring')
-  parser.add_argument('-Kdep', dest='Kdep', type=float, default=0.1, help='Sediment deposition constant')
-  parser.add_argument('-Ks', dest='Ks', type=float, default=0.1, help='Soil softness constant')
-  parser.add_argument('-Kc', dest='Kc', type=float, default=1.0, help='Sediment capacity')
-  parser.add_argument('-Ka', dest='Ka', type=float, default=2.0, help='Slope dependency of erosion')
-  parser.add_argument('-ri', action='store_true', dest='rawin', default=False, help='use Blender raw format for input')
-  parser.add_argument('-ro', action='store_true', dest='rawout', default=False, help='use Blender raw format for output')
-  parser.add_argument('-i',  action='store_true', dest='useinputfile', default=False, help='use an inputfile (instead of just a synthesized grid)')
-  parser.add_argument('-t',  action='store_true', dest='timingonly', default=False, help='do not write anything to an output file')
-  parser.add_argument('-infile', type=str, default="-", help='input filename')
-  parser.add_argument('-outfile', type=str, default="-", help='output filename')
-  parser.add_argument('-Gn', dest='gridsize', type=int, default=20, help='Gridsize (always square)')
-  parser.add_argument('-Gp', dest='gridpeak', type=float, default=0, help='Add peak with given height')
-  parser.add_argument('-Gs', dest='gridshelf', type=float, default=0, help='Add shelve with given height')
-  parser.add_argument('-Gm', dest='gridmesa', type=float, default=0, help='Add mesa with given height')
-  parser.add_argument('-Gr', dest='gridrandom', type=float, default=0, help='Add random values between 0 and given value')
-  parser.add_argument('-m', dest='threads', type=int, default=1, help='number of threads to use')
-  parser.add_argument('-u', action='store_true', dest='unittest', default=False, help='perfom unittests')
-  parser.add_argument('-a', action='store_true', dest='analyze', default=False, help='show some statistics of input and output meshes')
-  parser.add_argument('-d', action='store_true', dest='dump', default=False, help='show sediment and water meshes at end of run')
-  parser.add_argument('-n', action='store_true', dest='usenumexpr', default=False, help='use numexpr optimizations')
-
-  args = parser.parse_args()
-  print("\nInput arguments:")
-  print("\n".join("%-15s: %s"%t for t in sorted(vars(args).items())), file=sys.stderr)
-
-  if args.unittest:
-    unittest.main(argv=[sys.argv[0]])
-    sys.exit(0)
-
-  if args.useinputfile:
-    if args.rawin:
-      grid = Grid.fromRaw(args.infile)
-    else:
-      grid = Grid.fromFile(args.infile)
-  else:
-    grid = Grid(args.gridsize)
-
-  if args.gridpeak   > 0 : grid.peak(args.gridpeak)
-  if args.gridmesa   > 0 : grid.mesa(args.gridmesa)
-  if args.gridshelf  > 0 : grid.shelf(args.gridshelf)
-  if args.gridrandom > 0 : grid.random(args.gridrandom)
-
-  if args.analyze:
-    print('\nstatistics of the input grid:\n\n', grid.analyze(), file=sys.stderr, sep='' )
-  t = getptime()
-  for g in range(args.iterations):
-    if args.Kd > 0:
-      grid.diffuse(args.Kd, args.usenumexpr)
-    if args.Kh > 0 and args.Kp > rand():
-      grid.avalanche(args.Kh, args.usenumexpr)
-    if args.Kr > 0 or args.Kspring > 0:
-      grid.fluvial_erosion(args.Kr, args.Kc, args.Ks, args.Kdep, args.Ka, args.Kspring, args.Kspringx, args.Kspringy, args.Kspringr, args.usenumexpr)
-  t = getptime() - t
-  print("\nElapsed time: %.1f seconds, max memory %.1f Mb.\n"%(t,grid.maxrss), file=sys.stderr)
-  if args.analyze:
-    print('\nstatistics of the output grid:\n\n', grid.analyze(), file=sys.stderr, sep='')
-
-  if not args.timingonly:
-    if args.rawout:
-      grid.toRaw(args.outfile, vars(args))
+    import argparse
+
+    parser = argparse.ArgumentParser(description='Erode a terrain while assuming zero boundary conditions.')
+    parser.add_argument('-I', dest='iterations', type=int, default=1, help='the number of iterations')
+    parser.add_argument('-Kd', dest='Kd', type=float, default=0.01, help='Diffusion constant')
+    parser.add_argument('-Kh', dest='Kh', type=float, default=6, help='Maximum stable cliff height')
+    parser.add_argument('-Kp', dest='Kp', type=float, default=0.1, help='Avalanche probability for unstable cliffs')
+    parser.add_argument('-Kr', dest='Kr', type=float, default=0.1, help='Average amount of rain per iteration')
+    parser.add_argument('-Kspring', dest='Kspring', type=float, default=0.0, help='Average amount of wellwater per iteration')
+    parser.add_argument('-Kspringx', dest='Kspringx', type=float, default=0.5, help='relative x position of spring')
+    parser.add_argument('-Kspringy', dest='Kspringy', type=float, default=0.5, help='relative y position of spring')
+    parser.add_argument('-Kspringr', dest='Kspringr', type=float, default=0.02, help='radius of spring')
+    parser.add_argument('-Kdep', dest='Kdep', type=float, default=0.1, help='Sediment deposition constant')
+    parser.add_argument('-Ks', dest='Ks', type=float, default=0.1, help='Soil softness constant')
+    parser.add_argument('-Kc', dest='Kc', type=float, default=1.0, help='Sediment capacity')
+    parser.add_argument('-Ka', dest='Ka', type=float, default=2.0, help='Slope dependency of erosion')
+    parser.add_argument('-ri', action='store_true', dest='rawin', default=False, help='use Blender raw format for input')
+    parser.add_argument('-ro', action='store_true', dest='rawout', default=False, help='use Blender raw format for output')
+    parser.add_argument('-i',  action='store_true', dest='useinputfile', default=False, help='use an inputfile (instead of just a synthesized grid)')
+    parser.add_argument('-t',  action='store_true', dest='timingonly', default=False, help='do not write anything to an output file')
+    parser.add_argument('-infile', type=str, default="-", help='input filename')
+    parser.add_argument('-outfile', type=str, default="-", help='output filename')
+    parser.add_argument('-Gn', dest='gridsize', type=int, default=20, help='Gridsize (always square)')
+    parser.add_argument('-Gp', dest='gridpeak', type=float, default=0, help='Add peak with given height')
+    parser.add_argument('-Gs', dest='gridshelf', type=float, default=0, help='Add shelve with given height')
+    parser.add_argument('-Gm', dest='gridmesa', type=float, default=0, help='Add mesa with given height')
+    parser.add_argument('-Gr', dest='gridrandom', type=float, default=0, help='Add random values between 0 and given value')
+    parser.add_argument('-m', dest='threads', type=int, default=1, help='number of threads to use')
+    parser.add_argument('-u', action='store_true', dest='unittest', default=False, help='perfom unittests')
+    parser.add_argument('-a', action='store_true', dest='analyze', default=False, help='show some statistics of input and output meshes')
+    parser.add_argument('-d', action='store_true', dest='dump', default=False, help='show sediment and water meshes at end of run')
+    parser.add_argument('-n', action='store_true', dest='usenumexpr', default=False, help='use numexpr optimizations')
+
+    args = parser.parse_args()
+    print("\nInput arguments:")
+    print("\n".join("%-15s: %s"%t for t in sorted(vars(args).items())), file=sys.stderr)
+
+    if args.unittest:
+        unittest.main(argv=[sys.argv[0]])
+        sys.exit(0)
+
+    if args.useinputfile:
+        if args.rawin:
+            grid = Grid.fromRaw(args.infile)
+        else:
+            grid = Grid.fromFile(args.infile)
     else:
-      grid.toFile(args.outfile)
+        grid = Grid(args.gridsize)
+
+    if args.gridpeak > 0 : grid.peak(args.gridpeak)
+    if args.gridmesa > 0 : grid.mesa(args.gridmesa)
+    if args.gridshelf > 0 : grid.shelf(args.gridshelf)
+    if args.gridrandom > 0 : grid.random(args.gridrandom)
+
+    if args.analyze:
+        print('\nstatistics of the input grid:\n\n', grid.analyze(), file=sys.stderr, sep='' )
+    t = getptime()
+    for g in range(args.iterations):
+        if args.Kd > 0:
+            grid.diffuse(args.Kd, args.usenumexpr)
+        if args.Kh > 0 and args.Kp > rand():
+            grid.avalanche(args.Kh, args.usenumexpr)
+        if args.Kr > 0 or args.Kspring > 0:
+            grid.fluvial_erosion(args.Kr, args.Kc, args.Ks, args.Kdep, args.Ka, args.Kspring, args.Kspringx, args.Kspringy, args.Kspringr, args.usenumexpr)
+    t = getptime() - t
+    print("\nElapsed time: %.1f seconds, max memory %.1f Mb.\n"%(t,grid.maxrss), file=sys.stderr)
+    if args.analyze:
+        print('\nstatistics of the output grid:\n\n', grid.analyze(), file=sys.stderr, sep='')
+
+    if not args.timingonly:
+        if args.rawout:
+            grid.toRaw(args.outfile, vars(args))
+        else:
+            grid.toFile(args.outfile)
 
-  if args.dump:
-    print("sediment\n", np.array_str(grid.sediment,precision=3), file=sys.stderr)
-    print("water\n", np.array_str(grid.water,precision=3), file=sys.stderr)
-    print("sediment concentration\n", np.array_str(grid.sediment/grid.water,precision=3), file=sys.stderr)
+    if args.dump:
+        print("sediment\n", np.array_str(grid.sediment,precision=3), file=sys.stderr)
+        print("water\n", np.array_str(grid.water,precision=3), file=sys.stderr)
+        print("sediment concentration\n", np.array_str(grid.sediment/grid.water,precision=3), file=sys.stderr)
diff --git a/ant_landscape/mesh_ant_displace.py b/ant_landscape/mesh_ant_displace.py
index d40aa26d..a0ecebb4 100644
--- a/ant_landscape/mesh_ant_displace.py
+++ b/ant_landscape/mesh_ant_displace.py
@@ -17,7 +17,7 @@
 # ##### END GPL LICENSE BLOCK #####
 
 # Another Noise Tool - Mesh Displace
-# Jim Hazevoet
+# Jimmy Hazevoet
 
 # ------------------------------------------------------------
 # import modules
@@ -43,7 +43,7 @@ from .ant_functions import (
 class AntMeshDisplace(bpy.types.Operator):
     bl_idname = "mesh.ant_displace"
     bl_label = "Another Noise Tool - Displace"
-    bl_description = "A.N.T. Displace mesh vertices"
+    bl_description = "Displace mesh vertices"
     bl_options = {'REGISTER', 'UNDO', 'PRESET'}
 
     ant_terrain_name = StringProperty(
@@ -162,7 +162,7 @@ class AntMeshDisplace(bpy.types.Operator):
             )
     noise_size = FloatProperty(
             name="Noise Size",
-            default=1.0,
+            default=0.25,
             min=0.01,
             max=1000.0,
             description="Noise size"
@@ -187,9 +187,10 @@ class AntMeshDisplace(bpy.types.Operator):
                 ('vl_hTerrain', "vlNoise hTerrain", "A.N.T: vlNoise hTerrain", 12),
                 ('distorted_heteroTerrain', "Distorted hTerrain", "A.N.T distorted hTerrain", 13),
                 ('double_multiFractal', "Double MultiFractal", "A.N.T: double multiFractal", 14),
-                ('slick_rock', "Slick Rock", "A.N.T: slick rock", 15),
-                ('planet_noise', "Planet Noise", "Planet Noise by: Farsthary", 16),
-                ('blender_texture', "Blender Texture - Texture Nodes", "Blender texture data block", 17)]
+                ('rocks_noise', "Noise Rocks", "A.N.T: turbulence variation", 15),
+                ('slick_rock', "Slick Rock", "A.N.T: slick rock", 16),
+                ('planet_noise', "Planet Noise", "Planet Noise by: Farsthary", 17),
+                ('blender_texture', "Blender Texture - Texture Nodes", "Blender texture data block", 18)]
             )
     basis_type = EnumProperty(
             name="Noise Basis",
@@ -325,7 +326,7 @@ class AntMeshDisplace(bpy.types.Operator):
         )
     height = FloatProperty(
             name="Height",
-            default=0.5,
+            default=0.25,
             min=-10000.0,
             max=10000.0,
             description="Noise intensity scale"
@@ -408,8 +409,8 @@ class AntMeshDisplace(bpy.types.Operator):
                 ("1", "Smooth", "Smooth transitions", 1),
                 ("2", "Sharp Sub", "Sharp substract transitions", 2),
                 ("3", "Sharp Add", "Sharp add transitions", 3),
-                ("4", "Posterize", "Posterize", 4),
-                ("5", "Posterize Mix", "Posterize mixed", 5)]
+                ("4", "Quantize", "Quantize", 4),
+                ("5", "Quantize Mix", "Quantize mixed", 5)]
             )
     water_plane = BoolProperty(
             name="Water Plane",
@@ -547,24 +548,40 @@ class AntMeshDisplace(bpy.types.Operator):
             self.remove_double
             ]
 
+        '''
+        gi = ob.vertex_groups["Group"].index # get group index
+        for v in ob.data.vertices:
+          for g in v.groups:
+            if g.group == gi: # compare with index in VertexGroupElement
+              v.co[0] = 5
+        '''
         # do displace
         mesh = ob.data
 
         if self.use_vgroup is True:
             vertex_group = ob.vertex_groups.active
             if vertex_group:
+                gi = vertex_group.index
                 if self.direction == "X":
                     for v in mesh.vertices:
-                        v.co[0] += vertex_group.weight(v.index) * noise_gen(v.co, props)
+                        for g in v.groups:
+                            if g.group == gi:
+                                v.co[0] += vertex_group.weight(v.index) * noise_gen(v.co, props)
                 if self.direction == "Y":
                     for v in mesh.vertices:
-                        v.co[1] += vertex_group.weight(v.index) * noise_gen(v.co, props)
+                        for g in v.groups:
+                            if g.group == gi:
+                                v.co[1] += vertex_group.weight(v.index) * noise_gen(v.co, props)
                 if self.direction == "Z":
                     for v in mesh.vertices:
-                        v.co[2] += vertex_group.weight(v.index) * noise_gen(v.co, props)
+                        for g in v.groups:
+                            if g.group == gi:
+                                v.co[2] += vertex_group.weight(v.index) * noise_gen(v.co, props)
                 else:
                     for v in mesh.vertices:
-                        v.co += vertex_group.weight(v.index) * v.normal * noise_gen(v.co, props)
+                        for g in v.groups:
+                            if g.group == gi:
+                                v.co += vertex_group.weight(v.index) * v.normal * noise_gen(v.co, props)
 
         else:
             if self.direction == "X":
@@ -582,9 +599,6 @@ class AntMeshDisplace(bpy.types.Operator):
 
             mesh.update()
 
-        if bpy.ops.object.shade_smooth == True:
-            bpy.ops.object.shade_smooth()
-
         if self.auto_refresh is False:
             self.refresh = False
 
diff --git a/ant_landscape/stats.py b/ant_landscape/stats.py
index 983d9d94..77759804 100644
--- a/ant_landscape/stats.py
+++ b/ant_landscape/stats.py
@@ -1,52 +1,52 @@
 from time import time
 
 try:
-  import psutil
-  print('psutil available')
-  psutil_available=True
+    import psutil
+    print('psutil available')
+    psutil_available=True
 except ImportError:
-  psutil_available=False
+    psutil_available=False
 
 class Stats:
-  def __init__(self):
-    self.memstats_available = False
-    if psutil_available:
-      self.process=psutil.Process()
-      self.memstats_available = True
-    self.reset()
+    def __init__(self):
+        self.memstats_available = False
+        if psutil_available:
+            self.process=psutil.Process()
+            self.memstats_available = True
+        self.reset()
 
-  def reset(self):
-    self.lasttime = self._gettime()
-    self.lastmem = self._getmem()
-    self.basemem = self.lastmem
-    self.maxmem = 0
-    self.elapsedtime = 0
+    def reset(self):
+        self.lasttime = self._gettime()
+        self.lastmem = self._getmem()
+        self.basemem = self.lastmem
+        self.maxmem = 0
+        self.elapsedtime = 0
 
-  def _gettime(self):
-    """return the time in seconds used by the current process."""
-    if psutil_available:
-      m=self.process.get_cpu_times()
-      return m.user + m.system
-    return time()
+    def _gettime(self):
+        """return the time in seconds used by the current process."""
+        if psutil_available:
+            m=self.process.get_cpu_times()
+            return m.user + m.system
+        return time()
 
-  def _getmem(self):
-    """return the resident set size in bytes used by the current process."""
-    if psutil_available:
-      m = self.process.get_memory_info()
-      return m.rss
-    return 0
+    def _getmem(self):
+        """return the resident set size in bytes used by the current process."""
+        if psutil_available:
+            m = self.process.get_memory_info()
+            return m.rss
+        return 0
 
-  def time(self):
-    """return the time since the last call in seconds used by the current process."""
-    old = self.lasttime
-    self.lasttime = self._gettime()
-    self.elapsedtime = self.lasttime - old
-    return self.elapsedtime
+    def time(self):
+        """return the time since the last call in seconds used by the current process."""
+        old = self.lasttime
+        self.lasttime = self._gettime()
+        self.elapsedtime = self.lasttime - old
+        return self.elapsedtime
 
-  def memory(self):
-    """return the maximum resident set size since the first call in bytes used by the current process."""
-    self.lastmem = self._getmem()
-    d = self.lastmem - self.basemem
-    if d>self.maxmem:
-      self.maxmem = d
-    return self.maxmem
+    def memory(self):
+        """return the maximum resident set size since the first call in bytes used by the current process."""
+        self.lastmem = self._getmem()
+        d = self.lastmem - self.basemem
+        if d > self.maxmem:
+            self.maxmem = d
+        return self.maxmem
diff --git a/ant_landscape/test.py b/ant_landscape/test.py
index c32a496c..d56815f2 100644
--- a/ant_landscape/test.py
+++ b/ant_landscape/test.py
@@ -18,4 +18,4 @@ if __name__ == '__main__':
     print(stats.memory())
     a = cos(a)**2+sin(a)**2
     print(stats.time())
- print(stats.memory())
\ No newline at end of file
+    print(stats.memory())
diff --git a/ant_landscape/utils.py b/ant_landscape/utils.py
index 595d5250..be8f1f29 100644
--- a/ant_landscape/utils.py
+++ b/ant_landscape/utils.py
@@ -1,7 +1,7 @@
 numexpr_available=False
 try:
-  import numexpr
-  numexpr_available=True
+    import numexpr
+    numexpr_available=True
 except ImportError:
-  pass
+    pass