#include "dc_density.h"
#include "dc_noise.h"
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <limits>
// -------------------------------------------------------------------------
// Helpers
// -------------------------------------------------------------------------
static float clampf(float v, float lo, float hi) {
return (v < lo) ? lo : (v > hi) ? hi : v;
}
static float smoothstep(float t) {
return t * t * (3.0f - 2.0f * t);
}
// Seeded LCG for Voronoi cell generation
class SeededRNG {
public:
explicit SeededRNG(uint64_t seed) : state_(seed) {}
float randf_range(float lo, float hi) {
state_ = state_ * 6364136223846793005ULL + 1442695040888963407ULL;
uint32_t v = static_cast<uint32_t>(state_ >> 33);
float f = static_cast<float>(v) / static_cast<float>(0x7FFFFFFFu);
if (f > 1.0f) f = 1.0f;
return lo + f * (hi - lo);
}
private:
uint64_t state_;
};
// Hash-based procedural Voronoi: generates consistent jittered cell centers
// for any integer grid cell (ix, iz), tiling infinitely.
static uint32_t hash2d(int ix, int iz, int seed_offset) {
uint32_t h = static_cast<uint32_t>(ix * 374761393 + iz * 668265263 + seed_offset);
h = (h ^ (h >> 13)) * 1274126177u;
h = h ^ (h >> 16);
return h;
}
// Returns 0..1 float from hash
static float hash_float(int ix, int iz, int seed_offset) {
return static_cast<float>(hash2d(ix, iz, seed_offset) & 0xFFFFFFu) / 16777216.0f;
}
// -------------------------------------------------------------------------
// Noise pimpl — owned by DCDensity via unique_ptr
// -------------------------------------------------------------------------
struct DCDensityNoise {
SimplexNoise primary;
SimplexNoise tunnel_a;
SimplexNoise tunnel_b;
SimplexNoise ridge;
SimplexNoise warp;
SimplexNoise detail;
};
// -------------------------------------------------------------------------
// DCDensity lifecycle
// -------------------------------------------------------------------------
DCDensity::DCDensity() = default;
DCDensity::~DCDensity() = default;
DCDensity::DCDensity(DCDensity&&) noexcept = default;
DCDensity& DCDensity::operator=(DCDensity&&) noexcept = default;
// -------------------------------------------------------------------------
// DCDensity implementation
// -------------------------------------------------------------------------
// Compute max octaves for given base frequency before hitting Nyquist limit.
// With lacunarity=2, octave n has freq = base * 2^(n-1).
// Nyquist limit = 0.5 / voxel_size. At least 1 octave always.
static int nyquist_octaves(float base_freq, float voxel_size, int max_octaves) {
float nyquist = 0.5f / voxel_size;
int oct = 0;
float freq = base_freq;
while (oct < max_octaves && freq <= nyquist) {
oct++;
freq *= 2.0f;
}
return std::max(oct, 1);
}
void DCDensity::configure(const DensityConfig& config) {
config_ = config;
// Pre-compute Nyquist-clamped octave counts.
// Spherical mode needs many more octaves to span planet-scale (100km+)
// down to ground-scale (1m) detail. Each octave doubles frequency.
// The Nyquist clamp per-LOD ensures coarse chunks skip fine octaves
// (e.g., LOD 10 at 512m voxels only evaluates ~9 of 20 octaves).
float vs = config.voxel_size;
int max_primary = config.spherical ? 20 : 3;
int max_tunnel = config.spherical ? 12 : 2;
int max_ridge = config.spherical ? 16 : 2;
int max_warp = config.spherical ? 10 : 2;
int max_detail = config.spherical ? 12 : 2;
// User-specified octave cap overrides per-layer maximums
if (config.max_noise_octaves > 0) {
max_primary = std::min(max_primary, config.max_noise_octaves);
max_tunnel = std::min(max_tunnel, config.max_noise_octaves);
max_ridge = std::min(max_ridge, config.max_noise_octaves);
max_warp = std::min(max_warp, config.max_noise_octaves);
max_detail = std::min(max_detail, config.max_noise_octaves);
}
// Separate ridge octave cap (0 = use the general cap above)
if (config.max_ridge_octaves > 0) {
max_ridge = std::min(max_ridge, config.max_ridge_octaves);
}
primary_octaves_ = nyquist_octaves(config.primary_frequency, vs, max_primary);
tunnel_octaves_ = nyquist_octaves(config.tunnel_frequency, vs, max_tunnel);
ridge_octaves_ = nyquist_octaves(config.ridge_frequency, vs, max_ridge);
warp_octaves_ = nyquist_octaves(config.warp_frequency, vs, max_warp);
detail_octaves_ = nyquist_octaves(config.detail_frequency, vs, max_detail);
// Create or reset noise generators
noise_ = std::make_unique<DCDensityNoise>();
noise_->primary = SimplexNoise(config.seed);
noise_->tunnel_a = SimplexNoise(config.seed + 10);
noise_->tunnel_b = SimplexNoise(config.seed + 20);
noise_->ridge = SimplexNoise(config.seed + 30);
noise_->warp = SimplexNoise(config.seed + 40);
noise_->detail = SimplexNoise(config.seed + 50);
// Generate Voronoi cells (jittered grid with random floor heights)
voronoi_cells_.clear();
float half = config.terrain_size / 2.0f;
float cell_spacing = config.terrain_size / static_cast<float>(config.voronoi_cell_count);
SeededRNG rng(static_cast<uint64_t>(config.seed + 100));
for (int cy = 0; cy < config.voronoi_cell_count; cy++) {
for (int cx = 0; cx < config.voronoi_cell_count; cx++) {
float center_x = -half + (cx + rng.randf_range(0.2f, 0.8f)) * cell_spacing;
float center_z = -half + (cy + rng.randf_range(0.2f, 0.8f)) * cell_spacing;
float floor_h = rng.randf_range(config.floor_height_min, config.floor_height_max);
voronoi_cells_.push_back({center_x, center_z, floor_h});
}
}
}
std::vector<int> DCDensity::get_perm_tables() const {
std::vector<int> tables(6 * 512);
if (!noise_) return tables;
const SimplexNoise* noises[6] = {
&noise_->primary, &noise_->tunnel_a, &noise_->tunnel_b,
&noise_->ridge, &noise_->warp, &noise_->detail
};
for (int n = 0; n < 6; n++) {
const int* src = noises[n]->perm_data();
for (int i = 0; i < 512; i++) {
tables[n * 512 + i] = src[i];
}
}
return tables;
}
// -------------------------------------------------------------------------
// CSG primitive evaluation
// -------------------------------------------------------------------------
static float eval_csg_primitive(const CSGPrimitive& p, float x, float y, float z) {
float dx = x - p.center[0];
float dy = y - p.center[1];
float dz = z - p.center[2];
if (p.type == CSGPrimitive::CSG_SPHERE) {
float dist = std::sqrt(dx*dx + dy*dy + dz*dz);
return p.size[0] - dist; // positive inside
} else { // CSG_BOX
// Exact box SDF (inverted: positive inside)
float qx = std::fabs(dx) - p.size[0];
float qy = std::fabs(dy) - p.size[1];
float qz = std::fabs(dz) - p.size[2];
float outside = std::sqrt(
std::max(qx, 0.0f) * std::max(qx, 0.0f) +
std::max(qy, 0.0f) * std::max(qy, 0.0f) +
std::max(qz, 0.0f) * std::max(qz, 0.0f));
float inside = std::min(std::max(qx, std::max(qy, qz)), 0.0f);
return -(outside + inside); // negate: standard SDF is negative inside
}
}
// Evaluate CSG primitives grouped by group_id.
// Within a group: first primitive = base, rest apply their op.
// Each group result is UNION'd with the terrain density.
static float apply_csg_primitives(const std::vector<CSGPrimitive>& prims,
float x, float y, float z, float terrain_d) {
float result = terrain_d;
int i = 0;
int n = static_cast<int>(prims.size());
while (i < n) {
int group = prims[i].group;
// First primitive in group = base shape
float gd = eval_csg_primitive(prims[i], x, y, z);
i++;
// Subsequent primitives in same group apply their op
while (i < n && prims[i].group == group) {
float pd = eval_csg_primitive(prims[i], x, y, z);
switch (prims[i].op) {
case CSGPrimitive::CSG_UNION:
gd = std::max(gd, pd);
break;
case CSGPrimitive::CSG_SUBTRACT:
gd = std::min(gd, -pd);
break;
case CSGPrimitive::CSG_INTERSECT:
gd = std::min(gd, pd);
break;
}
i++;
}
// UNION group result with terrain
result = std::max(result, gd);
}
return result;
}
float DCDensity::evaluate(float x, float y, float z, float lod_voxel_size) const {
const auto& c = config_;
if (c.spherical) {
return evaluate_spherical_lod(x, y, z, lod_voxel_size);
}
if (c.floating_islands) {
return evaluate_floating_islands(x, y, z, lod_voxel_size);
}
// For flat terrain, per-LOD octave clamping is less critical
// (fewer total octaves), but still avoids aliasing at high LODs.
// Recompute octave counts for this LOD's voxel size.
int max_primary = c.spherical ? 20 : 3;
int max_tunnel = c.spherical ? 12 : 2;
int max_ridge = c.spherical ? 16 : 2;
int max_warp = c.spherical ? 10 : 2;
int max_detail = c.spherical ? 12 : 2;
if (c.max_noise_octaves > 0) {
max_primary = std::min(max_primary, c.max_noise_octaves);
max_tunnel = std::min(max_tunnel, c.max_noise_octaves);
max_ridge = std::min(max_ridge, c.max_noise_octaves);
max_warp = std::min(max_warp, c.max_noise_octaves);
max_detail = std::min(max_detail, c.max_noise_octaves);
}
if (c.max_ridge_octaves > 0) {
max_ridge = std::min(max_ridge, c.max_ridge_octaves);
}
int lod_primary = nyquist_octaves(c.primary_frequency, lod_voxel_size, max_primary);
int lod_tunnel = nyquist_octaves(c.tunnel_frequency, lod_voxel_size, max_tunnel);
int lod_ridge = nyquist_octaves(c.ridge_frequency, lod_voxel_size, max_ridge);
int lod_warp = nyquist_octaves(c.warp_frequency, lod_voxel_size, max_warp);
int lod_detail = nyquist_octaves(c.detail_frequency, lod_voxel_size, max_detail);
return evaluate_internal(x, y, z, lod_primary, lod_tunnel, lod_ridge, lod_warp, lod_detail);
}
float DCDensity::evaluate(float x, float y, float z) const {
const auto& c = config_;
if (c.spherical) {
return evaluate_spherical(x, y, z);
}
if (c.floating_islands) {
return evaluate_floating_islands(x, y, z, 1.0f);
}
return evaluate_internal(x, y, z,
primary_octaves_, tunnel_octaves_, ridge_octaves_,
warp_octaves_, detail_octaves_);
}
float DCDensity::evaluate_internal(float x, float y, float z,
int p_oct, int t_oct, int r_oct,
int w_oct, int d_oct) const {
const auto& c = config_;
// === Ground plane (always active) ===
// Positive below y=0 (solid), negative above (air).
float d = -y * c.ground_bias;
// === Noise layers (skipped when use_noise=false) ===
if (c.use_noise && noise_) {
const auto& n = *noise_;
auto sample = [](const SimplexNoise& sn, float freq, float px, float py, float pz, int oct) {
return sn.fbm3(px * freq, py * freq, pz * freq, oct);
};
// Domain warping
float wx = x + sample(n.warp, c.warp_frequency, x, y, z, w_oct) * c.warp_strength;
float wy = y + sample(n.warp, c.warp_frequency, x + 31.7f, y + 47.3f, z + 89.1f, w_oct) * c.warp_strength;
float wz = z + sample(n.warp, c.warp_frequency, x + 73.1f, y + 11.9f, z + 59.7f, w_oct) * c.warp_strength;
// Primary 3D noise
d += sample(n.primary, c.primary_frequency, wx, wy, wz, p_oct) * c.primary_strength;
// Worm tunnels
float ta = sample(n.tunnel_a, c.tunnel_frequency, wx, wy * 0.7f, wz, t_oct);
float tb = sample(n.tunnel_b, c.tunnel_frequency, wx + 200.0f, wy * 0.7f, wz + 200.0f, t_oct);
float tube_dist_sq = ta * ta + tb * tb;
float radius_sq = c.tunnel_radius * c.tunnel_radius;
float tube_carve = std::max(1.0f - tube_dist_sq / radius_sq, 0.0f);
d -= tube_carve * c.tunnel_strength;
// Ridged noise — smooth-abs avoids derivative discontinuity at ridge peaks.
// fabs(x) has a V-shaped kink at x=0 that the voxel grid can't resolve,
// causing stair-stepping on plateau ridges. sqrt(x²+ε²) is C∞ smooth.
float ridge_raw = sample(n.ridge, c.ridge_frequency, wx, wy, wz, r_oct);
float ridge_eps = 0.08f; // smoothing radius — larger = rounder ridges
float ridged = 1.0f - std::sqrt(ridge_raw * ridge_raw + ridge_eps * ridge_eps);
ridged = std::max(ridged, 0.0f);
ridged = ridged * ridged;
d += ridged * c.ridge_strength;
// Detail noise
d += sample(n.detail, c.detail_frequency, wx, wy, wz, d_oct) * c.detail_strength;
// Voronoi floors/ceilings
if (c.floor_strength > 0.0f) {
float cell_spacing = c.terrain_size / static_cast<float>(c.voronoi_cell_count);
int gx = static_cast<int>(std::floor(x / cell_spacing));
int gz = static_cast<int>(std::floor(z / cell_spacing));
float nearest_dist_sq = std::numeric_limits<float>::max();
float nearest_floor_h = 0.0f;
float second_dist_sq = std::numeric_limits<float>::max();
int seed_base = c.seed + 100;
for (int dz_i = -1; dz_i <= 1; dz_i++) {
for (int dx_i = -1; dx_i <= 1; dx_i++) {
int cx = gx + dx_i;
int cz = gz + dz_i;
float jx = hash_float(cx, cz, seed_base);
float jz = hash_float(cx, cz, seed_base + 7);
float cell_x = (cx + 0.2f + jx * 0.6f) * cell_spacing;
float cell_z = (cz + 0.2f + jz * 0.6f) * cell_spacing;
float fh = hash_float(cx, cz, seed_base + 13);
float floor_h = c.floor_height_min + fh * (c.floor_height_max - c.floor_height_min);
float ddx = x - cell_x;
float ddz = z - cell_z;
float vd = ddx * ddx + ddz * ddz;
if (vd < nearest_dist_sq) {
second_dist_sq = nearest_dist_sq;
nearest_dist_sq = vd;
nearest_floor_h = floor_h;
} else if (vd < second_dist_sq) {
second_dist_sq = vd;
}
}
}
float edge_ratio = 1.0f - nearest_dist_sq / (second_dist_sq + 0.001f);
float cell_blend = clampf(edge_ratio * 2.0f, 0.0f, 1.0f);
cell_blend = smoothstep(cell_blend);
float floor_dist = nearest_floor_h - y;
float floor_effect = 0.0f;
if (floor_dist > 0.0f) {
floor_effect = clampf(floor_dist * 0.5f, 0.0f, 2.0f);
} else if (floor_dist > -1.5f) {
floor_effect = floor_dist * 1.5f;
}
float ceiling_h = nearest_floor_h + c.room_height;
float ceil_dist = y - ceiling_h;
float ceil_effect = 0.0f;
if (ceil_dist > 0.0f && ceil_dist < 3.0f) {
ceil_effect = -clampf(ceil_dist * 0.5f, 0.0f, 2.0f);
}
d += (floor_effect + ceil_effect) * cell_blend * c.floor_strength;
}
// Height clamping (noise terrain only)
float floor_blend = std::max((-y - c.solid_floor_depth) * 0.1f, 0.0f);
d += floor_blend;
float sky_blend = std::max((y - c.sky_height) * 0.1f, 0.0f);
d -= sky_blend;
// Crater system: spawn crater + random craters across the island
struct CraterInfo { float cx, cz, radius; };
CraterInfo craters[6];
int num_craters = 0;
// Main spawn crater at origin
craters[num_craters++] = {0.0f, 0.0f, c.flat_center_radius};
// Random craters seeded from terrain seed
for (int ci = 0; ci < 5; ci++) {
float angle = hash_float(ci, c.seed, 200) * 6.2831853f;
float dist = 40.0f + hash_float(ci, c.seed, 201) * (c.island_radius * 0.5f);
float cr = 10.0f + hash_float(ci, c.seed, 202) * 20.0f;
craters[num_craters++] = {
std::cos(angle) * dist,
std::sin(angle) * dist,
cr
};
}
for (int ci = 0; ci < num_craters; ci++) {
float cdx = x - craters[ci].cx;
float cdz = z - craters[ci].cz;
float cdist = std::sqrt(cdx * cdx + cdz * cdz);
float cr = craters[ci].radius;
float ctransition = cr * 0.5f;
if (cdist < cr + ctransition) {
// Flat center
float blend = 1.0f - clampf(
(cdist - cr * 0.3f) / (cr * 0.7f + ctransition),
0.0f, 1.0f);
blend = smoothstep(blend);
float flat_d = -y * 0.8f;
d = d + (flat_d - d) * blend;
// Crater rim
float rim_start = cr * 0.7f;
float rim_end = cr + ctransition * 0.5f;
if (cdist > rim_start && cdist < rim_end) {
float rim_t = (cdist - rim_start) / (rim_end - rim_start);
rim_t = std::sin(rim_t * 3.14159f);
float rim_height = cr * 0.2f;
float rim_noise = noise_
? n.detail.fbm3(x * c.detail_frequency * 2.0f,
y * c.detail_frequency,
z * c.detail_frequency * 2.0f, 2) * rim_height * 0.4f
: 0.0f;
float rim_pull = (rim_height + rim_noise) * rim_t;
d += std::max(0.0f, (rim_pull - y) * 0.3f) * rim_t;
}
}
}
} // end use_noise
// Apply terrain edits (CSG sphere operations)
if (config_.edits && !config_.edits->empty()) {
d = DCEditList::apply(*config_.edits, x, y, z, d);
}
// Apply CSG primitives (group-based evaluation)
if (!c.csg_primitives.empty()) {
d = apply_csg_primitives(c.csg_primitives, x, y, z, d);
}
return d;
}
// -------------------------------------------------------------------------
// Spherical/planetoid density evaluation
// -------------------------------------------------------------------------
// Cepero's moonmaker approach: spherical SDF base + concentric noise shells
// for surface detail. Noise is sampled in world space (not surface-local)
// so the octree clipmap works unchanged.
//
// Spherical density evaluation.
//
// Layers (evaluated from cheapest to most expensive):
// 1. Sphere SDF (planet_radius - distance_from_center)
// 2. Domain warping (organic surface distortion, multi-octave)
// ---------------------------------------------------------------------------
// Floating island density evaluation
// ---------------------------------------------------------------------------
float DCDensity::evaluate_floating_islands(float x, float y, float z, float lod_voxel_size) const {
const auto& c = config_;
if (!noise_) return -1.0f;
const auto& n = *noise_;
float r = c.island_radius;
float depth = r * c.island_depth_ratio;
float xz_dist = std::sqrt(x*x + z*z);
// === ABOVE Y=0: Use the EXACT original terrain ===
// Same alien landscape as flat mode (evaluate_internal)
int lod_primary = nyquist_octaves(c.primary_frequency, lod_voxel_size, 3);
int lod_tunnel = nyquist_octaves(c.tunnel_frequency, lod_voxel_size, 2);
int lod_ridge = nyquist_octaves(c.ridge_frequency, lod_voxel_size, 2);
int lod_warp = nyquist_octaves(c.warp_frequency, lod_voxel_size, 2);
int lod_detail = nyquist_octaves(c.detail_frequency, lod_voxel_size, 2);
float terrain_d = evaluate_internal(x, y, z, lod_primary, lod_tunnel, lod_ridge, lod_warp, lod_detail);
// === Smooth blend: surface terrain → volumetric rock mass ===
// Starts above Y=0, fully volumetric well below. No sharp transition.
{
float vol_freq = c.primary_frequency * 0.3f;
int vol_oct = nyquist_octaves(vol_freq, lod_voxel_size, 3);
float vol_noise = n.primary.fbm3(x * vol_freq + 77.0f, y * vol_freq, z * vol_freq + 77.0f, vol_oct);
float volume_d = vol_noise * c.primary_strength * 1.5f + c.ground_bias * 3.0f;
// Blend zone: starts at Y=+20, fully volumetric at Y=-60
float blend_top = r * 0.1f; // 20 units above center
float blend_bottom = r * 0.3f; // 60 units below center
float blend_t = (blend_top - y) / (blend_top + blend_bottom);
blend_t = std::max(0.0f, std::min(1.0f, blend_t));
blend_t = blend_t * blend_t * (3.0f - 2.0f * blend_t); // smoothstep
terrain_d = terrain_d * (1.0f - blend_t) + volume_d * blend_t;
}
// === Island boundary: flat top, noise-deformed tapered bottom ===
float island_boundary;
if (y > 0.0f) {
// Top: gentle dome, barely clips the terrain
float y_top = r * 0.5f;
float ex = x / r;
float ey = y / y_top;
float ez = z / r;
float edist = std::sqrt(ex*ex + ey*ey + ez*ez);
island_boundary = (1.0f - edist) * r;
} else {
// Bottom: cone taper to a point, deformed by noise
// Cone narrows linearly: at y=0, radius=r. At y=-depth, radius=0.
float taper_t = std::min(1.0f, -y / depth); // 0 at top, 1 at point
float cone_radius = r * (1.0f - taper_t * taper_t); // quadratic taper (rounder top, sharper point)
// Noise deformation: warps the cone boundary for organic shape
float warp_freq = c.primary_frequency * 0.5f;
int warp_oct = nyquist_octaves(warp_freq, lod_voxel_size, 3);
float warp = n.warp.fbm3(x * warp_freq, y * warp_freq * 0.3f, z * warp_freq, warp_oct);
cone_radius += warp * r * 0.25f; // ±25% radius variation
cone_radius = std::max(0.0f, cone_radius);
float xz_cone = std::sqrt(x*x + z*z);
island_boundary = cone_radius - xz_cone;
// Rock texture displaces the surface
float rf = c.primary_frequency * 1.5f;
int rock_oct = nyquist_octaves(rf, lod_voxel_size, 3);
float rock = n.primary.fbm3(x * rf, y * rf, z * rf, rock_oct) * c.primary_strength * 0.4f;
if (c.ridge_strength > 0.0f && lod_ridge > 0) {
float rrf = c.ridge_frequency;
float ridge = 1.0f - std::fabs(n.ridge.fbm3(x * rrf, y * rrf, z * rrf, lod_ridge));
ridge *= ridge;
rock += ridge * c.ridge_strength * 0.5f;
}
island_boundary += rock;
// Bottom cap: close off below the point
if (y < -depth) {
island_boundary = -std::fabs(y + depth) * 2.0f;
}
}
// === COMBINE: intersection — terrain clipped by island shape ===
float d = std::min(terrain_d, island_boundary);
// === Cave systems: swiss cheese through the lower layers ===
// Low frequency, large caverns — different from the surface caves
if (c.tunnel_strength > 0.0f) {
float cave_freq = c.tunnel_frequency;
int cave_oct = nyquist_octaves(cave_freq, lod_voxel_size, 3);
// Surface tunnels (from evaluate_internal already applied)
// Additional deep cave network: lower frequency, bigger caverns
float deep_freq = cave_freq * 0.4f; // lower frequency = bigger caves
int deep_oct = nyquist_octaves(deep_freq, lod_voxel_size, 3);
float deep_a = n.tunnel_a.fbm3(x * deep_freq, y * deep_freq * 0.5f, z * deep_freq, deep_oct);
float deep_b = n.tunnel_b.fbm3(x * deep_freq * 0.7f + 300.0f, y * deep_freq * 0.4f, z * deep_freq * 0.7f + 300.0f, deep_oct);
// Third tunnel layer for more connectivity
float deep_c = n.tunnel_a.fbm3(x * deep_freq * 1.3f + 500.0f, y * deep_freq * 0.6f + 500.0f, z * deep_freq * 1.3f, deep_oct);
float cave = std::max({deep_a, deep_b, deep_c});
float cave_threshold = 0.25f; // lower threshold = more caves
// Caves get stronger below surface (swiss cheese the bottom)
float depth_boost = 1.0f;
if (y < 0.0f) {
depth_boost = 1.0f + std::min(3.0f, -y / (depth * 0.3f)); // up to 4x stronger below
}
if (cave > cave_threshold) {
float cs = (cave - cave_threshold) / (1.0f - cave_threshold);
d -= cs * c.tunnel_strength * 2.0f * depth_boost;
}
}
// Apply edits
if (c.edits && !c.edits->empty()) {
d = DCEditList::apply(*c.edits, x, y, z, d);
}
return d;
}
// 3. Primary noise (mountains/valleys — many octaves span planet to ground scale)
// 4. Ridged noise (mountain ranges, ridgelines, crags)
// 5. Detail noise (surface roughness, boulders)
// 6. Deep crust tunnels (cave systems that extend miles underground)
//
// The key insight for planet-scale terrain: primary_frequency might be 0.000002
// (500km wavelength) with up to 20 octaves via Nyquist clamping per LOD.
// Coarse LODs (1km voxels) evaluate ~9 octaves → continent shapes.
// Fine LODs (0.5m voxels) evaluate all 20 → ground-level rocks.
// This fills the entire frequency spectrum from planet to pebble.
float DCDensity::evaluate_spherical(float x, float y, float z) const {
const auto& n = *noise_;
const auto& c = config_;
auto sample = [](const SimplexNoise& sn, float freq, float px, float py, float pz, int oct) {
return sn.fbm3(px * freq, py * freq, pz * freq, oct);
};
// Distance from planet center
float dx = x - c.planet_center[0];
float dy = y - c.planet_center[1];
float dz = z - c.planet_center[2];
float dist = std::sqrt(dx * dx + dy * dy + dz * dz);
// === Base: Sphere SDF ===
// Positive inside sphere (solid), negative outside (air)
float d = c.planet_radius - dist;
// Early out: deep interior or far above surface — no noise needed.
// The surface band must encompass all possible displacement.
float max_displacement = c.primary_strength + c.ridge_strength +
c.detail_strength + c.warp_strength * 0.5f;
if (std::fabs(d) > max_displacement * 2.5f) {
return d;
}
// === Layer 1: Domain warping (multi-octave) ===
// Warps coordinates organically — prevents noise from looking gridded.
// With many octaves, warp creates both continent-scale distortion and
// local-scale organic variation.
float wx = x + sample(n.warp, c.warp_frequency, x, y, z, warp_octaves_) * c.warp_strength;
float wy = y + sample(n.warp, c.warp_frequency, x + 31.7f, y + 47.3f, z + 89.1f, warp_octaves_) * c.warp_strength;
float wz = z + sample(n.warp, c.warp_frequency, x + 73.1f, y + 11.9f, z + 59.7f, warp_octaves_) * c.warp_strength;
// === Layer 2: Primary noise (mountains/valleys) ===
// This is the main terrain shape. With 20 octaves spanning 500km→1m,
// it creates continents, mountain ranges, hills, and boulders all in one.
d += sample(n.primary, c.primary_frequency, wx, wy, wz, primary_octaves_) * c.primary_strength;
// === Layer 3: Ridged noise (mountain ranges, crags) ===
// Ridged noise creates sharp features: mountain ridges, cliff edges, arches.
// More octaves = finer ridges and crags down to meter scale.
float ridge_raw = sample(n.ridge, c.ridge_frequency, wx, wy, wz, ridge_octaves_);
float ridge_eps = 0.08f;
float ridged = 1.0f - std::sqrt(ridge_raw * ridge_raw + ridge_eps * ridge_eps);
ridged = std::max(ridged, 0.0f);
ridged = ridged * ridged; // sharpen ridges
d += ridged * c.ridge_strength;
// === Layer 4: Detail noise (surface roughness) ===
// High-frequency surface detail — rocks, bumps, texture.
d += sample(n.detail, c.detail_frequency, wx, wy, wz, detail_octaves_) * c.detail_strength;
// === Layer 5: Deep crust tunnels ===
// Two perpendicular noise fields form "worms" where both are near zero.
// On a planet, tunnels carve through the crust creating cave systems
// that can extend for miles. The depth range is proportional to tunnel_strength.
// d values near the surface band (where solid meets air) get carved.
float crust_depth = max_displacement * 2.0f;
if (d > -c.tunnel_strength && d < crust_depth) {
float ta = sample(n.tunnel_a, c.tunnel_frequency, wx, wy * 0.7f, wz, tunnel_octaves_);
float tb = sample(n.tunnel_b, c.tunnel_frequency, wx + 200.0f, wy * 0.7f, wz + 200.0f, tunnel_octaves_);
float tube_dist_sq = ta * ta + tb * tb;
float radius_sq = c.tunnel_radius * c.tunnel_radius;
float tube_carve = std::max(1.0f - tube_dist_sq / radius_sq, 0.0f);
// Fade tunnels out near the surface to avoid swiss-cheese terrain
float depth_below_surface = d; // positive = underground
float surface_fade = std::min(depth_below_surface / (c.tunnel_strength * 0.5f), 1.0f);
surface_fade = std::max(surface_fade, 0.0f);
d -= tube_carve * c.tunnel_strength * surface_fade;
}
// Apply terrain edits
if (config_.edits && !config_.edits->empty()) {
d = DCEditList::apply(*config_.edits, x, y, z, d);
}
// Apply CSG primitives
if (!c.csg_primitives.empty()) {
d = apply_csg_primitives(c.csg_primitives, x, y, z, d);
}
return d;
}
// Per-LOD spherical evaluation: recomputes octave counts for the actual
// LOD voxel size, preventing high-frequency aliasing at coarse LODs.
float DCDensity::evaluate_spherical_lod(float x, float y, float z, float lod_voxel_size) const {
const auto& n = *noise_;
const auto& c = config_;
// Recompute octave counts for this LOD's voxel size
int max_primary = 20, max_tunnel = 12, max_ridge = 16, max_warp = 10, max_detail = 12;
if (c.max_noise_octaves > 0) {
max_primary = std::min(max_primary, c.max_noise_octaves);
max_tunnel = std::min(max_tunnel, c.max_noise_octaves);
max_ridge = std::min(max_ridge, c.max_noise_octaves);
max_warp = std::min(max_warp, c.max_noise_octaves);
max_detail = std::min(max_detail, c.max_noise_octaves);
}
if (c.max_ridge_octaves > 0) {
max_ridge = std::min(max_ridge, c.max_ridge_octaves);
}
int lod_primary = nyquist_octaves(c.primary_frequency, lod_voxel_size, max_primary);
int lod_tunnel = nyquist_octaves(c.tunnel_frequency, lod_voxel_size, max_tunnel);
int lod_ridge = nyquist_octaves(c.ridge_frequency, lod_voxel_size, max_ridge);
int lod_warp = nyquist_octaves(c.warp_frequency, lod_voxel_size, max_warp);
int lod_detail = nyquist_octaves(c.detail_frequency, lod_voxel_size, max_detail);
auto sample = [](const SimplexNoise& sn, float freq, float px, float py, float pz, int oct) {
return sn.fbm3(px * freq, py * freq, pz * freq, oct);
};
float dx = x - c.planet_center[0];
float dy = y - c.planet_center[1];
float dz = z - c.planet_center[2];
float dist = std::sqrt(dx * dx + dy * dy + dz * dz);
float d = c.planet_radius - dist;
float max_displacement = c.primary_strength + c.ridge_strength +
c.detail_strength + c.warp_strength * 0.5f;
if (std::fabs(d) > max_displacement * 2.5f) {
return d;
}
float wx = x + sample(n.warp, c.warp_frequency, x, y, z, lod_warp) * c.warp_strength;
float wy = y + sample(n.warp, c.warp_frequency, x + 31.7f, y + 47.3f, z + 89.1f, lod_warp) * c.warp_strength;
float wz = z + sample(n.warp, c.warp_frequency, x + 73.1f, y + 11.9f, z + 59.7f, lod_warp) * c.warp_strength;
d += sample(n.primary, c.primary_frequency, wx, wy, wz, lod_primary) * c.primary_strength;
float ridge_raw = sample(n.ridge, c.ridge_frequency, wx, wy, wz, lod_ridge);
float ridge_eps = 0.08f;
float ridged = 1.0f - std::sqrt(ridge_raw * ridge_raw + ridge_eps * ridge_eps);
ridged = std::max(ridged, 0.0f);
ridged = ridged * ridged;
d += ridged * c.ridge_strength;
d += sample(n.detail, c.detail_frequency, wx, wy, wz, lod_detail) * c.detail_strength;
float crust_depth = max_displacement * 2.0f;
if (d > -c.tunnel_strength && d < crust_depth) {
float ta = sample(n.tunnel_a, c.tunnel_frequency, wx, wy * 0.7f, wz, lod_tunnel);
float tb = sample(n.tunnel_b, c.tunnel_frequency, wx + 200.0f, wy * 0.7f, wz + 200.0f, lod_tunnel);
float tube_dist_sq = ta * ta + tb * tb;
float radius_sq = c.tunnel_radius * c.tunnel_radius;
float tube_carve = std::max(1.0f - tube_dist_sq / radius_sq, 0.0f);
float depth_below_surface = d;
float surface_fade = std::min(depth_below_surface / (c.tunnel_strength * 0.5f), 1.0f);
surface_fade = std::max(surface_fade, 0.0f);
d -= tube_carve * c.tunnel_strength * surface_fade;
}
if (config_.edits && !config_.edits->empty()) {
d = DCEditList::apply(*config_.edits, x, y, z, d);
}
if (!c.csg_primitives.empty()) {
d = apply_csg_primitives(c.csg_primitives, x, y, z, d);
}
return d;
}