diff --git a/Saved/Research/intervals/REPORT.md b/Saved/Research/intervals/REPORT.md index 870af83..ba29e6e 100644 --- a/Saved/Research/intervals/REPORT.md +++ b/Saved/Research/intervals/REPORT.md @@ -114,3 +114,68 @@ All final runs exited 0. No Unreal process was launched, no `Source/` file was m ## Next step Extend the CPU corpus until it emits the intended ~60 MB / 2.4 million triangle workload, with page meshing parallelized deterministically. In parallel, replace the far sheet proxy with actual selected parent room/tunnel/strata shapes, radial horizon clipping, and a shared near/far transition. Re-run the same validation and gates; then decide whether the Unreal publication/collision test is justified. + +## Round 2: bounded blends, QEF mesh, and surface relief + +Date: 2026-10-02. The existing section above records the first-round measurements; this section records the changed field and mesh. + +### What changed + +- The corpus now has **363 shapes**: 18 rooms, 147 anisotropic tunnel segments with changing cross-sections, 40 pillars, 36 stalactite/stalagmite forms, 42 shelves, 24 floating masses, 12 craters, 40 crowded-junction edits, three sparse rock bands, and the isolated audit seam. The 12 repeated thin strata were removed. The three remaining bands have different orientations and are placed at a nearby air lattice node so the interval mesh sees their thin crossings. +- Smooth min/max is applied only in bounded overlap regions. At 70 cm blend width, the mask covers a median **3.08% of cell centers** (3.07–3.19% across five seeds) and **1.94% of lattice nodes** (1.94–1.98%). Everywhere else the original exact ordered interval field is used. The banded features keep their hard field so their crossings remain exact. +- The near mesh now uses per-cell QEF dual contouring with cached analytic Hermite crossings/normals, face-connectivity component splits, a deterministic ambiguous-face sign, and deterministic page ownership. Page-edge ownership is checked across all 64 pages. The 12 cm audit seam, which no lattice point can see, gets a closed 336-triangle analytic patch. +- A deterministic four-band value-noise displacement runs after meshing, with wavelengths **6.4, 3.2, 1.6, and 1.2 m**. The shortest band spans 4.8 voxel cells. The hard displacement bound is **0.60 m**; the five measured seeds reached 0.404–0.439 m. Vertex normals are rebuilt after displacement. No extra volume samples are taken for relief. +- `interval_viewer2.exe` renders the blended dual-contour mesh and keeps a separate displaced vertex buffer. Press **R** to switch displacement off/on. The world-space rock color mottling and D3D11 back-buffer PNG capture are part of this standalone viewer. Capture mode closes its own window after writing the screenshots. + +### Timing comparison + +The first-round values below are the medians from the original five-seed report. Round 2 reports the median/p95/worst of five fresh sequential processes per worker count; p95 is linearly interpolated at the 95th percentile. The blend-index row is the AABB overlap raster plus surface-pair filtering. Near total includes the listed stages; far-field time is separate. + +| Stage (ms) | Before medians, 1 / 6 / 8 workers | After, 1 worker med / p95 / worst | After, 6 workers med / p95 / worst | After, 8 workers med / p95 / worst | +|---|---:|---:|---:|---:| +| Corpus/hierarchy, excluding blend-zone index | 0.5 / 0.5 / 0.5 | 0.96 / 1.03 / 1.04 | 0.96 / 0.97 / 0.97 | 0.96 / 1.08 / 1.10 | +| Blend-zone index and mask | — / — / — | 186.68 / 203.10 / 205.61 | 186.07 / 203.27 / 205.51 | 186.97 / 202.43 / 204.94 | +| Interval solve, compose, and sign words | 154.2 / 48.2 / 43.2 | 154.88 / 156.72 / 157.12 | 48.89 / 50.47 / 50.65 | 46.10 / 46.97 / 47.16 | +| Blend-zone pointwise evaluation and crossings | — / — / — | 203.85 / 225.44 / 229.87 | 200.79 / 205.81 / 205.98 | 201.49 / 208.96 / 209.65 | +| Active-cell masks | 2.1 / 2.1 / 2.1 | 2.62 / 3.82 / 3.86 | 2.85 / 4.05 / 4.28 | 3.58 / 4.29 / 4.29 | +| Near meshing: old contour emitter / QEF dual contour | 188.4 / 191.2 / 192.8 | 315.62 / 324.11 / 325.76 | 322.14 / 344.69 / 349.24 | 330.38 / 331.55 / 331.57 | +| Surface displacement and normal rebuild | — / — / — | 61.43 / 62.29 / 62.45 | 61.04 / 61.80 / 61.91 | 61.26 / 63.41 / 63.69 | +| **Near CPU total** | **349.3 / 241.1 / 240.2** | **931.93 / 957.68 / 959.08** | **836.06 / 848.28 / 849.57** | **824.27 / 854.37 / 861.33** | +| Far-field proxy | 8.6 / 8.6 / 8.6 | 8.59 / 8.61 / 8.62 | 8.57 / 9.51 / 9.58 | 8.58 / 8.60 / 8.61 | +| **Near plus far** | **357.9 / 249.7 / 248.8** | **940.54 / 966.28 / 967.68** | **844.62 / 857.37 / 858.80** | **832.85 / 862.93 / 869.89** | + +That is a **2.67× / 3.47× / 3.43×** near-total cost against round 1 at 1/6/8 workers. Most of the new cost is visible in the table: roughly 187–190 ms to build the overlap mask, 201–204 ms to evaluate blended nodes/rebuild crossings, 316–330 ms for dual contouring, and about 61 ms for relief plus normal rebuild. The 6/8-worker interval speedup remains, but the zone work and mesher are single-threaded here. + +The median output across the five seeds is **265,508 near vertices**, **531,164 triangles**, and **12,747,440 near-buffer bytes**. Median working sets were 69.6 / 70.0 / 70.4 MB at 1 / 6 / 8 workers. The proposed 60 MB mesh workload is still not exercised. + +### Correctness, topology, and repeatability + +- All 15 final processes (seeds 1–5 at 1, 6, and 8 workers) exited 0. Each compared the interval/blend signs against the new pointwise blended-field reference at every lattice node: **zero mismatches** and zero blend-boundary sign disagreements. All page ownership checks passed, with zero missing analytic crossings, zero cell-contour faults, and zero open cell contours. +- The three lattice-snapped bands are present in the structural mesh. The 12 cm audit seam has no lattice sample and is retained by its closed 336-triangle patch. The thin-surface omission check still rejects a deliberately omitted patch. +- The edge-manifold diagnostic reports **zero boundary edges** on all five seeds, but **3–7 overused edge groups** remain. The JSON `non_manifold_edges_zero` field is therefore false. These are localized dual-contour singularities around ambiguous/multiple-crossing cells; they are visible in the diagnostics and were not waived as a zero count. The page-ownership/sign/thin-feature validation result is reported separately from this remaining topology count. +- The QEF clamp counter is 141,729–155,046 across seeds. It increments per clamped coordinate axis, not per unique vertex. There are also 1,654–2,035 retained near-zero-area triangles; keeping them preserves edge connectivity after QEF clamping. These facts, plus the remaining overused edges, mean the mesher still needs topology work before treating the output as production collision geometry. +- For each seed, raw mesh, displaced mesh, and combined hashes match at all three worker counts. A fresh repeated seed-1/one-worker process also matched: mesh `20d4d80f439c17f6`, displaced mesh `0592d64cb4f269ca`, combined `fe846e12af9fef67`. Per-seed combined hashes: `fe846e12af9fef67`, `8f0994fe2295a34e`, `ee13abc933f14a3f`, `129dcc0c9bac3b3e`, `66361d9daff2b1a6`. + +### Screenshots and visual judgement + +These PNGs were captured from the viewer’s own D3D11 swap-chain back buffer, at 1280×800: + +| View | Capture | +|---|---| +| Inside the first chamber, relief on | [`chamber.png`](views/chamber.png) | +| Same camera, relief off | [`chamber_relief_off.png`](views/chamber_relief_off.png) | +| Open point in the crowded junction | [`junction.png`](views/junction.png) | +| Close to the chamber wall | [`wall.png`](views/wall.png) | + +Compared with the reject image, the repeated thin slab stack is gone, the color has rock-like variation, and many primitive junctions read as rounded transitions. The relief-on chamber is visibly less smooth than the relief-off view. The result is still well short of the reference cave: the close-wall view shows an obvious rounded ellipsoid, the chamber has broad flat-looking facets, and the material lacks the target’s rich rock texture, crystals, and water. I would describe this as a measurable improvement over the rejection, not a finished natural-rock look. + +### Failed and retried Round 2 attempts + +- Two early compile invocations passed the `vcvars64.bat` path with broken shell quoting, so `cl.exe` never started. An intermediate C++ build also failed because `near` was used as an identifier; it was renamed. The first D3D11 viewer link failed with unresolved `CoInitializeEx`/`CoUninitialize`/`CoCreateInstance`; linking `ole32.lib` fixed it. A later corpus-snap build failed once because `eval_quad` was used before declaration; a forward declaration fixed it. The corrected builds succeeded. MSVC continued to print that `vswhere.exe` was unavailable from `vcvars64.bat`, but it still launched `cl.exe`. +- `round2_probe3` through `round2_probe12` were intermediate dual-contour runs that exited on validation failures while the Hermite cache, component topology, degenerate triangles, and thin patches were being corrected. At probe 12, missing crossings had reached zero, but the mesh still reported about 3.7k non-manifold edges. Skipping the zero-area triangles produced thousands of boundary edges; a one-QEF-vertex-per-cell experiment made the overused-edge count worse (798), so both approaches were reverted. The final five-seed check has zero boundary edges and 3–7 overused groups. +- The first two D3D capture runs failed at WIC pixel-format selection (`0x88982f80`, unsupported RGBA); they left zero-byte PNGs. The capture path now swaps the D3D11 RGBA readback into WIC BGRA and writes all four non-empty images above. The capture-mode viewer closed itself; no previously open viewer was touched. +- MSVC continued to print that `vswhere.exe` was unavailable from `vcvars64.bat`, but all final builds completed. All 15 final benchmark runs and the repeated determinism run exited 0. No Unreal build or launch was done, no plugin `Source/` file changed, and nothing was committed, pushed, or stashed. + +### Round 2 conclusion + +Smooth joins and broad relief improve the rejected world, with only about 3.1% of cell centers in blend zones. The price is substantial: the near build rises from 241 ms to 836 ms at six workers and from 240 ms to 824 ms at eight workers. QEF clamping and the remaining 3–7 overused edges also keep this from being production-ready geometry. The viewer and measurements are ready for owner review; the screenshots show honestly that more corpus shaping/material work is needed to approach the reference. diff --git a/Saved/Research/intervals/interval_spike.cpp b/Saved/Research/intervals/interval_spike.cpp index f1cb70c..fb873b9 100644 --- a/Saved/Research/intervals/interval_spike.cpp +++ b/Saved/Research/intervals/interval_spike.cpp @@ -19,6 +19,7 @@ #include #include #include +#include #include #include @@ -44,6 +45,8 @@ constexpr int64_t kQScale = int64_t{1} << 28; constexpr int64_t kSubPerCm = 16; // canonical interval/sign lattice: 0.625 mm constexpr int64_t kSurfaceSubPerCm = 32; // crossing positions: 0.3125 mm constexpr int kCellStepCm = 25; +constexpr int kBlendWidthCm = 70; +constexpr float kReliefBoundMeters = 0.60f; constexpr int kGridMinCm = -3225; constexpr int kGridMaxCm = kGridMinCm + (kNodes - 1) * kCellStepCm; constexpr int kWorldOriginCm[3] = {-12000, 6000, -4500}; @@ -87,6 +90,7 @@ struct Crossing { uint8_t beforeRock = 0; uint8_t afterRock = 0; }; +constexpr uint32_t kBlendCrossingShape = std::numeric_limits::max(); struct LineRange { uint32_t offset = 0; @@ -141,15 +145,34 @@ struct MeshOutput { uint64_t thinFeaturesHiddenInLattice = 0; uint64_t thinFeaturesMeshed = 0; uint64_t pageBorderVertices = 0; + uint64_t qefClampedVertices = 0; + uint64_t cellContourFaults = 0; + uint64_t cachedHermiteFallbacks = 0; + uint64_t nonManifoldEdges = 0; + uint64_t boundaryEdges = 0; + uint64_t boundaryThinEdges = 0; + uint64_t boundaryNearGridEdge = 0; + uint64_t degenerateTriangles = 0; + uint64_t overusedEdges = 0; + uint64_t displacedHash = 0; + double maxDisplacementMeters = 0.0; }; struct Bin3D { std::vector ids; }; +struct BlendBox { int32_t lo[3]{}, hi[3]{}; }; struct ProjectionBins { std::vector ids[kBucketCount * kBucketCount]; }; struct Corpus { std::vector shapes; Bin3D bins3d[kBucketCount * kBucketCount * kBucketCount]; + Bin3D blendBins3d[kBucketCount * kBucketCount * kBucketCount]; ProjectionBins projections[3]; + std::vector blendBoxes; + std::vector blendNodeMask; + std::vector blendNodeIndices; + uint64_t blendZoneNodes = 0; + uint64_t blendZoneCells = 0; + double blendZoneBuildMs = 0.0; uint64_t seed = 0; uint32_t maxMembershipAudit = 0; }; @@ -243,6 +266,77 @@ void add_sphere(Corpus& c, int x, int y, int z, int r, Operation op, FeatureKind add_shape(c, x, y, z, r, r, r, op, kind, 0); } +int64_t eval_quad(const Shape& s, const int64_t d[3]); + +void build_blend_zones(Corpus& c) { + const int pad = (kBlendWidthCm + 1) / 2 + kCellStepCm; + c.blendBoxes.clear(); + for (size_t i = 0; i < c.shapes.size(); ++i) { + const Shape& a = c.shapes[i]; + if (a.kind == FeatureKind::Stratum) continue; + for (size_t j = i + 1; j < c.shapes.size(); ++j) { + const Shape& b = c.shapes[j]; + if (b.kind == FeatureKind::Stratum) continue; + BlendBox box{}; + bool intersects = true; + for (int axis = 0; axis < 3; ++axis) { + box.lo[axis] = std::max(a.lo[axis] - pad, b.lo[axis] - pad); + box.hi[axis] = std::min(a.hi[axis] + pad, b.hi[axis] + pad); + box.lo[axis] = std::max(box.lo[axis], kGridMinCm); + box.hi[axis] = std::min(box.hi[axis], kGridMaxCm); + intersects &= box.lo[axis] <= box.hi[axis]; + } + if (intersects) c.blendBoxes.push_back(box); + } + } + + const size_t rowCount = static_cast(kNodes) * kNodes; + std::vector> nodeRows(rowCount); + auto setRange = [](auto& words, int lo, int hi) { + const int firstWord=lo>>6, lastWord=hi>>6; + for (int w=firstWord; w<=lastWord; ++w) { + const int firstBit=(w==firstWord)?(lo&63):0; + const int lastBit=(w==lastWord)?(hi&63):63; + const uint64_t low=~uint64_t{0}<(z) * kNodes + y], nlo[0], nhi[0]); + } + + c.blendNodeMask.assign(static_cast(kNodes) * kNodes * kNodes, 0); + c.blendNodeIndices.clear(); + c.blendNodeIndices.reserve(300000); + c.blendZoneNodes = 0; + for (size_t row = 0; row < rowCount; ++row) { + const size_t base=row*kNodes; + for (int w=0; w(std::countr_zero(bits)); bits&=bits-1; + if (x>=kNodes) continue; + const size_t index=base+x; + c.blendNodeMask[index]=1; + c.blendNodeIndices.push_back(static_cast(index)); + ++c.blendZoneNodes; + } + } + } + c.blendZoneCells = 0; +} + +void filter_blend_zone_mask(Corpus& c); +size_t node_mask_index(int x, int y, int z); + void build_corpus(Corpus& c, uint64_t seed) { c.seed = seed; Rng rng(seed ^ 0x0A57A20261002ull); @@ -278,8 +372,12 @@ void build_corpus(Corpus& c, uint64_t seed) { const int x = (a[0] * (6 - t) + b[0] * t) / 6 + (seg == 2 ? cx - (a[0] + b[0]) / 2 : 0); const int y = (a[1] * (6 - t) + b[1] * t) / 6 + (seg == 2 ? cy - (a[1] + b[1]) / 2 : 0); const int z = (a[2] * (6 - t) + b[2] * t) / 6 + (seg == 2 ? cz - (a[2] + b[2]) / 2 : 0); - const int r = rng.range(200, 300); - add_sphere(c, x, y, z, r, Operation::CarveAir, FeatureKind::Tunnel); + const int r = rng.range(210, 295); + // Adjacent ellipsoids use changing cross-sections and axis orientation to + // avoid a chain of identical spherical tunnel beads. + add_shape(c, x, y, z, r, r * rng.range(82, 115) / 100, + r * rng.range(78, 112) / 100, Operation::CarveAir, + FeatureKind::Tunnel, rng.next() % 8); } } } @@ -288,10 +386,12 @@ void build_corpus(Corpus& c, uint64_t seed) { for (int iy = 0; iy < 3; ++iy) for (int ix = 0; ix < 3; ++ix) { const auto a = centers[0][iy][ix], b = centers[1][iy][ix]; for (int seg = 1; seg <= 3; ++seg) { - add_sphere(c, (a[0] * (4 - seg) + b[0] * seg) / 4, + const int r = rng.range(195, 265); + add_shape(c, (a[0] * (4 - seg) + b[0] * seg) / 4, (a[1] * (4 - seg) + b[1] * seg) / 4, (a[2] * (4 - seg) + b[2] * seg) / 4, - rng.range(190, 270), Operation::CarveAir, FeatureKind::Tunnel); + r, r * rng.range(80, 115) / 100, r * rng.range(78, 112) / 100, + Operation::CarveAir, FeatureKind::Tunnel, rng.next() % 8); } } @@ -339,17 +439,56 @@ void build_corpus(Corpus& c, uint64_t seed) { add_shape(c, x, y, z, r, rng.range(70, 200), rng.range(75, 185), op, FeatureKind::CrowdedJunction, rng.next() % 8); } - // Thin strata and shelves, with several widths below the 25 cm structural cell size. - for (int i = 0; i < 12; ++i) { + // Three widely separated, differently oriented rock bands replace the repeated thin stack. + auto snapToLattice = [](int value) { + const int node = static_cast(std::llround((value - kGridMinCm) / static_cast(kCellStepCm))); + return kGridMinCm + node * kCellStepCm; + }; + auto isAirBeforeStratum = [&](int x, int y, int z) { + const int64_t p[3] = {static_cast(x)*kSubPerCm, + static_cast(y)*kSubPerCm, + static_cast(z)*kSubPerCm}; + bool rock = true; + for (const Shape& s : c.shapes) { + bool insideAabb = true; + for (int a=0;a<3;++a) insideAabb &= p[a]>=static_cast(s.lo[a])*kSubPerCm + && p[a]<=static_cast(s.hi[a])*kSubPerCm; + if (!insideAabb) continue; + const int64_t d[3] = {p[0]-static_cast(s.c[0])*kSubPerCm, + p[1]-static_cast(s.c[1])*kSubPerCm, + p[2]-static_cast(s.c[2])*kSubPerCm}; + if (eval_quad(s,d) <= kQScale*kSubPerCm*kSubPerCm) rock = s.op==Operation::AddRock; + } + return !rock; + }; + for (int i = 0; i < 3; ++i) { const auto center = centers[i % 2][(i / 2) % 3][(i / 6) % 3]; - add_shape(c, center[0] + rng.range(-180, 180), center[1] + rng.range(-180, 180), - center[2] + rng.range(-180, 180), rng.range(300, 580), rng.range(250, 520), - rng.range(4, 11), Operation::AddRock, FeatureKind::Stratum, rng.next() % 8); + const int desired[3] = {center[0] + rng.range(-180,180),center[1] + rng.range(-180,180), + center[2] + rng.range(-180,180)}; + int placed[3] = {snapToLattice(desired[0]),snapToLattice(desired[1]),snapToLattice(desired[2])}; + if (!isAirBeforeStratum(placed[0],placed[1],placed[2])) { + bool found=false; + for (int radius=25;radius<=300&&!found;radius+=25) { + const int offsets[6][3]={{radius,0,0},{-radius,0,0},{0,radius,0},{0,-radius,0},{0,0,radius},{0,0,-radius}}; + for(const auto& offset:offsets){ + int candidate[3]={snapToLattice(desired[0]+offset[0]),snapToLattice(desired[1]+offset[1]), + snapToLattice(desired[2]+offset[2])}; + if(isAirBeforeStratum(candidate[0],candidate[1],candidate[2])){std::copy(candidate,candidate+3,placed);found=true;break;} + } + } + } + add_shape(c, placed[0], placed[1], placed[2], rng.range(300, 580), rng.range(250, 520), + rng.range(7, 12), Operation::AddRock, FeatureKind::Stratum, + std::array{4,2,3}[i]); } // A guaranteed isolated sub-cell air seam in solid rock. No 25 cm lattice node lies inside it. add_shape(c, 2900, -2900, 2887, 500, 400, 6, Operation::CarveAir, FeatureKind::Stratum, 0, true); + const auto blendBuildBegin = Clock::now(); + build_blend_zones(c); + c.blendZoneBuildMs = ms(blendBuildBegin, Clock::now()); + // Build deterministic 3D and axis-projection bins in feature order. for (uint32_t id = 0; id < c.shapes.size(); ++id) { const Shape& s = c.shapes[id]; @@ -362,6 +501,15 @@ void build_corpus(Corpus& c, uint64_t seed) { for (int bz = bz0; bz <= bz1; ++bz) for (int by = by0; by <= by1; ++by) for (int bx = bx0; bx <= bx1; ++bx) c.bins3d[(bz * kBucketCount + by) * kBucketCount + bx].ids.push_back(id); + const int blendLoX = std::clamp(floor_div(s.lo[0] - kBlendWidthCm - kGridMinCm, kBucketCm), 0, kBucketCount - 1); + const int blendHiX = std::clamp(floor_div(s.hi[0] + kBlendWidthCm - kGridMinCm, kBucketCm), 0, kBucketCount - 1); + const int blendLoY = std::clamp(floor_div(s.lo[1] - kBlendWidthCm - kGridMinCm, kBucketCm), 0, kBucketCount - 1); + const int blendHiY = std::clamp(floor_div(s.hi[1] + kBlendWidthCm - kGridMinCm, kBucketCm), 0, kBucketCount - 1); + const int blendLoZ = std::clamp(floor_div(s.lo[2] - kBlendWidthCm - kGridMinCm, kBucketCm), 0, kBucketCount - 1); + const int blendHiZ = std::clamp(floor_div(s.hi[2] + kBlendWidthCm - kGridMinCm, kBucketCm), 0, kBucketCount - 1); + for (int bz = blendLoZ; bz <= blendHiZ; ++bz) for (int by = blendLoY; by <= blendHiY; ++by) + for (int bx = blendLoX; bx <= blendHiX; ++bx) + c.blendBins3d[(bz * kBucketCount + by) * kBucketCount + bx].ids.push_back(id); for (int axis = 0; axis < 3; ++axis) { int a = (axis + 1) % 3, b = (axis + 2) % 3; int a0 = std::clamp(floor_div(s.lo[a] - kGridMinCm, kBucketCm), 0, kBucketCount - 1); @@ -372,6 +520,9 @@ void build_corpus(Corpus& c, uint64_t seed) { c.projections[axis].ids[ib * kBucketCount + ia].push_back(id); } } + const auto zoneFilterBegin = Clock::now(); + filter_blend_zone_mask(c); + c.blendZoneBuildMs += ms(zoneFilterBegin, Clock::now()); } int bucket_of_sub(int64_t qSub) { @@ -390,6 +541,67 @@ int64_t eval_quad(const Shape& s, const int64_t d[3]) { return value; } +double shape_phi_approx(const Shape& s, const int64_t pSub[3]) { + int64_t d[3] = {pSub[0] - static_cast(s.c[0]) * kSubPerCm, + pSub[1] - static_cast(s.c[1]) * kSubPerCm, + pSub[2] - static_cast(s.c[2]) * kSubPerCm}; + const int64_t q = eval_quad(s, d); + const double q2 = std::max(0.0, static_cast(q) + / static_cast(kQScale * kSubPerCm * kSubPerCm)); + const double radius = static_cast(std::min({s.radius[0], s.radius[1], s.radius[2]})); + return (1.0 - std::sqrt(q2)) * radius; +} + +void filter_blend_zone_mask(Corpus& c) { + c.blendZoneNodes = 0; + std::vector candidates = std::move(c.blendNodeIndices); + c.blendNodeIndices.clear(); + c.blendNodeIndices.reserve(candidates.size() / 4 + 1024); + for (uint32_t rawIndex : candidates) { + const size_t nodeIndex = rawIndex; + const int x = static_cast(nodeIndex % kNodes); + const int y = static_cast((nodeIndex / kNodes) % kNodes); + const int z = static_cast(nodeIndex / (kNodes * kNodes)); + uint8_t& marked = c.blendNodeMask[nodeIndex]; + if (!marked) continue; + const int64_t p[3] = { + static_cast(kGridMinCm + x * kCellStepCm) * kSubPerCm, + static_cast(kGridMinCm + y * kCellStepCm) * kSubPerCm, + static_cast(kGridMinCm + z * kCellStepCm) * kSubPerCm}; + const int bx = std::clamp(floor_div(kGridMinCm + x * kCellStepCm - kGridMinCm, kBucketCm), 0, kBucketCount - 1); + const int by = std::clamp(floor_div(kGridMinCm + y * kCellStepCm - kGridMinCm, kBucketCm), 0, kBucketCount - 1); + const int bz = std::clamp(floor_div(kGridMinCm + z * kCellStepCm - kGridMinCm, kBucketCm), 0, kBucketCount - 1); + uint32_t nearSurfaces = 0; + const auto& candidates = c.blendBins3d[(bz * kBucketCount + by) * kBucketCount + bx].ids; + for (uint32_t id : candidates) { + const Shape& s = c.shapes[id]; + if (s.kind == FeatureKind::Stratum) continue; + bool insideExpandedAabb=true; + for (int axis=0;axis<3;++axis) + insideExpandedAabb &= p[axis]>=static_cast(s.lo[axis]-kBlendWidthCm)*kSubPerCm + && p[axis]<=static_cast(s.hi[axis]+kBlendWidthCm)*kSubPerCm; + if (!insideExpandedAabb) continue; + if (std::abs(shape_phi_approx(s, p)) <= 0.5 * kBlendWidthCm && ++nearSurfaces >= 2) break; + } + marked = nearSurfaces >= 2 ? 1 : 0; + c.blendZoneNodes += marked; + if (marked) c.blendNodeIndices.push_back(static_cast(nodeIndex)); + } + std::vector cellMask(static_cast(kCoreCells) * kCoreCells * kCoreCells, 0); + for (uint32_t nodeIndex : c.blendNodeIndices) { + const int x = static_cast(nodeIndex % kNodes); + const int y = static_cast((nodeIndex / kNodes) % kNodes); + const int z = static_cast(nodeIndex / (kNodes * kNodes)); + for (int cz = z - 1; cz <= z; ++cz) for (int cy = y - 1; cy <= y; ++cy) for (int cx = x - 1; cx <= x; ++cx) { + if (cx < kHalo || cy < kHalo || cz < kHalo + || cx >= kHalo + kCoreCells || cy >= kHalo + kCoreCells || cz >= kHalo + kCoreCells) continue; + const size_t index = (static_cast(cx - kHalo) * kCoreCells + (cy - kHalo)) * kCoreCells + (cz - kHalo); + cellMask[index] = 1; + } + } + c.blendZoneCells = std::count(cellMask.begin(), cellMask.end(), uint8_t{1}); +} + bool contains_sub(const Shape& s, const int64_t pSub[3]) { int64_t d[3] = {pSub[0] - static_cast(s.c[0]) * kSubPerCm, pSub[1] - static_cast(s.c[1]) * kSubPerCm, @@ -421,6 +633,159 @@ bool point_rock(const Corpus& c, const int64_t pSub[3]) { return rock; } +double smooth_min(double a, double b, double k) { + const double h = std::max(k - std::abs(a - b), 0.0) / k; + return std::min(a, b) - h * h * k * 0.25; +} + +double smooth_max(double a, double b, double k) { + return -smooth_min(-a, -b, k); +} + +double blended_field(const Corpus& c, const int64_t pSub[3]) { + const int xcm = static_cast(floor_div(static_cast(pSub[0]), static_cast(kSubPerCm))); + const int ycm = static_cast(floor_div(static_cast(pSub[1]), static_cast(kSubPerCm))); + const int zcm = static_cast(floor_div(static_cast(pSub[2]), static_cast(kSubPerCm))); + const int bx = std::clamp(floor_div(xcm - kGridMinCm, kBucketCm), 0, kBucketCount - 1); + const int by = std::clamp(floor_div(ycm - kGridMinCm, kBucketCm), 0, kBucketCount - 1); + const int bz = std::clamp(floor_div(zcm - kGridMinCm, kBucketCm), 0, kBucketCount - 1); + const auto& bucket = c.blendBins3d[(bz * kBucketCount + by) * kBucketCount + bx].ids; + double field = 100000.0; // positive means rock; the untouched world is solid rock + for (uint32_t id : bucket) { + const Shape& s = c.shapes[id]; + bool insideExpandedAabb = true; + for (int axis = 0; axis < 3; ++axis) + insideExpandedAabb &= pSub[axis] >= static_cast(s.lo[axis] - kBlendWidthCm) * kSubPerCm + && pSub[axis] <= static_cast(s.hi[axis] + kBlendWidthCm) * kSubPerCm; + if (!insideExpandedAabb) continue; + int64_t d[3] = {pSub[0] - static_cast(s.c[0]) * kSubPerCm, + pSub[1] - static_cast(s.c[1]) * kSubPerCm, + pSub[2] - static_cast(s.c[2]) * kSubPerCm}; + const double q2 = std::max(0.0, static_cast(eval_quad(s, d)) + / static_cast(kQScale * kSubPerCm * kSubPerCm)); + const double radius = static_cast(std::min({s.radius[0], s.radius[1], s.radius[2]})); + double phi = (1.0 - std::sqrt(q2)) * radius; + if (eval_quad(s, d) == kQScale * kSubPerCm * kSubPerCm) phi = 1.0e-8; + if (s.kind == FeatureKind::Stratum) { + // The few retained thin bands keep hard, exact CSG so their sub-cell roots survive. + field = s.op == Operation::AddRock ? std::max(field, phi) : std::min(field, -phi); + } else { + field = s.op == Operation::AddRock + ? smooth_max(field, phi, static_cast(kBlendWidthCm)) + : smooth_min(field, -phi, static_cast(kBlendWidthCm)); + } + } + return field; +} + +void blended_field_gradient(const Corpus& c, const int64_t pSub[3], double& field, double gradient[3]) { + const int xcm = static_cast(floor_div(static_cast(pSub[0]), static_cast(kSubPerCm))); + const int ycm = static_cast(floor_div(static_cast(pSub[1]), static_cast(kSubPerCm))); + const int zcm = static_cast(floor_div(static_cast(pSub[2]), static_cast(kSubPerCm))); + const int bx = std::clamp(floor_div(xcm-kGridMinCm,kBucketCm),0,kBucketCount-1); + const int by = std::clamp(floor_div(ycm-kGridMinCm,kBucketCm),0,kBucketCount-1); + const int bz = std::clamp(floor_div(zcm-kGridMinCm,kBucketCm),0,kBucketCount-1); + const auto& bucket=c.blendBins3d[(bz*kBucketCount+by)*kBucketCount+bx].ids; + field=100000.0; gradient[0]=gradient[1]=gradient[2]=0.0; + for (uint32_t id : bucket) { + const Shape& s=c.shapes[id]; + bool inside=true; + for (int a=0;a<3;++a) inside &= pSub[a]>=static_cast(s.lo[a]-kBlendWidthCm)*kSubPerCm + && pSub[a]<=static_cast(s.hi[a]+kBlendWidthCm)*kSubPerCm; + if (!inside) continue; + const double d[3] = {pSub[0]/static_cast(kSubPerCm)-s.c[0], + pSub[1]/static_cast(kSubPerCm)-s.c[1], + pSub[2]/static_cast(kSubPerCm)-s.c[2]}; + int64_t dSub[3] = {pSub[0]-static_cast(s.c[0])*kSubPerCm, + pSub[1]-static_cast(s.c[1])*kSubPerCm, + pSub[2]-static_cast(s.c[2])*kSubPerCm}; + const double q2=std::max(0.0,static_cast(eval_quad(s,dSub))/static_cast(kQScale*kSubPerCm*kSubPerCm)); + const double root=std::sqrt(q2); + const double radius=static_cast(std::min({s.radius[0],s.radius[1],s.radius[2]})); + double phi=(1.0-root)*radius; + double phiGrad[3]{}; + if (root>1.0e-12) for (int i=0;i<3;++i) { + double qd=0.0; + for (int j=0;j<3;++j) qd += (static_cast(s.q[i][j])/kQScale)*d[j]; + phiGrad[i]=-radius*qd/root; + } + if (eval_quad(s,dSub)==kQScale*kSubPerCm*kSubPerCm) phi=1.0e-8; + double b=phi, bg[3]={phiGrad[0],phiGrad[1],phiGrad[2]}; + const bool isMin=s.op==Operation::CarveAir; + if (isMin) { b=-b; for (double& v:bg) v=-v; } + if (s.kind==FeatureKind::Stratum) { + const bool take = isMin ? bfield; + if (take) { field=b; for (int a=0;a<3;++a) gradient[a]=bg[a]; } + continue; + } + const double k=static_cast(kBlendWidthCm); + const double weightA = isMin + ? std::clamp(0.5+(b-field)/(2.0*k),0.0,1.0) + : std::clamp(0.5+(field-b)/(2.0*k),0.0,1.0); + const double delta=std::abs(field-b); + if (delta>=k) { + const bool take = isMin ? bfield; + if (take) { field=b; for (int a=0;a<3;++a) gradient[a]=bg[a]; } + } else { + field=weightA*field+(1.0-weightA)*b-k*weightA*(1.0-weightA); + for (int a=0;a<3;++a) gradient[a]=weightA*gradient[a]+(1.0-weightA)*bg[a]; + } + } +} + +size_t node_mask_index(int x, int y, int z) { + return (static_cast(z) * kNodes + y) * kNodes + x; +} + +bool point_in_blend_zone(const Corpus& c, const int64_t pSub[3]) { + const int pad = (kBlendWidthCm + 1) / 2 + kCellStepCm; + uint32_t hits = 0; + for (const Shape& s : c.shapes) { + if (s.kind == FeatureKind::Stratum) continue; + bool inside = true; + for (int axis = 0; axis < 3; ++axis) + inside &= pSub[axis] >= static_cast(s.lo[axis] - pad) * kSubPerCm + && pSub[axis] <= static_cast(s.hi[axis] + pad) * kSubPerCm; + if (inside && ++hits >= 2) return true; + } + return false; +} + +bool point_rock_reference(const Corpus& c, const int64_t pSub[3], bool blendZone) { + return blendZone ? blended_field(c, pSub) >= 0.0 : point_rock(c, pSub); +} + +bool point_rock_reference(const Corpus& c, const int64_t pSub[3]) { + return point_rock_reference(c, pSub, point_in_blend_zone(c, pSub)); +} + +double hard_field(const Corpus& c, const int64_t pSub[3]) { + const int xcm = static_cast(floor_div(static_cast(pSub[0]), static_cast(kSubPerCm))); + const int ycm = static_cast(floor_div(static_cast(pSub[1]), static_cast(kSubPerCm))); + const int zcm = static_cast(floor_div(static_cast(pSub[2]), static_cast(kSubPerCm))); + const int bx = std::clamp(floor_div(xcm - kGridMinCm, kBucketCm), 0, kBucketCount - 1); + const int by = std::clamp(floor_div(ycm - kGridMinCm, kBucketCm), 0, kBucketCount - 1); + const int bz = std::clamp(floor_div(zcm - kGridMinCm, kBucketCm), 0, kBucketCount - 1); + const auto& bucket = c.blendBins3d[(bz * kBucketCount + by) * kBucketCount + bx].ids; + double field = 100000.0; + for (uint32_t id : bucket) { + const Shape& s = c.shapes[id]; + bool insideExpandedAabb = true; + for (int axis = 0; axis < 3; ++axis) + insideExpandedAabb &= pSub[axis] >= static_cast(s.lo[axis] - kBlendWidthCm) * kSubPerCm + && pSub[axis] <= static_cast(s.hi[axis] + kBlendWidthCm) * kSubPerCm; + if (!insideExpandedAabb) continue; + const double phi = shape_phi_approx(s, pSub); + field = s.op == Operation::AddRock ? std::max(field, phi) : std::min(field, -phi); + } + const bool exactRock = point_rock(c, pSub); + if ((field >= 0.0) != exactRock) { + const double magnitude = std::max(std::abs(field), 1.0e-8); + field = exactRock ? magnitude : -magnitude; + } + return field; +} + struct Interval { int64_t loSub, hiSub; int loSample, hiSample; }; int64_t eval_line(const Shape& s, int axis, int64_t variableSub, @@ -748,6 +1113,137 @@ uint64_t extract64(const AxisData& a, uint32_t line, int bit) { return v; } +struct BlendStats { + uint64_t nodeEvaluations = 0; + uint64_t edgeSegments = 0; + uint64_t blendedCrossings = 0; + uint64_t boundarySignDisagreements = 0; +}; + +void apply_blend_zones(const Corpus& c, AxisData axes[3], AxisStats stats[3], BlendStats& out) { + const int64_t baseSub = static_cast(kGridMinCm) * kSubPerCm; + const int64_t stepSub = static_cast(kCellStepCm) * kSubPerCm; + std::vector blendSigns(c.blendNodeMask.size(), 0); + std::array,3> affectedLines; + for (auto& v : affectedLines) v.assign(static_cast(kNodes) * kNodes, 0); + for (uint32_t rawIndex : c.blendNodeIndices) { + const size_t index = rawIndex; + const int x = static_cast(rawIndex % kNodes); + const int y = static_cast((rawIndex / kNodes) % kNodes); + const int z = static_cast(rawIndex / (kNodes * kNodes)); + const int64_t p[3] = { + static_cast(kGridMinCm + x * kCellStepCm) * kSubPerCm, + static_cast(kGridMinCm + y * kCellStepCm) * kSubPerCm, + static_cast(kGridMinCm + z * kCellStepCm) * kSubPerCm}; + blendSigns[index] = blended_field(c, p) >= 0.0; + for (int axis = 0; axis < 3; ++axis) affectedLines[axis][line_id(axis, x, y, z)] = 1; + ++out.nodeEvaluations; + } + for (int axis = 0; axis < 3; ++axis) { + AxisData& data = axes[axis]; + std::vector oldCrossings = std::move(data.crossings); + std::vector rebuilt; + rebuilt.reserve(oldCrossings.size() + 4096); + std::vector counts(data.lines.size(), 0); + for (uint32_t li = 0; li < data.lines.size(); ++li) { + LineRange& line = data.lines[li]; + const uint32_t oldOffset = line.offset, oldCount = line.count; + line.offset = static_cast(rebuilt.size()); + if (!affectedLines[axis][li]) { + rebuilt.insert(rebuilt.end(), oldCrossings.begin() + oldOffset, + oldCrossings.begin() + oldOffset + oldCount); + counts[li] = line.count; + continue; + } + int p0, p1; decode_line(axis, li, p0, p1); + const int a1 = (axis + 1) % 3, a2 = (axis + 2) % 3; + int node[3]{}; node[a1] = p0; node[a2] = p1; + bool zone[kNodes]{}, signs[kNodes]{}; + for (int sample = 0; sample < kNodes; ++sample) { + node[axis] = sample; + zone[sample] = c.blendNodeMask[node_mask_index(node[0], node[1], node[2])] != 0; + signs[sample] = get_sign(data, li, sample) != 0; + if (zone[sample]) { + signs[sample] = blendSigns[node_mask_index(node[0], node[1], node[2])] != 0; + } + uint64_t& word = data.signs[static_cast(li) * kWordsPerLine + (sample >> 6)]; + const uint64_t bit = uint64_t{1} << (sample & 63); + if (signs[sample]) word |= bit; else word &= ~bit; + } + uint32_t oldIndex = 0; + for (int edge = 0; edge + 1 < kNodes; ++edge) { + const int64_t startSub = baseSub + edge * stepSub; + const int64_t endSub = startSub + stepSub; + const int32_t start32 = static_cast(startSub * 2); + const int32_t end32 = static_cast(endSub * 2); + const bool edgeZone = zone[edge] || zone[edge + 1]; + if (edgeZone) ++out.edgeSegments; + while (oldIndex < oldCount + && oldCrossings[oldOffset + oldIndex].positionSub32 <= start32) { + ++oldIndex; + } + uint32_t scan = oldIndex; + for (; scan < oldCount && oldCrossings[oldOffset + scan].positionSub32 < end32; ++scan) { + const Crossing& old = oldCrossings[oldOffset + scan]; + if (!edgeZone) rebuilt.push_back(old); + else if (old.shapeId < c.shapes.size() + && c.shapes[old.shapeId].kind == FeatureKind::Stratum) rebuilt.push_back(old); + } + oldIndex = scan; + if (!edgeZone || signs[edge] == signs[edge + 1]) continue; + int64_t pLeft[3] = { + static_cast(kGridMinCm + node[0] * kCellStepCm) * kSubPerCm, + static_cast(kGridMinCm + node[1] * kCellStepCm) * kSubPerCm, + static_cast(kGridMinCm + node[2] * kCellStepCm) * kSubPerCm}; + int64_t pRight[3] = {pLeft[0], pLeft[1], pLeft[2]}; + pLeft[axis] = startSub; pRight[axis] = endSub; + const double leftValue = zone[edge] ? blended_field(c, pLeft) : hard_field(c, pLeft); + const double rightValue = zone[edge + 1] ? blended_field(c, pRight) : hard_field(c, pRight); + const bool leftFieldRock = leftValue >= 0.0; + const bool rightFieldRock = rightValue >= 0.0; + if (leftFieldRock != signs[edge] || rightFieldRock != signs[edge + 1]) + ++out.boundarySignDisagreements; + if (leftFieldRock == rightFieldRock) continue; + const bool beforeRock = leftFieldRock; + double rootSub = static_cast(startSub); + if (zone[edge] && zone[edge + 1]) { + int64_t lo = startSub, hi = endSub; + while (hi - lo > 4) { + const int64_t mid = lo + (hi - lo) / 2; + int64_t p[3] = {pLeft[0], pLeft[1], pLeft[2]}; p[axis] = mid; + const bool midRock = blended_field(c, p) >= 0.0; + if (midRock == beforeRock) lo = mid; else hi = mid; + } + rootSub = (lo + hi) * 0.5; + } else { + const double denominator = leftValue - rightValue; + const double t = std::clamp(std::abs(denominator) > 1.0e-12 + ? leftValue / denominator : 0.5, 0.0, 1.0); + rootSub = startSub + t * static_cast(stepSub); + } + int64_t p32 = std::llround(rootSub * 2.0); + p32=std::clamp(p32,static_cast(start32)+1,static_cast(end32)-1); + if (p32 < std::numeric_limits::min() || p32 > std::numeric_limits::max()) continue; + rebuilt.push_back({static_cast(p32), kBlendCrossingShape, + static_cast(beforeRock), static_cast(!beforeRock)}); + ++out.blendedCrossings; + } + std::sort(rebuilt.begin() + line.offset, rebuilt.end(), [](const Crossing& a, const Crossing& b) { + return a.positionSub32 < b.positionSub32 + || (a.positionSub32 == b.positionSub32 && a.shapeId < b.shapeId); + }); + line.count = static_cast(std::min(rebuilt.size() - line.offset, 65535)); + counts[li] = line.count; + } + data.crossings.swap(rebuilt); + std::sort(counts.begin(), counts.end()); + stats[axis].crossingTotal = std::accumulate(counts.begin(), counts.end(), uint64_t{0}); + stats[axis].crossingP50 = counts[counts.size() / 2]; + stats[axis].crossingP95 = counts[static_cast(0.95 * (counts.size() - 1))]; + stats[axis].crossingMax = counts.back(); + } +} + uint64_t index_hash_bytes(uint64_t h, const void* data, size_t bytes) { const auto* p = static_cast(data); for (size_t i = 0; i < bytes; ++i) { h ^= p[i]; h *= 1099511628211ull; } @@ -800,7 +1296,7 @@ bool find_edge_crossing(const AxisData& axisData, int axis, int x, int y, int z, for (uint32_t i = 0; i < lr.count; ++i) { const Crossing& c = axisData.crossings[lr.offset + i]; if (c.positionSub32 <= start32 || c.positionSub32 >= end32) continue; - if (c.beforeRock == leftRock && c.afterRock == rightRock) { + if ((c.beforeRock != 0) == leftRock && (c.afterRock != 0) == rightRock) { out = c; ++found; } } @@ -808,8 +1304,25 @@ bool find_edge_crossing(const AxisData& axisData, int axis, int x, int y, int z, return found != 0; } -void normal_for_crossing(const Shape& s, const Crossing& c, int axis, int local[3], +void normal_for_crossing(const Corpus& corpus, const Crossing& c, int axis, int local[3], double posCm, double outN[3]) { + if (c.shapeId == kBlendCrossingShape) { + int64_t p[3] = {static_cast(local[0]) * kSubPerCm, + static_cast(local[1]) * kSubPerCm, + static_cast(local[2]) * kSubPerCm}; + p[axis] = std::llround(posCm * kSubPerCm); + double field{}; + blended_field_gradient(corpus,p,field,outN); + const bool wantPositiveAxis = c.beforeRock && !c.afterRock; + const double sign = ((outN[axis] >= 0.0) == wantPositiveAxis) ? 1.0 : -1.0; + double len2 = 0.0; + for (double& v : std::span(outN,3)) { v *= sign; len2 += v * v; } + if (len2 <= 1.0e-28) { outN[0] = (axis == 0 ? 1.0 : 0.0); outN[1] = (axis == 1 ? 1.0 : 0.0); outN[2] = (axis == 2 ? 1.0 : 0.0); return; } + const double inv = 1.0 / std::sqrt(len2); + for (double& v : std::span(outN,3)) v *= inv; + return; + } + const Shape& s = corpus.shapes[c.shapeId]; double d[3] = {static_cast(local[0] - s.c[0]), static_cast(local[1] - s.c[1]), static_cast(local[2] - s.c[2])}; @@ -845,13 +1358,12 @@ MeshOutput emit_near_mesh(const Corpus& c, const AxisData axes[3], const AxisDat const uint64_t key = edge_key(e.axis, ex, ey, ez); auto it = edgeToVertex.find(key); if (it != edgeToVertex.end()) return it->second; - const Shape& s = c.shapes[crossing.shapeId]; int p[3] = {kGridMinCm + ex * kCellStepCm, kGridMinCm + ey * kCellStepCm, kGridMinCm + ez * kCellStepCm}; const double localPosCm = crossing.positionSub32 / static_cast(kSurfaceSubPerCm); p[e.axis] = static_cast(std::floor(localPosCm)); - double n[3]; normal_for_crossing(s, crossing, e.axis, p, localPosCm, n); + double n[3]; normal_for_crossing(c, crossing, e.axis, p, localPosCm, n); const double world[3] = { kWorldOriginCm[0] + (e.axis == 0 ? localPosCm : (kGridMinCm + ex * kCellStepCm)), kWorldOriginCm[1] + (e.axis == 1 ? localPosCm : (kGridMinCm + ey * kCellStepCm)), @@ -1109,6 +1621,432 @@ MeshOutput emit_near_mesh(const Corpus& c, const AxisData axes[3], const AxisDat return out; } +uint32_t cell_linear_key(int x, int y, int z) { + return static_cast((static_cast(x) * kCoreCells + y) * kCoreCells + z); +} + +uint32_t owner_page_for_cell(int x, int y, int z) { + return static_cast(((z / kPageCells) * kPagesPerAxis + y / kPageCells) * kPagesPerAxis + x / kPageCells); +} + +bool solve_qef3(double a[3][3], double b[3], double out[3]) { + double m[3][4]{}; + for (int r = 0; r < 3; ++r) { for (int c = 0; c < 3; ++c) m[r][c] = a[r][c]; m[r][3] = b[r]; } + for (int c = 0; c < 3; ++c) { + int pivot = c; + for (int r = c + 1; r < 3; ++r) if (std::abs(m[r][c]) > std::abs(m[pivot][c])) pivot = r; + if (std::abs(m[pivot][c]) < 1.0e-12) return false; + if (pivot != c) for (int j = c; j < 4; ++j) std::swap(m[pivot][j], m[c][j]); + const double inv = 1.0 / m[c][c]; + for (int j = c; j < 4; ++j) m[c][j] *= inv; + for (int r = 0; r < 3; ++r) if (r != c) { + const double f = m[r][c]; + for (int j = c; j < 4; ++j) m[r][j] -= f * m[c][j]; + } + } + for (int i = 0; i < 3; ++i) out[i] = m[i][3]; + return true; +} + +MeshOutput emit_dual_contour_mesh(const Corpus& c, const AxisData axes[3], const AxisData& zSigns, + const std::vector& active) { + MeshOutput out; + const size_t estimated = 300000; + out.vertices.reserve(estimated); + struct CellTopology { std::array edgeVertex; CellTopology() { edgeVertex.fill(std::numeric_limits::max()); } }; + std::unordered_map cellTopology; + cellTopology.reserve(estimated); + std::array, kPagesPerAxis * kPagesPerAxis * kPagesPerAxis> pageIndices; + uint64_t multiEdges = 0; + struct Hermite { Crossing crossing{}; double normal[3]{}; }; + std::unordered_map hermites; + hermites.reserve(axes[0].crossings.size()+axes[1].crossings.size()+axes[2].crossings.size()); + const int64_t domainStart32 = static_cast(kGridMinCm) * kSurfaceSubPerCm; + const int64_t edgeSpan32 = static_cast(kCellStepCm) * kSurfaceSubPerCm; + for (int axis = 0; axis < 3; ++axis) for (uint32_t li = 0; li < axes[axis].lines.size(); ++li) { + int p0, p1; decode_line(axis, li, p0, p1); + const int a1 = (axis + 1) % 3, a2 = (axis + 2) % 3; + for (uint32_t i = 0; i < axes[axis].lines[li].count; ++i) { + const Crossing& crossing = axes[axis].crossings[axes[axis].lines[li].offset+i]; + const int edge = static_cast(floor_div64(static_cast(crossing.positionSub32)-domainStart32, edgeSpan32)); + if (edge < 0 || edge >= kNodes-1) continue; + int coord[3]{}; coord[axis]=edge; coord[a1]=p0; coord[a2]=p1; + int localCm[3] = {kGridMinCm+coord[0]*kCellStepCm,kGridMinCm+coord[1]*kCellStepCm,kGridMinCm+coord[2]*kCellStepCm}; + const double posCm = crossing.positionSub32/static_cast(kSurfaceSubPerCm); + localCm[axis] = static_cast(std::floor(posCm)); + Hermite h{}; h.crossing=crossing; + normal_for_crossing(c,crossing,axis,localCm,posCm,h.normal); + const uint64_t key = edge_key(axis,coord[0],coord[1],coord[2]) + | (crossing.beforeRock ? (uint64_t{1}<<62) : 0); + auto [it,inserted] = hermites.emplace(key,h); + if (!inserted) { ++multiEdges; it->second=h; } + } + } + auto addQefComponent = [&](int ix,int iy,int iz,int cellX,int cellY,int cellZ,int component, + const std::vector& loop,const Hermite edgeHermite[12],CellTopology& topology) { + double ata[3][3]{},atb[3]{},avg[3]{},normalSum[3]{}; + for (int ei : loop) { + const Hermite& h=edgeHermite[ei]; + const EdgeDesc& e=kEdges[ei]; + const int ex=cellX+e.x,ey=cellY+e.y,ez=cellZ+e.z; + const double crossingCm=h.crossing.positionSub32/static_cast(kSurfaceSubPerCm); + const double* n=h.normal; + double hp[3]={static_cast(kGridMinCm+ex*kCellStepCm-(kGridMinCm+cellX*kCellStepCm)), + static_cast(kGridMinCm+ey*kCellStepCm-(kGridMinCm+cellY*kCellStepCm)), + static_cast(kGridMinCm+ez*kCellStepCm-(kGridMinCm+cellZ*kCellStepCm))}; + hp[e.axis]=crossingCm-(kGridMinCm+(e.axis==0?cellX:e.axis==1?cellY:cellZ)*kCellStepCm); + const double ndotp=n[0]*hp[0]+n[1]*hp[1]+n[2]*hp[2]; + for(int r=0;r<3;++r){atb[r]+=n[r]*ndotp;avg[r]+=hp[r];normalSum[r]+=n[r];for(int col=0;col<3;++col)ata[r][col]+=n[r]*n[col];} + } + if(loop.empty())return; + const int count=static_cast(loop.size()); + for(double& v:avg)v/=count; + const double regularizer=1.0e-5*count; + for(int a=0;a<3;++a){ata[a][a]+=regularizer;atb[a]+=regularizer*avg[a];} + double qef[3]{}; + if(!solve_qef3(ata,atb,qef))for(int a=0;a<3;++a)qef[a]=avg[a]; + for(int a=0;a<3;++a){if(qef[a]<0.0||qef[a]>kCellStepCm)++out.qefClampedVertices;qef[a]=std::clamp(qef[a],0.0,static_cast(kCellStepCm));} + double nlen=std::sqrt(normalSum[0]*normalSum[0]+normalSum[1]*normalSum[1]+normalSum[2]*normalSum[2]); + if(nlen<1.0e-20){normalSum[0]=0;normalSum[1]=0;normalSum[2]=1;nlen=1;} + const double wx=(kWorldOriginCm[0]+kGridMinCm+cellX*kCellStepCm+qef[0])*0.01; + const double wy=(kWorldOriginCm[1]+kGridMinCm+cellY*kCellStepCm+qef[1])*0.01; + const double wz=(kWorldOriginCm[2]+kGridMinCm+cellZ*kCellStepCm+qef[2])*0.01; + const uint32_t page=owner_page_for_cell(ix,iy,iz); + const uint64_t cellId=cell_linear_key(ix,iy,iz); + const uint64_t key=(uint64_t{4}<<60)|((cellId<<4)|static_cast(component&15)); + const uint32_t vertex=static_cast(out.vertices.size()); + out.vertices.push_back({key,page,{static_cast(wx),static_cast(wy),static_cast(wz), + static_cast(normalSum[0]/nlen),static_cast(normalSum[1]/nlen),static_cast(normalSum[2]/nlen)}}); + for(int ei:loop)topology.edgeVertex[ei]=vertex; + if((ix>0&&ix%kPageCells==0)||(ix+10&&iy%kPageCells==0)||(iy+10&&iz%kPageCells==0)||(iz+1(ix) * kCoreCells + iy) * 4 + block]; + while (mask) { + const unsigned bit=std::countr_zero(mask);mask&=mask-1; + const int iz=block*64+static_cast(bit),cellX=ix+kHalo,cellY=iy+kHalo,cellZ=iz+kHalo; + bool rock[8]{}; + for(int k=0;k<8;++k){const int x=cellX+kCornerOffset[k][0],y=cellY+kCornerOffset[k][1],z=cellZ+kCornerOffset[k][2];rock[k]=get_sign(zSigns,line_id(2,x,y,z),z)!=0;} + bool edgeCross[12]{}; + Hermite edgeHermite[12]{}; + for(int ei=0;ei<12;++ei){ + const EdgeDesc& e=kEdges[ei];if(rock[e.c0]==rock[e.c1])continue; + const int ex=cellX+e.x,ey=cellY+e.y,ez=cellZ+e.z; + const bool c0low=kCornerOffset[e.c0][e.axis]second;edgeCross[ei]=true;continue;} + Crossing crossing{}; + if(!find_edge_crossing(axes[e.axis],e.axis,ex,ey,ez,lowRock,highRock,crossing,multiEdges)){ + ++out.missingAnalyticCrossings; + continue; + } + const double pos=crossing.positionSub32/static_cast(kSurfaceSubPerCm); + int localCm[3]={kGridMinCm+ex*kCellStepCm,kGridMinCm+ey*kCellStepCm,kGridMinCm+ez*kCellStepCm};localCm[e.axis]=static_cast(std::floor(pos)); + edgeHermite[ei].crossing=crossing;normal_for_crossing(c,crossing,e.axis,localCm,pos,edgeHermite[ei].normal); + edgeCross[ei]=true;++out.cachedHermiteFallbacks; + } + std::array,12> adj{};std::array degree{}; + auto connect=[&](int a,int b){if(degree[a]<2&°ree[b]<2){adj[a][degree[a]++]=b;adj[b][degree[b]++]=a;}}; + for(int f=0;f<6;++f){ + int crossed[4],nCross=0;for(int k=0;k<4;++k)if(edgeCross[kFaceEdges[f][k]])crossed[nCross++]=k; + if(nCross==2)connect(kFaceEdges[f][crossed[0]],kFaceEdges[f][crossed[1]]); + else if(nCross==4){ + ++out.ambiguousFaces; + int64_t q[3]={static_cast(kGridMinCm+cellX*kCellStepCm)*kSubPerCm, + static_cast(kGridMinCm+cellY*kCellStepCm)*kSubPerCm, + static_cast(kGridMinCm+cellZ*kCellStepCm)*kSubPerCm}; + int sums[3]{};for(int k=0;k<4;++k)for(int a=0;a<3;++a)sums[a]+=kCornerOffset[kFaceCorners[f][k]][a]; + for(int a=0;a<3;++a)q[a]+=sums[a]*(kCellStepCm*kSubPerCm/4); + const bool centerRock=point_rock_reference(c,q),s0=rock[kFaceCorners[f][0]]; + if(centerRock==s0){connect(kFaceEdges[f][0],kFaceEdges[f][1]);connect(kFaceEdges[f][2],kFaceEdges[f][3]);} + else{connect(kFaceEdges[f][3],kFaceEdges[f][0]);connect(kFaceEdges[f][1],kFaceEdges[f][2]);} + } + } + CellTopology topology; + std::array visited{};int component=0; + for(int start=0;start<12;++start){ + if(!edgeCross[start]||visited[start])continue; + if(degree[start]!=2){++out.contourDegreeFaults;visited[start]=1;continue;} + std::vector loop;int prev=-1,cur=start; + for(int guard=0;guard<16;++guard){if(visited[cur])break;visited[cur]=1;loop.push_back(cur);const int next=adj[cur][0]==prev?adj[cur][1]:adj[cur][0];prev=cur;cur=next;if(cur==start)break;} + if(cur!=start){++out.openContourLoops;continue;} + if(loop.size()>=3)addQefComponent(ix,iy,iz,cellX,cellY,cellZ,component++,loop,edgeHermite,topology); + } + if(component>0)cellTopology.emplace(cell_linear_key(ix,iy,iz),std::move(topology)); + } + } + } + auto getCell = [&](int x,int y,int z,int localEdge)->uint32_t { + auto it=cellTopology.find(cell_linear_key(x,y,z)); + return it==cellTopology.end()?std::numeric_limits::max():it->second.edgeVertex[localEdge]; + }; + auto emitTriangle = [&](uint32_t page, uint32_t a, uint32_t b, uint32_t d, int axis, double sign) { + const Vertex& va = out.vertices[a].v; const Vertex& vb = out.vertices[b].v; const Vertex& vd = out.vertices[d].v; + const double ux = vb.x - va.x, uy = vb.y - va.y, uz = vb.z - va.z; + const double vx = vd.x - va.x, vy = vd.y - va.y, vz = vd.z - va.z; + const double n[3] = {uy*vz-uz*vy, uz*vx-ux*vz, ux*vy-uy*vx}; + if (std::abs(n[0]) + std::abs(n[1]) + std::abs(n[2]) < 1.0e-12) ++out.degenerateTriangles; + if (n[axis] * sign < 0.0) std::swap(b, d); + auto& indices = pageIndices[page]; + indices.push_back(a); indices.push_back(b); indices.push_back(d); + }; + auto emitQuad = [&](int axis, int xNode, int yNode, int zNode, std::array,4> cells, + std::array localEdges,bool lowRock) { + uint32_t v[4]; + for (int i = 0; i < 4; ++i) v[i] = getCell(cells[i][0], cells[i][1], cells[i][2],localEdges[i]); + if (std::any_of(std::begin(v), std::end(v), [](uint32_t q) { return q == std::numeric_limits::max(); })) { + ++out.contourDegreeFaults; return; + } + const uint32_t page = owner_page_for_edge(axis, xNode, yNode, zNode); + const double sign = lowRock ? 1.0 : -1.0; + emitTriangle(page, v[0], v[1], v[2], axis, sign); + emitTriangle(page, v[0], v[2], v[3], axis, sign); + }; + std::unordered_set emittedEdges; + emittedEdges.reserve(300000); + for (int axis = 0; axis < 3; ++axis) { + const int a1 = (axis+1)%3, a2 = (axis+2)%3; + for (uint32_t li = 0; li < axes[axis].lines.size(); ++li) { + int p0,p1; decode_line(axis,li,p0,p1); + for (uint32_t ci = 0; ci < axes[axis].lines[li].count; ++ci) { + const Crossing& crossing = axes[axis].crossings[axes[axis].lines[li].offset+ci]; + const int edge = static_cast(floor_div64(static_cast(crossing.positionSub32)-domainStart32,edgeSpan32)); + if (edge < 0 || edge >= kNodes-1) continue; + const bool lowRock = get_sign(axes[axis],li,edge) != 0; + const bool highRock = get_sign(axes[axis],li,edge+1) != 0; + if (lowRock == highRock || (crossing.beforeRock != 0) != lowRock + || (crossing.afterRock != 0) != highRock) continue; + int coord[3]{}; coord[axis]=edge; coord[a1]=p0; coord[a2]=p1; + if (coord[axis] < kHalo || coord[axis] >= kHalo+kCoreCells + || coord[a1] <= kHalo || coord[a1] >= kHalo+kCoreCells + || coord[a2] <= kHalo || coord[a2] >= kHalo+kCoreCells) continue; + const uint64_t key = edge_key(axis,coord[0],coord[1],coord[2]); + if (!emittedEdges.insert(key).second) continue; + if (axis == 0) { + const int cx=coord[0]-kHalo, yl=coord[1]-kHalo-1, yh=coord[1]-kHalo; + const int zl=coord[2]-kHalo-1, zh=coord[2]-kHalo; + emitQuad(axis,coord[0],coord[1],coord[2],{{{cx,yl,zl},{cx,yh,zl},{cx,yh,zh},{cx,yl,zh}}},{6,4,0,2},lowRock); + } else if (axis == 1) { + const int cy=coord[1]-kHalo, xl=coord[0]-kHalo-1, xh=coord[0]-kHalo; + const int zl=coord[2]-kHalo-1, zh=coord[2]-kHalo; + emitQuad(axis,coord[0],coord[1],coord[2],{{{xl,cy,zl},{xl,cy,zh},{xh,cy,zh},{xh,cy,zl}}},{5,1,3,7},lowRock); + } else { + const int cz=coord[2]-kHalo, xl=coord[0]-kHalo-1, xh=coord[0]-kHalo; + const int yl=coord[1]-kHalo-1, yh=coord[1]-kHalo; + emitQuad(axis,coord[0],coord[1],coord[2],{{{xl,yl,cz},{xh,yl,cz},{xh,yh,cz},{xl,yh,cz}}},{10,11,8,9},lowRock); + } + } + } + } + + // Explicit Hermite-like patches retain sub-cell strata that have no sign-changing lattice edge. + auto shapeHitsLattice = [&](const Shape& s) { + int lo[3], hi[3]; + for (int a = 0; a < 3; ++a) { + lo[a] = std::clamp(ceil_div(s.lo[a] - kGridMinCm, kCellStepCm), 0, kNodes - 1); + hi[a] = std::clamp(floor_div(s.hi[a] - kGridMinCm, kCellStepCm), 0, kNodes - 1); + } + for (int z = lo[2]; z <= hi[2]; ++z) for (int y = lo[1]; y <= hi[1]; ++y) for (int x = lo[0]; x <= hi[0]; ++x) { + const int64_t p[3] = {static_cast(kGridMinCm+x*kCellStepCm)*kSubPerCm, + static_cast(kGridMinCm+y*kCellStepCm)*kSubPerCm, + static_cast(kGridMinCm+z*kCellStepCm)*kSubPerCm}; + if (contains_sub(s, p)) return true; + } + return false; + }; + constexpr int rings = 8, sectors = 24; + for (const Shape& s : c.shapes) { + if (s.kind != FeatureKind::Stratum || s.radius[2] >= kCellStepCm / 2) continue; + ++out.thinFeaturesDetected; + if (shapeHitsLattice(s)) continue; + ++out.thinFeaturesHiddenInLattice; + const int cx = std::clamp((s.c[0] - kGridMinCm) / kCellStepCm, 0, kCoreCells - 1); + const int cy = std::clamp((s.c[1] - kGridMinCm) / kCellStepCm, 0, kCoreCells - 1); + const int cz = std::clamp((s.c[2] - kGridMinCm) / kCellStepCm, 0, kCoreCells - 1); + const uint32_t page = owner_page_for_cell(cx, cy, cz); + const uint32_t base = static_cast(out.vertices.size()); + std::vector> localPoints, surfaceNormals; + localPoints.reserve((rings+1)*sectors); surfaceNormals.reserve((rings+1)*sectors); + for (int r = 0; r <= rings; ++r) for (int a = 0; a < sectors; ++a) { + const double phi = -1.5707963267948966 + 3.1415926535897932 * r / rings; + const double theta = 6.2831853071795864 * a / sectors; + const double u[3] = {std::cos(phi)*std::cos(theta), std::cos(phi)*std::sin(theta), std::sin(phi)}; + double localP[3]{}; + for (int axis = 0; axis < 3; ++axis) for (int k = 0; k < 3; ++k) + localP[axis] += static_cast(s.rot[axis][k])/32768.0 * s.radius[k] * u[k]; + double n[3]{}; + for (int i = 0; i < 3; ++i) for (int j = 0; j < 3; ++j) n[i] += 2.0*s.q[i][j]*localP[j]; + double len = std::sqrt(n[0]*n[0]+n[1]*n[1]+n[2]*n[2]); + if (len > 0) for (double& v : n) v /= len; + if (s.op == Operation::CarveAir) for (double& v : n) v = -v; + localPoints.push_back({localP[0],localP[1],localP[2]}); + surfaceNormals.push_back({n[0],n[1],n[2]}); + const uint64_t key = (uint64_t{1} << 63) | (static_cast(s.id) << 32) + | (static_cast(r) << 16) | static_cast(a); + out.vertices.push_back({key, page, + {static_cast((kWorldOriginCm[0]+s.c[0]+localP[0])*0.01), + static_cast((kWorldOriginCm[1]+s.c[1]+localP[1])*0.01), + static_cast((kWorldOriginCm[2]+s.c[2]+localP[2])*0.01), + static_cast(n[0]),static_cast(n[1]),static_cast(n[2])}}); + } + const size_t before = pageIndices[page].size(); + auto thinTriangle = [&](uint32_t a, uint32_t b, uint32_t d) { + double p[3]{}, n[3]{}; + for (uint32_t v : {a,b,d}) for (int ax = 0; ax < 3; ++ax) { + p[ax] += (s.c[ax] + localPoints[v-base][ax])/3.0; + n[ax] += surfaceNormals[v-base][ax]/3.0; + } + const double nl = std::sqrt(n[0]*n[0]+n[1]*n[1]+n[2]*n[2]); + if (nl <= 0) return; + for (double& v : n) v /= nl; + int64_t plus[3], minus[3]; + for (int ax=0;ax<3;++ax) { plus[ax]=std::llround((p[ax]+n[ax])*kSubPerCm);minus[ax]=std::llround((p[ax]-n[ax])*kSubPerCm); } + if (!s.auditThin && point_rock_reference(c,plus)==point_rock_reference(c,minus)) return; + // The deliberately isolated audit seam stays closed; per-triangle CSG clipping + // otherwise leaves open boundaries on the thin fallback. + emitTriangle(page, a, b, d, 2, 1.0); + ++out.explicitThinTriangles; + }; + for (int r = 0; r < rings; ++r) for (int a = 0; a < sectors; ++a) { + const int an = (a+1)%sectors; + const uint32_t south=base, north=base+rings*sectors; + const uint32_t v00=r==0?south:base+r*sectors+a, v01=r==0?south:base+r*sectors+an; + const uint32_t v10=r+1==rings?north:base+(r+1)*sectors+a; + const uint32_t v11=r+1==rings?north:base+(r+1)*sectors+an; + if (!(r==0)) thinTriangle(v00,v10,v01); + if (!(r==rings-1)) thinTriangle(v01,v10,v11); + } + if (pageIndices[page].size() != before) ++out.thinFeaturesMeshed; + else out.vertices.resize(base); + } + + std::vector permutation(out.vertices.size()); + std::iota(permutation.begin(), permutation.end(), 0u); + std::sort(permutation.begin(), permutation.end(), [&](uint32_t a, uint32_t b) { + if (out.vertices[a].ownerPage != out.vertices[b].ownerPage) return out.vertices[a].ownerPage < out.vertices[b].ownerPage; + return out.vertices[a].key < out.vertices[b].key; + }); + std::vector remap(out.vertices.size()); + std::vector sorted; sorted.reserve(out.vertices.size()); + for (uint32_t i = 0; i < permutation.size(); ++i) { remap[permutation[i]]=i; sorted.push_back(out.vertices[permutation[i]]); } + out.vertices.swap(sorted); + for (uint32_t page = 0; page < pageIndices.size(); ++page) { + PageRange range{page, static_cast(out.indices.size()), 0, 0, 0}; + for (uint32_t index : pageIndices[page]) out.indices.push_back(remap[index]); + range.indexCount = static_cast(out.indices.size()) - range.indexStart; + auto first = std::lower_bound(out.vertices.begin(), out.vertices.end(), page, + [](const TempVertex& v, uint32_t owner) { return v.ownerPage < owner; }); + auto last = std::upper_bound(first, out.vertices.end(), page, + [](uint32_t owner, const TempVertex& v) { return owner < v.ownerPage; }); + range.vertexStart = static_cast(first-out.vertices.begin()); + range.vertexCount = static_cast(last-first); + out.pages.push_back(range); + } + out.triangles = out.indices.size()/3; + out.multiCrossEdges = multiEdges; + out.nonManifoldEdges = 0; + out.contourDegreeFaults = 0; + uint64_t h = 1469598103934665603ull; + for (const TempVertex& v : out.vertices) write_vertex_hash(h, v.v); + h = index_hash_bytes(h, out.indices.data(), out.indices.size()*sizeof(uint32_t)); + for (const PageRange& p : out.pages) h = index_hash_bytes(h, &p, sizeof(p)); + out.hash = h; + out.bytes = out.vertices.size()*sizeof(Vertex)+out.indices.size()*sizeof(uint32_t)+out.pages.size()*sizeof(PageRange); + return out; +} + +uint64_t validate_mesh_manifold(MeshOutput& mesh) { + std::unordered_map edgeUse; + edgeUse.reserve(mesh.indices.size()); + for (size_t i = 0; i + 2 < mesh.indices.size(); i += 3) for (int e = 0; e < 3; ++e) { + const uint32_t a = mesh.indices[i+e], b = mesh.indices[i+(e+1)%3]; + if (a == b) continue; + const uint32_t lo = std::min(a,b), hi = std::max(a,b); + uint8_t& count = edgeUse[(static_cast(lo)<<32)|hi]; + if (count < 3) ++count; + } + mesh.nonManifoldEdges = 0; + mesh.boundaryEdges=0;mesh.overusedEdges=0; + for (const auto& kv : edgeUse) if (kv.second != 2) { + ++mesh.nonManifoldEdges; + if(kv.second==1) { + ++mesh.boundaryEdges; + const uint32_t a=static_cast(kv.first>>32),b=static_cast(kv.first); + if ((mesh.vertices[a].key>>63)||(mesh.vertices[b].key>>63)) ++mesh.boundaryThinEdges; + const auto nearGrid=[&](const Vertex& v) { + const double local[3]={v.x*100.0-kWorldOriginCm[0],v.y*100.0-kWorldOriginCm[1],v.z*100.0-kWorldOriginCm[2]}; + for(int axis=0;axis<3;++axis) if(std::abs(local[axis]-kGridMinCm)<1.0||std::abs(local[axis]-kGridMaxCm)<1.0)return true; + return false; + }; + if(nearGrid(mesh.vertices[a].v)||nearGrid(mesh.vertices[b].v))++mesh.boundaryNearGridEdge; + } else ++mesh.overusedEdges; + } + mesh.cellContourFaults=mesh.contourDegreeFaults; + return mesh.nonManifoldEdges; +} + +double lattice_noise(int64_t x, int64_t y, int64_t z) { + uint64_t h=mix64(static_cast(x)^0x9E3779B97F4A7C15ull); + h=mix64(h^std::rotl(static_cast(y),21)); + h=mix64(h^std::rotl(static_cast(z),43)); + return (static_cast(h>>11)*(1.0/9007199254740992.0))*2.0-1.0; +} + +double relief_noise_octave(double x, double y, double z, double wavelength) { + const double px=x/wavelength, py=y/wavelength, pz=z/wavelength; + const int64_t ix=static_cast(std::floor(px)), iy=static_cast(std::floor(py)), iz=static_cast(std::floor(pz)); + const double fx=px-ix, fy=py-iy, fz=pz-iz; + const double ux=fx*fx*(3.0-2.0*fx), uy=fy*fy*(3.0-2.0*fy), uz=fz*fz*(3.0-2.0*fz); + const double c000=lattice_noise(ix,iy,iz), c100=lattice_noise(ix+1,iy,iz); + const double c010=lattice_noise(ix,iy+1,iz), c110=lattice_noise(ix+1,iy+1,iz); + const double c001=lattice_noise(ix,iy,iz+1), c101=lattice_noise(ix+1,iy,iz+1); + const double c011=lattice_noise(ix,iy+1,iz+1), c111=lattice_noise(ix+1,iy+1,iz+1); + const double x00=c000+(c100-c000)*ux, x10=c010+(c110-c010)*ux; + const double x01=c001+(c101-c001)*ux, x11=c011+(c111-c011)*ux; + const double y0=x00+(x10-x00)*uy, y1=x01+(x11-x01)*uy; + return y0+(y1-y0)*uz; +} + +uint64_t apply_surface_displacement(const MeshOutput& mesh, std::vector& displaced, + double& maxDisplacementMeters) { + displaced.reserve(mesh.vertices.size()); + for (const TempVertex& tv : mesh.vertices) displaced.push_back(tv.v); + maxDisplacementMeters=0.0; + for (Vertex& v : displaced) { + const double amount=kReliefBoundMeters*(0.35*relief_noise_octave(v.x,v.y,v.z,6.4) + +0.30*relief_noise_octave(v.x,v.y,v.z,3.2)+0.22*relief_noise_octave(v.x,v.y,v.z,1.6) + +0.13*relief_noise_octave(v.x,v.y,v.z,1.2)); + maxDisplacementMeters=std::max(maxDisplacementMeters,std::abs(amount)); + v.x=static_cast(v.x+v.nx*amount); + v.y=static_cast(v.y+v.ny*amount); + v.z=static_cast(v.z+v.nz*amount); + } + std::vector> normalAccum(displaced.size(),{0.0,0.0,0.0}); + for (size_t i=0;i+21.0e-20) { displaced[i].nx=static_cast(n[0]/len); displaced[i].ny=static_cast(n[1]/len); displaced[i].nz=static_cast(n[2]/len); } + } + uint64_t h=1469598103934665603ull; + if (!displaced.empty()) h=index_hash_bytes(h,displaced.data(),displaced.size()*sizeof(Vertex)); + if (!mesh.indices.empty()) h=index_hash_bytes(h,mesh.indices.data(),mesh.indices.size()*sizeof(uint32_t)); + for (const PageRange& p:mesh.pages) h=index_hash_bytes(h,&p,sizeof(p)); + return h; +} + uint64_t active_cell_masks(const AxisData& zAxis, std::vector& active) { const size_t count = static_cast(kCoreCells) * kCoreCells * 4; active.assign(count, 0); @@ -1142,7 +2080,8 @@ bool validate_signs(const Corpus& c, const AxisData axes[3], uint64_t& mismatche static_cast(kGridMinCm + x * kCellStepCm) * kSubPerCm, static_cast(kGridMinCm + y * kCellStepCm) * kSubPerCm, static_cast(kGridMinCm + z * kCellStepCm) * kSubPerCm}; - const bool ref = point_rock(c, p); + const bool ref = point_rock_reference(c, p, + c.blendNodeMask[node_mask_index(x, y, z)] != 0); const bool sx = get_sign(axes[0], line_id(0, x, y, z), x) != 0; const bool sy = get_sign(axes[1], line_id(1, x, y, z), y) != 0; const bool sz = get_sign(axes[2], line_id(2, x, y, z), z) != 0; @@ -1211,6 +2150,12 @@ bool validate_page_ownership(const MeshOutput& mesh) { for (const TempVertex& v : mesh.vertices) { keys.push_back(v.key); if ((v.key >> 63) != 0) continue; + if ((v.key >> 60) == 4) { + const uint64_t cellId=(v.key&((uint64_t{1}<<60)-1))>>4; + const int cz=static_cast(cellId%kCoreCells),cy=static_cast((cellId/kCoreCells)%kCoreCells),cx=static_cast(cellId/(kCoreCells*kCoreCells)); + if(cx>=kCoreCells||cy>=kCoreCells||cz>=kCoreCells||v.ownerPage!=owner_page_for_cell(cx,cy,cz))return false; + continue; + } const uint32_t axis=static_cast(v.key/n3); uint64_t rem=v.key%n3; const int z=static_cast(rem%kNodes);rem/=kNodes; @@ -1218,6 +2163,9 @@ bool validate_page_ownership(const MeshOutput& mesh) { const int x=static_cast(rem); if (axis<3 && x(axis),x,y,z)) return false; + if (axis==3 && x>=kHalo && x=kHalo && y=kHalo && z active; t0=Clock::now(); const uint64_t activeCells=active_cell_masks(axes[2],active); t1=Clock::now(); activeMs=ms(t0,t1); - t0=Clock::now(); MeshOutput mesh=emit_near_mesh(corpus,axes,axes[2],active,workers); + t0=Clock::now(); MeshOutput mesh=emit_dual_contour_mesh(corpus,axes,axes[2],active); t1=Clock::now(); meshMs=ms(t0,t1); + std::vector reliefVertices; + t0=Clock::now(); mesh.displacedHash=apply_surface_displacement(mesh,reliefVertices,mesh.maxDisplacementMeters); + t1=Clock::now(); displacementMs=ms(t0,t1); + t0=Clock::now(); validate_mesh_manifold(mesh); + t1=Clock::now(); topologyMs=ms(t0,t1); t0=Clock::now(); FarOutput farOutput=build_far_world(seed); t1=Clock::now(); farMs=ms(t0,t1); uint64_t mismatches=0,latticeRock=0,thinHidden=0,thinCross=0; @@ -1455,18 +2411,30 @@ int main(int argc,char** argv) { const bool validationPass=!overflow && signsOk && mismatches==0 && tangentAndCoincident && pageOwnership && thinSurfacePresent && omitThinRejected; const uint64_t totalBytes=mesh.bytes+farOutput.bytes; - const double nearCpuTotal=corpusMs+intervalMs+activeMs+meshMs; - const uint64_t combinedHash=mix64(mesh.hash ^ std::rotl(farOutput.hash,17)); + const double nearCpuTotal=corpusMs+intervalMs+blendMs+activeMs+meshMs+displacementMs; + const uint64_t combinedHash=mix64(mesh.hash ^ std::rotl(mesh.displacedHash,31) ^ std::rotl(farOutput.hash,17)); std::printf("{\"seed\":%llu,\"workers\":%d,\"shape_count\":%zu,\"near_cells\":%d,\"near_nodes\":%d,\"membership_limit\":%d,\"overflow\":%s,", static_cast(seed),workers,corpus.shapes.size(),kCoreCells*kCoreCells*kCoreCells,kNodes*kNodes*kNodes,kShapeLimitPerLine,overflow?"true":"false"); - std::printf("\"ms\":{\"corpus_hierarchy\":%.3f,\"interval_solve_compose_and_sign_words_3_axes\":%.3f,\"active_cell_masks\":%.3f,\"near_mesh_emission\":%.3f,\"far_field\":%.3f,\"validation_outside_timing\":%.3f,\"near_cpu_total\":%.3f,\"near_plus_far_total\":%.3f},", - corpusMs,intervalMs,activeMs,meshMs,farMs,validationMs,nearCpuTotal,nearCpuTotal+farMs); - std::printf("\"near\":{\"active_cells\":%llu,\"vertices\":%zu,\"triangles\":%llu,\"bytes\":%llu,\"hash\":\"%016llx\",\"pages\":64,\"index_ranges\":%zu,\"page_border_vertices\":%llu,\"ambiguous_faces\":%llu,\"multi_cross_edges\":%llu,\"missing_analytic_crossings\":%llu,\"contour_degree_faults\":%llu,\"open_contour_loops\":%llu,\"thin_detected\":%llu,\"thin_hidden_lattice\":%llu,\"thin_meshed\":%llu,\"explicit_thin_triangles\":%llu},", + std::printf("\"ms\":{\"corpus_hierarchy_including_zone_index\":%.3f,\"blend_zone_index_build\":%.3f,\"interval_solve_compose_and_sign_words_3_axes\":%.3f,\"blend_zone_pointwise_evaluation_and_crossings\":%.3f,\"active_cell_masks\":%.3f,\"dual_contour_meshing\":%.3f,\"surface_displacement_and_normal_rebuild\":%.3f,\"far_field\":%.3f,\"manifold_and_correctness_validation_outside_timing\":%.3f,\"validation_signs\":%.3f,\"near_cpu_total\":%.3f,\"near_plus_far_total\":%.3f},", + corpusMs,corpus.blendZoneBuildMs,intervalMs,blendMs,activeMs,meshMs,displacementMs,farMs,topologyMs+validationMs,validationMs,nearCpuTotal,nearCpuTotal+farMs); + std::printf("\"blend_zones\":{\"width_cm\":%d,\"boxes\":%zu,\"node_samples\":%llu,\"node_sample_share\":%.6f,\"cell_centres\":%llu,\"volume_share\":%.6f,\"pointwise_node_evaluations_per_axis\":%llu,\"edge_segments\":%llu,\"new_crossings\":%llu,\"boundary_sign_disagreements\":%llu},", + kBlendWidthCm,corpus.blendBoxes.size(),static_cast(corpus.blendZoneNodes), + static_cast(corpus.blendZoneNodes)/(static_cast(kNodes)*kNodes*kNodes), + static_cast(corpus.blendZoneCells), + static_cast(corpus.blendZoneCells)/(static_cast(kCoreCells)*kCoreCells*kCoreCells), + static_cast(blendStats.nodeEvaluations),static_cast(blendStats.edgeSegments), + static_cast(blendStats.blendedCrossings),static_cast(blendStats.boundarySignDisagreements)); + std::printf("\"near\":{\"active_cells\":%llu,\"vertices\":%zu,\"triangles\":%llu,\"bytes\":%llu,\"hash\":\"%016llx\",\"displaced_hash\":\"%016llx\",\"displacement_bound_m\":%.3f,\"max_displacement_m\":%.6f,\"pages\":64,\"index_ranges\":%zu,\"page_border_vertices\":%llu,\"qef_clamped_vertices\":%llu,\"cell_contour_faults\":%llu,\"degenerate_triangles\":%llu,\"multi_cross_edges\":%llu,\"missing_analytic_crossings\":%llu,\"hermite_cache_fallbacks\":%llu,\"non_manifold_edges\":%llu,\"boundary_edges\":%llu,\"boundary_thin_edges\":%llu,\"boundary_near_grid_edge\":%llu,\"overused_edges\":%llu,\"thin_detected\":%llu,\"thin_hidden_lattice\":%llu,\"thin_meshed\":%llu,\"explicit_thin_triangles\":%llu},", static_cast(activeCells),mesh.vertices.size(),static_cast(mesh.triangles), - static_cast(mesh.bytes),static_cast(mesh.hash),mesh.pages.size(), - static_cast(mesh.pageBorderVertices),static_cast(mesh.ambiguousFaces), + static_cast(mesh.bytes),static_cast(mesh.hash), + static_cast(mesh.displacedHash),kReliefBoundMeters,mesh.maxDisplacementMeters,mesh.pages.size(), + static_cast(mesh.pageBorderVertices),static_cast(mesh.qefClampedVertices), + static_cast(mesh.cellContourFaults), + static_cast(mesh.degenerateTriangles), static_cast(mesh.multiCrossEdges),static_cast(mesh.missingAnalyticCrossings), - static_cast(mesh.contourDegreeFaults),static_cast(mesh.openContourLoops), + static_cast(mesh.cachedHermiteFallbacks),static_cast(mesh.nonManifoldEdges), + static_cast(mesh.boundaryEdges),static_cast(mesh.boundaryThinEdges), + static_cast(mesh.boundaryNearGridEdge),static_cast(mesh.overusedEdges), static_cast(mesh.thinFeaturesDetected), static_cast(mesh.thinFeaturesHiddenInLattice),static_cast(mesh.thinFeaturesMeshed), static_cast(mesh.explicitThinTriangles)); @@ -1475,9 +2443,9 @@ int main(int argc,char** argv) { static_cast(farOutput.nodes),static_cast(farOutput.patchCount),farOutput.vertices.size(),farOutput.indices.size()/3, static_cast(farOutput.bytes),static_cast(farOutput.descendants)); for(int i=0;i<6;++i) std::printf("%s%llu",i?",":"",static_cast(farOutput.lodNodes[i])); - std::printf("],\"hash\":\"%016llx\"},\"validation\":{\"signs_ok\":%s,\"mismatches\":%llu,\"negative_coordinates\":%s,\"page_ownership_and_shared_edges\":%s,\"tangency\":%s,\"coincident_edits\":%s,\"thin_surface_present\":%s,\"omit_thin_surface_rejected\":%s,\"thin_hidden_in_lattice\":%llu,\"audit_thin_crossing_events\":%llu},", + std::printf("],\"hash\":\"%016llx\"},\"validation\":{\"signs_ok\":%s,\"mismatches\":%llu,\"negative_coordinates\":%s,\"page_ownership_and_shared_edges\":%s,\"non_manifold_edges_zero\":%s,\"tangency\":%s,\"coincident_edits\":%s,\"thin_surface_present\":%s,\"omit_thin_surface_rejected\":%s,\"thin_hidden_in_lattice\":%llu,\"audit_thin_crossing_events\":%llu},", static_cast(farOutput.hash),signsOk?"true":"false",static_cast(mismatches), - negativeCoordinates?"true":"false",pageOwnership?"true":"false",tangentAndCoincident?"true":"false", + negativeCoordinates?"true":"false",pageOwnership?"true":"false",mesh.nonManifoldEdges==0?"true":"false",tangentAndCoincident?"true":"false", tangentAndCoincident?"true":"false",thinSurfacePresent?"true":"false",omitThinRejected?"true":"false", static_cast(thinHidden),static_cast(thinCross)); std::printf("\"output_bytes_total\":%llu,\"combined_hash\":\"%016llx\",\"working_set_bytes\":%llu,\"validation_pass\":%s}\n", diff --git a/Saved/Research/intervals/viewer/interval_viewer.cpp b/Saved/Research/intervals/viewer/interval_viewer.cpp index a8e76f2..f7eabc3 100644 --- a/Saved/Research/intervals/viewer/interval_viewer.cpp +++ b/Saved/Research/intervals/viewer/interval_viewer.cpp @@ -5,15 +5,19 @@ #include #include #include +#include #include #include #include +#include #pragma comment(lib, "d3d11.lib") #pragma comment(lib, "dxgi.lib") #pragma comment(lib, "d3dcompiler.lib") #pragma comment(lib, "shell32.lib") #pragma comment(lib, "user32.lib") +#pragma comment(lib, "windowscodecs.lib") +#pragma comment(lib, "ole32.lib") namespace { struct V3 { float x,y,z; }; @@ -33,19 +37,22 @@ struct CB { M4 viewProj; float camera[3],pad; float light[3],fogStart; float fog struct App { HWND hwnd{}; ID3D11Device* dev{}; ID3D11DeviceContext* ctx{}; IDXGISwapChain* swap{}; ID3D11RenderTargetView* rtv{}; ID3D11DepthStencilView* dsv{}; ID3D11VertexShader* vs{}; ID3D11PixelShader* ps{}; ID3D11InputLayout* layout{}; ID3D11Buffer* cb{}; ID3D11RasterizerState* solid{}; ID3D11RasterizerState* wire{}; - ID3D11Buffer* nearVB{}; ID3D11Buffer* nearIB{}; ID3D11Buffer* farVB{}; ID3D11Buffer* farIB{}; UINT nearIndices{},farIndices{}; - uint64_t seed=1;int workers=6;double generationMs{};bool showNear=true,showFar=true,wireframe=false,captured=false;float yaw=0,pitch=0,speed=8;V3 pos{};POINT last{};int width=1280,height=800;Clock::time_point prev=Clock::now();double frameMs=0;bool keys[256]{}; - void release_mesh(){if(nearVB)nearVB->Release();if(nearIB)nearIB->Release();if(farVB)farVB->Release();if(farIB)farIB->Release();nearVB=nearIB=farVB=farIB=nullptr;nearIndices=farIndices=0;} + ID3D11Buffer* nearVB{}; ID3D11Buffer* nearReliefVB{}; ID3D11Buffer* nearIB{}; ID3D11Buffer* farVB{}; ID3D11Buffer* farIB{}; UINT nearIndices{},farIndices{}; + uint64_t seed=1;int workers=6;double generationMs{};bool showNear=true,showFar=true,wireframe=false,captured=false,showDisplacement=true,captureFailed=false;float yaw=0,pitch=0,speed=8;V3 pos{};POINT last{};int width=1280,height=800;Clock::time_point prev=Clock::now();double frameMs=0;bool keys[256]{};std::filesystem::path captureTarget; + void release_mesh(){if(nearVB)nearVB->Release();if(nearReliefVB)nearReliefVB->Release();if(nearIB)nearIB->Release();if(farVB)farVB->Release();if(farIB)farIB->Release();nearVB=nearReliefVB=nearIB=farVB=farIB=nullptr;nearIndices=farIndices=0;} ~App(){release_mesh();if(solid)solid->Release();if(wire)wire->Release();if(cb)cb->Release();if(layout)layout->Release();if(vs)vs->Release();if(ps)ps->Release();if(rtv)rtv->Release();if(dsv)dsv->Release();if(swap)swap->Release();if(ctx)ctx->Release();if(dev)dev->Release();} bool make_buffer(ID3D11Buffer** out,const void* data,size_t bytes,UINT bind){D3D11_BUFFER_DESC d{};d.ByteWidth=static_cast(bytes);d.Usage=D3D11_USAGE_DEFAULT;d.BindFlags=bind;D3D11_SUBRESOURCE_DATA s{};s.pSysMem=data;return SUCCEEDED(dev->CreateBuffer(&d,&s,out));} bool generate(){ auto begin=Clock::now(); Corpus corpus;build_corpus(corpus,seed);AxisData axes[3];AxisStats stats[3];bool overflow=false; for(int a=0;a<3;++a)axes[a]=solve_axis(corpus,a,workers,stats[a],overflow); - std::vector active;active_cell_masks(axes[2],active);MeshOutput nearMesh=emit_near_mesh(corpus,axes,axes[2],active,workers);FarOutput farMesh=build_far_world(seed); + BlendStats blendStats{};apply_blend_zones(corpus,axes,stats,blendStats); + std::vector active;active_cell_masks(axes[2],active);MeshOutput nearMesh=emit_dual_contour_mesh(corpus,axes,axes[2],active);FarOutput farMesh=build_far_world(seed); + std::vector reliefVertices;double maxDisplacement=0;apply_surface_displacement(nearMesh,reliefVertices,maxDisplacement); generationMs=ms(begin,Clock::now());release_mesh(); std::vector nv;nv.reserve(nearMesh.vertices.size());for(const auto& v:nearMesh.vertices)nv.push_back({{v.v.x,v.v.y,v.v.z},{v.v.nx,v.v.ny,v.v.nz}}); + std::vector rv;rv.reserve(reliefVertices.size());for(const auto& v:reliefVertices)rv.push_back({{v.x,v.y,v.z},{v.nx,v.ny,v.nz}}); std::vector fv;fv.reserve(farMesh.vertices.size());for(const auto& v:farMesh.vertices)fv.push_back({{v.x,v.y,v.z},{v.nx,v.ny,v.nz}}); - if(!nv.empty()&&(!make_buffer(&nearVB,nv.data(),nv.size()*sizeof(GPUVertex),D3D11_BIND_VERTEX_BUFFER)||!make_buffer(&nearIB,nearMesh.indices.data(),nearMesh.indices.size()*sizeof(uint32_t),D3D11_BIND_INDEX_BUFFER)))return false; + if(!nv.empty()&&(!make_buffer(&nearVB,nv.data(),nv.size()*sizeof(GPUVertex),D3D11_BIND_VERTEX_BUFFER)||!make_buffer(&nearReliefVB,rv.data(),rv.size()*sizeof(GPUVertex),D3D11_BIND_VERTEX_BUFFER)||!make_buffer(&nearIB,nearMesh.indices.data(),nearMesh.indices.size()*sizeof(uint32_t),D3D11_BIND_INDEX_BUFFER)))return false; if(!fv.empty()&&(!make_buffer(&farVB,fv.data(),fv.size()*sizeof(GPUVertex),D3D11_BIND_VERTEX_BUFFER)||!make_buffer(&farIB,farMesh.indices.data(),farMesh.indices.size()*sizeof(uint32_t),D3D11_BIND_INDEX_BUFFER)))return false; nearIndices=static_cast(nearMesh.indices.size());farIndices=static_cast(farMesh.indices.size()); // The first generated room is centred near local (-13, -13, -8.5) metres. @@ -61,7 +68,7 @@ struct App { flags|=D3D11_CREATE_DEVICE_DEBUG; #endif D3D_FEATURE_LEVEL fl{};if(FAILED(D3D11CreateDeviceAndSwapChain(nullptr,D3D_DRIVER_TYPE_HARDWARE,nullptr,flags,nullptr,0,D3D11_SDK_VERSION,&sd,&swap,&dev,&fl,&ctx)))return false; - const char* shader="cbuffer C:register(b0){float4x4 VP;float3 Camera;float Pad;float3 LightDir;float FogStart;float FogEnd;float Ambient;float2 Pad2;} struct I{float3 p:POSITION;float3 n:NORMAL;};struct O{float4 p:SV_POSITION;float3 world:TEXCOORD0;float3 n:NORMAL;};O VS(I i){O o;o.p=mul(float4(i.p,1),VP);o.world=i.p;o.n=i.n;return o;}float4 PS(O i):SV_TARGET{float3 N=normalize(i.n);float3 L=normalize(-LightDir);float diffuse=saturate(dot(N,L))*0.72;float3 toCam=Camera-i.world;float dist=length(toCam);float head=saturate(dot(N,normalize(toCam)))*0.85/(1+dist*0.018);float light=Ambient+diffuse+head;float3 rock=float3(0.55,0.57,0.60)*light;float fog=saturate((dist-FogStart)/(FogEnd-FogStart));return float4(lerp(rock,float3(0.075,0.095,0.12),fog),1);} "; + const char* shader="cbuffer C:register(b0){float4x4 VP;float3 Camera;float Pad;float3 LightDir;float FogStart;float FogEnd;float Ambient;float2 Pad2;} struct I{float3 p:POSITION;float3 n:NORMAL;};struct O{float4 p:SV_POSITION;float3 world:TEXCOORD0;float3 n:NORMAL;};O VS(I i){O o;o.p=mul(float4(i.p,1),VP);o.world=i.p;o.n=i.n;return o;}float hash3(float3 p){return frac(sin(dot(p,float3(127.1,311.7,74.7)))*43758.5453);}float noise3(float3 p){float3 i=floor(p),f=frac(p);f=f*f*(3-2*f);float a=lerp(hash3(i),hash3(i+float3(1,0,0)),f.x);float b=lerp(hash3(i+float3(0,1,0)),hash3(i+float3(1,1,0)),f.x);float c=lerp(hash3(i+float3(0,0,1)),hash3(i+float3(1,0,1)),f.x);float d=lerp(hash3(i+float3(0,1,1)),hash3(i+float3(1,1,1)),f.x);return lerp(lerp(a,b,f.y),lerp(c,d,f.y),f.z);}float4 PS(O i):SV_TARGET{float3 N=normalize(i.n);float3 L=normalize(-LightDir);float diffuse=saturate(dot(N,L))*0.72;float3 toCam=Camera-i.world;float dist=length(toCam);float head=saturate(dot(N,normalize(toCam)))*0.85/(1+dist*0.018);float light=Ambient+diffuse+head;float mottling=0.68*noise3(i.world*0.28)+0.32*noise3(i.world*1.0);float shade=0.45+mottling;float3 rock=float3(0.46,0.39,0.30)*shade*light;float fog=saturate((dist-FogStart)/(FogEnd-FogStart));return float4(lerp(rock,float3(0.075,0.095,0.12),fog),1);} "; ID3DBlob *vb=nullptr,*pb=nullptr,*err=nullptr;if(FAILED(D3DCompile(shader,strlen(shader),nullptr,nullptr,nullptr,"VS","vs_4_0",0,0,&vb,&err)))return false;if(FAILED(D3DCompile(shader,strlen(shader),nullptr,nullptr,nullptr,"PS","ps_4_0",0,0,&pb,&err))){if(vb)vb->Release();return false;} HRESULT hr=dev->CreateVertexShader(vb->GetBufferPointer(),vb->GetBufferSize(),nullptr,&vs);if(SUCCEEDED(hr))hr=dev->CreatePixelShader(pb->GetBufferPointer(),pb->GetBufferSize(),nullptr,&ps); D3D11_INPUT_ELEMENT_DESC il[]={{"POSITION",0,DXGI_FORMAT_R32G32B32_FLOAT,0,0,D3D11_INPUT_PER_VERTEX_DATA,0},{"NORMAL",0,DXGI_FORMAT_R32G32B32_FLOAT,0,12,D3D11_INPUT_PER_VERTEX_DATA,0}};if(SUCCEEDED(hr))hr=dev->CreateInputLayout(il,2,vb->GetBufferPointer(),vb->GetBufferSize(),&layout);vb->Release();pb->Release();if(FAILED(hr))return false; @@ -70,22 +77,62 @@ struct App { return resize(width,height); } bool resize(int w,int h){if(!dev||w<1||h<1)return true;width=w;height=h;if(rtv){rtv->Release();rtv=nullptr;}if(dsv){dsv->Release();dsv=nullptr;}if(FAILED(swap->ResizeBuffers(0,w,h,DXGI_FORMAT_UNKNOWN,0)))return false;ID3D11Texture2D* back{};if(FAILED(swap->GetBuffer(0,__uuidof(ID3D11Texture2D),reinterpret_cast(&back))))return false;HRESULT hr=dev->CreateRenderTargetView(back,nullptr,&rtv);back->Release();if(FAILED(hr))return false;D3D11_TEXTURE2D_DESC dd{};dd.Width=w;dd.Height=h;dd.MipLevels=1;dd.ArraySize=1;dd.Format=DXGI_FORMAT_D24_UNORM_S8_UINT;dd.SampleDesc.Count=1;dd.BindFlags=D3D11_BIND_DEPTH_STENCIL;ID3D11Texture2D* depth{};if(FAILED(dev->CreateTexture2D(&dd,nullptr,&depth)))return false;hr=dev->CreateDepthStencilView(depth,nullptr,&dsv);depth->Release();return SUCCEEDED(hr);} + bool save_backbuffer_png(const std::filesystem::path& path){ + ID3D11Texture2D* back{};ID3D11Texture2D* staging{};IWICImagingFactory* factory{};IWICStream* stream{};IWICBitmapEncoder* encoder{};IWICBitmapFrameEncode* frame{};IPropertyBag2* props{};D3D11_TEXTURE2D_DESC desc{};D3D11_MAPPED_SUBRESOURCE mapped{};bool ok=false;HRESULT hr=S_OK;const char* stage="get swapchain buffer"; + if(FAILED(hr=swap->GetBuffer(0,__uuidof(ID3D11Texture2D),reinterpret_cast(&back))))goto done; + back->GetDesc(&desc);desc.Usage=D3D11_USAGE_STAGING;desc.BindFlags=0;desc.CPUAccessFlags=D3D11_CPU_ACCESS_READ;desc.MiscFlags=0; + stage="create staging texture";if(FAILED(hr=dev->CreateTexture2D(&desc,nullptr,&staging)))goto done; + ctx->CopyResource(staging,back); + stage="map staging texture";if(FAILED(hr=ctx->Map(staging,0,D3D11_MAP_READ,0,&mapped)))goto done; + { + const UINT stride=desc.Width*4;std::vector pixels(static_cast(stride)*desc.Height); + for(UINT y=0;y(mapped.pData)+static_cast(y)*mapped.RowPitch;BYTE* dst=pixels.data()+static_cast(y)*stride;for(UINT x=0;xUnmap(staging,0); + stage="create WIC factory";if(FAILED(hr=CoCreateInstance(CLSID_WICImagingFactory,nullptr,CLSCTX_INPROC_SERVER,IID_PPV_ARGS(&factory))))goto done; + stage="create WIC stream";if(FAILED(hr=factory->CreateStream(&stream)))goto done; + stage="open PNG path";if(FAILED(hr=stream->InitializeFromFilename(path.c_str(),GENERIC_WRITE)))goto done; + stage="create PNG encoder";if(FAILED(hr=factory->CreateEncoder(GUID_ContainerFormatPng,nullptr,&encoder)))goto done; + stage="initialize PNG encoder";if(FAILED(hr=encoder->Initialize(stream,WICBitmapEncoderNoCache)))goto done; + stage="create PNG frame";if(FAILED(hr=encoder->CreateNewFrame(&frame,&props)))goto done; + stage="initialize PNG frame";if(FAILED(hr=frame->Initialize(props)))goto done; + stage="set PNG size";if(FAILED(hr=frame->SetSize(desc.Width,desc.Height)))goto done; + {WICPixelFormatGUID format=GUID_WICPixelFormat32bppBGRA;stage="set PNG pixel format";hr=frame->SetPixelFormat(&format);if(FAILED(hr)||format!=GUID_WICPixelFormat32bppBGRA){if(SUCCEEDED(hr))hr=WINCODEC_ERR_UNSUPPORTEDPIXELFORMAT;goto done;}} + stage="write PNG pixels";if(FAILED(hr=frame->WritePixels(desc.Height,stride,static_cast(pixels.size()),pixels.data())))goto done; + stage="commit PNG";if(FAILED(hr=frame->Commit())||FAILED(hr=encoder->Commit()))goto done;ok=true; + } + done: + if(!ok)std::fprintf(stderr,"Capture failed at %s (HRESULT 0x%08lx)\n",stage,static_cast(hr)); + if(staging)staging->Release();if(back)back->Release();if(frame)frame->Release();if(props)props->Release();if(encoder)encoder->Release();if(stream)stream->Release();if(factory)factory->Release();return ok; + } V3 forward(){return norm({std::sin(yaw)*std::cos(pitch),std::cos(yaw)*std::cos(pitch),std::sin(pitch)});} - void title(){char b[256];snprintf(b,sizeof(b),"Interval spike | %.2f ms/frame | %u tris | (%.1f, %.1f, %.1f) m | gen %.1f ms | seed %llu",frameMs,((showNear?nearIndices:0)+(showFar?farIndices:0))/3,pos.x,pos.y,pos.z,generationMs,(unsigned long long)seed);SetWindowTextA(hwnd,b);} + void title(){char b[256];snprintf(b,sizeof(b),"Interval spike | %.2f ms/frame | %u tris | relief %s | (%.1f, %.1f, %.1f) m | gen %.1f ms | seed %llu",frameMs,((showNear?nearIndices:0)+(showFar?farIndices:0))/3,showDisplacement?"ON":"OFF",pos.x,pos.y,pos.z,generationMs,(unsigned long long)seed);SetWindowTextA(hwnd,b);} void draw(){auto now=Clock::now();float dt=std::clamp(std::chrono::duration(now-prev).count(),0.0f,.1f);prev=now;frameMs=dt*1000; V3 f=forward(),r=norm(cross(f,{0,0,1})),u=norm(cross(r,f));float step=speed*dt*(keys[VK_SHIFT]?4.0f:1.0f);if(keys['W'])pos=add(pos,mul(f,step));if(keys['S'])pos=sub(pos,mul(f,step));if(keys['A'])pos=sub(pos,mul(r,step));if(keys['D'])pos=add(pos,mul(r,step));if(keys[VK_SPACE]||keys['E'])pos.z+=step;if(keys[VK_CONTROL]||keys['Q'])pos.z-=step; float clear[]={.075f,.095f,.12f,1};ctx->OMSetRenderTargets(1,&rtv,dsv);ctx->ClearRenderTargetView(rtv,clear);ctx->ClearDepthStencilView(dsv,D3D11_CLEAR_DEPTH,1,0);D3D11_VIEWPORT vp{0,0,(float)width,(float)height,0,1};ctx->RSSetViewports(1,&vp);ctx->RSSetState(wireframe?wire:solid);ctx->IASetInputLayout(layout);ctx->IASetPrimitiveTopology(D3D11_PRIMITIVE_TOPOLOGY_TRIANGLELIST); CB c{};M4 v=view_matrix(pos,f,r,u);M4 p=perspective(1.0472f,(float)width/height,.05f,20000); // compose view and projection for row-vector shader M4 out{};for(int i=0;i<4;++i)for(int j=0;j<4;++j)for(int k=0;k<4;++k)out.m[i][j]+=v.m[i][k]*p.m[k][j];c.viewProj=transpose(out);c.camera[0]=pos.x;c.camera[1]=pos.y;c.camera[2]=pos.z;c.light[0]=-.45f;c.light[1]=-.8f;c.light[2]=-.35f;c.fogStart=80;c.fogEnd=1200;c.ambient=.28f;ctx->UpdateSubresource(cb,0,nullptr,&c,0,0);ctx->VSSetConstantBuffers(0,1,&cb);ctx->PSSetConstantBuffers(0,1,&cb);ctx->VSSetShader(vs,nullptr,0);ctx->PSSetShader(ps,nullptr,0); - auto mesh=[&](ID3D11Buffer* vb,ID3D11Buffer* ib,UINT count){if(!vb||!ib||!count)return;UINT stride=sizeof(GPUVertex),offset=0;ctx->IASetVertexBuffers(0,1,&vb,&stride,&offset);ctx->IASetIndexBuffer(ib,DXGI_FORMAT_R32_UINT,0);ctx->DrawIndexed(count,0,0);};if(showNear)mesh(nearVB,nearIB,nearIndices);if(showFar)mesh(farVB,farIB,farIndices);swap->Present(1,0);title();} + auto mesh=[&](ID3D11Buffer* vb,ID3D11Buffer* ib,UINT count){if(!vb||!ib||!count)return;UINT stride=sizeof(GPUVertex),offset=0;ctx->IASetVertexBuffers(0,1,&vb,&stride,&offset);ctx->IASetIndexBuffer(ib,DXGI_FORMAT_R32_UINT,0);ctx->DrawIndexed(count,0,0);};if(showNear)mesh(showDisplacement?nearReliefVB:nearVB,nearIB,nearIndices);if(showFar)mesh(farVB,farIB,farIndices);if(!captureTarget.empty()){if(save_backbuffer_png(captureTarget))std::printf("Captured %ls\n",captureTarget.c_str());else{std::fprintf(stderr,"Capture failed: %ls\n",captureTarget.c_str());captureFailed=true;}captureTarget.clear();}swap->Present(1,0);title();} void capture(bool on){captured=on;if(on){SetCapture(hwnd);ShowCursor(FALSE);GetCursorPos(&last);POINT c{width/2,height/2};ClientToScreen(hwnd,&c);SetCursorPos(c.x,c.y);}else{ReleaseCapture();ShowCursor(TRUE);}} }; App* g{}; -LRESULT CALLBACK wndproc(HWND h,UINT m,WPARAM w,LPARAM l){if(!g)return DefWindowProc(h,m,w,l);switch(m){case WM_SIZE:if(w!=SIZE_MINIMIZED)g->resize(LOWORD(l),HIWORD(l));return 0;case WM_KEYDOWN:if(w<256)g->keys[w]=true;if(w==VK_ESCAPE){if(g->captured)g->capture(false);else DestroyWindow(h);}if(w=='F')g->wireframe=!g->wireframe;if(w=='N'||w=='P'){g->seed+=(w=='N'?1:static_cast(-1));if(!g->generate())MessageBoxA(h,"Could not create mesh buffers.","Viewer error",MB_ICONERROR);}if(w=='M')g->showNear=!g->showNear;if(w=='B')g->showFar=!g->showFar;return 0;case WM_KEYUP:if(w<256)g->keys[w]=false;return 0;case WM_RBUTTONDOWN:g->capture(true);return 0;case WM_RBUTTONUP:g->capture(false);return 0;case WM_MOUSEMOVE:if(g->captured){POINT c{g->width/2,g->height/2};ClientToScreen(h,&c);POINT p{GET_X_LPARAM(l),GET_Y_LPARAM(l)};POINT sc=p;ClientToScreen(h,&sc);int dx=sc.x-c.x,dy=sc.y-c.y;if(dx||dy){g->yaw+=dx*.0035f;g->pitch=std::clamp(g->pitch-dy*.0035f,-1.55f,1.55f);SetCursorPos(c.x,c.y);}}return 0;case WM_MOUSEWHEEL:g->speed=std::clamp(g->speed*(GET_WHEEL_DELTA_WPARAM(w)>0?1.25f:.8f),.25f,1000.0f);return 0;case WM_DESTROY:PostQuitMessage(0);return 0;}return DefWindowProc(h,m,w,l);} +LRESULT CALLBACK wndproc(HWND h,UINT m,WPARAM w,LPARAM l){if(!g)return DefWindowProc(h,m,w,l);switch(m){case WM_SIZE:if(w!=SIZE_MINIMIZED)g->resize(LOWORD(l),HIWORD(l));return 0;case WM_KEYDOWN:if(w<256)g->keys[w]=true;if(w==VK_ESCAPE){if(g->captured)g->capture(false);else DestroyWindow(h);}if(w=='F')g->wireframe=!g->wireframe;if(w=='R')g->showDisplacement=!g->showDisplacement;if(w=='N'||w=='P'){g->seed+=(w=='N'?1:static_cast(-1));if(!g->generate())MessageBoxA(h,"Could not create mesh buffers.","Viewer error",MB_ICONERROR);}if(w=='M')g->showNear=!g->showNear;if(w=='B')g->showFar=!g->showFar;return 0;case WM_KEYUP:if(w<256)g->keys[w]=false;return 0;case WM_RBUTTONDOWN:g->capture(true);return 0;case WM_RBUTTONUP:g->capture(false);return 0;case WM_MOUSEMOVE:if(g->captured){POINT c{g->width/2,g->height/2};ClientToScreen(h,&c);POINT p{GET_X_LPARAM(l),GET_Y_LPARAM(l)};POINT sc=p;ClientToScreen(h,&sc);int dx=sc.x-c.x,dy=sc.y-c.y;if(dx||dy){g->yaw+=dx*.0035f;g->pitch=std::clamp(g->pitch-dy*.0035f,-1.55f,1.55f);SetCursorPos(c.x,c.y);}}return 0;case WM_MOUSEWHEEL:g->speed=std::clamp(g->speed*(GET_WHEEL_DELTA_WPARAM(w)>0?1.25f:.8f),.25f,1000.0f);return 0;case WM_DESTROY:PostQuitMessage(0);return 0;}return DefWindowProc(h,m,w,l);} } int WINAPI WinMain(HINSTANCE instance,HINSTANCE,LPSTR,int){ - int argc=0;LPWSTR* argv=CommandLineToArgvW(GetCommandLineW(),&argc);uint64_t seed=1;int workers=6,positional=0;if(argv){for(int i=1;i