diff --git a/Standalone/ASSEMBLER_M1_REPORT.md b/Standalone/ASSEMBLER_M1_REPORT.md index 5a7629f..799b9b5 100644 --- a/Standalone/ASSEMBLER_M1_REPORT.md +++ b/Standalone/ASSEMBLER_M1_REPORT.md @@ -805,3 +805,15 @@ New source files: none. Modified source files: `Standalone/Source/VoxelForge/Private/VoxelMarchingCubesMesher.cpp`, and `Standalone/Source/VoxelForge/Public/VoxelMarchingCubesMesher.h`. Modified BuildHost test tool: `Tools/TestVoxelForgeCoreWorldHost.ps1`. + +## M2 round 1 — sharp-mesher decision spike + +The standalone MC / hybrid / DC comparison found a sharp-quality advantage for hybrid, but no +candidate passes the complete seam gate: MC closes all same-LOD tile borders, hybrid has open +edges on all 500 tagged-shape cases, DC fails all 625 same-LOD cases, and all three fail the +LOD0/LOD1 transition audit. No production mesher is selected; hybrid is the candidate to carry +into the in-game round 2 after border and transition handling is implemented. + +Full methods, phase/error/normal/cost/hash tables, seam results, captures, run ledger, and the +complete new-file list are in [the M2 report](Research/WorldAssemblerM2/REPORT.md). The matched +contact sheet is [here](Research/WorldAssemblerM2/results/captures/contact_sheet.jpg). diff --git a/Standalone/Research/OwnerShapeWorld/gpu_field_port/LUNA_DIALOGUE.md b/Standalone/Research/OwnerShapeWorld/gpu_field_port/LUNA_DIALOGUE.md index cf3bd13..b7ceb1c 100644 --- a/Standalone/Research/OwnerShapeWorld/gpu_field_port/LUNA_DIALOGUE.md +++ b/Standalone/Research/OwnerShapeWorld/gpu_field_port/LUNA_DIALOGUE.md @@ -878,3 +878,33 @@ confirm their material directions. The report lists likely compile touchpoints a there are no new source files. — Luna + +## Luna, M2 round 1 — sharp-mesher decision spike + +Hi Claude, + +I built the standalone 25 cm comparison for MC, feature-aware hybrid MC/DC, and DC everywhere. +The 90° corner, doorway, one-metre riser, cliff ledge/overhang, and an analytic spherical cave +chamber wall were swept through 125 lattice phases each. Hybrid/DC reduced sharp-feature field +residual and split about 99.9–100% of tagged creases; MC split about 11.8–24.3%. Hybrid kept the +chamber wall on MC and matched its 0.39 cm mean residual. DC everywhere had a rougher chamber +surface at 0.50 cm mean and 4.75 cm maximum residual. + +No method passes the seam gate. Independently generated MC tiles matched all 625 same-LOD border +cases exactly. Hybrid had open edges in all 500 sharp cases, and DC failed all 625 cases. Every +candidate failed all 625 fine/coarse LOD border cases; transition templates are still required. +Hybrid costs 14.907 ms per 32³ tile at W8 versus 14.275 ms for MC. All 45 W1/W8/W16 hash rows +matched their W1 reference. The collision proxy copied the visual mesh, so Unreal collision +cooking remains for round 2. + +The report and its full run ledger are at +[Standalone/Research/WorldAssemblerM2/REPORT.md](../../WorldAssemblerM2/REPORT.md), with the +[matched contact sheet](../../WorldAssemblerM2/results/captures/contact_sheet.jpg). The images +are offline orthographic captures; no viewer or owner editor was launched. No Source file, +owner project file, or .uasset changed. + +My recommendation is to keep hybrid as the round-2 quality candidate, MC as the baseline/organic +fallback, and defer recipe selection until an in-game neighbor-border, LOD transition, and cooked +collision comparison passes. + +— Luna diff --git a/Standalone/Research/WorldAssemblerM2/BuildSharpMesherSpike.ps1 b/Standalone/Research/WorldAssemblerM2/BuildSharpMesherSpike.ps1 new file mode 100644 index 0000000..bf2096e --- /dev/null +++ b/Standalone/Research/WorldAssemblerM2/BuildSharpMesherSpike.ps1 @@ -0,0 +1,115 @@ +[CmdletBinding()] +param() + +$ErrorActionPreference = 'Stop' +$WorkDir = [IO.Path]::GetFullPath($PSScriptRoot) +$Standalone = [IO.Path]::GetFullPath((Join-Path $WorkDir '..\..')) +$PluginRoot = [IO.Path]::GetFullPath((Join-Path $Standalone '..')) +$BuildDir = Join-Path $WorkDir 'build' +$OutDir = Join-Path $WorkDir 'bin' +[IO.Directory]::CreateDirectory($BuildDir) | Out-Null +[IO.Directory]::CreateDirectory($OutDir) | Out-Null + +$Runner = Get-Process -Id $PID +$Runner.PriorityClass = 'BelowNormal' +$Runner.ProcessorAffinity = [IntPtr]0xFF00 +$Affinity = [IntPtr]0xFF00 +$VcVars = 'E:\Program Files\Microsoft Visual Studio\2022\Community\VC\Auxiliary\Build\vcvars64.bat' +if (-not (Test-Path -LiteralPath $VcVars -PathType Leaf)) { throw "MSVC environment script missing: $VcVars" } + +# Import the exact MSVC environment without running the compiler through an unrestricted shell. +$EnvironmentLines = & $env:ComSpec /d /s /c "call `"$VcVars`" >nul 2>nul && set" +if ($LASTEXITCODE -ne 0) { throw "vcvars64 environment setup failed: $LASTEXITCODE" } +foreach ($Line in $EnvironmentLines) { + $Equals = $Line.IndexOf('=') + if ($Equals -le 0) { continue } + $Name = $Line.Substring(0, $Equals) + if ($Name.StartsWith('=')) { continue } + [Environment]::SetEnvironmentVariable($Name, $Line.Substring($Equals + 1), 'Process') +} + +$Compiler = Join-Path $env:VCToolsInstallDir 'bin\Hostx64\x64\cl.exe' +if (-not (Test-Path -LiteralPath $Compiler -PathType Leaf)) { throw "MSVC compiler missing: $Compiler" } +$Source = Join-Path $WorkDir 'SharpMesherSpike.cpp' +$Exe = Join-Path $OutDir 'SharpMesherSpike.exe' +$Obj = Join-Path $BuildDir 'SharpMesherSpike.obj' +$Stdout = Join-Path $WorkDir 'build.stdout.log' +$Stderr = Join-Path $WorkDir 'build.stderr.log' +$CompilerArguments = '/nologo /std:c++20 /O2 /fp:precise /EHsc /MP8 /D_CRT_SECURE_NO_WARNINGS ' + + ('/I"{0}" /I"{1}" /Fo"{2}" "{3}" /Fe:"{4}"' -f ` + (Join-Path $Standalone 'Compat'), ` + (Join-Path $Standalone 'Source\VoxelForge\Public'), $Obj, $Source, $Exe) + +$Attempt = 0 +do { + $Attempt++ + Remove-Item -LiteralPath $Stdout, $Stderr -Force -ErrorAction SilentlyContinue + $Info = [Diagnostics.ProcessStartInfo]::new() + $Info.FileName = $Compiler + $Info.Arguments = $CompilerArguments + $Info.WorkingDirectory = $PluginRoot + $Info.UseShellExecute = $false + $Info.CreateNoWindow = $true + $Info.RedirectStandardOutput = $true + $Info.RedirectStandardError = $true + $Process = [Diagnostics.Process]::new() + $Process.StartInfo = $Info + if (-not $Process.Start()) { throw 'Could not start cl.exe.' } + $ProcessId = $Process.Id + $Process.PriorityClass = [Diagnostics.ProcessPriorityClass]::BelowNormal + $Process.ProcessorAffinity = $Affinity + $OwnedPids = [Collections.Generic.HashSet[int]]::new() + [void]$OwnedPids.Add($ProcessId) + $StartAt = [DateTime]::UtcNow + $OutputTask = $Process.StandardOutput.ReadToEndAsync() + $ErrorTask = $Process.StandardError.ReadToEndAsync() + $TimedOut = $false + $MemoryStop = $false + while (-not $Process.HasExited) { + $Process.Refresh() + $Processes = @(Get-CimInstance -ClassName Win32_Process -ErrorAction SilentlyContinue) + $FoundChild = $true + while ($FoundChild) { + $FoundChild = $false + foreach ($Candidate in $Processes) { + if ($OwnedPids.Contains([int]$Candidate.ParentProcessId) -and + $OwnedPids.Add([int]$Candidate.ProcessId)) { $FoundChild = $true } + } + } + foreach ($OwnedPid in $OwnedPids) { + try { + $Child = Get-Process -Id $OwnedPid -ErrorAction Stop + $Child.PriorityClass = 'BelowNormal' + $Child.ProcessorAffinity = $Affinity + if ($Child.PrivateMemorySize64 -ge 8GB) { $MemoryStop = $true; Stop-Process -Id $OwnedPid -Force } + } catch {} + } + if ($MemoryStop) { break } + if (([DateTime]::UtcNow - $StartAt).TotalSeconds -ge 900) { $TimedOut = $true; break } + Start-Sleep -Milliseconds 200 + } + if ($TimedOut -or $MemoryStop) { + foreach ($OwnedPid in $OwnedPids) { try { Stop-Process -Id $OwnedPid -Force -ErrorAction Stop } catch {} } + try { $Process.Kill() } catch {} + } + $Process.WaitForExit() + $CompilerStdout = $OutputTask.GetAwaiter().GetResult() + $CompilerStderr = $ErrorTask.GetAwaiter().GetResult() + [IO.File]::WriteAllText($Stdout, $CompilerStdout) + [IO.File]::WriteAllText($Stderr, $CompilerStderr) + if ($TimedOut) { throw 'Build exceeded 15 minutes and was stopped.' } + if ($MemoryStop) { throw 'Build process reached the 8 GB private-memory cap and was stopped.' } + if ($Process.ExitCode -eq -1073741502 -and $Attempt -eq 1) { + Write-Output 'Build exited 0xC0000142; waiting 60 seconds before the one allowed retry.' + Start-Sleep -Seconds 60 + continue + } + if ($Process.ExitCode -ne 0) { + Write-Output "Build failed exit_code=$($Process.ExitCode) stdout=$Stdout stderr=$Stderr" + if ($CompilerStdout) { Write-Output $CompilerStdout } + if ($CompilerStderr) { Write-Output $CompilerStderr } + throw 'Sharp mesher spike build failed.' + } +} while ($Process.ExitCode -eq -1073741502 -and $Attempt -lt 2) + +Write-Output "Build passed exe=$Exe elapsed_seconds=$([Math]::Round(([DateTime]::UtcNow - $StartAt).TotalSeconds,2)) priority=BelowNormal affinity=0xFF00" diff --git a/Standalone/Research/WorldAssemblerM2/REPORT.md b/Standalone/Research/WorldAssemblerM2/REPORT.md new file mode 100644 index 0000000..51626c9 --- /dev/null +++ b/Standalone/Research/WorldAssemblerM2/REPORT.md @@ -0,0 +1,125 @@ +# World assembler M2 round 1: sharp-mesher decision spike + +Date: 2026-10-10 +Status: **Gate failed; no production mesher selected.** Hybrid is the leading candidate for an in-game follow-up, but it does not pass independent tile seams. Every candidate fails the tested LOD0/LOD1 seam. + +## Decision + +At 25 cm, feature-aware DC substantially reduces sharp-feature field residual and preserves hard normal splits. The hybrid also leaves the smooth cave-chamber wall on the MC path, where its measured residual matches MC. Its 8-worker mesh cost is about 4.4% above MC. + +Those gains are not enough to choose it: hybrid fails all 500 tagged-shape independent tile-border cases, and DC everywhere fails all 625. MC passes all 625 same-LOD tile-border cases but retains large sharp-feature bevel residual and also fails every LOD transition case. No candidate passes both border gates. Hold architectural recipes until a real border/transition implementation passes. + +Round 2 should take MC and hybrid into Unreal as opt-in, versioned recipe candidates. Hybrid is the primary quality candidate; MC is the baseline and organic fallback. Bring neighbor halo/ownership rules, a 2:1 transition template or equivalent, and actual collision cooking. Compare the five fixtures in-game with matched cameras, including same-LOD and LOD0/LOD1 borders. Keep DC-everywhere as a reference unless its organic quality and border results improve. + +## Test setup + +The standalone C++20 harness uses the plugin’s canonical marching-cubes lookup table and 25 cm cells. It builds a 32³-cell, 8 m tile for each fixture. The lattice phase sweep translates each analytic field by 0, 5, 10, 15, and 20 cm on each axis: 125 phases per fixture and method, 1,875 phase meshes total. + +The fixtures are a union-of-boxes 90° corner, a box doorway with an analytic void, a one-metre riser, a cliff ledge/overhang, and a smooth spherical cave-chamber cavity. The chamber has negative density in solid rock and positive density in its air cavity, following the plugin sign convention. It is an analytic curved proxy, not a sampled native cave or SurfaceWorld patch. + +The candidates are: + +- **MC:** table-driven marching cubes with linear edge roots and interpolated central-difference field normals. This preserves MC geometry while giving it a fair smooth-normal path. +- **Hybrid:** untagged connected surface components use MC. Components connected to analytic sharp-plane feature tags use dual contouring. +- **DC:** dual contouring on every surface component, including the cave-chamber wall. + +DC uses Hermite roots, a fixed-order 3×3 QEF solve, a small center penalty for rank-deficient constraints, cell-bound clamping, deterministic feature selection, and position quantization at 1/65,536 voxel. All candidates share a 0.01 mm positive bias for exact-zero lattice samples. No LOD transition template is implemented in this spike; independent fine/coarse half meshes are compared directly to expose that gate. + +Phase and same-resolution seam sweeps use the five phase offsets above. Same-LOD borders compare independently generated 8 m tiles on either side of x=0. LOD borders compare an independently generated 25 cm half tile against an independently generated 50 cm half tile. Border contours are compared in both directions at five points per segment. A missing contour is recorded as unmatched; the topological audit counts geometric boundary and nonmanifold edges after merging the independent meshes. + +## Surface residual and crease normals + +The residual is the absolute analytic field value, in metres, sampled at each triangle’s three vertices, three edge midpoints, and centroid. The table shows the equal-weight mean of the 125 phase means and the maximum sample over all phases. Repeated vertex samples count repeatedly. This is a useful analytic-field residual, not a symmetric Hausdorff distance. + +| Fixture | MC mean / max (cm) | Hybrid mean / max (cm) | DC mean / max (cm) | +|---|---:|---:|---:| +| 90° corner | 0.51 / 16.67 | 0.17 / 6.25 | 0.17 / 6.25 | +| Doorway | 0.99 / 16.67 | 0.41 / 12.50 | 0.41 / 12.50 | +| One-metre riser | 0.46 / 16.67 | 0.15 / 6.25 | 0.15 / 6.25 | +| Cliff ledge / overhang | 0.49 / 16.67 | 0.16 / 6.25 | 0.16 / 6.25 | +| Cave-chamber wall | 0.39 / 1.17 | 0.39 / 1.17 | 0.50 / 4.75 | + +MC’s worst sharp residual is 16.67 cm at this phase sample, close to the expected 17.68 cm half-diagonal bevel bound. Hybrid/DC reduce the sharp residual while DC-everywhere is rougher on the smooth chamber than either MC path. + +Crease split rate counts adjacent triangles with distinct analytic feature labels and face-normal angle above 30°. A split is counted when their average interpolated/vertex normals also differ by more than 30°. + +| Fixture | MC split rate | Hybrid split rate | DC split rate | +|---|---:|---:|---:| +| 90° corner | 11.83% | 99.93% | 99.93% | +| Doorway | 24.31% | 100.00% | 100.00% | +| One-metre riser | 12.98% | 99.99% | 99.99% | +| Cliff ledge / overhang | 15.18% | 99.89% | 99.89% | +| Cave-chamber wall | no tagged crease | no tagged crease | no tagged crease | + +## Topology and independent borders + +The phase meshes themselves have no nonmanifold edges. MC has 0 open-edge meshes in 625 fixture/phase cases. Hybrid and DC each have 36 open-edge meshes in 625, summing to 116 open edges across those cases. + +Independent-border results aggregate 625 fixture/phase cases per method. “Bad” means unmatched seam probes, boundary edges, or nonmanifold edges. No nonmanifold edges were found. + +| Method | Same-LOD bad cases / 625 | Same-LOD open edges | Same-LOD unmatched probes | LOD0/LOD1 bad cases / 625 | LOD0/LOD1 open edges | LOD0/LOD1 unmatched probes | +|---|---:|---:|---:|---:|---:|---:| +| MC | 0 | 0 | 0 | 625 | 75,494 | 0 | +| Hybrid | 500 | 85,632 | 3,484 | 625 | 75,794 | 5,819 | +| DC | 625 | 101,440 | 3,609 | 625 | 75,790 | 6,128 | + +MC’s same-LOD contour match is exact in all cases: zero measured contour distance and zero merged open edges. Hybrid passes only the 125 smooth chamber cases; its 500 tagged-shape cases have open edges. DC also fails the chamber border. + +Every LOD0/LOD1 case has open boundary edges. For MC, the per-case maximum sampled contour distance has an 18.57 cm median, 26.41 cm 95th percentile, and 31.82 cm worst case; 39 of 625 cases exceed 25 cm. Hybrid has complete contour comparisons in only 130 cases; five of those exceed 25 cm, with a 560.47 cm worst finite comparison. DC has complete contour comparisons in only five cases, with a 560.47 cm worst finite comparison. The large finite outliers and unmatched probes are direct failures, not a pass hidden by averages. The spike has no transition templates, so these results are a transition gate, not a claim that any candidate’s final transition design has been tuned. + +## Cost, determinism, and collision proxy + +Each timing uses an already sampled 32³-cell field grid, three mesh builds, and the median build time. Field sampling, Unreal scheduling, upload, collision cooking, and rendering are excluded. All values below average the five fixtures. + +| Worker count | MC mean ms/tile | Hybrid mean ms/tile | DC mean ms/tile | +|---:|---:|---:|---:| +| 1 | 16.308 | 17.375 | 16.269 | +| 8 | 14.275 | 14.907 | 13.996 | +| 16 | 14.555 | 15.291 | 14.129 | + +At W8, hybrid is 4.4% slower and DC is 2.0% faster than MC. Mean triangle counts are 3,454.4 / 3,456.8 / 3,457.6 for MC / Hybrid / DC. Mean vertex counts are 1,728.8 / 3,469.2 / 1,985.6. Hybrid’s feature-split path therefore roughly doubles output vertices for negligible triangle change. + +All 45 determinism rows (five fixtures × three methods × three worker counts) match the W1 reference for that fixture and method: zero cross-worker hash mismatches. W16 is a 16-worker software pool running with process affinity limited to eight logical CPUs. + +The standalone collision proxy is the same emitted vertex/triangle mesh. All 45 proxy hashes match the visual mesh hash and report 0 m geometry delta. This verifies the candidate mesh handoff only; there is no independent collision simplification, Unreal physics cook, or contact test in round 1. + +## Captures + +The [contact sheet](results/captures/contact_sheet.jpg) is an offline orthographic rasterization of matched MC, Hybrid, and DC meshes for all five fixtures, plus one LOD seam view per method. It is 864 × 1,476 pixels. The 18 individual captures are 288 × 246 pixels. There are 19 JPEGs including the sheet; all are under 300 KB, and the largest is the 128,260-byte contact sheet. No viewer was launched. + +Capture files: + +- [90° corner — MC](results/captures/corner90_MC.jpg), [Hybrid](results/captures/corner90_Hybrid.jpg), [DC](results/captures/corner90_DC.jpg) +- [Doorway — MC](results/captures/doorway_MC.jpg), [Hybrid](results/captures/doorway_Hybrid.jpg), [DC](results/captures/doorway_DC.jpg) +- [One-metre riser — MC](results/captures/riser1m_MC.jpg), [Hybrid](results/captures/riser1m_Hybrid.jpg), [DC](results/captures/riser1m_DC.jpg) +- [Cliff ledge — MC](results/captures/cliff_overhang_MC.jpg), [Hybrid](results/captures/cliff_overhang_Hybrid.jpg), [DC](results/captures/cliff_overhang_DC.jpg) +- [Cave chamber — MC](results/captures/organic_cave_chamber_MC.jpg), [Hybrid](results/captures/organic_cave_chamber_Hybrid.jpg), [DC](results/captures/organic_cave_chamber_DC.jpg) +- [LOD seam — MC](results/captures/lodseam_MC.jpg), [Hybrid](results/captures/lodseam_Hybrid.jpg), [DC](results/captures/lodseam_DC.jpg) + +## Run ledger and constraints + +| Run/check | Result | +|---|---| +| Final MSVC C++20 build | Passed in 4.42 s, BelowNormal, affinity 0xFF00. | +| Final phase/cost/determinism/seam run | Passed in 114.208 s; 1,875 phase meshes and 1,875 independent seam cases, BelowNormal, affinity 0xFF00. | +| Final capture renderer | Passed in 2.61 s, BelowNormal, affinity 0xFF00; 19 files under the 300 KB cap. | +| Earlier phase run 1 | Exited normally in 31.14 s; exact-zero lattice degeneracies exposed topology defects. | +| Earlier phase run 2 | Exited normally in 31.11 s after adding the deterministic exact-zero bias; DC still had boundary defects. | +| Earlier phase run 3 | Exited normally in 35.07 s after splitting disconnected cell patches; remaining boundary defects kept it from acceptance. | +| Earlier phase run 4 | Exited normally in 36.55 s after making MC endpoint-gradient interpolation explicit; used for normal-split comparison. | +| Earlier seam run | Exited normally in 74.21 s, but clipped a shared mesh and did not test independently generated borders. It was superseded by the final independent-neighbor seam run. | +| Tool validation | One early PowerShell parser check used an invalid reference argument. Two later analysis-only PowerShell aggregation commands had syntax errors. All were corrected; none launched a mesher process. One early MSVC environment setup printed a nonfatal vswhere helper warning; compilation succeeded, and the final build emitted no warning. | + +Every final build/test/render process used BelowNormal priority and affinity 0xFF00, limiting it to eight of the 16 logical processors. The build/test/render scripts enforce an 8 GB private-memory guard and 15-minute timeout. No process hit either limit, crashed, or needed a 0xC0000142 retry. The earlier phase runs above exited normally; their topology/seam measurements failed or were superseded as noted. No viewer or owner editor was opened. No owner project was built and no .uasset was touched. + +## Files added for the spike + +New implementation and launch files: + +- [SharpMesherSpike.cpp](SharpMesherSpike.cpp) +- [BuildSharpMesherSpike.ps1](BuildSharpMesherSpike.ps1) +- [RunSharpMesherSpike.ps1](RunSharpMesherSpike.ps1) +- [RenderSharpMesherCaptures.ps1](RenderSharpMesherCaptures.ps1) +- [RenderSharpMesherCaptures.py](RenderSharpMesherCaptures.py) + +The report, CSV measurements, exported OBJ meshes, offline captures, executable/object build outputs, and run logs are also in this WorldAssemblerM2 directory. No files under Source were changed. diff --git a/Standalone/Research/WorldAssemblerM2/RenderSharpMesherCaptures.ps1 b/Standalone/Research/WorldAssemblerM2/RenderSharpMesherCaptures.ps1 new file mode 100644 index 0000000..4588275 --- /dev/null +++ b/Standalone/Research/WorldAssemblerM2/RenderSharpMesherCaptures.ps1 @@ -0,0 +1,53 @@ +[CmdletBinding()] +param() + +$ErrorActionPreference = 'Stop' +$WorkDir = [IO.Path]::GetFullPath($PSScriptRoot) +$Python = (Get-Command python.exe -ErrorAction Stop).Source +$Script = Join-Path $WorkDir 'RenderSharpMesherCaptures.py' +$Stdout = Join-Path $WorkDir 'render.stdout.log' +$Stderr = Join-Path $WorkDir 'render.stderr.log' +$Runner = Get-Process -Id $PID +$Runner.PriorityClass = 'BelowNormal' +$Runner.ProcessorAffinity = [IntPtr]0xFF00 + +$Info = [Diagnostics.ProcessStartInfo]::new() +$Info.FileName = $Python +$Info.Arguments = ('"{0}"' -f $Script) +$Info.WorkingDirectory = $WorkDir +$Info.UseShellExecute = $false +$Info.CreateNoWindow = $true +$Info.RedirectStandardOutput = $true +$Info.RedirectStandardError = $true +$Process = [Diagnostics.Process]::new() +$Process.StartInfo = $Info +if (-not $Process.Start()) { throw 'Could not start capture renderer.' } +$Process.PriorityClass = [Diagnostics.ProcessPriorityClass]::BelowNormal +$Process.ProcessorAffinity = [IntPtr]0xFF00 +$StartAt = [DateTime]::UtcNow +$OutputTask = $Process.StandardOutput.ReadToEndAsync() +$ErrorTask = $Process.StandardError.ReadToEndAsync() +$TimedOut = $false +$OverMemory = $false +$EightGB = 8L * 1024L * 1024L * 1024L +while (-not $Process.HasExited) { + $Process.Refresh() + if ($Process.PrivateMemorySize64 -ge $EightGB) { $OverMemory = $true; break } + if (([DateTime]::UtcNow - $StartAt).TotalSeconds -ge 900) { $TimedOut = $true; break } + Start-Sleep -Milliseconds 250 +} +if ($TimedOut -or $OverMemory) { try { $Process.Kill() } catch {} } +$Process.WaitForExit() +$RenderStdout = $OutputTask.GetAwaiter().GetResult() +$RenderStderr = $ErrorTask.GetAwaiter().GetResult() +[IO.File]::WriteAllText($Stdout, $RenderStdout) +[IO.File]::WriteAllText($Stderr, $RenderStderr) +if ($TimedOut) { throw "Capture rendering exceeded 15 minutes. Logs: $Stdout / $Stderr" } +if ($OverMemory) { throw "Capture renderer reached the 8 GB private-memory cap. Logs: $Stdout / $Stderr" } +if ($Process.ExitCode -ne 0) { + Write-Output "Capture rendering failed exit_code=$($Process.ExitCode) stdout=$Stdout stderr=$Stderr" + if ($RenderStdout) { Write-Output $RenderStdout } + if ($RenderStderr) { Write-Output $RenderStderr } + throw 'Matched capture rendering failed.' +} +Write-Output "Capture render passed elapsed_seconds=$([Math]::Round(([DateTime]::UtcNow - $StartAt).TotalSeconds,2)) priority=BelowNormal affinity=0xFF00 stdout=$Stdout stderr=$Stderr" diff --git a/Standalone/Research/WorldAssemblerM2/RenderSharpMesherCaptures.py b/Standalone/Research/WorldAssemblerM2/RenderSharpMesherCaptures.py new file mode 100644 index 0000000..6391d59 --- /dev/null +++ b/Standalone/Research/WorldAssemblerM2/RenderSharpMesherCaptures.py @@ -0,0 +1,205 @@ +#!/usr/bin/env python3 +"""Render the standalone spike's exported meshes from matched orthographic views.""" +from __future__ import annotations + +import math +from pathlib import Path + +from PIL import Image, ImageDraw + + +ROOT = Path(__file__).resolve().parent +RESULTS = ROOT / "results" +CAPTURES = RESULTS / "captures" +CAPTURES.mkdir(parents=True, exist_ok=True) +METHODS = ("MC", "Hybrid", "DC") +PIECES = ( + ("corner90", "90 degree corner"), + ("doorway", "doorway"), + ("riser1m", "one metre riser"), + ("cliff_overhang", "cliff ledge / overhang"), + ("organic_cave_chamber", "organic cave chamber wall"), +) +WIDTH, HEIGHT = 288, 216 +BG = (239, 242, 244) +INK = (36, 44, 51) +SEAM_FINE = (73, 133, 174) +SEAM_COARSE = (188, 112, 68) + + +def normalize(v): + n = math.sqrt(sum(c * c for c in v)) + return tuple(c / n for c in v) if n > 1.0e-12 else (0.0, 0.0, 1.0) + + +def dot(a, b): + return sum(x * y for x, y in zip(a, b)) + + +def cross(a, b): + return (a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]) + + +def read_obj(path: Path): + vertices, normals, faces = [], [], [] + for line in path.read_text(encoding="ascii").splitlines(): + if line.startswith("v "): + vertices.append(tuple(map(float, line.split()[1:4]))) + elif line.startswith("vn "): + normals.append(tuple(map(float, line.split()[1:4]))) + elif line.startswith("f "): + face = [] + for item in line.split()[1:4]: + fields = item.split("//") + face.append((int(fields[0]) - 1, int(fields[1]) - 1 if len(fields) > 1 else int(fields[0]) - 1)) + faces.append(tuple(face)) + return vertices, normals, faces + + +def project_mesh(vertices, faces, seam=False): + # One view is reused for every candidate in a row: azimuth 38°, elevation 27°. + target = (0.0, 0.0, 0.0) + view_dir = normalize((8.0, -10.0, 7.0)) + forward = tuple(-c for c in view_dir) + right = normalize(cross(forward, (0.0, 0.0, 1.0))) + up = normalize(cross(right, forward)) + scale = 20.2 + projected = [] + for p in vertices: + rel = tuple(p[i] - target[i] for i in range(3)) + projected.append((WIDTH * 0.5 + dot(rel, right) * scale, + HEIGHT * 0.51 - dot(rel, up) * scale, + dot(rel, forward))) + return projected + + +def raster(path: Path, color_by_x=False): + vertices, normals, faces = read_obj(path) + if color_by_x: + clipped_vertices, clipped_normals, clipped_faces = [], [], [] + for face in faces: + center_x = sum(vertices[i][0] for i, _ in face) / 3.0 + left = center_x < 0.0 + polygon = [(vertices[vi], normals[ni] if ni < len(normals) else (0.0, 0.0, 1.0)) for vi, ni in face] + clipped = [] + for i, (pa, na) in enumerate(polygon): + pb, nb = polygon[(i + 1) % len(polygon)] + inside_a = pa[0] <= 0.0 if left else pa[0] >= 0.0 + inside_b = pb[0] <= 0.0 if left else pb[0] >= 0.0 + if inside_a: + clipped.append((pa, na)) + if inside_a != inside_b: + t = (0.0 - pa[0]) / (pb[0] - pa[0]) + point = tuple(pa[k] + t * (pb[k] - pa[k]) for k in range(3)) + normal = normalize(tuple(na[k] + t * (nb[k] - na[k]) for k in range(3))) + clipped.append((point, normal)) + if len(clipped) < 3: + continue + first = len(clipped_vertices) + for point, normal in clipped: + clipped_vertices.append(point) + clipped_normals.append(normal) + for i in range(1, len(clipped) - 1): + clipped_faces.append(((first, first), (first + i, first + i), (first + i + 1, first + i + 1))) + vertices, normals, faces = clipped_vertices, clipped_normals, clipped_faces + projected = project_mesh(vertices, faces, color_by_x) + pixels = Image.new("RGB", (WIDTH, HEIGHT), BG) + zbuf = [float("inf")] * (WIDTH * HEIGHT) + px = pixels.load() + light = normalize((-0.45, -0.55, 0.75)) + for face in faces: + ids = [v for v, _ in face] + points = [projected[i] for i in ids] + x0 = max(0, int(math.floor(min(p[0] for p in points)))) + x1 = min(WIDTH - 1, int(math.ceil(max(p[0] for p in points)))) + y0 = max(0, int(math.floor(min(p[1] for p in points)))) + y1 = min(HEIGHT - 1, int(math.ceil(max(p[1] for p in points)))) + if x0 > x1 or y0 > y1: + continue + (ax, ay, az), (bx, by, bz), (cx, cy, cz) = points + denom = (by - cy) * (ax - cx) + (cx - bx) * (ay - cy) + if abs(denom) < 1.0e-9: + continue + face_normals = [normals[ni] if ni < len(normals) else (0.0, 0.0, 1.0) for _, ni in face] + for y in range(y0, y1 + 1): + for x in range(x0, x1 + 1): + sx, sy = x + 0.5, y + 0.5 + u = ((by - cy) * (sx - cx) + (cx - bx) * (sy - cy)) / denom + v = ((cy - ay) * (sx - cx) + (ax - cx) * (sy - cy)) / denom + w = 1.0 - u - v + if u < -1.0e-8 or v < -1.0e-8 or w < -1.0e-8: + continue + depth = u * az + v * bz + w * cz + idx = y * WIDTH + x + if depth >= zbuf[idx]: + continue + zbuf[idx] = depth + n = normalize(tuple(u * face_normals[0][k] + v * face_normals[1][k] + w * face_normals[2][k] for k in range(3))) + shade = 0.36 + 0.64 * max(0.0, dot(n, light)) + if color_by_x: + center_x = sum(vertices[i][0] for i in ids) / 3.0 + base = SEAM_FINE if center_x < 0.0 else SEAM_COARSE + else: + base = (130, 145, 151) + px[x, y] = tuple(max(0, min(255, int(c * shade))) for c in base) + if color_by_x: + draw = ImageDraw.Draw(pixels) + # Project a visible segment on the LOD boundary plane, x=0. + seam_world = [(-0.0, y, z) for y, z in ((-4, -4), (4, -4), (4, 4), (-4, 4), (-4, -4))] + seam_points = project_mesh(seam_world, [], True) + draw.line([(p[0], p[1]) for p in seam_points], fill=(181, 55, 49), width=2) + return pixels + + +def capture_title(image: Image.Image, title: str, subtitle: str): + canvas = Image.new("RGB", (WIDTH, HEIGHT + 30), BG) + canvas.paste(image, (0, 30)) + draw = ImageDraw.Draw(canvas) + draw.rectangle((0, 0, WIDTH, 29), fill=(31, 42, 50)) + draw.text((8, 4), title, fill=(255, 255, 255)) + draw.text((8, 16), subtitle, fill=(196, 210, 218)) + return canvas + + +def main(): + captures = [] + for key, label in PIECES: + row = [] + for method in METHODS: + obj = RESULTS / "meshes" / f"{key}_{method}.obj" + image = capture_title(raster(obj), label, f"{method} | matched orthographic view") + path = CAPTURES / f"{key}_{method}.jpg" + image.save(path, quality=86, optimize=True) + row.append(image) + captures.append(row) + seam_row = [] + for method in METHODS: + obj = RESULTS / "meshes" / f"lodseam_cliff_overhang_{method}.obj" + image = capture_title(raster(obj, color_by_x=True), "LOD0 to LOD1 cliff seam", f"{method} | fine left, coarse right, x=0 seam") + path = CAPTURES / f"lodseam_{method}.jpg" + image.save(path, quality=86, optimize=True) + seam_row.append(image) + captures.append(seam_row) + + sheet = Image.new("RGB", (WIDTH * 3, (HEIGHT + 30) * 6), (221, 227, 230)) + draw = ImageDraw.Draw(sheet) + for row_index, row in enumerate(captures): + for col_index, image in enumerate(row): + x, y = col_index * WIDTH, row_index * (HEIGHT + 30) + sheet.paste(image, (x, y)) + draw.rectangle((x, y, x + WIDTH - 1, y + HEIGHT + 29), outline=(73, 83, 89), width=1) + sheet_path = CAPTURES / "contact_sheet.jpg" + sheet.save(sheet_path, quality=86, optimize=True) + oversized = [] + for path in sorted(CAPTURES.glob("*.jpg")): + size = path.stat().st_size + print(f"capture file={path.name} bytes={size}") + if size >= 300_000: + oversized.append((path.name, size)) + if oversized: + raise RuntimeError(f"captures exceed 300 KB: {oversized}") + print(f"contact_sheet={sheet_path} bytes={sheet_path.stat().st_size} image_count={len(list(CAPTURES.glob('*.jpg')))}") + + +if __name__ == "__main__": + main() diff --git a/Standalone/Research/WorldAssemblerM2/RunSharpMesherSpike.ps1 b/Standalone/Research/WorldAssemblerM2/RunSharpMesherSpike.ps1 new file mode 100644 index 0000000..1c63234 --- /dev/null +++ b/Standalone/Research/WorldAssemblerM2/RunSharpMesherSpike.ps1 @@ -0,0 +1,69 @@ +[CmdletBinding()] +param( + [ValidateRange(2,5)][int]$PhaseCount = 5 +) + +$ErrorActionPreference = 'Stop' +$WorkDir = [IO.Path]::GetFullPath($PSScriptRoot) +$Exe = Join-Path $WorkDir 'bin\SharpMesherSpike.exe' +$OutDir = Join-Path $WorkDir 'results' +if (-not (Test-Path -LiteralPath $Exe -PathType Leaf)) { throw "Spike executable missing: $Exe" } +[IO.Directory]::CreateDirectory($OutDir) | Out-Null + +$Runner = Get-Process -Id $PID +$Runner.PriorityClass = 'BelowNormal' +$Runner.ProcessorAffinity = [IntPtr]0xFF00 +$Affinity = [IntPtr]0xFF00 +$Stdout = Join-Path $WorkDir 'run.stdout.log' +$Stderr = Join-Path $WorkDir 'run.stderr.log' +$EightGB = 8L * 1024L * 1024L * 1024L +$Attempt = 0 +do { + $Attempt++ + Remove-Item -LiteralPath $Stdout, $Stderr -Force -ErrorAction SilentlyContinue + $Info = [Diagnostics.ProcessStartInfo]::new() + $Info.FileName = $Exe + $Info.Arguments = ('--out="{0}" --phase-count={1}' -f $OutDir, $PhaseCount) + $Info.WorkingDirectory = [IO.Path]::GetFullPath((Join-Path $WorkDir '..\..')) + $Info.UseShellExecute = $false + $Info.CreateNoWindow = $true + $Info.RedirectStandardOutput = $true + $Info.RedirectStandardError = $true + $Process = [Diagnostics.Process]::new() + $Process.StartInfo = $Info + if (-not $Process.Start()) { throw 'Could not start SharpMesherSpike.exe.' } + $Process.PriorityClass = [Diagnostics.ProcessPriorityClass]::BelowNormal + $Process.ProcessorAffinity = $Affinity + $StartAt = [DateTime]::UtcNow + $OutputTask = $Process.StandardOutput.ReadToEndAsync() + $ErrorTask = $Process.StandardError.ReadToEndAsync() + $TimedOut = $false + $OverMemory = $false + while (-not $Process.HasExited) { + $Process.Refresh() + if ($Process.PrivateMemorySize64 -ge $EightGB) { $OverMemory = $true; break } + if (([DateTime]::UtcNow - $StartAt).TotalSeconds -ge 900) { $TimedOut = $true; break } + Start-Sleep -Milliseconds 250 + } + if ($TimedOut -or $OverMemory) { try { $Process.Kill() } catch {} } + $Process.WaitForExit() + $RunStdout = $OutputTask.GetAwaiter().GetResult() + $RunStderr = $ErrorTask.GetAwaiter().GetResult() + [IO.File]::WriteAllText($Stdout, $RunStdout) + [IO.File]::WriteAllText($Stderr, $RunStderr) + if ($TimedOut) { throw "Spike exceeded 15 minutes and was stopped. Logs: $Stdout / $Stderr" } + if ($OverMemory) { throw "Spike reached the 8 GB private-memory cap and was stopped. Logs: $Stdout / $Stderr" } + if ($Process.ExitCode -eq -1073741502 -and $Attempt -eq 1) { + Write-Output 'Spike exited 0xC0000142; waiting 60 seconds before the one allowed retry.' + Start-Sleep -Seconds 60 + continue + } + if ($Process.ExitCode -ne 0) { + Write-Output "Spike failed exit_code=$($Process.ExitCode) stdout=$Stdout stderr=$Stderr" + if ($RunStdout) { Write-Output $RunStdout } + if ($RunStderr) { Write-Output $RunStderr } + throw 'Sharp mesher spike failed.' + } +} while ($Process.ExitCode -eq -1073741502 -and $Attempt -lt 2) + +Write-Output "Spike passed elapsed_seconds=$([Math]::Round(([DateTime]::UtcNow - $StartAt).TotalSeconds,2)) phase_count=$PhaseCount priority=BelowNormal affinity=0xFF00 stdout=$Stdout stderr=$Stderr" diff --git a/Standalone/Research/WorldAssemblerM2/SharpMesherSpike.cpp b/Standalone/Research/WorldAssemblerM2/SharpMesherSpike.cpp new file mode 100644 index 0000000..fe9a707 --- /dev/null +++ b/Standalone/Research/WorldAssemblerM2/SharpMesherSpike.cpp @@ -0,0 +1,1174 @@ +#include "CoreMinimal.h" +#include "MarchingCubesTables.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace +{ + constexpr double VoxelMeters = 0.25; + constexpr double PositionQuantum = VoxelMeters / 65536.0; + constexpr double NormalSplitCos = 0.86602540378443864676; // 30 degrees + constexpr int CornerCount = 8; + constexpr int EdgeCount = 12; + constexpr std::array, CornerCount> CornerOffsets{{ + {{0,0,0}}, {{1,0,0}}, {{1,1,0}}, {{0,1,0}}, + {{0,0,1}}, {{1,0,1}}, {{1,1,1}}, {{0,1,1}} + }}; + constexpr std::array, EdgeCount> EdgeCorners{{ + {{0,1}}, {{1,2}}, {{2,3}}, {{3,0}}, {{4,5}}, {{5,6}}, + {{6,7}}, {{7,4}}, {{0,4}}, {{1,5}}, {{2,6}}, {{3,7}} + }}; + + struct V3 + { + double x=0, y=0, z=0; + V3 operator+(V3 b) const { return {x+b.x,y+b.y,z+b.z}; } + V3 operator-(V3 b) const { return {x-b.x,y-b.y,z-b.z}; } + V3 operator*(double s) const { return {x*s,y*s,z*s}; } + V3 operator/(double s) const { return {x/s,y/s,z/s}; } + V3& operator+=(V3 b) { x+=b.x; y+=b.y; z+=b.z; return *this; } + }; + double Dot(V3 a,V3 b){return a.x*b.x+a.y*b.y+a.z*b.z;} + V3 Cross(V3 a,V3 b){return {a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x};} + double Length(V3 a){return std::sqrt(Dot(a,a));} + V3 Normalize(V3 a) + { + const double n=Length(a); + return n>1.0e-20?a/n:V3{0,0,1}; + } + V3 Quantize(V3 p) + { + const auto q=[](double x){return std::nearbyint(x/PositionQuantum)*PositionQuantum;}; + return {q(p.x),q(p.y),q(p.z)}; + } + + template void ParallelFor(std::size_t count,int workers,F&& function) + { + workers=std::clamp(workers,1,16); + if(workers==1||count<512) + { + for(std::size_t i=0;i(workers,static_cast(count)); + std::vector threads; + threads.reserve(actual); + for(int worker=0;worker(worker)/actual; + const std::size_t end=count*static_cast(worker+1)/actual; + threads.emplace_back([begin,end,&function](){for(std::size_t i=begin;i planes{}; + std::uint8_t planeCount=0; + V3 gradient{}; + bool sharp() const { return planeCount!=0; } + std::uint64_t primaryFeature() const + { + if(!planeCount)return 0; + std::uint64_t id=planes[0].id; + for(std::uint8_t i=1;i planes{}; + std::uint8_t planeCount=0; + V3 gradient{}; + }; + + BoxSample EvalBox(const Box& box,V3 p) + { + const V3 center=(box.low+box.high)*0.5; + const V3 half=(box.high-box.low)*0.5; + const V3 delta=p-center; + const V3 q{std::abs(delta.x)-half.x,std::abs(delta.y)-half.y,std::abs(delta.z)-half.z}; + const V3 outside{std::max(q.x,0.0),std::max(q.y,0.0),std::max(q.z,0.0)}; + BoxSample out; + out.density=Length(outside)+std::min(std::max(q.x,std::max(q.y,q.z)),0.0); + const double maxQ=std::max(q.x,std::max(q.y,q.z)); + V3 gradient{}; + if(Length(outside)>1.0e-14) + { + gradient=Normalize({outside.x*(delta.x<0?-1:1),outside.y*(delta.y<0?-1:1),outside.z*(delta.z<0?-1:1)}); + } + for(int axis=0;axis<3;++axis) + { + const double qAxis=axis==0?q.x:(axis==1?q.y:q.z); + if(std::abs(qAxis-maxQ)>1.0e-10)continue; + const double dAxis=axis==0?delta.x:(axis==1?delta.y:delta.z); + V3 normal{}; + if(axis==0)normal.x=dAxis<0?-1:1; + if(axis==1)normal.y=dAxis<0?-1:1; + if(axis==2)normal.z=dAxis<0?-1:1; + const std::uint64_t side=dAxis<0?0ull:1ull; + const std::uint64_t feature=(static_cast(box.id)<<4) + |static_cast(axis*2)+side+1ull; + if(box.inverted)normal=normal*-1.0; + out.planes[out.planeCount++]={feature,normal}; + gradient+=normal; + } + if(box.inverted) + { + out.density=-out.density; + gradient=gradient*-1.0; + } + out.gradient=Length(gradient)>1.0e-14?Normalize(gradient):Normalize(outside); + return out; + } + + enum class Piece : int { Corner=0, Doorway=1, Step=2, Ledge=3, Organic=4 }; + const char* PieceName(Piece p) + { + switch(p) + { + case Piece::Corner:return "corner90"; + case Piece::Doorway:return "doorway"; + case Piece::Step:return "riser1m"; + case Piece::Ledge:return "cliff_overhang"; + default:return "organic_cave_chamber"; + } + } + struct Field + { + Piece piece=Piece::Corner; + V3 phaseShift{}; + std::array boxes{}; + std::uint8_t boxCount=0; + Box doorVoid{}; + void AddBox(V3 low,V3 high,std::uint32_t id) + { boxes[boxCount++]={low,high,id,false}; } + explicit Field(Piece p,V3 phase={}):piece(p),phaseShift(phase) + { + switch(p) + { + case Piece::Corner: + AddBox({-2.5,-2.5,-1.5},{0.25,2.5,1.5},1); + AddBox({-2.5,-0.25,-1.5},{2.5,2.5,1.5},2); + break; + case Piece::Doorway: + AddBox({-0.5,-2.5,-2.0},{0.5,2.5,3.0},3); + doorVoid={{-0.8,-0.85,-1.75},{0.8,0.85,2.0},4,true}; + break; + case Piece::Step: + AddBox({-3.0,-3.0,-2.0},{3.0,0.25,0.0},5); + AddBox({-3.0,-0.25,-2.0},{3.0,3.0,1.0},6); + break; + case Piece::Ledge: + AddBox({-3.0,-3.0,-2.0},{0.25,3.0,3.0},7); + AddBox({-0.25,-3.0,1.5},{2.2,3.0,2.0},8); + break; + case Piece::Organic:break; + } + } + + FieldSample Eval(V3 p) const + { + p=p+phaseShift; + if(piece==Piece::Organic) + { + const V3 d=p-V3{0,0,0}; + FieldSample result; + // Exact spherical chamber cavity: negative solid outside, positive air inside. + result.density=2.0-Length(d); + result.gradient=Normalize(d)*-1.0; + return result; + } + FieldSample result; + result.density=std::numeric_limits::infinity(); + for(std::uint8_t i=0;iresult.density+1.0e-10) + { + result.density=cutDensity; + result.planeCount=cut.planeCount; + for(std::uint8_t k=0;k1.0e-14?Normalize(sum):V3{0,0,1}; + return result; + } + }; + + struct Grid + { + int nx=0,ny=0,nz=0,step=1; + V3 origin{}; + std::vector density; + std::size_t NodeIndex(int x,int y,int z) const + { return (static_cast(z)*(ny+1)+y)*(nx+1)+x; } + V3 Position(int x,int y,int z) const + { return origin+V3{x*step*VoxelMeters,y*step*VoxelMeters,z*step*VoxelMeters}; } + std::size_t CellIndex(int x,int y,int z) const + { return (static_cast(z)*ny+y)*nx+x; } + std::size_t CellCount() const { return static_cast(nx)*ny*nz; } + }; + Grid SampleGrid(const Field& field,int nx,int ny,int nz,int step,V3 origin,int workers) + { + Grid grid;grid.nx=nx;grid.ny=ny;grid.nz=nz;grid.step=step;grid.origin=origin; + const std::size_t count=static_cast(nx+1)*(ny+1)*(nz+1); + grid.density.resize(count); + ParallelFor(count,workers,[&](std::size_t index) + { + const int x=static_cast(index%(nx+1)); + const int y=static_cast((index/(nx+1))%(ny+1)); + const int z=static_cast(index/((nx+1)*(ny+1))); + const double density=field.Eval(grid.Position(x,y,z)).density; + // Stable iso-level tie rule: zero remains air, matching the Core sign convention, but + // receives a sub-cell positive bias so edge roots do not collapse onto one lattice + // vertex from several unrelated edges. This is only 0.01 mm and is below the output + // surface-error precision; it is shared by all three candidates. + grid.density[index]=density==0.0?1.0e-5:density; + }); + return grid; + } + + struct Root + { + bool active=false; + V3 position{},gradient{}; + std::array planes{}; + std::uint8_t planeCount=0; + std::uint64_t feature=0; + }; + struct EdgeArrays + { + std::array,3> roots; + std::array,3> vertexIds; + std::array counts{}; + std::size_t Total() const{return counts[0]+counts[1]+counts[2];} + std::size_t Index(int axis,int x,int y,int z,const Grid& g) const + { + if(axis==0)return (static_cast(z)*(g.ny+1)+y)*g.nx+x; + if(axis==1)return (static_cast(z)*g.ny+y)*(g.nx+1)+x; + return (static_cast(z)*(g.ny+1)+y)*(g.nx+1)+x; + } + Root& Get(int axis,int x,int y,int z,const Grid& g){return roots[axis][Index(axis,x,y,z,g)];} + const Root& Get(int axis,int x,int y,int z,const Grid& g) const{return roots[axis][Index(axis,x,y,z,g)];} + int& Vertex(int axis,int x,int y,int z,const Grid& g){return vertexIds[axis][Index(axis,x,y,z,g)];} + int Vertex(int axis,int x,int y,int z,const Grid& g) const{return vertexIds[axis][Index(axis,x,y,z,g)];} + }; + + EdgeArrays BuildRoots(const Grid& grid,const Field& field,int workers) + { + EdgeArrays edges; + edges.counts={static_cast(grid.nx)*(grid.ny+1)*(grid.nz+1), + static_cast(grid.nx+1)*grid.ny*(grid.nz+1), + static_cast(grid.nx+1)*(grid.ny+1)*grid.nz}; + for(int axis=0;axis<3;++axis) + { + edges.roots[axis].resize(edges.counts[axis]); + edges.vertexIds[axis].assign(edges.counts[axis],-1); + } + const std::size_t total=edges.Total(); + ParallelFor(total,workers,[&](std::size_t flat) + { + int axis=0;std::size_t index=flat; + while(axis<2&&index>=edges.counts[axis]){index-=edges.counts[axis];++axis;} + int x=0,y=0,z=0; + if(axis==0) + { + x=static_cast(index%grid.nx);y=static_cast((index/grid.nx)%(grid.ny+1)); + z=static_cast(index/(grid.nx*(grid.ny+1))); + } + else if(axis==1) + { + x=static_cast(index%(grid.nx+1));y=static_cast((index/(grid.nx+1))%grid.ny); + z=static_cast(index/((grid.nx+1)*grid.ny)); + } + else + { + x=static_cast(index%(grid.nx+1));y=static_cast((index/(grid.nx+1))%(grid.ny+1)); + z=static_cast(index/((grid.nx+1)*(grid.ny+1))); + } + int x2=x,y2=y,z2=z; + if(axis==0)++x2;else if(axis==1)++y2;else ++z2; + const double d0=grid.density[grid.NodeIndex(x,y,z)]; + const double d1=grid.density[grid.NodeIndex(x2,y2,z2)]; + if((d0<0.0)==(d1<0.0))return; + double t=d0/(d0-d1);t=std::clamp(t,0.0,1.0); + Root root;root.active=true; + root.position=Quantize(grid.Position(x,y,z)+(grid.Position(x2,y2,z2)-grid.Position(x,y,z))*t); + const FieldSample sample=field.Eval(root.position); + const double h=grid.step*VoxelMeters; + const auto gradientAt=[&](V3 point) + { + const double dx=field.Eval(point+V3{h,0,0}).density-field.Eval(point-V3{h,0,0}).density; + const double dy=field.Eval(point+V3{0,h,0}).density-field.Eval(point-V3{0,h,0}).density; + const double dz=field.Eval(point+V3{0,0,h}).density-field.Eval(point-V3{0,0,h}).density; + return Normalize({dx,dy,dz}); + }; + const V3 g0=gradientAt(grid.Position(x,y,z)); + const V3 g1=gradientAt(grid.Position(x2,y2,z2)); + root.gradient=Normalize(g0*(1.0-t)+g1*t); + root.planeCount=sample.planeCount; + for(std::uint8_t i=0;i index{}; + std::uint64_t feature=0; + }; + struct Mesh + { + std::vector vertices; + std::vector triangles; + std::uint64_t hash=0; + std::uint64_t skippedDualEdges=0; + }; + enum class Method : int { MC=0, Hybrid=1, DC=2 }; + const char* MethodName(Method m){return m==Method::MC?"MC":m==Method::Hybrid?"Hybrid": "DC";} + + struct Cell + { + std::uint8_t cubeCase=0; + bool active=false; + std::uint8_t componentCount=0; + std::array edgeComponent{}; + std::array componentSharp{},componentUseDC{}; + std::array dual{}; + }; + + int FindLocal(std::array& parent,int value) + { + while(parent[value]!=value) + { + parent[value]=parent[parent[value]]; + value=parent[value]; + } + return value; + } + void UniteLocal(std::array& parent,int a,int b) + { + a=FindLocal(parent,a);b=FindLocal(parent,b);if(a==b)return; + if(a>b)std::swap(a,b); + parent[b]=static_cast(a); + } + + struct ComponentSets + { + std::vector parent; + std::vector rank,sharp; + explicit ComponentSets(std::size_t size):parent(size),rank(size,0),sharp(size,0) + {std::iota(parent.begin(),parent.end(),0u);} + std::uint32_t Find(std::uint32_t value) + { + while(parent[value]!=value){parent[value]=parent[parent[value]];value=parent[value];} + return value; + } + void Unite(std::uint32_t a,std::uint32_t b) + { + a=Find(a);b=Find(b);if(a==b)return; + if(rank[a](sharp[a]|sharp[b]); + if(rank[a]==rank[b])++rank[a]; + } + }; + + std::array LocalEdge(int x,int y,int z,int edge) + { + const auto c=EdgeCorners[edge];const auto a=CornerOffsets[c[0]],b=CornerOffsets[c[1]]; + const int axis=a[0]!=b[0]?0:a[1]!=b[1]?1:2; + return {axis,x+std::min(a[0],b[0]),y+std::min(a[1],b[1]),z+std::min(a[2],b[2])}; + } + int FindLocalEdge(int cellX,int cellY,int cellZ,int axis,int edgeX,int edgeY,int edgeZ) + { + for(int edge=0;edge>& constraints,V3 fallback,V3 low,V3 high) + { + if(constraints.empty())return fallback; + double a[3][4]{}; + // Fixed constraint order; a small center penalty resolves rank-deficient planes without + // moving a two/three-plane intersection at measurable scale. + constexpr double lambda=1.0e-10; + const double c[3]={fallback.x,fallback.y,fallback.z}; + for(int i=0;i<3;++i){a[i][i]=lambda;a[i][3]=lambda*c[i];} + for(const auto& constraint:constraints) + { + const V3 n=Normalize(constraint.second);const double b=Dot(n,constraint.first); + const double nv[3]={n.x,n.y,n.z}; + for(int r=0;r<3;++r){for(int k=0;k<3;++k)a[r][k]+=nv[r]*nv[k];a[r][3]+=nv[r]*b;} + } + bool singular=false; + for(int column=0;column<3;++column) + { + int pivot=column; + for(int row=column+1;row<3;++row)if(std::abs(a[row][column])>std::abs(a[pivot][column]))pivot=row; + if(std::abs(a[pivot][column])<1.0e-14){singular=true;break;} + if(pivot!=column)for(int k=column;k<4;++k)std::swap(a[pivot][k],a[column][k]); + const double divisor=a[column][column]; + for(int k=column;k<4;++k)a[column][k]/=divisor; + for(int row=0;row<3;++row)if(row!=column) + { + const double factor=a[row][column]; + for(int k=column;k<4;++k)a[row][k]-=factor*a[column][k]; + } + } + if(singular)return fallback; + return Quantize({std::clamp(a[0][3],low.x,high.x), + std::clamp(a[1][3],low.y,high.y),std::clamp(a[2][3],low.z,high.z)}); + } + + std::uint64_t SelectFeature(const Root& root,int edgeAxis) + { + if(!root.planeCount)return 0; + std::uint8_t selected=0; + for(std::uint8_t i=1;ib+1.0e-12||(std::abs(a-b)<=1.0e-12&&root.planes[i].id(word>>(i*8));h*=1099511628211ull;} + }; + add(mesh.vertices.size(),8); + for(const Vertex& v:mesh.vertices) + { + for(double x:{v.position.x,v.position.y,v.position.z,v.normal.x,v.normal.y,v.normal.z}) + add(std::bit_cast(x),8); + } + add(mesh.triangles.size(),8); + for(const Triangle& t:mesh.triangles) + {add(t.index[0],4);add(t.index[1],4);add(t.index[2],4);add(t.feature,8);} + return h; + } + + Mesh BuildMesh(const Grid& grid,const Field& field,Method method,int workers) + { + EdgeArrays edges=BuildRoots(grid,field,workers); + std::vector cells(grid.CellCount()); + ParallelFor(cells.size(),workers,[&](std::size_t ci) + { + const int x=static_cast(ci%grid.nx); + const int y=static_cast((ci/grid.nx)%grid.ny); + const int z=static_cast(ci/(grid.nx*grid.ny)); + Cell& cell=cells[ci];cell.edgeComponent.fill(0xff); + double values[CornerCount]; + for(int corner=0;corner(1u< localParent{}; + for(int edge=0;edge(edge); + for(int k=0;k<16&&MCTriTable[cell.cubeCase][k]>=0;k+=3) + { + const int a=MCTriTable[cell.cubeCase][k],b=MCTriTable[cell.cubeCase][k+1],c=MCTriTable[cell.cubeCase][k+2]; + UniteLocal(localParent,a,b);UniteLocal(localParent,b,c); + } + std::array rootToComponent{};rootToComponent.fill(0xff); + for(int edge=0;edge patches; + if(method==Method::Hybrid) + { + patches=std::make_unique(cells.size()*EdgeCount); + for(std::size_t cellIndex=0;cellIndexsharp[cellIndex*EdgeCount+component]=cells[cellIndex].componentSharp[component]?1:0; + for(int axis=0;axis<3;++axis) + { + for(std::size_t i=0;i(i%grid.nx);y=static_cast((i/grid.nx)%(grid.ny+1));z=static_cast(i/(grid.nx*(grid.ny+1)));} + else if(axis==1){x=static_cast(i%(grid.nx+1));y=static_cast((i/(grid.nx+1))%grid.ny);z=static_cast(i/((grid.nx+1)*grid.ny));} + else{x=static_cast(i%(grid.nx+1));y=static_cast((i/(grid.nx+1))%(grid.ny+1));z=static_cast(i/((grid.nx+1)*(grid.ny+1)));} + std::array,4> around{}; + if(axis==0)around={{{x,y-1,z-1},{x,y,z-1},{x,y,z},{x,y-1,z}}}; + else if(axis==1)around={{{x-1,y,z-1},{x-1,y,z},{x,y,z},{x,y,z-1}}}; + else around={{{x-1,y-1,z},{x,y-1,z},{x,y,z},{x-1,y,z}}}; + std::array patchIds{};bool valid=true; + for(int corner=0;corner<4;++corner) + { + const auto& c=around[corner]; + if(c[0]<0||c[1]<0||c[2]<0||c[0]>=grid.nx||c[1]>=grid.ny||c[2]>=grid.nz){valid=false;break;} + const std::size_t cellIndex=grid.CellIndex(c[0],c[1],c[2]); + const int localEdge=FindLocalEdge(c[0],c[1],c[2],axis,x,y,z); + if(localEdge<0||cells[cellIndex].edgeComponent[localEdge]==0xff){valid=false;break;} + patchIds[corner]=static_cast(cellIndex*EdgeCount+cells[cellIndex].edgeComponent[localEdge]); + } + if(valid)for(int corner=1;corner<4;++corner)patches->Unite(patchIds[0],patchIds[corner]); + } + } + } + for(std::size_t cellIndex=0;cellIndex(cellIndex*EdgeCount+component); + cells[cellIndex].componentUseDC[component]=method==Method::DC|| + (method==Method::Hybrid&&patches->sharp[patches->Find(patch)]!=0); + } + + ParallelFor(cells.size(),workers,[&](std::size_t ci) + { + Cell& cell=cells[ci];if(cell.componentCount==0)return; + const int x=static_cast(ci%grid.nx); + const int y=static_cast((ci/grid.nx)%grid.ny); + const int z=static_cast(ci/(grid.nx*grid.ny)); + for(std::uint8_t component=0;component> constraints; + V3 rootSum{};int rootCount=0; + for(int edge=0;edge(rootCount):(low+high)*0.5; + cell.dual[component]=SolveQEF(constraints,center,low,high); + } + }); + + Mesh mesh; + std::unordered_map dcVertexMap; + if(method==Method::MC||method==Method::Hybrid) + { + std::map,std::uint32_t> positionMap; + for(int axis=0;axis<3;++axis)for(std::size_t index=0;index(std::llround(x/PositionQuantum));}; + const V3 p=root.position; + const auto key=std::tuple{q(p.x),q(p.y),q(p.z)}; + auto found=positionMap.find(key); + std::uint32_t vertex=0; + if(found==positionMap.end()) + { + vertex=static_cast(mesh.vertices.size()); + mesh.vertices.push_back({p,root.gradient});positionMap.emplace(key,vertex); + } + else vertex=found->second; + edges.vertexIds[axis][index]=static_cast(vertex); + } + // Coincident zero-at-node roots use the deterministic first edge normal. The shared + // vertex position, matching the owner's 0.01 cm-style MC welding contract. + } + const auto getDCVertex=[&](int cellId,std::uint8_t component,const Root& root,int axis)->std::uint32_t + { + const std::uint64_t feature=SelectFeature(root,axis); + const std::uint64_t patch=static_cast(cellId)*EdgeCount+component; + const std::uint64_t key=(patch<<32)^(feature&0xffffffffull); + const auto existing=dcVertexMap.find(key); + if(existing!=dcVertexMap.end())return existing->second; + const V3 normal=feature?NormalForFeature(root,feature):root.gradient; + const std::uint32_t index=static_cast(mesh.vertices.size()); + mesh.vertices.push_back({cells[cellId].dual[component],Normalize(normal)}); + dcVertexMap.emplace(key,index); + return index; + }; + for(int z=0;z localVertices{};localVertices.fill(-1); + for(int edge=0;edge=0;k+=3) + { + const int e0=MCTriTable[cell.cubeCase][k],e1=MCTriTable[cell.cubeCase][k+1],e2=MCTriTable[cell.cubeCase][k+2]; + if(localVertices[e0]<0||localVertices[e1]<0||localVertices[e2]<0)continue; + const std::uint8_t component=cell.edgeComponent[e0]; + if(component==0xff||cell.componentUseDC[component])continue; + std::array labels{}; + const std::array triEdges{{e0,e1,e2}}; + for(int edgeIndex=0;edgeIndex<3;++edgeIndex) + { + const auto le=LocalEdge(x,y,z,triEdges[edgeIndex]); + const Root& root=edges.Get(le[0],le[1],le[2],le[3],grid); + labels[edgeIndex]=root.feature; + } + std::uint64_t label=0;int bestVotes=0; + for(std::uint64_t candidate:labels)if(candidate) + { + const int votes=static_cast(std::count(labels.begin(),labels.end(),candidate)); + if(votes>bestVotes||(votes==bestVotes&&(!label||candidate(localVertices[e0]), + static_cast(localVertices[e1]),static_cast(localVertices[e2]), + label,field.Eval(center).gradient); + } + } + for(int axis=0;axis<3;++axis) + { + for(std::size_t index=0;index(index%grid.nx);y=static_cast((index/grid.nx)%(grid.ny+1));z=static_cast(index/(grid.nx*(grid.ny+1)));} + else if(axis==1){x=static_cast(index%(grid.nx+1));y=static_cast((index/(grid.nx+1))%grid.ny);z=static_cast(index/((grid.nx+1)*grid.ny));} + else{x=static_cast(index%(grid.nx+1));y=static_cast((index/(grid.nx+1))%(grid.ny+1));z=static_cast(index/((grid.nx+1)*(grid.ny+1)));} + std::array,4> around{}; + if(axis==0)around={{{x,y-1,z-1},{x,y,z-1},{x,y,z},{x,y-1,z}}}; + else if(axis==1)around={{{x-1,y,z-1},{x-1,y,z},{x,y,z},{x,y,z-1}}}; + else around={{{x-1,y-1,z},{x,y-1,z},{x,y,z},{x-1,y,z}}}; + std::array incident{};bool valid=true;bool anyDC=false; + std::array incidentComponents{}; + for(int i=0;i<4;++i) + { + const auto& c=around[i]; + if(c[0]<0||c[1]<0||c[2]<0||c[0]>=grid.nx||c[1]>=grid.ny||c[2]>=grid.nz) + {valid=false;break;} + incident[i]=static_cast(grid.CellIndex(c[0],c[1],c[2])); + const int localEdge=FindLocalEdge(c[0],c[1],c[2],axis,x,y,z); + if(localEdge<0){valid=false;break;} + incidentComponents[i]=cells[incident[i]].edgeComponent[localEdge]; + if(incidentComponents[i]==0xff){valid=false;break;} + anyDC|=cells[incident[i]].componentUseDC[incidentComponents[i]]; + } + if(!valid||!anyDC||!cells[incident[0]].componentUseDC[incidentComponents[0]])continue; + std::array v{}; + for(int i=0;i<4;++i) + { + if(!cells[incident[i]].componentUseDC[incidentComponents[i]]){valid=false;break;} + v[i]=getDCVertex(incident[i],incidentComponents[i],root,axis); + } + if(!valid){++mesh.skippedDualEdges;continue;} + if(Dot(Cross(mesh.vertices[v[1]].position-mesh.vertices[v[0]].position, + mesh.vertices[v[2]].position-mesh.vertices[v[0]].position),root.gradient)<0.0) + std::swap(v[1],v[3]); + const std::uint64_t label=SelectFeature(root,axis); + const std::uint32_t edgeTie=static_cast(x+y+z+axis)&1u; + if(edgeTie==0) + { + AddTriangle(mesh,v[0],v[1],v[2],label,root.gradient); + AddTriangle(mesh,v[0],v[2],v[3],label,root.gradient); + } + else + { + AddTriangle(mesh,v[0],v[1],v[3],label,root.gradient); + AddTriangle(mesh,v[1],v[2],v[3],label,root.gradient); + } + } + } + mesh.hash=HashMesh(mesh); + return mesh; + } + + struct ErrorMetrics{double mean=0,maximum=0;std::uint64_t samples=0;}; + ErrorMetrics SurfaceError(const Mesh& mesh,const Field& field) + { + ErrorMetrics out;double total=0; + const auto sample=[&](V3 p) + { + const double error=std::abs(field.Eval(p).density); + total+=error;out.maximum=std::max(out.maximum,error);++out.samples; + }; + for(const Triangle& t:mesh.triangles) + { + const V3 a=mesh.vertices[t.index[0]].position,b=mesh.vertices[t.index[1]].position,c=mesh.vertices[t.index[2]].position; + sample(a);sample(b);sample(c);sample((a+b)*0.5);sample((b+c)*0.5);sample((c+a)*0.5);sample((a+b+c)/3.0); + } + if(out.samples)out.mean=total/out.samples; + return out; + } + + struct GeomKey + { + std::int64_t x=0,y=0,z=0; + friend bool operator<(const GeomKey& a,const GeomKey& b) + {return std::tie(a.x,a.y,a.z)(std::llround(p.x/PositionQuantum)), + static_cast(std::llround(p.y/PositionQuantum)), + static_cast(std::llround(p.z/PositionQuantum))};} + struct EdgeInfo{int count=0;std::vector> faces;}; + std::map,EdgeInfo> GeomEdges(const Mesh& mesh) + { + std::map,EdgeInfo> result; + for(std::size_t ti=0;ti,3> edges{{{t.index[0],t.index[1]},{t.index[1],t.index[2]},{t.index[2],t.index[0]}}}; + for(const auto& e:edges) + { + GeomKey ka=Key(mesh.vertices[e.first].position),kb=Key(mesh.vertices[e.second].position); + if(kb(ti),normal); + } + } + return result; + } + struct TopologyMetrics + { + std::uint64_t boundaryEdges=0,nonManifoldEdges=0,creaseEdgePairs=0,splitPairs=0; + }; + TopologyMetrics AnalyzeTopology(const Mesh& mesh) + { + TopologyMetrics out; + const auto edges=GeomEdges(mesh); + for(const auto& [key,info]:edges) + { + if(info.count==1)++out.boundaryEdges; + else if(info.count!=2)++out.nonManifoldEdges; + if(info.faces.size()==2) + { + const Triangle& a=mesh.triangles[info.faces[0].first]; + const Triangle& b=mesh.triangles[info.faces[1].first]; + if(a.feature&&b.feature&&a.feature!=b.feature&& + Dot(Normalize(info.faces[0].second),Normalize(info.faces[1].second))sideFilter)continue; + if(side>0&¢er.x SeamSegments(const Mesh& mesh,double xPlane,int side) + { + std::vector segments; + for(const Triangle& triangle:mesh.triangles) + { + std::array input{{mesh.vertices[triangle.index[0]].position, + mesh.vertices[triangle.index[1]].position,mesh.vertices[triangle.index[2]].position}}; + std::vector output; + for(int i=0;i<3;++i) + { + const V3 a=input[i],b=input[(i+1)%3]; + const bool ina=side<0?a.x<=xPlane+PositionQuantum:a.x>=xPlane-PositionQuantum; + const bool inb=side<0?b.x<=xPlane+PositionQuantum:b.x>=xPlane-PositionQuantum; + if(ina)output.push_back(a); + if(ina!=inb) + { + const double t=(xPlane-a.x)/(b.x-a.x); + output.push_back(Quantize(a+(b-a)*t)); + } + } + std::vector onPlane; + for(V3 p:output)if(std::abs(p.x-xPlane)<=2.0*PositionQuantum) + { + bool found=false;for(V3 q:onPlane)if(Length(p-q)PositionQuantum) + segments.push_back({onPlane[0],onPlane[1]}); + } + return segments; + } + double PointSegmentDistance(V3 p,const Segment& s) + { + const V3 d=s.b-s.a;const double den=Dot(d,d); + const double t=den>0?std::clamp(Dot(p-s.a,d)/den,0.0,1.0):0.0; + return Length(p-(s.a+d*t)); + } + struct SeamMetric{double mean=0,maximum=0;std::uint64_t probes=0,unmatched=0;}; + SeamMetric CompareSegments(const std::vector& a,const std::vector& b) + { + SeamMetric out;double total=0; + const auto probe=[&](const std::vector& source,const std::vector& target) + { + for(const Segment& s:source)for(double t:{0.0,0.25,0.5,0.75,1.0}) + { + const V3 p=s.a+(s.b-s.a)*t;double nearest=std::numeric_limits::infinity(); + for(const Segment& candidate:target)nearest=std::min(nearest,PointSegmentDistance(p,candidate)); + if(!std::isfinite(nearest)){++out.unmatched;++out.probes;continue;} + total+=nearest;out.maximum=std::max(out.maximum,nearest);++out.probes; + } + }; + probe(a,b);probe(b,a); + if(out.unmatched){out.mean=-1.0;out.maximum=-1.0;} + else if(out.probes)out.mean=total/out.probes; + return out; + } + bool HasIsoContourOnX(const Field& field,double xPlane,double low,double high,double step) + { + const int cells=static_cast(std::llround((high-low)/step)); + for(int y=0;y(result.vertices.size()); + result.vertices.insert(result.vertices.end(),b.vertices.begin(),b.vertices.end()); + for(Triangle triangle:b.triangles) + { + for(std::uint32_t& index:triangle.index)index+=offset; + result.triangles.push_back(triangle); + } + result.hash=HashMesh(result);return result; + } + + bool WritePhaseCsv(const std::filesystem::path& path,int phaseCount) + { + std::ofstream csv(path,std::ios::binary); + csv<<"piece,method,phase_x_m,phase_y_m,phase_z_m,mean_abs_field_residual_m,max_abs_field_residual_m,probe_count,triangles,boundary_edges,nonmanifold_edges,crease_edge_pairs,normal_split_pairs,normal_split_rate,hash\n"; + const std::array pieces{{Piece::Corner,Piece::Doorway,Piece::Step,Piece::Ledge,Piece::Organic}}; + const std::array methods{{Method::MC,Method::Hybrid,Method::DC}}; + const std::array phase{{0.0,0.05,0.10,0.15,0.20}}; + constexpr int n=32; + std::uint64_t completed=0; + for(Piece piece:pieces) + { + for(int px=0;px(topology.splitPairs)/topology.creaseEdgePairs:0.0; + csv< pieces{{Piece::Corner,Piece::Doorway,Piece::Step,Piece::Ledge,Piece::Organic}}; + const std::array methods{{Method::MC,Method::Hybrid,Method::DC}}; + const std::array workerCounts{{1,8,16}}; + const V3 phase{0.10,0.10,0.10};constexpr int n=32; + for(Piece piece:pieces) + { + const Field field(piece,phase); + Grid grid=SampleGrid(field,n,n,n,1,{-4.0,-4.0,-4.0},8); + for(Method method:methods) + { + std::uint64_t reference=0; + for(int workers:workerCounts) + { + std::vector times;Mesh last; + for(int repeat=0;repeat<3;++repeat) + { + const auto start=std::chrono::steady_clock::now(); + last=BuildMesh(grid,field,method,workers); + const auto finish=std::chrono::steady_clock::now(); + times.push_back(std::chrono::duration(finish-start).count()); + } + std::sort(times.begin(),times.end()); + if(workers==1)reference=last.hash; + const bool same=last.hash==reference; + det< pieces{{Piece::Corner,Piece::Doorway,Piece::Step,Piece::Ledge,Piece::Organic}}; + const std::array methods{{Method::MC,Method::Hybrid,Method::DC}}; + const std::array phaseValues{{0.0,0.05,0.10,0.15,0.20}}; + std::uint64_t completed=0; + for(Piece piece:pieces) + { + for(double px:phaseValues)for(double py:phaseValues)for(double pz:phaseValues) + { + const V3 phase{px,py,pz}; + const Field field(piece,phase); + const bool seamExpected=HasIsoContourOnX(field,0.0,-4.0,4.0,VoxelMeters); + const Grid tileLeft=SampleGrid(field,32,32,32,1,{-8.0,-4.0,-4.0},8); + const Grid tileRight=SampleGrid(field,32,32,32,1,{0.0,-4.0,-4.0},8); + const Grid lodFineLeft=SampleGrid(field,16,32,32,1,{-4.0,-4.0,-4.0},8); + const Grid lodCoarseRight=SampleGrid(field,8,16,16,2,{0.0,-4.0,-4.0},8); + for(Method method:methods) + { + const Mesh leftMesh=BuildMesh(tileLeft,field,method,8); + const Mesh rightMesh=BuildMesh(tileRight,field,method,8); + const Mesh lodFineMesh=BuildMesh(lodFineLeft,field,method,8); + const Mesh lodCoarseMesh=BuildMesh(lodCoarseRight,field,method,8); + const auto leftSegments=SeamSegments(leftMesh,0.0,-1); + const auto rightSegments=SeamSegments(rightMesh,0.0,1); + const auto lodFineSegments=SeamSegments(lodFineMesh,0.0,-1); + const auto lodCoarseSegments=SeamSegments(lodCoarseMesh,0.0,1); + SeamMetric tile=CompareSegments(leftSegments,rightSegments); + SeamMetric lod=CompareSegments(lodFineSegments,lodCoarseSegments); + if(seamExpected&&leftSegments.empty()&&rightSegments.empty()) + {tile.unmatched=1;tile.mean=-1;tile.maximum=-1;} + if(seamExpected&&lodFineSegments.empty()&&lodCoarseSegments.empty()) + {lod.unmatched=1;lod.mean=-1;lod.maximum=-1;} + const TopologyMetrics tileTopology=AnalyzeTopology(MergeMeshes(leftMesh,rightMesh)); + const TopologyMetrics lodTopology=AnalyzeTopology(MergeMeshes(lodFineMesh,lodCoarseMesh)); + csv< pieces{{Piece::Corner,Piece::Doorway,Piece::Step,Piece::Ledge,Piece::Organic}}; + const std::array methods{{Method::MC,Method::Hybrid,Method::DC}}; + const V3 phase{0.10,0.10,0.10};constexpr int n=32; + for(Piece piece:pieces) + { + const Field field(piece,phase);Grid grid=SampleGrid(field,n,n,n,1,{-4.0,-4.0,-4.0},8); + for(Method method:methods) + { + const Mesh mesh=BuildMesh(grid,field,method,8); + WriteObj(outDir/"meshes"/(std::string(PieceName(piece))+"_"+MethodName(method)+".obj"),mesh); + } + } + return true; + } +} + +int main(int argc,char** argv) +{ + try + { + std::filesystem::path outDir="Standalone/Research/WorldAssemblerM2/results"; + int phaseCount=5; + for(int i=1;i(std::chrono::steady_clock::now()-start).count(); + std::cout<<"[SharpMesherSpike] complete seconds="<