Adopting new changes from plate-tectonics

.gitignore

@@ -1,3 +1,5 @@
+venv
+
 *.py[cod]
 
 # C extensions

platec_src/heightmap.hpp

@@ -10,7 +10,7 @@ using namespace std;
 template <typename Value>
 class Matrix
 {
-	public:
+public:
 
     Matrix(unsigned int width, unsigned int height)
         : _width(width), _height(height)
@@ -39,9 +39,10 @@ class Matrix
 
     const void set_all(const Value& value)
     {
-        for (int x=0; x<_width;x++){
-            for (int y=0; y<_height;y++){
-                set(x,y,value);
+        // we cannot use memset to make it very general
+        for (int x=0; x<_width; x++){
+            for (int y=0; y<_height; y++){
+                set(x, y, value);
             }
         }
     }
@@ -63,17 +64,17 @@ class Matrix
         return _data[y * _width + x];
     }
 
-    Matrix<Value> & operator=(const Matrix<Value>& other)
+    Matrix<Value>& operator=(const Matrix<Value>& other)
     {
         if (this != &other) // prevent self-assignment
         {
             _width  = other._width;
             _height = other._height;
             delete[] _data;
             _data = new Value[_width * _height];
-            for (int x=0; x<_width;x++){
-                for (int y=0; y<_height;y++){
-                    set(x,y,other.get(x,y));
+            for (int x=0; x<_width; x++){
+                for (int y=0; y<_height; y++){
+                    set(x, y, other.get(x,y));
                 }
             }
         }
@@ -106,8 +107,8 @@ class Matrix
     }
 
     bool equals(Value* other)
-    {
-        for (int i=0; i<(_width*_height);i++){
+    {                
+        for (int i=0; i< (_width *_height); i++){
             if (_data[i] != other[i]){
                 return false;
             }
@@ -127,15 +128,21 @@ class Matrix
 
     int xMod(int x) const
     {
-        return (x + _width) % _width;
+        while (x < 0) {
+            x += _width;
+        }
+        return x % _width;
     }
 
     int yMod(int y) const
     {
-        return (y + _height) % _height;
+        while (y < 0) {
+            y += _height;
+        }
+        return y % _height;
     }    
 
-	private:
+private:
 
     Value* _data;
     unsigned int _width;

platec_src/lithosphere.cpp

@@ -20,6 +20,7 @@
 #include "plate.hpp"
 #include "sqrdmd.hpp"
 #include "simplexnoise.hpp"
+#include "noise.hpp"
 
 #include <cfloat>
 #include <cmath>
@@ -49,7 +50,6 @@ class plateArea
     size_t hgt; ///< Height of area in pixels.
 };
 
-static const float SQRDMD_ROUGHNESS = 0.5f;
 static const float SUBDUCT_RATIO = 0.5f;
 
 static const float BUOYANCY_BONUS_X = 3;
@@ -64,45 +64,9 @@ size_t findBound(const size_t* map, size_t length, size_t x0, size_t y0,
                  int dx, int dy);
 size_t findPlate(plate** plates, float x, float y, size_t num_plates);
 
-static uint32_t nearest_pow(uint32_t num)
-{
-    uint32_t n = 1;
-
-    while (n < num){
-        n <<= 1;
-    }
-
-    return n;
-}
-
 void lithosphere::createNoise(float* tmp, const WorldDimension& tmpDim, bool useSimplex)
 {
-try {
-    if (useSimplex) {
-        simplexnoise(_randsource(), tmp, 
-            tmpDim.getWidth(), 
-            tmpDim.getHeight(), 
-            SQRDMD_ROUGHNESS);
-    } else {        
-        size_t side = tmpDim.getMax();
-        side = nearest_pow(side)+1;
-        float* squareTmp = new float[side*side];
-        memset(squareTmp, 0, sizeof(float)*side*side);
-        for (int y=0; y<tmpDim.getHeight(); y++){
-            memcpy(&squareTmp[y*side],&tmp[y*tmpDim.getWidth()],sizeof(float)*tmpDim.getWidth());
-        }
-        sqrdmd(_randsource(), squareTmp, side, SQRDMD_ROUGHNESS);
-        for (int y=0; y<tmpDim.getHeight(); y++){
-            memcpy(&tmp[y*tmpDim.getWidth()],&squareTmp[y*side],sizeof(float)*tmpDim.getWidth());
-        }
-        delete[] squareTmp;
-    }    
-} catch (const exception& e){
-    std::string msg = "Problem during lithosphere::createNoise, tmpDim+=";
-    msg = msg + to_string(tmpDim.getWidth()) + "x" + to_string(tmpDim.getHeight()) + " ";
-    msg = msg + e.what();
-    throw runtime_error(msg.c_str());
-}
+    ::createNoise(tmp, tmpDim, _randsource, useSimplex);
 }
 
 lithosphere::lithosphere(long seed, size_t width, size_t height, float sea_level,
@@ -121,7 +85,8 @@ lithosphere::lithosphere(long seed, size_t width, size_t height, float sea_level
     max_plates(0), 
     num_plates(0),
     _worldDimension(width, height),
-    _randsource(seed)
+    _randsource(seed),
+    _steps(0)
 {
     WorldDimension tmpDim = WorldDimension(width+1, height+1);
     const size_t A = tmpDim.getArea();
@@ -366,9 +331,15 @@ float* lithosphere::getTopography() const throw()
     return hmap.raw_data();
 }
 
+bool lithosphere::isFinished() const 
+{
+    return getPlateCount() == 0;
+}
+
 void lithosphere::update()
 {
 try {
+    _steps++;
     float totalVelocity = 0;
     float systemKineticEnergy = 0;
 
@@ -749,7 +720,7 @@ try {
     cycle_count += max_cycles > 0; // No increment if running for ever.
     if (cycle_count > max_cycles)
         return;
-
+    
     // Update height map to include all recent changes.
     hmap.set_all(0);
     for (size_t i = 0; i < num_plates; ++i)
@@ -852,7 +823,7 @@ try {
     }
 
     for (size_t y = 0; y < _worldDimension.getHeight(); y += 8)
-    {
+    {        
         for (size_t x = 0; x < _worldDimension.getWidth(); x += 8)
         {
             size_t i = _worldDimension.indexOf(x, y);
@@ -958,6 +929,10 @@ try {
         h_highest = h_highest > hmap[i] ? h_highest : hmap[i];
     }
 
+    /*
+    Removing this piece of code because is causing a regular grid of dots
+    polluting square maps not using a side which is a power of 2
+    -----------------------------------------------
     for (size_t y = 0; y < _worldDimension.getHeight(); y += 4)
     {
         for (size_t x = 0; x < _worldDimension.getWidth(); x += 4)
@@ -991,7 +966,7 @@ try {
         }
 
         memset(&hmap[_worldDimension.lineIndex((y+3) % _worldDimension.getHeight())], 0, line_size);
-    }
+    }*/
 
     for (size_t y = 0; y < _worldDimension.getHeight(); ++y) // Copy map into fractal buffer.
     {
@@ -1020,6 +995,7 @@ try {
                 hmap[_worldDimension.indexOf(x, y)] = tmp[tmpDim.indexOf(x, y)];
             else
                 hmap[_worldDimension.indexOf(x, y)] = original[_worldDimension.indexOf(x, y)];
+
 } catch (const exception& e){
     std::string msg = "Problem during restart: ";
     msg = msg + e.what();

platec_src/lithosphere.hpp

@@ -88,6 +88,7 @@ class lithosphere
 	void update(); ///< Simulate one step of plate tectonics.	
 	size_t getWidth() const;
 	size_t getHeight() const;
+	bool isFinished() const;
 
   protected:
   private:
@@ -146,6 +147,7 @@ class lithosphere
 
 	const WorldDimension _worldDimension;
 	mt19937 _randsource;
+	int _steps;
 };
 
 

platec_src/noise.cpp

@@ -0,0 +1,83 @@
+#include "noise.hpp"
+#include "sqrdmd.hpp"
+#include "simplexnoise.hpp"
+
+static const float SQRDMD_ROUGHNESS = 0.35f;
+
+static uint32_t nearest_pow(uint32_t num)
+{
+    uint32_t n = 1;
+
+    while (n < num){
+        n <<= 1;
+    }
+
+    return n;
+}
+
+void createNoise(float* tmp, const WorldDimension& tmpDim, mt19937 randsource, bool useSimplex)
+{
+try {
+    if (useSimplex) {
+        simplexnoise(randsource(), tmp, 
+            tmpDim.getWidth(), 
+            tmpDim.getHeight(), 
+            SQRDMD_ROUGHNESS);
+    } else {        
+        size_t side = tmpDim.getMax();
+        side = nearest_pow(side)+1;
+        float* squareTmp = new float[side*side];
+        memset(squareTmp, 0, sizeof(float)*side*side);
+        for (int y=0; y<tmpDim.getHeight(); y++){
+            memcpy(&squareTmp[y*side],&tmp[y*tmpDim.getWidth()],sizeof(float)*tmpDim.getWidth());
+        }
+        // to make it tileable we need to insert proper values in the padding area
+        // 1) on the right of the valid area        
+        for (int y=0; y<tmpDim.getHeight(); y++){
+        	for (int x=tmpDim.getWidth(); x<side; x++){
+        		// we simply put it as a mix between the east and west border (they should be fairly
+        		// similar because it is a toroidal world)
+        		squareTmp[y*side+x] = (squareTmp[y*side+0] + squareTmp[y*side+(tmpDim.getWidth()-1)])/2;
+        	}
+        }
+        // 2) below the valid area
+        for (int y=tmpDim.getHeight(); y<side; y++){
+        	for (int x=0; x<side; x++){
+        		// we simply put it as a mix between the north and south border (they should be fairly
+        		// similar because it is a toroidal world)
+        		squareTmp[y*side+x] = (squareTmp[(0)*side+x] + squareTmp[(tmpDim.getHeight()-1)*side+x])/2;
+        	}
+        }        
+
+        sqrdmd(randsource(), squareTmp, side, SQRDMD_ROUGHNESS);
+
+        // Calcuate deltas (noise introduced)
+        float* deltas = new float[tmpDim.getWidth()*tmpDim.getHeight()];
+        for (int y=0; y<tmpDim.getHeight(); y++){
+        	for (int x=0; x<tmpDim.getWidth(); x++){
+        		deltas[y*tmpDim.getWidth()+x] = squareTmp[y*side+x]-tmp[y*tmpDim.getWidth()+x];
+        	}
+        }
+
+        // make it tileable
+        for (int y=0; y<tmpDim.getHeight(); y++){
+        	for (int x=0; x<tmpDim.getWidth(); x++){
+        		int specularX = tmpDim.getWidth() - 1 - x;
+        		int specularY = tmpDim.getHeight() -1 - y;
+        		float myDelta = deltas[y*tmpDim.getWidth()+x];
+        		float specularWidthDelta = deltas[y*tmpDim.getWidth()+specularX];
+        		float specularHeightDelta = deltas[specularY*tmpDim.getWidth()+x];
+        		float oppositeDelta = deltas[specularY*tmpDim.getWidth()+specularX];
+        		tmp[y*tmpDim.getWidth()+x] += (myDelta + specularWidthDelta + specularHeightDelta + oppositeDelta)/4;
+        	}
+        }
+
+        delete[] squareTmp;
+    }    
+} catch (const exception& e){
+    std::string msg = "Problem during lithosphere::createNoise, tmpDim+=";
+    msg = msg + to_string(tmpDim.getWidth()) + "x" + to_string(tmpDim.getHeight()) + " ";
+    msg = msg + e.what();
+    throw runtime_error(msg.c_str());
+}
+}

platec_src/noise.hpp

@@ -0,0 +1,9 @@
+#ifndef NOISE_HPP
+#define NOISE_HPP
+
+#include <random>
+#include "rectangle.hpp"
+
+void createNoise(float* tmp, const WorldDimension& tmpDim, mt19937 _randsource, bool useSimplex = false);
+
+#endif

platec_src/platecapi.cpp

@@ -102,8 +102,11 @@ lithosphere* platec_api_get_lithosphere(size_t id)
 size_t platec_api_is_finished(void *pointer)
 {
 	lithosphere* litho = (lithosphere*)pointer;
-
-	return litho->getPlateCount() == 0;
+	if (litho->isFinished()) {
+		return 1;
+	} else {
+		return 0;
+	}
 }
 
 void platec_api_step(void *pointer)