#ifndef DC_DENSITY_GRID_H
#define DC_DENSITY_GRID_H
// Pre-computed density grid with trilinear interpolation.
// Replaces millions of expensive noise evaluations with cheap grid lookups.
//
// Usage:
// 1. Fill grid from density function: fill()
// 2. Use sample() for trilinear interpolated density at any world position
// 3. Use gradient() for central-difference normal at any world position
//
// The grid stores (size+1)^3 density values at voxel corners.
// Positions between corners are trilinear interpolated.
#include <vector>
#include <cmath>
#include <algorithm>
#include "dc_edit.h"
#include "dc_generator.h"
struct DensityGrid {
int grid_size; // number of cells per axis (e.g., 64)
int stride; // grid_size + 1
float voxel_scale;
float origin[3];
std::vector<float> data; // (stride)^3 density values
std::vector<float> normal_x; // stride^3 pre-computed normal X components
std::vector<float> normal_y; // stride^3 pre-computed normal Y components
std::vector<float> normal_z; // stride^3 pre-computed normal Z components
void allocate(int size, float vs, const float org[3]) {
grid_size = size;
stride = size + 1;
voxel_scale = vs;
origin[0] = org[0]; origin[1] = org[1]; origin[2] = org[2];
data.resize(stride * stride * stride);
int n = stride * stride * stride;
normal_x.resize(n);
normal_y.resize(n);
normal_z.resize(n);
}
// Fill grid from density function. This is the only expensive step.
// Parallelized with OpenMP — each Z-slice is independent.
template<typename DFunc>
void fill(DFunc density_fn, void* userdata) {
#pragma omp parallel for schedule(static)
for (int z = 0; z < stride; z++) {
for (int y = 0; y < stride; y++) {
for (int x = 0; x < stride; x++) {
float wx = origin[0] + x * voxel_scale;
float wy = origin[1] + y * voxel_scale;
float wz = origin[2] + z * voxel_scale;
data[z * stride * stride + y * stride + x] =
density_fn(wx, wy, wz, userdata);
}
}
}
}
// Fill from a GENERATOR, one Z-slice per batch.
//
// Same grid, same order, same values as fill() with a lambda that calls
// evaluate() -- DCGenerator::evaluate_batch's default IS that loop. What it
// buys is that a generator which resolves per-grid state can do it once per
// slice instead of once per corner. Measured on the proc-nodes bridge, whose
// evaluate() scanned a vector of Godot Strings per point.
//
// A slice rather than the whole grid so the scratch buffers stay small and
// the OpenMP parallelism over z is unchanged.
void fill_generator(DCGenerator* gen, float lod_voxel_size) {
#pragma omp parallel for schedule(static)
for (int z = 0; z < stride; z++) {
std::vector<float> pos((size_t)stride * stride * 3);
std::vector<float> out((size_t)stride * stride);
float wz = origin[2] + z * voxel_scale;
size_t k = 0;
for (int y = 0; y < stride; y++) {
float wy = origin[1] + y * voxel_scale;
for (int x = 0; x < stride; x++) {
pos[k * 3 + 0] = origin[0] + x * voxel_scale;
pos[k * 3 + 1] = wy;
pos[k * 3 + 2] = wz;
k++;
}
}
gen->evaluate_batch(pos.data(), (int)k, lod_voxel_size, out.data());
std::copy(out.begin(), out.begin() + k,
data.begin() + (size_t)z * stride * stride);
}
}
// Direct grid lookup at integer voxel coordinates.
inline float at(int x, int y, int z) const {
x = std::max(0, std::min(x, grid_size));
y = std::max(0, std::min(y, grid_size));
z = std::max(0, std::min(z, grid_size));
return data[z * stride * stride + y * stride + x];
}
// Trilinear interpolation at world position.
float sample(float wx, float wy, float wz) const {
// Convert world → grid coordinates
float gx = (wx - origin[0]) / voxel_scale;
float gy = (wy - origin[1]) / voxel_scale;
float gz = (wz - origin[2]) / voxel_scale;
// Clamp to grid
gx = std::max(0.0f, std::min(gx, (float)grid_size));
gy = std::max(0.0f, std::min(gy, (float)grid_size));
gz = std::max(0.0f, std::min(gz, (float)grid_size));
int x0 = (int)gx, y0 = (int)gy, z0 = (int)gz;
if (x0 >= grid_size) x0 = grid_size - 1;
if (y0 >= grid_size) y0 = grid_size - 1;
if (z0 >= grid_size) z0 = grid_size - 1;
float fx = gx - x0, fy = gy - y0, fz = gz - z0;
// Trilinear interpolation
float c000 = at(x0, y0, z0), c100 = at(x0+1, y0, z0);
float c010 = at(x0, y0+1, z0), c110 = at(x0+1, y0+1, z0);
float c001 = at(x0, y0, z0+1), c101 = at(x0+1, y0, z0+1);
float c011 = at(x0, y0+1, z0+1), c111 = at(x0+1, y0+1, z0+1);
float c00 = c000 + fx * (c100 - c000);
float c10 = c010 + fx * (c110 - c010);
float c01 = c001 + fx * (c101 - c001);
float c11 = c011 + fx * (c111 - c011);
float c0 = c00 + fy * (c10 - c00);
float c1 = c01 + fy * (c11 - c01);
return c0 + fz * (c1 - c0);
}
// Pre-compute normals at every grid point using the gradient() method.
// 3x cheaper than computing gradient on-the-fly during hermite data collection:
// sample_normal() does 1 trilinear interp per component (3 total) vs
// gradient() which does 4 trilinear interps of density.
// Parallelized with OpenMP — each Z-slice is independent.
void compute_normals() {
#pragma omp parallel for schedule(static)
for (int z = 0; z < stride; z++) {
for (int y = 0; y < stride; y++) {
for (int x = 0; x < stride; x++) {
float wx = origin[0] + x * voxel_scale;
float wy = origin[1] + y * voxel_scale;
float wz = origin[2] + z * voxel_scale;
float n[3];
gradient(wx, wy, wz, n);
int idx = z * stride * stride + y * stride + x;
normal_x[idx] = n[0];
normal_y[idx] = n[1];
normal_z[idx] = n[2];
}
}
}
}
// Trilinear interpolation of pre-computed normals at world position.
// Much cheaper than gradient() which requires 4 trilinear density samples.
void sample_normal(float wx, float wy, float wz, float normal[3]) const {
float gx = (wx - origin[0]) / voxel_scale;
float gy = (wy - origin[1]) / voxel_scale;
float gz = (wz - origin[2]) / voxel_scale;
gx = std::max(0.0f, std::min(gx, (float)grid_size));
gy = std::max(0.0f, std::min(gy, (float)grid_size));
gz = std::max(0.0f, std::min(gz, (float)grid_size));
int x0 = (int)gx, y0 = (int)gy, z0 = (int)gz;
if (x0 >= grid_size) x0 = grid_size - 1;
if (y0 >= grid_size) y0 = grid_size - 1;
if (z0 >= grid_size) z0 = grid_size - 1;
float fx = gx - x0, fy = gy - y0, fz = gz - z0;
int s = stride, s2 = stride * stride;
auto interp_component = [&](const std::vector<float>& arr) -> float {
float c000 = arr[z0*s2 + y0*s + x0], c100 = arr[z0*s2 + y0*s + x0+1];
float c010 = arr[z0*s2 + (y0+1)*s + x0], c110 = arr[z0*s2 + (y0+1)*s + x0+1];
float c001 = arr[(z0+1)*s2 + y0*s + x0], c101 = arr[(z0+1)*s2 + y0*s + x0+1];
float c011 = arr[(z0+1)*s2 + (y0+1)*s + x0], c111 = arr[(z0+1)*s2 + (y0+1)*s + x0+1];
float c00 = c000 + fx*(c100-c000), c10 = c010 + fx*(c110-c010);
float c01 = c001 + fx*(c101-c001), c11 = c011 + fx*(c111-c011);
float c_0 = c00 + fy*(c10-c00), c_1 = c01 + fy*(c11-c01);
return c_0 + fz*(c_1-c_0);
};
normal[0] = interp_component(normal_x);
normal[1] = interp_component(normal_y);
normal[2] = interp_component(normal_z);
float len = std::sqrt(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]);
if (len > 1e-10f) { normal[0]/=len; normal[1]/=len; normal[2]/=len; }
}
// Gradient via Quilez tetrahedron technique (4 samples instead of 6).
// Negated so normals point outward (away from solid/positive density).
// Reference: https://iquilezles.org/articles/normalsSDF/
void gradient(float wx, float wy, float wz, float normal[3]) const {
float h = voxel_scale * 0.25f; // quarter cell: captures local gradient accurately (density already smoothed)
// Tetrahedron vertices: k0=(1,-1,-1), k1=(-1,-1,1), k2=(-1,1,-1), k3=(1,1,1)
float d0 = sample(wx + h, wy - h, wz - h);
float d1 = sample(wx - h, wy - h, wz + h);
float d2 = sample(wx - h, wy + h, wz - h);
float d3 = sample(wx + h, wy + h, wz + h);
// Negate for outward-facing normals (density positive = solid)
normal[0] = -(d0 - d1 - d2 + d3);
normal[1] = -(-d0 - d1 + d2 + d3);
normal[2] = -(-d0 + d1 - d2 + d3);
float len = std::sqrt(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]);
if (len > 1e-10f) {
normal[0] /= len; normal[1] /= len; normal[2] /= len;
}
}
// Raw (unnormalized) gradient for Newton-step vertex polish.
// Returns the actual gradient vector so |grad|² is meaningful.
void gradient_raw(float wx, float wy, float wz, float grad[3]) const {
float h = voxel_scale * 0.25f;
float d0 = sample(wx + h, wy - h, wz - h);
float d1 = sample(wx - h, wy - h, wz + h);
float d2 = sample(wx - h, wy + h, wz - h);
float d3 = sample(wx + h, wy + h, wz + h);
grad[0] = -(d0 - d1 - d2 + d3);
grad[1] = -(-d0 - d1 + d2 + d3);
grad[2] = -(-d0 + d1 - d2 + d3);
}
// Plain Laplacian smooth (Cepero's "Just Relax").
void smooth(int iterations = 1, float lambda = 0.5f) {
std::vector<float> tmp(data.size());
for (int iter = 0; iter < iterations; iter++) {
#pragma omp parallel for schedule(static)
for (int z = 0; z < stride; z++) {
for (int y = 0; y < stride; y++) {
for (int x = 0; x < stride; x++) {
int idx = z * stride * stride + y * stride + x;
float center = data[idx];
if (x > 0 && x < grid_size && y > 0 && y < grid_size && z > 0 && z < grid_size) {
float avg = (
data[idx - 1] + data[idx + 1] +
data[idx - stride] + data[idx + stride] +
data[idx - stride*stride] + data[idx + stride*stride]
) / 6.0f;
tmp[idx] = center + lambda * (avg - center);
} else {
tmp[idx] = center;
}
}
}
}
std::swap(data, tmp);
}
}
// Two-band octave-aware smoothing.
// Separates the density field into low-frequency (terrain shape) and
// high-frequency (noise detail that aliases into stair-steps) bands.
// Only the high-frequency band is smoothed, preserving terrain shape
// exactly — no LOD height offset regardless of smoothing strength.
//
// low_freq_fn: density function evaluated with fewer noise octaves
// (e.g., 2x the voxel_size to drop the highest octave).
// iterations: smoothing passes on the high-frequency band
// lambda: blend factor per pass
template<typename DFunc>
void smooth_octave_aware(DFunc low_freq_fn, void* userdata,
int iterations = 4, float lambda = 0.5f) {
int n = stride * stride * stride;
// 1. Sample low-frequency density at every grid point
std::vector<float> low_freq(n);
#pragma omp parallel for schedule(static)
for (int z = 0; z < stride; z++) {
for (int y = 0; y < stride; y++) {
for (int x = 0; x < stride; x++) {
float wx = origin[0] + x * voxel_scale;
float wy = origin[1] + y * voxel_scale;
float wz = origin[2] + z * voxel_scale;
low_freq[z * stride * stride + y * stride + x] =
low_freq_fn(wx, wy, wz, userdata);
}
}
}
// 2. Extract high-frequency band: high = full - low
std::vector<float> high_freq(n);
for (int i = 0; i < n; i++)
high_freq[i] = data[i] - low_freq[i];
// 3. Smooth ONLY the high-frequency band (plain Laplacian is fine —
// high-freq component has near-zero mean, no volume to preserve)
std::vector<float> tmp(n);
for (int iter = 0; iter < iterations; iter++) {
#pragma omp parallel for schedule(static)
for (int z = 0; z < stride; z++) {
for (int y = 0; y < stride; y++) {
for (int x = 0; x < stride; x++) {
int idx = z * stride * stride + y * stride + x;
float center = high_freq[idx];
if (x > 0 && x < grid_size && y > 0 && y < grid_size && z > 0 && z < grid_size) {
float avg = (
high_freq[idx - 1] + high_freq[idx + 1] +
high_freq[idx - stride] + high_freq[idx + stride] +
high_freq[idx - stride*stride] + high_freq[idx + stride*stride]
) / 6.0f;
tmp[idx] = center + lambda * (avg - center);
} else {
tmp[idx] = center;
}
}
}
}
std::swap(high_freq, tmp);
}
// 4. Recombine: result = low_freq + smoothed_high_freq
for (int i = 0; i < n; i++)
data[i] = low_freq[i] + high_freq[i];
}
// Apply terrain edits onto an existing density grid.
// Used after loading a cached base grid to overlay edit effects
// without re-evaluating expensive noise functions.
void apply_edits(const std::vector<TerrainEdit>& edits) {
if (edits.empty()) return;
#pragma omp parallel for schedule(static)
for (int z = 0; z < stride; z++) {
for (int y = 0; y < stride; y++) {
for (int x = 0; x < stride; x++) {
float wx = origin[0] + x * voxel_scale;
float wy = origin[1] + y * voxel_scale;
float wz = origin[2] + z * voxel_scale;
int idx = z * stride * stride + y * stride + x;
data[idx] = DCEditList::apply(edits, wx, wy, wz, data[idx]);
}
}
}
}
// Copy only density data (not normals) from another grid.
// Used to restore cached base density before applying new edits.
void copy_density_from(const DensityGrid& src) {
data = src.data;
}
};
#endif // DC_DENSITY_GRID_H