Render sediment differently.

This commit is contained in:
Joshua Yanovski 2019-11-20 11:33:53 +01:00
parent e289849c87
commit e01bb8a8c9
3 changed files with 168 additions and 34 deletions

View File

@ -436,13 +436,14 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
(temp < CONFIG.snow_temp ||
downhill_alt.sub(CONFIG.sea_level) >= CONFIG.mountain_scale * 0.25)*/
is_rocky {
sim.get_interpolated_monotone(wpos, |chunk| chunk.alt)?
// sim.get_interpolated_monotone(wpos, |chunk| chunk.alt)?
// sim.get_interpolated_bilinear(wpos, |chunk| chunk.alt)?
// sim.get_interpolated(wpos, |chunk| chunk.alt)?
sim.get_interpolated(wpos, |chunk| chunk.alt)?
} else {
sim.get_interpolated_monotone(wpos, |chunk| chunk.alt)?
// sim.get_interpolated(wpos, |chunk| chunk.alt)?
// sim.get_interpolated_monotone(wpos, |chunk| chunk.alt)?
sim.get_interpolated(wpos, |chunk| chunk.alt)?
};
let basement = sim.get_interpolated/*get_interpolated_monotone*/(wpos, |chunk| chunk.basement)?;
// Find the average distance to each neighboring body of water.
let mut river_count = 0.0f64;
@ -828,6 +829,7 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
let riverless_alt_delta = Lerp::lerp(0.0, riverless_alt_delta, warp_factor);
let alt = alt_ + riverless_alt_delta;
let basement = basement.min(alt);
let rock = (sim.gen_ctx.small_nz.get(
Vec3::new(wposf.x, wposf.y, alt as f64)
@ -917,6 +919,7 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
// For below desert humidity, we are always sand or rock, depending on altitude and
// temperature.
let ground = Rgb::lerp(
cliff,
Rgb::lerp(
dead_tundra,
sand,
@ -924,10 +927,11 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
.div(CONFIG.desert_temp.sub(CONFIG.snow_temp))
.mul(/*4.5*/ 0.5),
),
cliff,
alt.sub(CONFIG.sea_level)
alt.sub(basement)
.mul(0.25)
/* alt.sub(CONFIG.sea_level)
.sub(CONFIG.mountain_scale * 0.25)
.div(CONFIG.mountain_scale * 0.125),
.div(CONFIG.mountain_scale * 0.125), */
);
// From desert to forest humidity, we go from tundra to dirt to grass to moss to sand,
// depending on temperature.
@ -1091,14 +1095,15 @@ impl<'a> Sampler<'a> for ColumnGen<'a> {
water_level,
warp_factor,
surface_color: Rgb::lerp(
sand,
Rgb::lerp(cliff, sand, alt.sub(basement).mul(0.25)),
// Land
ground,
// Beach
((ocean_level - 1.0) / 2.0).max(0.0),
),
sub_surface_color: /*warm_grass*/dirt,
tree_density,
sub_surface_color: /*warm_grass*/Lerp::lerp(cliff, dirt, alt.sub(basement).mul(0.25)),
// No growing directly on bedrock.
tree_density: if alt - basement < 2.0 { 0.0 } else { tree_density },
forest_kind: sim_chunk.forest_kind,
close_structures: self.gen_close_structures(wpos),
cave_xy,

View File

@ -521,7 +521,7 @@ fn get_max_slope(h: &[f32], rock_strength_nz: &(impl NoiseFn<Point3<f64>> + Sync
// sure anyone would be able to usefully tweak them on a per-map basis? Plus it's just a
// hacky heuristic anyway.
let center = 0.25;
let dmin = center - 0.15;//0.05;
let dmin = center - /*0.15;//0.05*/0.05;
let dmax = center + 0.05;//0.05;
let log_odds = |x: f64| logit(x) - logit(center);
let rock_strength = logistic_cdf(
@ -592,6 +592,7 @@ fn get_max_slope(h: &[f32], rock_strength_nz: &(impl NoiseFn<Point3<f64>> + Sync
fn erode(
h: &mut [f32],
b: &mut [f32],
wh: &mut [f32],
is_done: &mut BitBox,
done_val: bool,
erosion_base: f32,
@ -625,7 +626,7 @@ fn erode(
let dt = max_uplift as f64 / height_scale /* * CONFIG.mountain_scale as f64*/ / 5.010e-4;
println!("dt={:?}", dt);
// Landslide constant: ideally scaled to 10e-2 m / y^-1
let l = /*200.0 * max_uplift as f64;*/10.0e-2 /*/ CONFIG.mountain_scale as f64*/ * height_scale * dt;
let l = /*200.0 * max_uplift as f64;*/1.0e-2 /*/ CONFIG.mountain_scale as f64*/ * height_scale * dt;
// Net precipitation rate (m / year)
let p = 1.0 * height_scale;
// Stream power erosion constant (bedrock), in m^(1-2m) / year (times dt).
@ -638,12 +639,13 @@ fn erode(
// 8e-6 * dt
// 2e-5 * dt;
// 2.244/2048*5*32/(250000/4)*10^6
erosion_base as f64 + 2.244 / mmaxh as f64 * 5.0 * max_uplift as f64;
// ln(tan(30/360*2*pi))-ln(tan(6/360*2*pi))*1500 = 3378
erosion_base as f64 + 2.244 / mmaxh as f64 * /*10.0*//*5.0*//*9.0*//*7.5*/5.0/*3.75*/ * max_uplift as f64;
// 1e-5 * dt;
// (2.244*(5.010e-4)/512)/(2.244*(5.010e-4)/2500) = 4.88...
// 2.444 * 5
// Stream power erosion constant (sediment), in m^(1-2m) / year (times dt).
let k_fs = k_fb * 1.0/* 2.0*//*4.0*/;
let k_fs = k_fb * /*1.0*//*2.0*/1.0/*4.0*/;
let ((dh, indirection, newh, area), mut max_slopes) = rayon::join(
|| {
let mut dh = downhill(h, |posi| is_ocean(posi));
@ -677,11 +679,12 @@ fn erode(
if posj == -1 {
panic!("Disconnected lake!");
}
// wh for oceans is always 0.
wh[posi] = 0.0;
// max_slopes[posi] = kd(posi);
// Egress with no outgoing flows.
} else {
// *is_done.at(posi) = done_val;
nland += 1;
let posj = posj as usize;
let dxy = (uniform_idx_as_vec2(posi) - uniform_idx_as_vec2(posj)).map(|e| e as f64);
@ -694,8 +697,9 @@ fn erode(
// height_scale^2, to make sure we were working in the correct units for dz
// (which has height_scale height unit = 1.0 meters).
let old_h_i = h[posi] as f64;
let old_b_i = b[posi] as f64;
let uplift_i = uplift(posi) as f64;
let k = if (old_h_i - b[posi] as f64) > 1.0e-6 {
let k = if (old_h_i - old_b_i as f64) > 1.0e-6 {
// Sediment
k_fs
} else {
@ -706,13 +710,30 @@ fn erode(
let flux = k * (p * chunk_area * area[posi] as f64).sqrt() / neighbor_distance;
// h[i](t + dt) = (h[i](t) + δt * (uplift[i] + flux(i) * h[j](t + δt))) / (1 + flux(i) * δt).
// NOTE: posj has already been computed since it's downhill from us.
// Therefore, we can rely on wh being set to the water height for that node.
let h_j = h[posj] as f64;
let wh_j = wh[posj] as f64;
// hi(t + ∂t) = (hi(t) + ∂t(ui + kp^mAi^m(hj(t + ∂t)/||pi - pj||))) / (1 + ∂t * kp^mAi^m / ||pi - pj||)
let mut new_h_i = (old_h_i + (uplift_i + flux * h_j)) / (1.0 + flux);
let mut new_h_i = old_h_i + uplift_i;
// Only perform erosion if we are above the water level of the previous node.
if new_h_i > wh_j {
new_h_i = (new_h_i + flux * h_j) / (1.0 + flux);
// If we dipped below the receiver's water level, set our height to the receiver's
// water level.
if new_h_i <= wh_j {
new_h_i = wh_j;
} else {
nland += 1;
sumh += new_h_i;
}
}
// Set max_slope to this node's water height (max of receiver's water height and
// this node's height).
wh[posi] = wh_j.max(new_h_i) as f32;
// Prevent erosion from dropping us below our receiver, unless we were already below it.
// new_h_i = h_j.min(old_h_i + uplift_i).max(new_h_i);
// Find out if this is a lake bottom.
let indirection_idx = indirection[posi];
/* let indirection_idx = indirection[posi];
let is_lake_bottom = indirection_idx < 0;
let _fake_neighbor = is_lake_bottom && dxy.x.abs() > 1.0 && dxy.y.abs() > 1.0;
// Test the slope.
@ -736,12 +757,109 @@ fn erode(
let dxy = (uniform_idx_as_vec2(posi) - uniform_idx_as_vec2(posk)).map(|e| e as f64);
let neighbor_distance = (neighbor_coef * dxy).magnitude();
let dz = (new_h_i - /*h_j*/h_k).max(0.0) / height_scale/* * CONFIG.mountain_scale as f64*/;
let mag_slope = dz/*.abs()*/ / neighbor_distance;
let mag_slope = dz/*.abs()*/ / neighbor_distance; */
// If you're on the lake bottom and not right next to your neighbor, don't compute a
// slope.
if
/* !is_lake_bottom */ /* !fake_neighbor */
true {
/* if
/* !is_lake_bottom && */
mag_slope > max_slope {
// println!("old slope: {:?}, new slope: {:?}, dz: {:?}, h_j: {:?}, new_h_i: {:?}", mag_slope, max_slope, dz, h_j, new_h_i);
// Thermal erosion says this can't happen, so we reduce dh_i to make the slope
// exactly max_slope.
// max_slope = (old_h_i + dh - h_j) / height_scale/* * CONFIG.mountain_scale */ / NEIGHBOR_DISTANCE
// dh = max_slope * NEIGHBOR_DISTANCE * height_scale/* / CONFIG.mountain_scale */ + h_j - old_h_i.
let dh = max_slope * neighbor_distance * height_scale/* / CONFIG.mountain_scale as f64*/;
new_h_i = /*h_j.max*/(h_k + dh).max(new_h_i - l * (mag_slope - max_slope));
let dz = (new_h_i - /*h_j*/h_k).max(0.0) / height_scale/* * CONFIG.mountain_scale as f64*/;
let slope = dz/*.abs()*/ / neighbor_distance;
sums += slope;
/* max_slopes[posi] = /*(mag_slope - max_slope) * */kd(posi);
sums += mag_slope; */
// let slope = dz.signum() * max_slope;
// new_h_i = slope * neighbor_distance * height_scale /* / CONFIG.mountain_scale as f64 */ + h_j;
// sums += max_slope;
} else {
// max_slopes[posi] = 0.0;
sums += mag_slope;
// Just use the computed rate.
} */
h[posi] = new_h_i as f32;
// Make sure to update the basement as well!
b[posi] = (old_b_i + uplift_i).min(new_h_i) as f32;
}
}
// *is_done.at(posi) = done_val;
maxh = h[posi].max(maxh);
}
log::info!(
"Done applying stream power (max height: {:?}) (avg height: {:?}) (avg slope: {:?})",
maxh,
if nland == 0 {
f64::INFINITY
} else {
sumh / nland as f64
},
if nland == 0 {
f64::INFINITY
} else {
sums / nland as f64
},
);
// Apply thermal erosion.
maxh = 0.0;
sumh = 0.0;
for &posi in &*newh {
let posi = posi as usize;
let posj = dh[posi];
if posj < 0 {
if posj == -1 {
panic!("Disconnected lake!");
}
// max_slopes[posi] = kd(posi);
// Egress with no outgoing flows.
} else {
let posj = posj as usize;
let old_h_i = h[posi] as f64;
let old_b_i = b[posi] as f64;
// Find the water height for this chunk's receiver; we only apply thermal erosion
// for chunks above water.
let wh_j = wh[posj] as f64;
// If you're on the lake bottom and not right next to your neighbor, don't compute a
// slope.
if
/* !is_lake_bottom */ /* !fake_neighbor */
wh_j < old_h_i {
let mut new_h_i = old_h_i;
let h_j = h[posj] as f64;
/* let indirection_idx = indirection[posi];
let is_lake_bottom = indirection_idx < 0;
let _fake_neighbor = is_lake_bottom && dxy.x.abs() > 1.0 && dxy.y.abs() > 1.0; */
// Test the slope.
let max_slope = max_slopes[posi] as f64;
// Hacky version of thermal erosion: only consider lowest neighbor, don't redistribute
// uplift to other neighbors.
let (posk, h_k) = /* neighbors(posi)
.filter(|&posk| *is_done.at(posk) == done_val)
// .filter(|&posk| *is_done.at(posk) == done_val || is_ocean(posk))
.map(|posk| (posk, h[posk] as f64))
// .filter(|&(posk, h_k)| *is_done.at(posk) == done_val || h_k < 0.0)
.min_by(|&(_, a), &(_, b)| a.partial_cmp(&b).unwrap())
.unwrap_or((posj, h_j)); */
(posj, h_j);
// .max(h_j);
let (posk, h_k) = if h_k < h_j {
(posk, h_k)
} else {
(posj, h_j)
};
let dxy = (uniform_idx_as_vec2(posi) - uniform_idx_as_vec2(posk)).map(|e| e as f64);
let neighbor_distance = (neighbor_coef * dxy).magnitude();
let dz = (new_h_i - /*h_j*/h_k).max(0.0) / height_scale/* * CONFIG.mountain_scale as f64*/;
let mag_slope = dz/*.abs()*/ / neighbor_distance;
if
/* !is_lake_bottom && */
mag_slope > max_slope {
@ -767,7 +885,7 @@ fn erode(
}
h[posi] = new_h_i as f32;
// Make sure to update the basement as well!
b[posi] = (old_h_i + uplift_i).min(new_h_i) as f32;
b[posi] = old_b_i.min(new_h_i) as f32;
sumh += new_h_i;
}
}
@ -775,7 +893,7 @@ fn erode(
maxh = h[posi].max(maxh);
}
log::debug!(
"Done applying stream power (max height: {:?}) (avg height: {:?}) (avg slope: {:?})",
"Done applying thermal erosion (max height: {:?}) (avg height: {:?}) (avg slope: {:?})",
maxh,
if nland == 0 {
f64::INFINITY
@ -1255,6 +1373,7 @@ pub fn do_erosion(
.map(|posi| oldb(posi))
.collect::<Vec<_>>()
.into_boxed_slice();
let mut wh = vec![0.0; WORLD_SIZE.x * WORLD_SIZE.y].into_boxed_slice();
// TODO: Don't do this, maybe?
// (To elaborate, maybe we should have varying uplift or compute it some other way).
let uplift = (0..oldh_.len())
@ -1282,8 +1401,15 @@ pub fn do_erosion(
// on land, but in theory we probably want to at least differentiate between soil, bedrock, and
// sediment.
let height_scale = 1.0; // 1.0 / CONFIG.mountain_scale as f64;
let kdsed = /*1.5e-2*//*1e-4*/1e-2 * height_scale * height_scale/* / (CONFIG.mountain_scale as f64 * CONFIG.mountain_scale as f64) */;
let kd_bedrock = /*1e-2*//*0.25e-2*/1e-2 * height_scale * height_scale/* / (CONFIG.mountain_scale as f64 * CONFIG.mountain_scale as f64) */;
let mmaxh = CONFIG.mountain_scale as f64 * height_scale;
let dt = max_uplift as f64 / height_scale /* * CONFIG.mountain_scale as f64*/ / 5.010e-4;
let k_fb = (erosion_base as f64 + 2.244 / mmaxh as f64 * /*10.0*//*5.0*//*9.0*//*7.5*/5.0/*3.75*/ * max_uplift as f64) / dt;
let kd_bedrock =
/*1e-2*//*0.25e-2*/1e-2 * height_scale * height_scale/* / (CONFIG.mountain_scale as f64 * CONFIG.mountain_scale as f64) */
/* * k_fb / 2e-5 */;
let kdsed =
/*1.5e-2*//*1e-4*//*1.25e-2*/1.5e-2 * height_scale * height_scale/* / (CONFIG.mountain_scale as f64 * CONFIG.mountain_scale as f64) */
/* * k_fb / 2e-5 */;
let kd = |posi: usize| if is_ocean(posi) { /*0.0*/kd_bedrock } else { kd_bedrock };
// Hillslope diffusion coefficient for sediment.
let mut is_done = bitbox![0; WORLD_SIZE.x * WORLD_SIZE.y];
@ -1292,6 +1418,7 @@ pub fn do_erosion(
erode(
&mut h,
&mut b,
&mut wh,
&mut is_done,
// The value to use to indicate that erosion is complete on a chunk. Should toggle
// once per iteration, to avoid having to reset the bits, and start at true, since

View File

@ -190,7 +190,7 @@ impl WorldSim {
let river_seed = RandomField::new(rng.gen());
let rock_lacunarity = 0.5/*HybridMulti::DEFAULT_LACUNARITY*/;
let rock_strength_nz = /*Fbm*/HybridMulti_/*BasicMulti*//*Fbm*/::new()
.set_octaves(/*6*//*5*//*4*//*5*/4)
.set_octaves(/*6*//*5*//*4*//*5*//*4*/6)
// persistence = lacunarity^(-(1.0 - fractal increment))
.set_persistence(/*0.9*/ /*2.0*//*1.5*//*HybridMulti::DEFAULT_LACUNARITY*/rock_lacunarity.powf(-(1.0 - 0.9)))
// 256*32/2^4
@ -199,7 +199,7 @@ impl WorldSim {
// (0.5^(-(1.0-0.9)))^1/256/32*2^4*256*4
// (2^(-(1.0-0.9)))^4
// 16.0
.set_frequency(/*0.9*/ /*Fbm*/HybridMulti_::DEFAULT_FREQUENCY / (/*64.0*/256.0/*1.0*//*16.0*/ * 32.0/* TerrainChunkSize::RECT_SIZE.x as f64 */))
.set_frequency(/*0.9*/ /*Fbm*/HybridMulti_::DEFAULT_FREQUENCY / (64.0/*256.0*//*1.0*//*16.0*/ * 32.0/* TerrainChunkSize::RECT_SIZE.x as f64 */))
// .set_persistence(/*0.9*/ /*2.0*/0.67)
// .set_frequency(/*0.9*/ Fbm::DEFAULT_FREQUENCY / (2.0 * 32.0))
// .set_lacunarity(0.5)
@ -216,12 +216,12 @@ impl WorldSim {
.set_scale(1.0 / rock_strength_scale); */
let height_scale = 1.0f64; // 1.0 / CONFIG.mountain_scale as f64;
let max_erosion_per_delta_t = 32.0 * height_scale;
let max_erosion_per_delta_t = /*32.0*//*128.0*/128.0 * height_scale;
let erosion_pow_low = /*0.25*//*1.5*//*2.0*//*0.5*//*4.0*//*0.25*//*1.0*//*2.0*//*1.5*//*1.5*//*0.35*//*0.43*//*0.5*//*0.45*//*0.37*/1.002;
let erosion_pow_high = /*1.5*//*1.0*//*0.55*//*0.51*//*2.0*/1.002;
let erosion_center = /*0.45*//*0.75*//*0.75*//*0.5*//*0.75*/0.5;
let n_steps = 150;//37/*100*/;//50;//50;//37;//50;//37; // /*37*//*29*//*40*//*150*/37; //150;//200;
let n_small_steps = 0;//8;//8; // 8
let n_steps = /*100*/25/*100*//*37*/;//150;//37/*100*/;//50;//50;//37;//50;//37; // /*37*//*29*//*40*//*150*/37; //150;//200;
let n_small_steps = 8;//8;//8; // 8
// fractal dimension should be between 0 and 0.9999...
// (lacunarity^octaves)^(-H) = persistence^(octaves)
@ -525,7 +525,7 @@ impl WorldSim {
/* let erosion_factor = |x: f64| logistic_cdf(logistic_base * if x <= /*erosion_center*/alt_old_center_uniform/*alt_old_center*/ { erosion_pow_low.ln() } else { erosion_pow_high.ln() } * log_odds(x))/*0.5 + (x - 0.5).signum() * ((x - 0.5).mul(2.0).abs(
).powf(erosion_pow).mul(0.5))*/; */
let erosion_factor = |x: f64| (/*if x <= /*erosion_center*/alt_old_center_uniform/*alt_old_center*/ { erosion_pow_low.ln() } else { erosion_pow_high.ln() } * */(/*exp_inverse_cdf*//*logit*/inv_func(x) - alt_exp_min_uniform) / (alt_exp_max_uniform - alt_exp_min_uniform))/*0.5 + (x - 0.5).signum() * ((x - 0.5).mul(2.0).abs(
).powf(erosion_pow).mul(0.5))*//*.powf(0.5)*//*.powf(1.5)*//*.powf(2.0)*/;
).powf(erosion_pow).mul(0.5))*//*.powf(0.5)*//*.powf(1.5)*/.powf(2.0);
let uplift_fn =
|posi| {
if is_ocean_fn(posi) {
@ -606,7 +606,7 @@ impl WorldSim {
// assert!(alt_main >= 0.0);
let (bump_factor, bump_max) = if
/*height < f32::EPSILON as f64 * 0.5/*false*/*/
false {
true {
(
/*(alt_main./*to_le_bytes()[7]*/to_bits() & 1) as f64*/
(alt_main / CONFIG.mountain_scale as f64 * 128.0).mul(0.1).powf(1.2) * /*(1.0 / CONFIG.mountain_scale as f64)*/(f32::EPSILON * 0.5) as f64,
@ -615,11 +615,13 @@ impl WorldSim {
} else {
(0.0, 0.0)
};
let height = 1.0f64;
// let height = 1.0f64;
// let height = 1.0 / 7.0f64;
let height = height
/* .mul(15.0 / 16.0)
.add(1.0 / 16.0) */
/* .mul(5.0 / 8.0)
.add(3.0 / 8.0) */
.mul(7.0 / 8.0)
.add(1.0 / 8.0)
.mul(max_erosion_per_delta_t)
@ -670,7 +672,7 @@ impl WorldSim {
/* + spring(alt_main.abs().powf(0.5).min(0.75).mul(60.0).sin(), 4.0)
.mul(0.045)*/)
};
old_height_uniform(posi)
old_height_uniform(posi)/*.powf(2.0)*/
/* // 0.0
// -/*CONFIG.sea_level / CONFIG.mountain_scale*//* 0.75 */1.0
// ((old_height(posi) - alt_old_min) as f64 / (alt_old_max - alt_old_min) as f64) as f32
@ -1485,8 +1487,8 @@ impl SimChunk {
_ => {}
}
// No trees in the ocean or with zero humidity (currently)
let tree_density = if is_underwater {
// No trees in the ocean, with zero humidity (currently), or directly on bedrock.
let tree_density = if is_underwater/* || alt - basement.min(alt) < 2.0 */ {
0.0
} else {
let tree_density = (gen_ctx.tree_nz.get((wposf.div(1024.0)).into_array()))