diff --git a/HighMap/include/highmap/authoring.hpp b/HighMap/include/highmap/authoring.hpp index 928558c11..a1ce41de8 100644 --- a/HighMap/include/highmap/authoring.hpp +++ b/HighMap/include/highmap/authoring.hpp @@ -24,6 +24,62 @@ enum StampingBlendMethod : int SUBSTRACT, ///< substract }; +/** + * @brief Penalty type of a deformation constraint (see @ref + * DeformationConstraint). + */ +enum DeformationConstraintType : int +{ + MATCH, ///< penalize any deviation from the target height + ABOVE, ///< penalize only heights below the target (target acts as a floor) + BELOW, ///< penalize only heights above the target (target acts as a + ///< ceiling) +}; + +/** + * @brief Per-vertex penalty term used by @ref sls_deformation. + * + * Each constraint contributes `scale * weight(x, y) * (z(x, y) - target(x, + * y))^2` to the fitness of a vertex, gated by its type: for `MATCH` the + * penalty always applies, for `ABOVE` only where `z < target`, for `BELOW` + * only where `z > target`. Vertices with a zero weight are unconstrained by + * that term. Several constraints are summed, so the `scale` factor can be + * used to normalize competing terms. + * + * Typical uses (after Stachniak & Stuerzlinger, 2005): + * - match a reference height map: `target = ref, weight = 1, MATCH`; + * - shape mask (island): `target = water_level, weight = mask, ABOVE` plus + * `target = water_level, weight = 1 - mask, BELOW`; + * - flat road: `target = road_height, weight = path, MATCH` plus + * `target = original, weight = 1 - path, MATCH` (preserve the rest); + * - edge matching: `target` = neighbour terrain edge values, `weight` non-zero + * along the edge only, `MATCH`. + */ +struct DeformationConstraint +{ + Array target; ///< desired height field + Array weight; ///< per-vertex weight (0 = free) + DeformationConstraintType type = MATCH; ///< penalty type + float scale = 1.f; ///< global weight of the term +}; + +/** + * @brief Truncated-Gaussian push operation, the elementary deformation used + * by @ref sls_deformation. + * + * The push adds `amplitude * G(d)` to every vertex within `ir` pixels of the + * center `(i, j)`, where `G` is a truncated Gaussian kernel equal to 1 at the + * center and to 0 at distance `ir`. The center vertex is therefore displaced + * by exactly `amplitude`. + */ +struct GaussianPush +{ + int i; ///< center index (x) + int j; ///< center index (y) + int ir; ///< kernel radius in pixels + float amplitude; ///< height displacement at the center +}; + /** * @brief Point-wise alteration: locally enforce a new elevation value while * maintaining the 'shape' of the heightmap. @@ -72,6 +128,21 @@ void alter_elevation(Array &array, glm::vec2 shift = {0.f, 0.f}, glm::vec2 scale = {1.f, 1.f}); +/** + * @brief Apply a sequence of truncated-Gaussian pushes to a heightmap. + * + * Replays, in order, the deformations recorded by @ref sls_deformation. Since + * the deformed terrain is fully described by the original terrain and the + * push sequence, this also provides a compact storage of the deformation. + * + * @param array Heightmap to deform (modified in place). + * @param pushes Push operations to apply, in order. + * + * @see sls_deformation + */ +void apply_gaussian_pushes(Array &array, + const std::vector &pushes); + /** * @brief Generate a heightmap from a coarse grid of control points with defined * elevation values. @@ -326,6 +397,78 @@ Array ridgelines_bezier(glm::ivec2 shape, const Array *p_stretching = nullptr, glm::vec4 bbox_array = {0.f, 1.f, 0.f, 1.f}); +/** + * @brief Deform a heightmap to satisfy a set of constraints using stochastic + * local search over truncated-Gaussian push operations. + * + * Implements the constraint-based terrain deformation of S. Stachniak and W. + * Stuerzlinger, "An Algorithm for Automated Fractal Terrain Deformation", + * WSCG 2005. The terrain `T` is deformed into `T'` by searching for a + * sequence of local push operations `(location, amplitude, radius)` (see @ref + * GaussianPush) that minimize a fitness function `F(T) = sum F(x, y)`, where + * the per-vertex penalty `F(x, y)` is the sum of the supplied constraint + * terms (see @ref DeformationConstraint). Because the penalties are summed, + * several constraints (shape masks, fixed paths, edge matching, reference + * heights...) can be satisfied simultaneously. + * + * At every iteration, a set of candidate vertices is drawn: half on a + * jittered uniform grid, half sampled proportionally to the current penalty + * map so that the search concentrates where constraints are violated. For + * each candidate and each radius, the push amplitude is the least-squares + * optimum of the (locally quadratic) penalty over the kernel footprint, + * clamped by the slope limit `talus_max * radius` which suppresses + * high-frequency spikes (frequency limitation of the paper). The best + * deformation is applied with probability `p_best`, otherwise one of the + * `top_fraction` best deformations is applied (the stochastic noise of the + * search, which prevents stalling in local minima). The search stops when the + * fitness drops below `tolerance` times its initial value, when no improving + * deformation can be found for a number of consecutive iterations, or after + * `iterations` pushes. + * + * @param z Input heightmap. + * @param constraints Constraint terms defining the fitness function. Every + * target and weight array must have the shape of `z`. + * @param seed Random seed number. + * @param iterations Maximum number of iterations (one push per iteration). + * @param ir_min Smallest push radius, in pixels. + * @param ir_max Largest push radius, in pixels. + * @param n_radii Number of radii, geometrically spaced in [ir_min, + * ir_max]. + * @param talus_max Maximum push amplitude per pixel of radius (slope + * limit). If zero or negative, defaults to the elevation + * range of the inputs divided by `ir_max`. + * @param n_candidates Number of candidate vertices evaluated per iteration. + * @param p_best Probability of applying the best deformation found + * (0.65 in the original paper). + * @param top_fraction Fraction of the best deformations a sub-optimal choice + * is drawn from. + * @param tolerance Relative fitness tolerance for early termination. + * @param p_pushes Optional pointer to a vector receiving the applied + * pushes, in order (see @ref apply_gaussian_pushes). + * @return Array Deformed heightmap. + * + * **Example** + * @include ex_sls_deformation.cpp + * + * **Result** + * @image html ex_sls_deformation.png + * + * @see apply_gaussian_pushes, DeformationConstraint + */ +Array sls_deformation(const Array &z, + const std::vector &constraints, + std::uint32_t seed, + int iterations = 200, + int ir_min = 4, + int ir_max = 64, + int n_radii = 7, + float talus_max = 0.f, + int n_candidates = 512, + float p_best = 0.65f, + float top_fraction = 0.1f, + float tolerance = 1e-3f, + std::vector *p_pushes = nullptr); + /** * @brief Generate a heightmap by stamping a kernel at predefined locations. * diff --git a/HighMap/src/authoring/sls_deformation.cpp b/HighMap/src/authoring/sls_deformation.cpp new file mode 100644 index 000000000..5ee9e69d9 --- /dev/null +++ b/HighMap/src/authoring/sls_deformation.cpp @@ -0,0 +1,443 @@ +/* Copyright (c) 2026 Otto Link. Distributed under the terms of the GNU General + * Public License. The full license is in the file LICENSE, distributed with + * this software. */ +#include +#include +#include +#include +#include +#include +#include + +#include "highmap/array.hpp" +#include "highmap/authoring.hpp" +#include "highmap/internal/validation.hpp" +#include "highmap/logger.hpp" + +namespace hmap +{ + +namespace +{ + +// Truncated Gaussian kernel of radius ir, equal to 1 at the center and to 0 at +// distance ir (and beyond), stored on a (2 ir + 1)^2 square window +struct PushKernel +{ + int ir = 0; + int n = 0; + std::vector values; + + float at(int di, int dj) const + { + return this->values[(dj + this->ir) * this->n + (di + this->ir)]; + } +}; + +PushKernel helper_make_kernel(int ir) +{ + PushKernel kernel; + kernel.ir = ir; + kernel.n = 2 * ir + 1; + kernel.values.assign(static_cast(kernel.n) * kernel.n, 0.f); + + float sigma = static_cast(ir) / 3.f; + float inv_2s2 = 1.f / (2.f * sigma * sigma); + float r2 = static_cast(ir * ir); + float e_r = std::exp(-r2 * inv_2s2); + + for (int dj = -ir; dj <= ir; ++dj) + for (int di = -ir; di <= ir; ++di) + { + float d2 = static_cast(di * di + dj * dj); + if (d2 >= r2) continue; + + kernel.values[(dj + ir) * kernel.n + (di + ir)] = (std::exp(-d2 * + inv_2s2) - + e_r) / + (1.f - e_r); + } + + return kernel; +} + +// kernel window clipped to the array, [i0, i1[ x [j0, j1[ +struct Window +{ + int i0, i1, j0, j1; +}; + +Window helper_window(glm::ivec2 shape, int ic, int jc, int ir) +{ + return {std::max(0, ic - ir), + std::min(shape.x, ic + ir + 1), + std::max(0, jc - ir), + std::min(shape.y, jc + ir + 1)}; +} + +void helper_apply_push(Array &z, + const PushKernel &kernel, + int ic, + int jc, + float amplitude) +{ + Window w = helper_window(z.shape, ic, jc, kernel.ir); + + for (int j = w.j0; j < w.j1; ++j) + for (int i = w.i0; i < w.i1; ++i) + { + float g = kernel.at(i - ic, j - jc); + if (g > 0.f) z(i, j) += amplitude * g; + } +} + +// lightweight view of a constraint +struct Term +{ + const float *target; + const float *weight; + DeformationConstraintType type; + float scale; +}; + +bool helper_is_active(DeformationConstraintType type, float d) +{ + switch (type) + { + case ABOVE: + return d < 0.f; + case BELOW: + return d > 0.f; + case MATCH: + default: + return true; + } +} + +float helper_penalty(const std::vector &terms, size_t k, float z) +{ + float p = 0.f; + for (const Term &t : terms) + { + float w = t.weight[k]; + if (w == 0.f) continue; + + float d = z - t.target[k]; + if (helper_is_active(t.type, d)) p += t.scale * w * d * d; + } + return p; +} + +struct Candidate +{ + int i = 0; + int j = 0; + int ik = 0; // kernel index + float amplitude = 0.f; + float delta = 0.f; // fitness change, negative means improvement +}; + +// Evaluate the best push of a given kernel at (ic, jc). The amplitude is the +// least-squares optimum of the local quadratic model of the penalty (exact +// for MATCH terms, using the currently active vertices for one-sided terms), +// clamped by the slope limit; the fitness change is then computed exactly. +bool helper_evaluate(const Array &z, + const std::vector &penalty, + const std::vector &terms, + const PushKernel &kernel, + int ic, + int jc, + float amplitude_max, + Candidate &out) +{ + Window w = helper_window(z.shape, ic, jc, kernel.ir); + + float s1 = 0.f; + float s2 = 0.f; + + for (int j = w.j0; j < w.j1; ++j) + for (int i = w.i0; i < w.i1; ++i) + { + float g = kernel.at(i - ic, j - jc); + if (g <= 0.f) continue; + + size_t k = static_cast(j) * z.shape.x + i; + float zv = z.vector[k]; + + for (const Term &t : terms) + { + float wt = t.weight[k]; + if (wt == 0.f) continue; + + float d = zv - t.target[k]; + if (!helper_is_active(t.type, d)) continue; + + float sw = t.scale * wt; + s1 += sw * g * d; + s2 += sw * g * g; + } + } + + if (s2 <= 0.f) return false; + + float amplitude = std::clamp(-s1 / s2, -amplitude_max, amplitude_max); + if (std::abs(amplitude) <= std::numeric_limits::epsilon()) + return false; + + float delta = 0.f; + + for (int j = w.j0; j < w.j1; ++j) + for (int i = w.i0; i < w.i1; ++i) + { + float g = kernel.at(i - ic, j - jc); + if (g <= 0.f) continue; + + size_t k = static_cast(j) * z.shape.x + i; + delta += helper_penalty(terms, k, z.vector[k] + amplitude * g) - + penalty[k]; + } + + if (delta >= 0.f) return false; + + out.i = ic; + out.j = jc; + out.amplitude = amplitude; + out.delta = delta; + return true; +} + +} // namespace + +void apply_gaussian_pushes(Array &array, + const std::vector &pushes) +{ + if (!validate_non_empty(array)) return; + + std::map kernels; + + for (const GaussianPush &p : pushes) + { + if (p.ir < 1) continue; + + auto it = kernels.find(p.ir); + if (it == kernels.end()) + it = kernels.emplace(p.ir, helper_make_kernel(p.ir)).first; + + helper_apply_push(array, it->second, p.i, p.j, p.amplitude); + } +} + +Array sls_deformation(const Array &z, + const std::vector &constraints, + std::uint32_t seed, + int iterations, + int ir_min, + int ir_max, + int n_radii, + float talus_max, + int n_candidates, + float p_best, + float top_fraction, + float tolerance, + std::vector *p_pushes) +{ + Array z_out = z; + + if (p_pushes) p_pushes->clear(); + + if (!validate_non_empty(z)) return z_out; + if (constraints.empty() || iterations <= 0) return z_out; + + for (const DeformationConstraint &c : constraints) + { + if (!validate_same_shape(z, c.target)) return z_out; + if (!validate_same_shape(z, c.weight)) return z_out; + } + + // --- parameters + + ir_min = std::max(1, ir_min); + ir_max = std::max(ir_min, ir_max); + n_radii = std::max(1, n_radii); + n_candidates = std::max(1, n_candidates); + p_best = std::clamp(p_best, 0.f, 1.f); + top_fraction = std::clamp(top_fraction, 0.f, 1.f); + + std::vector terms; + terms.reserve(constraints.size()); + for (const DeformationConstraint &c : constraints) + terms.push_back( + {c.target.vector.data(), c.weight.vector.data(), c.type, c.scale}); + + // geometrically spaced radii + std::vector kernels; + for (int k = 0; k < n_radii; ++k) + { + float t = n_radii > 1 ? static_cast(k) / (n_radii - 1) : 0.f; + int ir = static_cast( + std::round(ir_min * std::pow(static_cast(ir_max) / ir_min, t))); + if (kernels.empty() || kernels.back().ir != ir) + kernels.push_back(helper_make_kernel(ir)); + } + + // default slope limit: cross the whole elevation range with the largest + // radius + if (talus_max <= 0.f) + { + float vmin = z.min(); + float vmax = z.max(); + for (const DeformationConstraint &c : constraints) + { + vmin = std::min(vmin, c.target.min()); + vmax = std::max(vmax, c.target.max()); + } + talus_max = (vmax - vmin) / static_cast(ir_max); + if (talus_max <= 0.f) return z_out; + } + + // --- penalty map and initial fitness + + size_t n = z.vector.size(); + std::vector penalty(n); + + for (size_t k = 0; k < n; ++k) + penalty[k] = helper_penalty(terms, k, z_out.vector[k]); + + double fitness = std::accumulate(penalty.begin(), penalty.end(), 0.0); + if (fitness <= 0.0) return z_out; + + double fitness_stop = static_cast(tolerance) * fitness; + + // --- search + + std::mt19937 gen(seed); + std::uniform_real_distribution dis(0.f, 1.f); + + int n_uniform = n_candidates / 2; + int n_adaptive = n_candidates - n_uniform; + int nx = std::max( + 1, + static_cast(std::round( + std::sqrt(static_cast(n_uniform) * z.shape.x / z.shape.y)))); + int ny = std::max(1, n_uniform / nx); + + const int stall_max = 20; + int stall = 0; + + std::vector cand_ij; + std::vector cand; + std::vector cand_valid; + std::vector improving; + + for (int it = 0; it < iterations; ++it) + { + if (fitness <= fitness_stop) break; + + // candidate vertices: jittered uniform grid + penalty-weighted samples + cand_ij.clear(); + + if (n_uniform > 0) + for (int gj = 0; gj < ny; ++gj) + for (int gi = 0; gi < nx; ++gi) + { + int i = std::min(z.shape.x - 1, + static_cast((gi + dis(gen)) * z.shape.x / nx)); + int j = std::min(z.shape.y - 1, + static_cast((gj + dis(gen)) * z.shape.y / ny)); + cand_ij.push_back({i, j}); + } + + if (n_adaptive > 0) + { + std::discrete_distribution dd(penalty.begin(), penalty.end()); + for (int s = 0; s < n_adaptive; ++s) + { + int k = dd(gen); + cand_ij.push_back({k % z.shape.x, k / z.shape.x}); + } + } + + // evaluate the best push per candidate over all radii + int nc = static_cast(cand_ij.size()); + cand.assign(nc, Candidate()); + cand_valid.assign(nc, 0); + +#pragma omp parallel for schedule(dynamic) + for (int c = 0; c < nc; ++c) + { + Candidate best; + bool found = false; + + for (int ik = 0; ik < static_cast(kernels.size()); ++ik) + { + Candidate tmp; + float amplitude_max = talus_max * kernels[ik].ir; + + if (helper_evaluate(z_out, + penalty, + terms, + kernels[ik], + cand_ij[c].x, + cand_ij[c].y, + amplitude_max, + tmp) && + tmp.delta < best.delta) + { + tmp.ik = ik; + best = tmp; + found = true; + } + } + + cand[c] = best; + cand_valid[c] = found ? 1 : 0; + } + + improving.clear(); + for (int c = 0; c < nc; ++c) + if (cand_valid[c]) improving.push_back(c); + + if (improving.empty()) + { + if (++stall >= stall_max) break; + continue; + } + stall = 0; + + // stochastic local search: best choice with probability p_best, + // otherwise one of the top_fraction best choices + std::sort(improving.begin(), + improving.end(), + [&cand](int a, int b) { return cand[a].delta < cand[b].delta; }); + + int pick = 0; + if (dis(gen) > p_best) + { + int n_top = std::max( + 1, + static_cast(std::ceil(top_fraction * improving.size()))); + pick = std::min(n_top - 1, static_cast(dis(gen) * n_top)); + } + + const Candidate &sel = cand[improving[pick]]; + const PushKernel &kernel = kernels[sel.ik]; + + helper_apply_push(z_out, kernel, sel.i, sel.j, sel.amplitude); + + Window w = helper_window(z.shape, sel.i, sel.j, kernel.ir); + for (int j = w.j0; j < w.j1; ++j) + for (int i = w.i0; i < w.i1; ++i) + { + size_t k = static_cast(j) * z.shape.x + i; + penalty[k] = helper_penalty(terms, k, z_out.vector[k]); + } + + fitness = std::accumulate(penalty.begin(), penalty.end(), 0.0); + + if (p_pushes) p_pushes->push_back({sel.i, sel.j, kernel.ir, sel.amplitude}); + } + + return z_out; +} + +} // namespace hmap diff --git a/docs/images/ex_sls_deformation.png b/docs/images/ex_sls_deformation.png new file mode 100644 index 000000000..12cba3733 Binary files /dev/null and b/docs/images/ex_sls_deformation.png differ diff --git a/examples/ex_sls_deformation/CMakeLists.txt b/examples/ex_sls_deformation/CMakeLists.txt new file mode 100644 index 000000000..74397f072 --- /dev/null +++ b/examples/ex_sls_deformation/CMakeLists.txt @@ -0,0 +1,2 @@ +add_executable(ex_sls_deformation ex_sls_deformation.cpp) +target_link_libraries(ex_sls_deformation highmap) diff --git a/examples/ex_sls_deformation/ex_sls_deformation.cpp b/examples/ex_sls_deformation/ex_sls_deformation.cpp new file mode 100644 index 000000000..a97bd94a2 --- /dev/null +++ b/examples/ex_sls_deformation/ex_sls_deformation.cpp @@ -0,0 +1,131 @@ +#include +#include + +#include "highmap.hpp" + +int main(void) +{ + const glm::ivec2 shape = {512, 512}; + const uint32_t seed = 1; + + const float water_level = 0.5f; + const float coast_margin = 0.1f; + const float road_level = 0.65f; + + hmap::Array z = hmap::noise_fbm(hmap::NoiseType::PERLIN, + shape, + {4.f, 4.f}, + seed); + hmap::remap(z); + + // 1. Shape constraint: a five-pointed star island. Vertices inside the + // mask must be above the water level, vertices outside must be below (with + // a small margin on each side to get a clear coastline). + + hmap::Path star; + for (int k = 0; k < 10; ++k) + { + float angle = 0.5f * M_PI + 2.f * M_PI * static_cast(k) / 10.f; + float r = (k % 2 == 0) ? 0.42f : 0.18f; + star.add_point( + {0.5f + r * std::cos(angle), 0.5f + r * std::sin(angle), 1.f}); + } + star.set_closed(true); + + hmap::Array mask(shape, 0.f); + star.to_array(mask, {0.f, 1.f, 0.f, 1.f}, true); + + hmap::Array sea(shape, 1.f); + sea -= mask; + + hmap::DeformationConstraint land; + land.target = hmap::Array(shape, water_level + coast_margin); + land.weight = mask; + land.type = hmap::DeformationConstraintType::ABOVE; + + hmap::DeformationConstraint ocean; + ocean.target = hmap::Array(shape, water_level - coast_margin); + ocean.weight = sea; + ocean.type = hmap::DeformationConstraintType::BELOW; + + // 2. Path constraint: a flat S-shaped road at a fixed elevation, the rest + // of the terrain being preserved as much as possible (to avoid a global + // flattening of the terrain) + + hmap::Path road; + road.add_point({0.08f, 0.15f, 1.f}); + road.add_point({0.35f, 0.25f, 1.f}); + road.add_point({0.50f, 0.50f, 1.f}); + road.add_point({0.65f, 0.75f, 1.f}); + road.add_point({0.92f, 0.85f, 1.f}); + road = hmap::catmullrom(road); + + hmap::Array path(shape, 0.f); + road.to_array(path); + path = hmap::dilation(path, 4); + + hmap::Array not_path(shape, 1.f); + not_path -= path; + + hmap::DeformationConstraint flat_road; + flat_road.target = hmap::Array(shape, road_level); + flat_road.weight = path; + flat_road.type = hmap::DeformationConstraintType::MATCH; + + hmap::DeformationConstraint preserve; + preserve.target = z; + preserve.weight = not_path; + preserve.type = hmap::DeformationConstraintType::MATCH; + + // 3. Run the search for each set of constraints, and for the combination + // (the road is allowed to cross the sea as a causeway: sea vertices under + // the road are released from the ocean constraint) + + auto run = [&](const std::string &label, + const std::vector &constraints, + int iterations) + { + std::vector pushes; + + auto t0 = std::chrono::steady_clock::now(); + hmap::Array out = hmap::sls_deformation(z, + constraints, + seed, + iterations, + 4, + 64, + 7, + 0.f, + 512, + 0.65f, + 0.1f, + 1e-3f, + &pushes); + auto t1 = std::chrono::steady_clock::now(); + + std::cout << label << ": " << pushes.size() << " pushes in " + << std::chrono::duration(t1 - t0).count() << " s\n"; + return out; + }; + + hmap::Array z_star = run("star island", {land, ocean}, 600); + hmap::Array z_road = run("flat road", {flat_road, preserve}, 800); + + hmap::DeformationConstraint ocean_causeway = ocean; + ocean_causeway.weight *= not_path; + + hmap::Array z_both = run("star island + road", + {land, ocean_causeway, flat_road}, + 1000); + + // render with the sea flattened at the water level + for (hmap::Array *p : {&z_star, &z_road, &z_both}) + hmap::clamp_min(*p, water_level); + + hmap::export_banner_png("ex_sls_deformation.png", + {z, z_star, z_road, z_both}, + hmap::Cmap::TERRAIN, + true); + + return 0; +} diff --git a/tests/src/test_sls_deformation.cpp b/tests/src/test_sls_deformation.cpp new file mode 100644 index 000000000..c401f8025 --- /dev/null +++ b/tests/src/test_sls_deformation.cpp @@ -0,0 +1,333 @@ +#include "highmap/authoring.hpp" +#include "highmap/primitives.hpp" +#include "highmap/range.hpp" + +#include + +using namespace hmap; + +namespace +{ + +float helper_mse(const Array &a, const Array &b) +{ + float sum = 0.f; + for (size_t k = 0; k < a.vector.size(); ++k) + { + float d = a.vector[k] - b.vector[k]; + sum += d * d; + } + return sum / static_cast(a.vector.size()); +} + +// fraction of vertices on the right side of the water level w.r.t. the mask +float helper_mask_agreement(const Array &z, const Array &mask, float level) +{ + int count = 0; + for (size_t k = 0; k < z.vector.size(); ++k) + { + bool land = z.vector[k] > level; + bool want = mask.vector[k] > 0.5f; + if (land == want) ++count; + } + return static_cast(count) / static_cast(z.vector.size()); +} + +Array helper_disk_mask(glm::ivec2 shape, float radius_ratio) +{ + Array mask(shape, 0.f); + float cx = 0.5f * static_cast(shape.x - 1); + float cy = 0.5f * static_cast(shape.y - 1); + float r = radius_ratio * static_cast(std::min(shape.x, shape.y)); + for (int j = 0; j < shape.y; ++j) + for (int i = 0; i < shape.x; ++i) + { + float dx = static_cast(i) - cx; + float dy = static_cast(j) - cy; + if (dx * dx + dy * dy < r * r) mask(i, j) = 1.f; + } + return mask; +} + +} // namespace + +TEST(SlsDeformation, GaussianPushSetsCenterAndVanishesOutsideRadius) +{ + Array z({32, 32}, 0.f); + + std::vector pushes = {{10, 10, 4, 0.5f}}; + apply_gaussian_pushes(z, pushes); + + EXPECT_FLOAT_EQ(z(10, 10), 0.5f); + EXPECT_GT(z(10, 12), 0.f); + EXPECT_LT(z(10, 12), 0.5f); + EXPECT_FLOAT_EQ(z(10, 15), 0.f); + EXPECT_FLOAT_EQ(z(15, 10), 0.f); + EXPECT_FLOAT_EQ(z(0, 0), 0.f); +} + +TEST(SlsDeformation, GaussianPushIsClippedAtArrayBorder) +{ + Array z({16, 16}, 0.f); + + std::vector pushes = {{0, 0, 8, 1.f}}; + apply_gaussian_pushes(z, pushes); + + EXPECT_FLOAT_EQ(z(0, 0), 1.f); + EXPECT_GT(z(3, 3), 0.f); +} + +TEST(SlsDeformation, MatchConstraintReducesDistanceToTarget) +{ + glm::ivec2 shape = {64, 64}; + + Array z(shape, 0.f); + Array target = gaussian_pulse(shape, 12.f); + + DeformationConstraint c; + c.target = target; + c.weight = Array(shape, 1.f); + c.type = DeformationConstraintType::MATCH; + + float mse_before = helper_mse(z, target); + + Array z_out = sls_deformation(z, {c}, 1u, 100, 2, 16, 4); + + ASSERT_EQ(z_out.shape, shape); + float mse_after = helper_mse(z_out, target); + + EXPECT_LT(mse_after, 0.25f * mse_before); +} + +TEST(SlsDeformation, ShapeMaskConstraintsEnforceLandAndSea) +{ + glm::ivec2 shape = {64, 64}; + float level = 0.5f; + + Array z = noise(NoiseType::PERLIN, shape, {4.f, 4.f}, 3u); + remap(z, 0.f, 1.f); + + Array mask = helper_disk_mask(shape, 0.3f); + Array sea(shape, 1.f); + sea -= mask; + + DeformationConstraint land_c; + land_c.target = Array(shape, level); + land_c.weight = mask; + land_c.type = DeformationConstraintType::ABOVE; + + DeformationConstraint sea_c; + sea_c.target = Array(shape, level); + sea_c.weight = sea; + sea_c.type = DeformationConstraintType::BELOW; + + float agreement_before = helper_mask_agreement(z, mask, level); + + Array z_out = sls_deformation(z, {land_c, sea_c}, 7u, 200, 2, 16, 4); + + float agreement_after = helper_mask_agreement(z_out, mask, level); + + EXPECT_GT(agreement_after, agreement_before); + EXPECT_GT(agreement_after, 0.9f); +} + +TEST(SlsDeformation, RecordedPushesReplayExactly) +{ + glm::ivec2 shape = {48, 48}; + + Array z = noise(NoiseType::PERLIN, shape, {3.f, 3.f}, 5u); + Array target = gaussian_pulse(shape, 10.f); + + DeformationConstraint c; + c.target = target; + c.weight = Array(shape, 1.f); + + std::vector pushes; + + Array z_out = sls_deformation(z, + {c}, + 11u, + 50, + 2, + 12, + 3, + 0.f, + 256, + 0.65f, + 0.1f, + 1e-3f, + &pushes); + + EXPECT_FALSE(pushes.empty()); + + Array z_replay = z; + apply_gaussian_pushes(z_replay, pushes); + + ASSERT_EQ(z_replay.vector.size(), z_out.vector.size()); + for (size_t k = 0; k < z_out.vector.size(); ++k) + EXPECT_FLOAT_EQ(z_replay.vector[k], z_out.vector[k]); +} + +TEST(SlsDeformation, SameSeedIsDeterministic) +{ + glm::ivec2 shape = {48, 48}; + + Array z = noise(NoiseType::PERLIN, shape, {3.f, 3.f}, 5u); + + DeformationConstraint c; + c.target = gaussian_pulse(shape, 10.f); + c.weight = Array(shape, 1.f); + + Array z1 = sls_deformation(z, {c}, 42u, 30, 2, 12, 3); + Array z2 = sls_deformation(z, {c}, 42u, 30, 2, 12, 3); + + for (size_t k = 0; k < z1.vector.size(); ++k) + EXPECT_FLOAT_EQ(z1.vector[k], z2.vector[k]); +} + +TEST(SlsDeformation, SatisfiedConstraintLeavesTerrainUnchanged) +{ + glm::ivec2 shape = {32, 32}; + + Array z = noise(NoiseType::PERLIN, shape, {2.f, 2.f}, 9u); + + DeformationConstraint c; + c.target = z; + c.weight = Array(shape, 1.f); + + std::vector pushes; + + Array z_out = sls_deformation(z, + {c}, + 1u, + 20, + 2, + 8, + 3, + 0.f, + 128, + 0.65f, + 0.1f, + 1e-3f, + &pushes); + + EXPECT_TRUE(pushes.empty()); + for (size_t k = 0; k < z.vector.size(); ++k) + EXPECT_FLOAT_EQ(z_out.vector[k], z.vector[k]); +} + +TEST(SlsDeformation, SlopeLimitBoundsPushAmplitude) +{ + glm::ivec2 shape = {48, 48}; + float talus_max = 0.01f; + + Array z(shape, 0.f); + + DeformationConstraint c; + c.target = gaussian_pulse(shape, 10.f); + c.weight = Array(shape, 1.f); + + std::vector pushes; + sls_deformation(z, + {c}, + 3u, + 40, + 2, + 12, + 3, + talus_max, + 256, + 0.65f, + 0.1f, + 1e-3f, + &pushes); + + EXPECT_FALSE(pushes.empty()); + for (auto &p : pushes) + EXPECT_LE(std::abs(p.amplitude), + talus_max * static_cast(p.ir) + 1e-6f); +} + +TEST(SlsDeformation, EmptyConstraintsReturnInput) +{ + glm::ivec2 shape = {16, 16}; + Array z = noise(NoiseType::PERLIN, shape, {2.f, 2.f}, 1u); + + std::vector pushes; + + Array z_out = sls_deformation(z, + {}, + 1u, + 10, + 2, + 4, + 2, + 0.f, + 64, + 0.65f, + 0.1f, + 1e-3f, + &pushes); + + EXPECT_TRUE(pushes.empty()); + for (size_t k = 0; k < z.vector.size(); ++k) + EXPECT_FLOAT_EQ(z_out.vector[k], z.vector[k]); +} + +TEST(SlsDeformation, MismatchedConstraintShapeReturnsInput) +{ + glm::ivec2 shape = {16, 16}; + Array z = noise(NoiseType::PERLIN, shape, {2.f, 2.f}, 1u); + + DeformationConstraint c; + c.target = Array({8, 8}, 1.f); + c.weight = Array({8, 8}, 1.f); + + std::vector pushes; + + Array z_out = sls_deformation(z, + {c}, + 1u, + 10, + 2, + 4, + 2, + 0.f, + 64, + 0.65f, + 0.1f, + 1e-3f, + &pushes); + + EXPECT_TRUE(pushes.empty()); + for (size_t k = 0; k < z.vector.size(); ++k) + EXPECT_FLOAT_EQ(z_out.vector[k], z.vector[k]); +} + +TEST(SlsDeformation, ConstraintScaleWeightsCompetingTerms) +{ + glm::ivec2 shape = {48, 48}; + + Array z(shape, 0.f); + + // two competing full-weight targets: a bump and a flat plane + DeformationConstraint bump; + bump.target = gaussian_pulse(shape, 10.f); + bump.weight = Array(shape, 1.f); + + DeformationConstraint flat; + flat.target = Array(shape, 0.f); + flat.weight = Array(shape, 1.f); + + // heavily favour the bump + bump.scale = 100.f; + flat.scale = 1.f; + Array z_bump = sls_deformation(z, {bump, flat}, 1u, 60, 2, 12, 3); + + // heavily favour the plane + bump.scale = 1.f; + flat.scale = 100.f; + Array z_flat = sls_deformation(z, {bump, flat}, 1u, 60, 2, 12, 3); + + EXPECT_LT(helper_mse(z_bump, bump.target), helper_mse(z_flat, bump.target)); +}