2020-12-19 16:55:54 -08:00
# include <cstdint>
# include <vector>
# include <unordered_map>
# include <stack>
# include <ctime>
# include <limits>
# include <cmath>
# include <iostream>
# include "rivermapper.h"
# include "random.h"
using namespace std ;
RiverMapper : : RiverMapper ( Map < double > * dem_map ) : width ( dem_map - > width ) , height ( dem_map - > height ) , size ( dem_map - > size )
{
dem = dem_map - > data ;
water = new double [ size ] ;
dirs_ref = new size_t [ size ] ;
lakes = new double [ size ] ;
dirs = new uint8_t [ size ] ;
ndonors = new uint8_t [ size ] ;
basin_id = new uint32_t [ size ] ;
2020-12-27 09:41:36 -08:00
erosion_time = new double [ size ] ;
2020-12-19 16:55:54 -08:00
pcg = new PcgRandom ( ) ;
}
RiverMapper : : ~ RiverMapper ( )
{
delete [ ] water ;
delete [ ] dirs_ref ;
delete [ ] lakes ;
delete [ ] dirs ;
delete [ ] ndonors ;
delete [ ] basin_id ;
2020-12-27 09:41:36 -08:00
delete [ ] erosion_time ;
2020-12-19 16:55:54 -08:00
delete [ ] pcg ;
}
2020-12-27 09:41:36 -08:00
inline double diff_or_zero ( const double a , const double b ) {
2020-12-19 16:55:54 -08:00
return ( a > b ) ? a - b : 0 ;
}
2020-12-27 09:41:36 -08:00
inline double min ( const double a , const double b ) {
return ( a < b ) ? a : b ;
}
inline double max ( const double a , const double b ) {
2020-12-19 16:55:54 -08:00
return ( a > b ) ? a : b ;
}
2020-12-27 09:41:36 -08:00
inline double max3 ( const double a , const double b , const double c ) {
return ( a > b ) ? ( ( a > c ) ? a : c ) : ( ( b > c ) ? b : c ) ;
}
2020-12-19 16:55:54 -08:00
inline int RiverMapper : : flow_local ( const double zdiffs [ ] , const int nz )
{
double sum = 0.0f ;
for ( int i = 0 ; i < nz ; i + + ) {
sum + = zdiffs [ i ] ;
}
if ( sum < = 0.0f )
return 0 ;
double rn = pcg - > rand ( ) * sum ;
for ( int i = 0 ; i < nz ; i + + ) {
if ( zdiffs [ i ] > rn ) {
return i + 1 ;
}
rn - = zdiffs [ i ] ;
}
return 0 ;
}
inline uint64_t get_key ( uint32_t k1 , uint32_t k2 ) {
return ( ( uint64_t ) k1 < < 32 ) | ( uint64_t ) k2 ;
}
struct Link { uint32_t basins [ 2 ] ; size_t x ; size_t y ; double elev ; bool isY ; } ;
struct LinkManager {
Link * linklist ;
unordered_map < uint64_t , Link * > linkmap ;
size_t n = 0 ;
LinkManager ( size_t num_basins ) {
size_t max_num_links = 1 ;
if ( num_basins > 2 )
max_num_links = ( size_t ) num_basins * 3 - 6 ; // From Euler's characteristics, a planar graph of n vertices can have at most 3n-6 edges.
linklist = new Link [ max_num_links ] ;
linkmap . reserve ( max_num_links ) ;
} ;
~ LinkManager ( ) {
delete [ ] linklist ;
} ;
void add_link ( uint32_t b1 , uint32_t b2 , size_t x , size_t y , double elev , bool isY ) {
if ( b1 > b2 ) {
uint32_t b3 = b1 ;
b1 = b2 ; b2 = b3 ;
}
uint64_t key = get_key ( b1 , b2 ) ;
auto item = linkmap . find ( get_key ( b1 , b2 ) ) ;
if ( item = = linkmap . end ( ) ) {
Link * ln = & linklist [ n + + ] ;
* ln = { b1 , b2 , x , y , elev , isY } ;
linkmap . insert ( { key , ln } ) ;
} else {
Link * ln = item - > second ;
if ( elev < ln - > elev ) {
ln - > x = x ;
ln - > y = y ;
ln - > elev = elev ;
ln - > isY = isY ;
}
}
}
} ;
struct Adjacency { size_t next ; size_t nbasin ; } ;
struct Basin { size_t begin ; size_t size ; } ;
void planar_boruvka ( Link * linklist , Link * * graph , const size_t nlinks , const size_t nbasins )
{
//cout << "Building adjacency list" << endl;
/*
adjacencylist : List that stores all links of linklist , in the same order , but duplicated ( one for both directions ) . The ' . next ' parameters points to another link coming from the same node .
basinlist : first edge of every node , and number of edges .
*/
Adjacency * adjacencylist = new Adjacency [ nlinks * 2 ] ;
Basin * basinlist = new Basin [ nbasins ] ;
bool * active_edges = new bool [ nlinks ] ;
size_t * bucket = new size_t [ nbasins ] ;
for ( int i = 0 ; i < nbasins ; i + + )
basinlist [ i ] = { ( size_t ) - 1 , 0 } ;
for ( size_t i = 0 ; i < nlinks ; i + + ) { // Fill the adjacency list
active_edges [ i ] = true ;
Link & link = linklist [ i ] ;
// Add link from basin 1 to basin 2
Basin * basin = & basinlist [ link . basins [ 0 ] ] ;
adjacencylist [ i * 2 ] = { basin - > begin , link . basins [ 1 ] } ;
basin - > begin = i * 2 ;
basin - > size + + ;
// Add link from basin 2 to basin 1
basin = & basinlist [ link . basins [ 1 ] ] ;
adjacencylist [ i * 2 + 1 ] = { basin - > begin , link . basins [ 0 ] } ;
basin - > begin = i * 2 + 1 ;
basin - > size + + ;
}
vector < size_t > lowlevel ;
vector < size_t > highlevel ;
for ( size_t i = 0 ; i < nbasins ; i + + ) // Select basins based on their number of links
if ( basinlist [ i ] . size < = 8 )
lowlevel . push_back ( i ) ;
else {
highlevel . push_back ( i ) ;
}
for ( size_t i = 0 ; i < nbasins ; i + + )
bucket [ i ] = ( size_t ) - 1 ;
vector < size_t > in_bucket ;
size_t g = 0 ;
while ( ! lowlevel . empty ( ) ) {
for ( const size_t b1 : lowlevel ) { // Process all basins with 8 edges or less
if ( basinlist [ b1 ] . size > 8 ) { // Check that the basin has less than 8 edges. It may have changed.
highlevel . push_back ( b1 ) ;
continue ;
}
size_t lowest_last_n ;
size_t lowest_n ;
double lowest_elev = numeric_limits < double > : : infinity ( ) ;
size_t n = basinlist [ b1 ] . begin ;
size_t nlast = ( size_t ) - 1 ;
size_t count = 0 ;
while ( n ! = ( size_t ) - 1 ) { // Loop over b1's edges
const size_t nlink = n > > 1 ; // Corresponding index in linklist
if ( ! active_edges [ nlink ] ) {
// If edge is inactive, remove it and skip
// Removing the nth edge: making the n-1th edge pointing to the n+1th one.
if ( nlast = = ( size_t ) - 1 )
basinlist [ b1 ] . begin = adjacencylist [ n ] . next ;
else
adjacencylist [ nlast ] . next = adjacencylist [ n ] . next ;
n = adjacencylist [ n ] . next ; // Move to next link
continue ;
}
// If it is the lowest edge yet
if ( linklist [ nlink ] . elev < lowest_elev ) {
lowest_elev = linklist [ nlink ] . elev ;
lowest_last_n = nlast ;
lowest_n = n ;
}
count + + ;
nlast = n ;
n = adjacencylist [ n ] . next ; // Move to next link
}
if ( count ! = basinlist [ b1 ] . size )
cout < < " Inconsistent size for node B " < < b1 < < " ;: " < < count < < " links instead of " < < basinlist [ b1 ] . size < < endl ;
if ( ! isfinite ( lowest_elev ) ) {
// If node is empty! (Should happen only for last node)
continue ;
}
// Add link to the graph
graph [ g + + ] = & linklist [ lowest_n > > 1 ] ;
const size_t b2 = adjacencylist [ lowest_n ] . nbasin ;
// We are going to merge basin b1 into b2.
// First, edges that were pointing to basin b1 should be redirected to b2
n = basinlist [ b1 ] . begin ;
while ( n ! = ( size_t ) - 1 ) { // Loop again over the b1's edges
if ( adjacencylist [ n ] . nbasin = = b2 ) { // Links that point to b2 should be removed
active_edges [ n > > 1 ] = false ; // Mark the link as disused, it will be removed when encountered
basinlist [ b1 ] . size - - ; // Decrement basin sizes (applies to both basins)
basinlist [ b2 ] . size - - ;
}
// Redirect reciprocal edges to b2
adjacencylist [ n ^ ( size_t ) 1 ] . nbasin = b2 ; // n^1 corresponds to the reciprocal edge (because pairs of edges are stored next to each other)
n = adjacencylist [ n ] . next ;
}
adjacencylist [ nlast ] . next = basinlist [ b2 ] . begin ; // Append b2 after b1 by referencing last b1 member on first b2 one
basinlist [ b2 ] . begin = basinlist [ b1 ] . begin ; // Merge b2 and b1
basinlist [ b2 ] . size + = basinlist [ b1 ] . size ; // Increment b2's size
basinlist [ b1 ] = { ( size_t ) - 1 , 0 } ; // Delete b1
}
// Clean the graph, to remove many duplicate or self-loop edges that may appear
size_t cur_highlevel = 0 ;
lowlevel . clear ( ) ;
if ( highlevel . size ( ) = = 0 )
break ;
for ( const size_t b1 : highlevel ) { // High level nodes are the only ones remaining
in_bucket . clear ( ) ;
size_t n = basinlist [ b1 ] . begin ;
size_t j = 0 ;
while ( n ! = ( size_t ) - 1 ) { // Loop over edges
const size_t b2 = adjacencylist [ n ] . nbasin ;
if ( ! active_edges [ n > > 1 ] ) { // Skip disused edges
n = adjacencylist [ n ] . next ;
continue ;
}
size_t jref = bucket [ b2 ] ; // Check the content of the bucket for the opposite basin. If that basin has already been seen, it contains an index to the corresponding item in 'in_bucket'.
if ( jref = = ( size_t ) - 1 ) { // First time a link to this basin is encountered, keep it
bucket [ b2 ] = j + + ;
in_bucket . push_back ( n ) ;
} else { // There is already a link to this basin
if ( linklist [ n > > 1 ] . elev < linklist [ in_bucket [ jref ] > > 1 ] . elev ) { // Check elevation: if the new link found is lower, replace the older, otherwise, ignore it.
active_edges [ in_bucket [ jref ] > > 1 ] = false ;
in_bucket [ jref ] = n ;
} else {
active_edges [ n > > 1 ] = false ;
}
}
n = adjacencylist [ n ] . next ; // Move to next edge
}
size_t nlast = ( size_t ) - 1 ;
if ( in_bucket . size ( ) = = 0 ) {
basinlist [ b1 ] = { ( size_t ) - 1 , 0 } ;
continue ;
}
// Reconstruct adjacency of b1
for ( const size_t na : in_bucket ) { // na to avoid confusion with n
bucket [ adjacencylist [ na ] . nbasin ] = ( size_t ) - 1 ;
if ( nlast = = ( size_t ) - 1 ) {
basinlist [ b1 ] . begin = na ;
} else
adjacencylist [ nlast ] . next = na ;
nlast = na ;
}
adjacencylist [ nlast ] . next = ( size_t ) - 1 ;
// Recompute size and push into appropriate list
basinlist [ b1 ] . size = in_bucket . size ( ) ;
if ( ( basinlist [ b1 ] . size < = 8 ) & & ( basinlist [ b1 ] . size > 0 ) )
lowlevel . push_back ( b1 ) ;
else {
highlevel [ cur_highlevel + + ] = b1 ;
}
}
highlevel . resize ( cur_highlevel ) ;
}
delete [ ] adjacencylist ;
delete [ ] basinlist ;
delete [ ] active_edges ;
delete [ ] bucket ;
}
inline double disptime ( clock_t delta ) {
double time = ( double ) delta / ( double ) CLOCKS_PER_SEC ;
return time ;
}
void RiverMapper : : flow ( )
{
clock_t t0 , t1 ;
clock_t t_total = 0 ;
cout < < " Calculating local flow directions... " < < endl ;
t0 = clock ( ) ;
// Reinitialize array
for ( size_t i = 0 ; i < size ; i + + ) {
dirs [ i ] = 0u ;
}
const size_t Xmax = width - 1 ;
const size_t Ymax = height - 1 ;
vector < size_t > singular ;
// Determine local flow dirs
size_t i = 0 ;
double zdiffs [ 4 ] ;
for ( size_t y = 0 ; y < height ; y + + ) {
for ( size_t x = 0 ; x < width ; x + + , i + + ) {
double z = dem [ i ] ;
zdiffs [ 0 ] = ( x < Xmax ) ? diff_or_zero ( z , dem [ i + 1 ] ) : 0 ;
zdiffs [ 1 ] = ( y < Ymax ) ? diff_or_zero ( z , dem [ i + width ] ) : 0 ;
zdiffs [ 2 ] = ( x > 0 ) ? diff_or_zero ( z , dem [ i - 1 ] ) : 0 ;
zdiffs [ 3 ] = ( y > 0 ) ? diff_or_zero ( z , dem [ i - width ] ) : 0 ;
uint8_t d = flow_local ( zdiffs , 4u ) ;
switch ( d ) {
case 1 : {
dirs [ i + 1 ] | = ( uint8_t ) 1u ;
dirs_ref [ i ] = i + 1 ;
break ;
}
case 2 : {
dirs [ i + width ] | = ( uint8_t ) 2u ;
dirs_ref [ i ] = i + width ;
break ;
}
case 3 : {
dirs [ i - 1 ] | = ( uint8_t ) 4u ;
dirs_ref [ i ] = i - 1 ;
break ;
}
case 4 : {
dirs [ i - width ] | = ( uint8_t ) 8u ;
dirs_ref [ i ] = i - width ;
break ;
}
default : {
singular . push_back ( i ) ;
dirs_ref [ i ] = SIZE_MAX ;
}
}
}
}
t1 = clock ( ) ;
t_total + = ( t1 - t0 ) ;
cout < < " \t Completed in " < < disptime ( t1 - t0 ) < < " s " < < endl ;
cout < < " Percolating drainage basins ( " < < singular . size ( ) < < " basins)... " < < endl ;
t0 = clock ( ) ;
uint32_t nb = 0 ;
std : : stack < size_t > s ;
for ( size_t i : singular ) {
s . push ( i ) ;
while ( ! s . empty ( ) ) {
i = s . top ( ) ;
s . pop ( ) ;
uint8_t d = dirs [ i ] ;
if ( d & ( uint8_t ) 1u )
s . push ( i - 1 ) ;
if ( d & ( uint8_t ) 2u )
s . push ( i - width ) ;
if ( d & ( uint8_t ) 4u )
s . push ( i + 1 ) ;
if ( d & ( uint8_t ) 8u )
s . push ( i + width ) ;
basin_id [ i ] = nb ;
}
nb + + ;
}
const size_t basin_max = nb ;
t1 = clock ( ) ;
t_total + = ( t1 - t0 ) ;
cout < < " \t Completed in " < < disptime ( t1 - t0 ) < < " s " < < endl ;
cout < < " Finding lowest passes between basins... " < < endl ;
t0 = clock ( ) ;
LinkManager linkmgr ( basin_max + 1 ) ;
// Link basins
i = 0 ;
for ( size_t y = 0 ; y < height ; y + + ) {
uint32_t b1 = basin_id [ i ] ;
linkmgr . add_link ( b1 , basin_max , 0 , y , dem [ i ] , false ) ;
i + + ;
for ( size_t x = 1 ; x < width ; x + + , i + + ) {
const uint32_t b2 = basin_id [ i ] ;
if ( b1 ! = b2 )
linkmgr . add_link ( b1 , b2 , x , y , max ( dem [ i ] , dem [ i - 1 ] ) , false ) ;
b1 = b2 ;
}
linkmgr . add_link ( b1 , basin_max , width , y , dem [ i - 1 ] , false ) ;
}
for ( size_t x = 0 ; x < width ; x + + ) {
i = x ;
uint32_t b1 = basin_id [ i ] ;
linkmgr . add_link ( b1 , basin_max , x , 0 , dem [ i ] , true ) ;
i + = width ;
for ( size_t y = 1 ; y < height ; y + + , i + = width ) {
const uint32_t b2 = basin_id [ i ] ;
if ( b1 ! = b2 )
linkmgr . add_link ( b1 , b2 , x , y , max ( dem [ i ] , dem [ i - width ] ) , true ) ;
b1 = b2 ;
}
linkmgr . add_link ( b1 , basin_max , x , height , dem [ i - width ] , true ) ;
}
t1 = clock ( ) ;
t_total + = ( t1 - t0 ) ;
cout < < " \t Completed in " < < disptime ( t1 - t0 ) < < " s " < < endl ;
cout < < " Building basin flow tree... " < < endl ;
t0 = clock ( ) ;
Link * * basin_tree = new Link * [ basin_max ] ;
planar_boruvka ( linkmgr . linklist , basin_tree , linkmgr . n , basin_max + 1 ) ; // This computes basin tree
t1 = clock ( ) ;
t_total + = ( t1 - t0 ) ;
cout < < " \t Completed in " < < disptime ( t1 - t0 ) < < " s " < < endl ;
cout < < " Orienting basin tree and updating flow directions... " < < endl ;
t0 = clock ( ) ;
// Compute adjacency list of the basin tree
size_t * adjacencylist = new size_t [ basin_max * 2 ] ;
size_t * basinlist = new size_t [ basin_max + 1 ] ;
for ( size_t i = 0 ; i < basin_max + 1 ; i + + )
basinlist [ i ] = ( size_t ) - 1 ;
for ( size_t i = 0 ; i < basin_max ; i + + ) {
Link & link = * basin_tree [ i ] ;
adjacencylist [ i * 2 ] = basinlist [ link . basins [ 0 ] ] ;
basinlist [ link . basins [ 0 ] ] = i * 2 ;
adjacencylist [ i * 2 + 1 ] = basinlist [ link . basins [ 1 ] ] ;
basinlist [ link . basins [ 1 ] ] = i * 2 + 1 ;
}
// Resolve links
double * lake_elev = new double [ basin_max + 1 ] ;
lake_elev [ basin_max ] = - numeric_limits < double > : : infinity ( ) ;
stack < size_t > basin_stack ;
basin_stack . push ( basin_max ) ;
while ( ! basin_stack . empty ( ) ) {
const size_t b1 = basin_stack . top ( ) ;
const double b1_lake_elev = lake_elev [ b1 ] ;
basin_stack . pop ( ) ;
size_t n = basinlist [ b1 ] ;
while ( n ! = ( size_t ) - 1 ) {
if ( basin_tree [ n > > 1 ] = = nullptr ) {
n = adjacencylist [ n ] ;
continue ;
}
Link & link = * basin_tree [ n > > 1 ] ;
const size_t b2 = ( link . basins [ 0 ] = = b1 ) ? link . basins [ 1 ] : link . basins [ 0 ] ;
basin_stack . push ( b2 ) ;
lake_elev [ b2 ] = max ( b1_lake_elev , link . elev ) ;
// Resolve link
size_t offset = link . isY ? width : 1 ;
size_t i1 = link . y * width + link . x ;
bool forward = ( link . x = = width ) | | ( link . y = = height ) | | ( basin_id [ i1 ] ! = b2 ) ;
bool border = link . isY ? ( ( link . y = = 0 ) | | ( link . y = = height ) ) : ( ( link . x = = 0 ) | | ( link . x = = width ) ) ;
size_t i2 ;
if ( forward ) {
i2 = i1 ;
i1 - = offset ;
} else {
i2 = i1 - offset ;
}
if ( border ) {
i2 = SIZE_MAX ;
}
while ( i1 ! = SIZE_MAX ) {
if ( basin_id [ i1 ] ! = b2 )
cout < < " Basin " < < basin_id [ i1 ] < < " instead of " < < b2 < < " ! " < < endl ;
const size_t temp = dirs_ref [ i1 ] ;
dirs_ref [ i1 ] = i2 ;
i2 = i1 ;
i1 = temp ;
}
basin_tree [ n > > 1 ] = nullptr ; // To avoid looping twice on the same link
n = adjacencylist [ n ] ;
}
}
for ( size_t i = 0 ; i < size ; i + + ) {
lakes [ i ] = lake_elev [ basin_id [ i ] ] ;
const int dir_diff = dirs_ref [ i ] - i ;
uint8_t d = 0u ;
if ( dir_diff = = - width )
d = 1u ;
else if ( dir_diff = = 1 )
d = 2u ;
else if ( dir_diff = = width )
d = 3u ;
else if ( dir_diff = = - 1 )
d = 4u ;
dirs [ i ] = d ;
}
delete [ ] basin_tree ;
delete [ ] adjacencylist ;
delete [ ] basinlist ;
delete [ ] lake_elev ;
t1 = clock ( ) ;
t_total + = ( t1 - t0 ) ;
cout < < " \t Completed in " < < disptime ( t1 - t0 ) < < " s " < < endl ;
cout < < " Accumulating water... " < < endl ;
t0 = clock ( ) ;
for ( size_t i = 0 ; i < size ; i + + )
water [ i ] = 1 ;
accumulate ( ) ;
t1 = clock ( ) ;
t_total + = ( t1 - t0 ) ;
cout < < " \t Completed in " < < disptime ( t1 - t0 ) < < " s " < < endl ;
cout < < " Full flow calculation completed in " < < disptime ( t_total ) < < " s " < < endl ;
}
void RiverMapper : : accumulate ( )
{
for ( size_t i = 0 ; i < size ; i + + )
ndonors [ i ] = 0 ;
for ( size_t i = 0 ; i < size ; i + + ) {
size_t d = dirs_ref [ i ] ;
if ( d = = SIZE_MAX )
continue ;
ndonors [ d ] + + ;
}
for ( size_t i = 0 ; i < size ; i + + ) {
if ( ndonors [ i ] > 0 )
continue ;
double * w1 = & water [ i ] ;
size_t iw = i ;
while ( true ) {
iw = dirs_ref [ iw ] ;
if ( iw = = SIZE_MAX ) {
break ;
}
double * w2 = & water [ iw ] ;
* w2 + = * w1 ; // Flow w1 into w2
w1 = w2 ; // Now we will follow w2
if ( ndonors [ iw ] > 1 ) {
ndonors [ iw ] - - ;
break ;
}
}
}
}
2020-12-27 09:41:36 -08:00
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 < < " \t Completed in " < < disptime ( t1 - t0 ) < < " s " < < endl ;
}