Added erosion, using the same algorithm than the Python version.

main.cpp

@@ -13,21 +13,26 @@ int main(int argc, char** argv)
 	const char* file_to_read = (argc > 1) ? argv[1] : "dem";
 	if (file_exists(file_to_read)) {
 		std::cout << "Reading file '" << file_to_read << "' ..." << std::endl;
-		Map<double>* map = read_file<double>(file_to_read);
-		RiverMapper rm(map);
+		Map<double>* dem = read_file<double>(file_to_read);
+		RiverMapper rm(dem);
+		rm.flow();
+		rm.erode(1.0f);
 		rm.flow();
 
 		std::cout << "Writing files..." << std::endl;
-		Map<uint8_t> dirs(map->width, map->height, rm.dirs);
-		const char* file_dirs = (argc > 2) ? argv[2] : "dirs";
+		const char* file_dem = (argc > 2) ? argv[2] : "dem_new";
+		dem->write_file(file_dem);
+
+		Map<uint8_t> dirs(dem->width, dem->height, rm.dirs);
+		const char* file_dirs = (argc > 3) ? argv[3] : "dirs";
 		dirs.write_file(file_dirs);
 
-		Map<double> water(map->width, map->height, rm.water);
-		const char* file_rivers = (argc > 3) ? argv[3] : "rivers";
+		Map<double> water(dem->width, dem->height, rm.water);
+		const char* file_rivers = (argc > 4) ? argv[4] : "rivers";
 		water.write_file(file_rivers);
 
-		Map<double> lakes(map->width, map->height, rm.lakes);
-		const char* file_lakes = (argc > 4) ? argv[4] : "lakes";
+		Map<double> lakes(dem->width, dem->height, rm.lakes);
+		const char* file_lakes = (argc > 5) ? argv[5] : "lakes";
 		lakes.write_file(file_lakes);
 		std::cout << "Done." << std::endl;
 	} else {

rivermapper.cpp

@@ -20,6 +20,7 @@ RiverMapper::RiverMapper(Map<double>* dem_map) : width(dem_map->width), height(d
 	dirs = new uint8_t[size];
 	ndonors = new uint8_t[size];
 	basin_id = new uint32_t[size];
+	erosion_time = new double[size];
 	pcg = new PcgRandom();
 }
 
@@ -31,17 +32,26 @@ RiverMapper::~RiverMapper()
 	delete[] dirs;
 	delete[] ndonors;
 	delete[] basin_id;
+	delete[] erosion_time;
 	delete[] pcg;
 }
 
-inline double diff_or_zero(double a, double b) {
+inline double diff_or_zero(const double a, const double b) {
 	return (a > b) ? a-b : 0;
 }
 
-inline double max(double a, double b) {
+inline double min(const double a, const double b) {
+	return (a < b) ? a : b;
+}
+
+inline double max(const double a, const double b) {
 	return (a > b) ? a : b;
 }
 
+inline double max3(const double a, const double b, const double c) {
+	return (a > b) ? ((a > c) ? a : c) : ((b > c) ? b : c);
+}
+
 inline int RiverMapper::flow_local(const double zdiffs[], const int nz)
 {
 	double sum = 0.0f;
@@ -602,3 +612,38 @@ void RiverMapper::accumulate()
 		}
 	}
 }
+
+void RiverMapper::erode(const double time, const double K, const double m, const double sea_level) {
+	cout << "Eroding landscape..." << endl;
+	clock_t t0 = clock();
+
+	for (size_t i=1; i<size; i++) {
+		erosion_time[i] = 1.0f / (K * std::pow(water[i], m));
+		lakes[i] = max3(dem[i], lakes[i], sea_level); 
+	}
+
+	for (size_t i=1; i<size; i++) {
+		size_t iw = i;
+		double remaining = time;
+		double new_elev;
+		while (true) {
+			if (dirs_ref[iw] == SIZE_MAX) {
+				new_elev = lakes[iw];
+				break;
+			}
+			if (remaining <= erosion_time[iw]) {
+				double c = remaining / erosion_time[iw];
+				new_elev = (1-c) * lakes[iw] + c * lakes[dirs_ref[iw]];
+				break;
+			}
+
+			remaining -= erosion_time[iw];
+			iw = dirs_ref[iw];
+		}
+
+		dem[i] = min(dem[i], new_elev);
+	}
+
+	clock_t t1 = clock();
+	cout << "\tCompleted in " << disptime(t1-t0) << " s" << endl;
+}

rivermapper.h

@@ -17,6 +17,7 @@ public:
 
 	void flow();
 	void accumulate();
+	void erode(const double time, const double K=1.0f, const double m=0.5f, const double sea_level=0.0f);
 
 	const std::size_t width;
 	const std::size_t height;
@@ -27,5 +28,6 @@ private:
 	size_t* dirs_ref;
 	uint8_t* ndonors;
 	uint32_t* basin_id;
+	double* erosion_time;
 	PcgRandom* pcg;
 };