diff --git a/README.md b/README.md index 1228736..29eb115 100644 --- a/README.md +++ b/README.md @@ -30,6 +30,12 @@ shader semantics. - `Voronoi.SampleTones`, a C# batch surface that calls the Rust `cultmath-core` native kernel when `cultmath_core` is available and falls back to the managed parity path otherwise. +- `PlanetaryField`, a scale-aware CPU/HLSL spherical field contract with + deterministic identity, cube-sphere topology, residual pages, LOD residency, + ray/path/region queries, coarse patches, and projected map-tile baking. See + `docs/planetary.md`. +- `CultMath.Unity`, an engine-facing assembly containing only Unity mesh and + page-upload adapters for the canonical planetary contracts. ## Native Core diff --git a/docs/design.md b/docs/design.md index 8250ba7..f4ca5b4 100644 --- a/docs/design.md +++ b/docs/design.md @@ -19,6 +19,12 @@ Rust hot kernels, C# host code, and HLSL shader code. - Keep the surface boring and exact. Noise primitives are allowed when they are source-grounded, shader-shaped, and shared by more than one project. Matrices, quaternions, and packing helpers still need an owner before they enter. +- Planetary geometry is a scale-aware field query, not a renderer asset. CPU + simulation is the reference consumer; pages, patches, and maps are derived + lowerings carrying field identity and unresolved error. +- Application base fields and gameplay meaning enter through narrow inputs. + CultMath does not own biomes, resources, traversability, render materials, or + committed game consequences. ## Non-Goals diff --git a/docs/planetary.md b/docs/planetary.md new file mode 100644 index 0000000..8196c1c --- /dev/null +++ b/docs/planetary.md @@ -0,0 +1,117 @@ +# Planetary Field Runtime + +CultMath owns the deterministic, renderer-independent planetary field. The same +definition can be queried by a CPU simulation, baked into residual pages, +sampled through a map projection, or evaluated by HLSL. Render resources and +gameplay interpretation remain consumer responsibilities. + +## Authority + +- `PlanetaryFieldDefinition` identifies radius, seed, erosion parameters, and a + version that includes the consumer-owned base-field version. +- `IPlanetaryBaseField` supplies deterministic application-authored base + evidence without making CultMath depend on an application. +- `PlanetaryQueryScale` makes the physical footprint and accepted unresolved + error explicit. There is deliberately no scale-free authoritative + `HeightAt` call. +- `PlanetarySurfaceSample` carries displacement, tangent gradient, normal, + slope, ridge/gully evidence, resolved wavelength, and unresolved bound. +- pages, residency, patches, and map tiles are disposable lowerings of that + contract. + +The C# surface targets .NET and Unity-compatible `netstandard2.1`. The matching +HLSL lives in `shaders/Planetary.hlsl` and is included by `CultMath.hlsl`. + +## Point and batch queries + +```csharp +var erosion = AdvancedErosionParameters.Default with +{ + Scale = radius * 0.075f, + Strength = 0.12f, + Octaves = 7, +}; +var field = PlanetaryFieldDefinition.Create( + baseFieldVersion: 17, + radius, + seed: 42, + erosion); +var scale = PlanetaryQueryScale.AtFootprint( + footprintMeters: 10, + maximumUnresolvedHeight: 0.5f); + +PlanetaryBaseFieldSample basis = source.Sample(direction); +PlanetarySurfaceSample sample = PlanetaryField.Sample( + field, direction, basis, scale); +``` + +`PlanetaryField.SampleBatch` writes into a caller-owned span. `SamplePosition` +accepts `double3`, reduces the astronomical position to a stable unit direction, +and then evaluates the same float field contract. + +Higher-level CPU queries include bounded ray intersection, segment clearance, +great-circle profiles, and region summaries. Their answers remain field queries; +only a consuming simulation can turn them into committed game consequences. + +## Cube-sphere topology and pages + +`PlanetaryTopology` owns the tangent/QSC cube-sphere mapping, inverse face +coordinates, tile selection, local page coordinates, and radial-graph surface +normal. `PlanetaryPageSampling` owns bordered directions, nominal spacing, +bilinear sampling, and summaries. + +`PlanetaryPageBaker` evaluates a root page at its physical sample spacing. A +child stores `child filtered field - parent filtered field`. A +`PlanetaryPageSetSampling` query adds every containing ancestor contribution +and applies transition weights only to residuals. + +`PlanetaryLodSelector` selects six roots plus one camera-facing ancestor chain +from the unresolved-height error. `PlanetaryResidualResidency` owns arrival and +departure presentation weights, content and presentation versions, and nothing +about field truth. + +## Patches + +`PlanetaryPatch.CreateFace` and `CreateCubeSphere` produce renderer-neutral +coarse meshes from the canonical topology. HLSL calls +`cultmath_planetary_face_direction` for procedural vertices and +`cultmath_planetary_radial_refinement_step` for bounded pixel refinement. + +Renderers still own cameras, clip space, depth, materials, buffers, dispatch, +and draw submission. + +The optional `CultMath.Unity` assembly converts a canonical face patch into a +Unity `Mesh` and converts a residency snapshot into Unity-native page input and +metadata payloads. It owns no field parameters, page selection, or simulation +state. Unity render code remains responsible for allocating/binding graphics +buffers and dispatching application-specific base-field entry points. + +## Maps + +`PlanetaryProjection` supplies forward and inverse transforms for: + +- equirectangular; +- Web Mercator with an explicit polar cutoff; +- Equal Earth; +- orthographic; +- azimuthal equidistant and equal-area; +- cube atlas; +- local tangent/gnomonic views. + +`PlanetaryMapTileBaker` maps each valid tile sample through the inverse +projection and canonical field query. `PlanetaryMapTileKey` includes field, +projection, layer, layout, and query-scale identity. Palettes, contours, labels, +and strategic overlays do not enter CultMath. + +## Numeric contract + +Integer hashes and published field versions are exact. Basic arithmetic follows +CultMath's matched C#/HLSL semantics. Hardware transcendental functions such as +`atan` and `tan` use measured tolerances rather than fictional bit identity. +Every filtered terrain query carries an unresolved-height bound so consumers can +reject insufficient fidelity explicitly. + +Current D3D12 parity evidence covers the advanced erosion kernel and the QSC +forward/inverse topology. The remaining production page and patch shaders use +the common functions and retain their seam, residual, summary, and radial-hit +integration probes in Fensalir. diff --git a/shaders/CultMath.hlsl b/shaders/CultMath.hlsl index e7281c9..cf76c8c 100644 --- a/shaders/CultMath.hlsl +++ b/shaders/CultMath.hlsl @@ -264,6 +264,7 @@ float cultmath_value_noise_texture_bicubic(Texture2D noise_texture, Sampl } #include "AdvancedErosionFilter.hlsl" +#include "Planetary.hlsl" struct CultMathSphericalErosionParameters { diff --git a/shaders/Planetary.hlsl b/shaders/Planetary.hlsl new file mode 100644 index 0000000..729c633 --- /dev/null +++ b/shaders/Planetary.hlsl @@ -0,0 +1,260 @@ +#ifndef CULTMATH_PLANETARY_HLSL +#define CULTMATH_PLANETARY_HLSL + +static const float CULTMATH_PLANETARY_QUARTER_PI = 0.78539816339744830962; + +float3 cultmath_planetary_face_direction(int face, float2 coordinate) +{ + float2 tangent = tan(coordinate * CULTMATH_PLANETARY_QUARTER_PI); + float3 cube = face == 0 ? float3(1.0, tangent.y, -tangent.x) + : face == 1 ? float3(-1.0, tangent.y, tangent.x) + : face == 2 ? float3(tangent.x, 1.0, -tangent.y) + : face == 3 ? float3(tangent.x, -1.0, tangent.y) + : face == 4 ? float3(tangent.x, tangent.y, 1.0) + : float3(-tangent.x, tangent.y, -1.0); + return cultmath_normalize(cube); +} + +float2 cultmath_planetary_face_coordinate(float3 direction, out int face) +{ + float3 d = cultmath_normalize(direction); + float3 a = abs(d); + float2 cube_uv; + if (a.x >= a.y && a.x >= a.z) + { + if (d.x >= 0.0) { face = 0; cube_uv = float2(-d.z, d.y) / a.x; } + else { face = 1; cube_uv = float2(d.z, d.y) / a.x; } + } + else if (a.y >= a.z) + { + if (d.y >= 0.0) { face = 2; cube_uv = float2(d.x, -d.z) / a.y; } + else { face = 3; cube_uv = float2(d.x, d.z) / a.y; } + } + else + { + if (d.z >= 0.0) { face = 4; cube_uv = d.xy / a.z; } + else { face = 5; cube_uv = float2(-d.x, d.y) / a.z; } + } + return atan(cube_uv) / CULTMATH_PLANETARY_QUARTER_PI; +} + +bool cultmath_planetary_page_local(float3 direction, float4 tile_address, float active, out float2 local) +{ + int face; + float2 uv = cultmath_planetary_face_coordinate(direction, face); + local = 0.0; + if (active < 0.5 || abs((float)face - tile_address.x) > 0.25) return false; + float axis_tiles = exp2(tile_address.y); + float2 scaled = (uv * 0.5 + 0.5) * axis_tiles; + float2 canonical = floor(min(scaled, axis_tiles - 1.0e-5)); + if (any(abs(canonical - tile_address.zw) > 0.25)) return false; + local = scaled - tile_address.zw; + return all(local >= 0.0) && all(local <= 1.0); +} + +float3 cultmath_planetary_surface_normal(float3 direction, float3 world_distance_tangent_gradient) +{ + direction = cultmath_normalize(direction); + float3 tangent_gradient = world_distance_tangent_gradient - direction * dot(world_distance_tangent_gradient, direction); + return cultmath_normalize(direction - tangent_gradient); +} + +float2 cultmath_planetary_equirectangular_forward(float3 direction) +{ + direction = cultmath_normalize(direction); + return float2(atan2(direction.y, direction.x) / CULTMATH_PI, asin(clamp(direction.z, -1.0, 1.0)) / CULTMATH_HALF_PI); +} + +float3 cultmath_planetary_equirectangular_inverse(float2 coordinate) +{ + float longitude = coordinate.x * CULTMATH_PI; + float latitude = coordinate.y * CULTMATH_HALF_PI; + float cos_latitude = cos(latitude); + return cultmath_normalize(float3(cos_latitude * cos(longitude), cos_latitude * sin(longitude), sin(latitude))); +} + +static const float CULTMATH_EQUAL_EARTH_A1 = 1.340264; +static const float CULTMATH_EQUAL_EARTH_A2 = -0.081106; +static const float CULTMATH_EQUAL_EARTH_A3 = 0.000893; +static const float CULTMATH_EQUAL_EARTH_A4 = 0.003796; +static const float CULTMATH_EQUAL_EARTH_M = 0.8660254037844386; +static const float CULTMATH_EQUAL_EARTH_X_MAX = 2.706629; +static const float CULTMATH_EQUAL_EARTH_Y_MAX = 1.3173628; + +float cultmath_equal_earth_y(float theta) +{ + float theta2 = theta * theta; + float theta6 = theta2 * theta2 * theta2; + return theta * (CULTMATH_EQUAL_EARTH_A1 + CULTMATH_EQUAL_EARTH_A2 * theta2 + theta6 * (CULTMATH_EQUAL_EARTH_A3 + CULTMATH_EQUAL_EARTH_A4 * theta2)); +} + +float2 cultmath_planetary_equal_earth_forward(float3 direction) +{ + direction = cultmath_normalize(direction); + float longitude = atan2(direction.y, direction.x); + float latitude = asin(clamp(direction.z, -1.0, 1.0)); + float theta = asin(CULTMATH_EQUAL_EARTH_M * sin(latitude)); + float theta2 = theta * theta; + float theta6 = theta2 * theta2 * theta2; + float denominator = CULTMATH_EQUAL_EARTH_M * (CULTMATH_EQUAL_EARTH_A1 + 3.0 * CULTMATH_EQUAL_EARTH_A2 * theta2 + theta6 * (7.0 * CULTMATH_EQUAL_EARTH_A3 + 9.0 * CULTMATH_EQUAL_EARTH_A4 * theta2)); + return float2(longitude * cos(theta) / denominator / CULTMATH_EQUAL_EARTH_X_MAX, cultmath_equal_earth_y(theta) / CULTMATH_EQUAL_EARTH_Y_MAX); +} + +float3 cultmath_planetary_equal_earth_inverse(float2 coordinate) +{ + float target_y = coordinate.y * CULTMATH_EQUAL_EARTH_Y_MAX; + float theta = target_y / CULTMATH_EQUAL_EARTH_A1; + [unroll] for (int iteration = 0; iteration < 8; iteration++) + { + float theta2 = theta * theta; + float theta6 = theta2 * theta2 * theta2; + float derivative = CULTMATH_EQUAL_EARTH_A1 + 3.0 * CULTMATH_EQUAL_EARTH_A2 * theta2 + theta6 * (7.0 * CULTMATH_EQUAL_EARTH_A3 + 9.0 * CULTMATH_EQUAL_EARTH_A4 * theta2); + theta -= (cultmath_equal_earth_y(theta) - target_y) / derivative; + } + float latitude = asin(clamp(sin(theta) / CULTMATH_EQUAL_EARTH_M, -1.0, 1.0)); + float theta2 = theta * theta; + float theta6 = theta2 * theta2 * theta2; + float denominator = CULTMATH_EQUAL_EARTH_M * (CULTMATH_EQUAL_EARTH_A1 + 3.0 * CULTMATH_EQUAL_EARTH_A2 * theta2 + theta6 * (7.0 * CULTMATH_EQUAL_EARTH_A3 + 9.0 * CULTMATH_EQUAL_EARTH_A4 * theta2)); + float longitude = coordinate.x * CULTMATH_EQUAL_EARTH_X_MAX * denominator / cos(theta); + float cos_latitude = cos(latitude); + return cultmath_normalize(float3(cos_latitude * cos(longitude), cos_latitude * sin(longitude), sin(latitude))); +} + +struct CultMathPlanetaryFieldDefinition +{ + float radius; + int seed; + CultMathAdvancedErosionParameters erosion; +}; + +struct CultMathPlanetaryBaseFieldSample +{ + float radial_displacement; + float3 radial_gradient; + float field_value; + float3 field_gradient; + float fade_target; +}; + +struct CultMathPlanetarySurfaceSample +{ + float radial_displacement; + float3 tangent_gradient; + float3 surface_normal; + float slope; + float ridge; + float gully; + float finest_resolved_wavelength; + float unresolved_height_bound; +}; + +struct CultMathPlanetaryPageSample +{ + float4 height_gradient; + float2 masks; +}; + +CultMathPlanetaryPageSample cultmath_planetary_page_lerp( + CultMathPlanetaryPageSample a, + CultMathPlanetaryPageSample b, + CultMathPlanetaryPageSample c, + CultMathPlanetaryPageSample d, + float2 fraction) +{ + CultMathPlanetaryPageSample result; + result.height_gradient = lerp(lerp(a.height_gradient, b.height_gradient, fraction.x), lerp(c.height_gradient, d.height_gradient, fraction.x), fraction.y); + result.masks = lerp(lerp(a.masks, b.masks, fraction.x), lerp(c.masks, d.masks, fraction.x), fraction.y); + return result; +} + +CultMathPlanetaryPageSample cultmath_planetary_residual_sample( + CultMathPlanetarySurfaceSample child, + CultMathPlanetarySurfaceSample parent, + bool has_parent) +{ + CultMathPlanetaryPageSample result; + result.height_gradient = float4(child.radial_displacement, child.tangent_gradient); + if (has_parent) result.height_gradient -= float4(parent.radial_displacement, parent.tangent_gradient); + result.masks = float2(child.ridge, child.gully); + return result; +} + +float4 cultmath_planetary_summary_empty() +{ + return float4(3.402823466e+38, -3.402823466e+38, 0.0, 0.0); +} + +float4 cultmath_planetary_summary_accumulate(float4 summary, float4 height_gradient, float unresolved_height) +{ + summary.x = min(summary.x, height_gradient.x); + summary.y = max(summary.y, height_gradient.x); + summary.z = max(summary.z, length(height_gradient.yzw)); + summary.w = max(summary.w, max(unresolved_height, 0.0)); + return summary; +} + +float cultmath_planetary_pow4(float value) { float square = value * value; return square * square; } + +CultMathPlanetarySurfaceSample cultmath_planetary_field_sample( + CultMathPlanetaryFieldDefinition definition, + float3 unit_direction, + CultMathPlanetaryBaseFieldSample base_sample, + float footprint_meters) +{ + float3 direction = cultmath_normalize(unit_direction); + float3 radial_gradient = base_sample.radial_gradient - direction * dot(base_sample.radial_gradient, direction); + float3 field_gradient = base_sample.field_gradient - direction * dot(base_sample.field_gradient, direction); + CultMathAdvancedErosionParameters p = definition.erosion; + CultMathErosionBandSelection band = cultmath_select_erosion_bands( + p.scale * p.cell_scale, footprint_meters, p.octaves, p.lacunarity, + p.strength * p.scale, p.gain, 2.0); + float3 powers = float3( + cultmath_planetary_pow4(abs(direction.x)), + cultmath_planetary_pow4(abs(direction.y)), + cultmath_planetary_pow4(abs(direction.z))); + float power_sum = max(powers.x + powers.y + powers.z, 1.0e-6); + float3 weights = powers / power_sum; + float3 world = direction * definition.radius; + float seed = (float)definition.seed; + CultMathAdvancedErosionResult xy = cultmath_advanced_erosion_filter_banded( + world.xy + float2(713.0 + seed * 0.754877666, -291.0 + seed * 0.569840296), + float3(base_sample.field_value, field_gradient.x, field_gradient.y), base_sample.fade_target, p, band); + CultMathAdvancedErosionResult yz = cultmath_advanced_erosion_filter_banded( + world.yz + float2(-431.0 + seed * 0.438289027, 887.0 + seed * 0.328438163), + float3(base_sample.field_value, field_gradient.y, field_gradient.z), base_sample.fade_target, p, band); + CultMathAdvancedErosionResult zx = cultmath_advanced_erosion_filter_banded( + world.zx + float2(197.0 + seed * 0.219783071, 557.0 + seed * 0.145898034), + float3(base_sample.field_value, field_gradient.z, field_gradient.x), base_sample.fade_target, p, band); + float erosion_height = xy.delta.x * weights.z + yz.delta.x * weights.x + zx.delta.x * weights.y; + float3 erosion_gradient = float3( + xy.delta.y * weights.z + zx.delta.z * weights.y, + xy.delta.z * weights.z + yz.delta.y * weights.x, + yz.delta.z * weights.x + zx.delta.y * weights.y); + float3 power_gradient_x = float3(4.0 * direction.x * direction.x * direction.x / definition.radius, 0.0, 0.0); + float3 power_gradient_y = float3(0.0, 4.0 * direction.y * direction.y * direction.y / definition.radius, 0.0); + float3 power_gradient_z = float3(0.0, 0.0, 4.0 * direction.z * direction.z * direction.z / definition.radius); + power_gradient_x -= direction * dot(power_gradient_x, direction); + power_gradient_y -= direction * dot(power_gradient_y, direction); + power_gradient_z -= direction * dot(power_gradient_z, direction); + float3 power_sum_gradient = power_gradient_x + power_gradient_y + power_gradient_z; + float inverse_power_sum_squared = 1.0 / (power_sum * power_sum); + float3 weight_gradient_x = (power_gradient_x * power_sum - power_sum_gradient * powers.x) * inverse_power_sum_squared; + float3 weight_gradient_y = (power_gradient_y * power_sum - power_sum_gradient * powers.y) * inverse_power_sum_squared; + float3 weight_gradient_z = (power_gradient_z * power_sum - power_sum_gradient * powers.z) * inverse_power_sum_squared; + erosion_gradient += yz.delta.x * weight_gradient_x + zx.delta.x * weight_gradient_y + xy.delta.x * weight_gradient_z; + float3 gradient = radial_gradient + erosion_gradient; + gradient -= direction * dot(gradient, direction); + float ridge_evidence = xy.ridge_map * weights.z + yz.ridge_map * weights.x + zx.ridge_map * weights.y; + CultMathPlanetarySurfaceSample result; + result.radial_displacement = base_sample.radial_displacement + erosion_height; + result.tangent_gradient = gradient; + result.surface_normal = cultmath_planetary_surface_normal(direction, gradient); + result.slope = length(gradient); + result.ridge = saturate(ridge_evidence * 0.5 + 0.5); + result.gully = saturate(0.5 - ridge_evidence * 0.5); + result.finest_resolved_wavelength = band.finest_included_wavelength; + result.unresolved_height_bound = band.unresolved_height_bound; + return result; +} + +#endif diff --git a/src/CultMath/CultMath.csproj b/src/CultMath/CultMath.csproj index 01901b6..61b57c9 100644 --- a/src/CultMath/CultMath.csproj +++ b/src/CultMath/CultMath.csproj @@ -19,6 +19,8 @@ + + diff --git a/src/CultMath/PlanetaryField.cs b/src/CultMath/PlanetaryField.cs new file mode 100644 index 0000000..13f36c0 --- /dev/null +++ b/src/CultMath/PlanetaryField.cs @@ -0,0 +1,213 @@ +namespace CultMath; + +public readonly record struct PlanetaryQueryScale(float FootprintMeters, float MaximumUnresolvedHeight) +{ + public static PlanetaryQueryScale AtFootprint(float footprintMeters, float maximumUnresolvedHeight = float.PositiveInfinity) + => new PlanetaryQueryScale(footprintMeters, maximumUnresolvedHeight).Validate(); + + public PlanetaryQueryScale Validate() + { + if (!float.IsFinite(FootprintMeters) || FootprintMeters < 0) throw new ArgumentOutOfRangeException(nameof(FootprintMeters)); + if (float.IsNaN(MaximumUnresolvedHeight) || MaximumUnresolvedHeight < 0) throw new ArgumentOutOfRangeException(nameof(MaximumUnresolvedHeight)); + return this; + } +} + +public readonly record struct PlanetaryFieldDefinition( + ulong FieldVersion, + float Radius, + int Seed, + AdvancedErosionParameters Erosion) +{ + public static PlanetaryFieldDefinition Create( + ulong baseFieldVersion, + float radius, + int seed, + AdvancedErosionParameters erosion) + => new PlanetaryFieldDefinition(PlanetaryFieldIdentity.Compute(baseFieldVersion, radius, seed, erosion), radius, seed, erosion).Validate(); + + public PlanetaryFieldDefinition Validate() + { + if (FieldVersion == 0) throw new ArgumentOutOfRangeException(nameof(FieldVersion)); + if (!float.IsFinite(Radius) || Radius <= 0) throw new ArgumentOutOfRangeException(nameof(Radius)); + if (!float.IsFinite(Erosion.Scale) || Erosion.Scale <= 0) throw new ArgumentOutOfRangeException(nameof(Erosion)); + return this; + } +} + +public static class PlanetaryFieldIdentity +{ + public static ulong Compute(ulong baseFieldVersion, float radius, int seed, in AdvancedErosionParameters erosion) + { + if (baseFieldVersion == 0) throw new ArgumentOutOfRangeException(nameof(baseFieldVersion)); + const ulong offset = 14695981039346656037UL; + const ulong prime = 1099511628211UL; + var hash = offset; + Add(baseFieldVersion); + Add((uint)BitConverter.SingleToInt32Bits(radius)); Add((uint)seed); + Add((uint)BitConverter.SingleToInt32Bits(erosion.Scale)); Add((uint)BitConverter.SingleToInt32Bits(erosion.Strength)); + Add((uint)BitConverter.SingleToInt32Bits(erosion.GullyWeight)); Add((uint)BitConverter.SingleToInt32Bits(erosion.Detail)); + AddFloat4(erosion.Rounding); AddFloat4(erosion.Onset); + Add((uint)BitConverter.SingleToInt32Bits(erosion.AssumedSlope.x)); Add((uint)BitConverter.SingleToInt32Bits(erosion.AssumedSlope.y)); + Add((uint)BitConverter.SingleToInt32Bits(erosion.CellScale)); Add((uint)BitConverter.SingleToInt32Bits(erosion.Normalization)); + Add((uint)erosion.Octaves); Add((uint)BitConverter.SingleToInt32Bits(erosion.Lacunarity)); Add((uint)BitConverter.SingleToInt32Bits(erosion.Gain)); + return hash == 0 ? 1 : hash; + + void Add(ulong value) => hash = unchecked((hash ^ value) * prime); + void AddFloat4(float4 value) + { + Add((uint)BitConverter.SingleToInt32Bits(value.x)); Add((uint)BitConverter.SingleToInt32Bits(value.y)); + Add((uint)BitConverter.SingleToInt32Bits(value.z)); Add((uint)BitConverter.SingleToInt32Bits(value.w)); + } + } +} + +/// +/// Application-supplied base evidence at one unit spherical direction. +/// FieldValue and FieldGradient drive erosion. RadialDisplacement and +/// RadialGradient describe the uneroded surface in metres. +/// +public readonly record struct PlanetaryBaseFieldSample( + float RadialDisplacement, + float3 RadialGradient, + float FieldValue, + float3 FieldGradient, + float FadeTarget); + +public readonly record struct PlanetarySurfaceSample( + ulong FieldVersion, + float3 UnitDirection, + float Radius, + float RadialDisplacement, + float3 TangentGradient, + float3 SurfaceNormal, + float Slope, + float Ridge, + float Gully, + float FinestResolvedWavelength, + float UnresolvedHeightBound) +{ + public float SurfaceRadius => Radius + RadialDisplacement; + public float3 PlanetLocalPosition => UnitDirection * SurfaceRadius; + public bool Meets(PlanetaryQueryScale scale) => UnresolvedHeightBound <= scale.MaximumUnresolvedHeight; +} + +public static class PlanetaryField +{ + public static PlanetarySurfaceSample Sample( + in PlanetaryFieldDefinition definition, + float3 unitDirection, + in PlanetaryBaseFieldSample baseSample, + in PlanetaryQueryScale scale) + { + definition.Validate(); + scale.Validate(); + PlanetaryTopology.ValidateDirection(unitDirection); + PlanetaryTopology.ValidateFinite(baseSample.RadialGradient, nameof(baseSample)); + PlanetaryTopology.ValidateFinite(baseSample.FieldGradient, nameof(baseSample)); + if (!float.IsFinite(baseSample.RadialDisplacement) || !float.IsFinite(baseSample.FieldValue) || !float.IsFinite(baseSample.FadeTarget)) + throw new ArgumentOutOfRangeException(nameof(baseSample)); + + var direction = math.normalize(unitDirection); + var radialGradient = Tangent(direction, baseSample.RadialGradient); + var fieldGradient = Tangent(direction, baseSample.FieldGradient); + var p = definition.Erosion; + var band = ErosionFrequencyBands.Select( + p.Scale * p.CellScale, + scale.FootprintMeters, + p.Octaves, + p.Lacunarity, + p.Strength * p.Scale, + p.Gain); + + var powers = new float3( + Pow4(MathF.Abs(direction.x)), + Pow4(MathF.Abs(direction.y)), + Pow4(MathF.Abs(direction.z))); + var powerSum = math.max(powers.x + powers.y + powers.z, 1.0e-6f); + var weights = powers / powerSum; + var world = direction * definition.Radius; + var seed = definition.Seed; + var xy = AdvancedErosionFilter.Sample( + world.xy + new float2(713 + seed * 0.754877666f, -291 + seed * 0.569840296f), + new float3(baseSample.FieldValue, fieldGradient.x, fieldGradient.y), baseSample.FadeTarget, p, band); + var yz = AdvancedErosionFilter.Sample( + world.yz + new float2(-431 + seed * 0.438289027f, 887 + seed * 0.328438163f), + new float3(baseSample.FieldValue, fieldGradient.y, fieldGradient.z), baseSample.FadeTarget, p, band); + var zx = AdvancedErosionFilter.Sample( + new float2(world.z, world.x) + new float2(197 + seed * 0.219783071f, 557 + seed * 0.145898034f), + new float3(baseSample.FieldValue, fieldGradient.z, fieldGradient.x), baseSample.FadeTarget, p, band); + + var erosionHeight = xy.Delta.x * weights.z + yz.Delta.x * weights.x + zx.Delta.x * weights.y; + var erosionGradient = new float3( + xy.Delta.y * weights.z + zx.Delta.z * weights.y, + xy.Delta.z * weights.z + yz.Delta.y * weights.x, + yz.Delta.z * weights.x + zx.Delta.y * weights.y); + var powerGradientX = Tangent(direction, new float3(4 * direction.x * direction.x * direction.x / definition.Radius, 0, 0)); + var powerGradientY = Tangent(direction, new float3(0, 4 * direction.y * direction.y * direction.y / definition.Radius, 0)); + var powerGradientZ = Tangent(direction, new float3(0, 0, 4 * direction.z * direction.z * direction.z / definition.Radius)); + var sumGradient = powerGradientX + powerGradientY + powerGradientZ; + var weightGradientX = (powerGradientX * powerSum - sumGradient * powers.x) / (powerSum * powerSum); + var weightGradientY = (powerGradientY * powerSum - sumGradient * powers.y) / (powerSum * powerSum); + var weightGradientZ = (powerGradientZ * powerSum - sumGradient * powers.z) / (powerSum * powerSum); + erosionGradient += yz.Delta.x * weightGradientX + zx.Delta.x * weightGradientY + xy.Delta.x * weightGradientZ; + var gradient = Tangent(direction, radialGradient + erosionGradient); + var ridgeEvidence = xy.RidgeMap * weights.z + yz.RidgeMap * weights.x + zx.RidgeMap * weights.y; + var ridge = math.saturate(ridgeEvidence * 0.5f + 0.5f); + var gully = math.saturate(0.5f - ridgeEvidence * 0.5f); + var displacement = baseSample.RadialDisplacement + erosionHeight; + + return new( + definition.FieldVersion, + direction, + definition.Radius, + displacement, + gradient, + PlanetaryTopology.SurfaceNormal(direction, gradient), + math.length(gradient), + ridge, + gully, + band.FinestIncludedWavelength, + band.UnresolvedHeightBound); + } + + public static void SampleBatch( + in PlanetaryFieldDefinition definition, + ReadOnlySpan directions, + in PlanetaryQueryScale scale, + TSource source, + Span destination) + where TSource : IPlanetaryBaseField + { + if (destination.Length < directions.Length) throw new ArgumentException("Destination is shorter than directions.", nameof(destination)); + for (var i = 0; i < directions.Length; i++) + { + var direction = math.normalize(directions[i]); + var baseSample = source.Sample(direction); + destination[i] = Sample(definition, direction, baseSample, scale); + } + } + + public static PlanetarySurfaceSample SamplePosition( + in PlanetaryFieldDefinition definition, + double3 planetLocalPosition, + in PlanetaryQueryScale scale, + TSource source) + where TSource : IPlanetaryBaseField + { + if (!double.IsFinite(planetLocalPosition.x) || !double.IsFinite(planetLocalPosition.y) || !double.IsFinite(planetLocalPosition.z)) + throw new ArgumentOutOfRangeException(nameof(planetLocalPosition)); + var length = Math.Sqrt(planetLocalPosition.x * planetLocalPosition.x + planetLocalPosition.y * planetLocalPosition.y + planetLocalPosition.z * planetLocalPosition.z); + if (length < 1.0e-12) throw new ArgumentOutOfRangeException(nameof(planetLocalPosition)); + var direction = new float3((float)(planetLocalPosition.x / length), (float)(planetLocalPosition.y / length), (float)(planetLocalPosition.z / length)); + return Sample(definition, direction, source.Sample(direction), scale); + } + + private static float3 Tangent(float3 direction, float3 gradient) => gradient - direction * math.dot(gradient, direction); + private static float Pow4(float value) { var square = value * value; return square * square; } +} + +public interface IPlanetaryBaseField +{ + PlanetaryBaseFieldSample Sample(float3 unitDirection); +} diff --git a/src/CultMath/PlanetaryPages.cs b/src/CultMath/PlanetaryPages.cs new file mode 100644 index 0000000..c9cff51 --- /dev/null +++ b/src/CultMath/PlanetaryPages.cs @@ -0,0 +1,185 @@ +namespace CultMath; + +public readonly record struct PlanetaryPageLayout(PlanetaryTileAddress Tile, int InteriorSize, int BorderSize) +{ + public int StorageSize => checked(InteriorSize + BorderSize * 2); + + public PlanetaryPageLayout Validate() + { + if (InteriorSize < 2) throw new ArgumentOutOfRangeException(nameof(InteriorSize)); + if (BorderSize < 0 || BorderSize >= InteriorSize) throw new ArgumentOutOfRangeException(nameof(BorderSize)); + return this; + } +} + +public readonly record struct PlanetaryPageSummary( + float MinimumHeight, + float MaximumHeight, + float MaximumSlope, + float UnresolvedHeightBound, + float AngularTexelSize); + +public readonly record struct PlanetaryPageSample(float Height, float3 TangentGradient, float Ridge, float Gully); + +public static class PlanetaryPageSampling +{ + public static float3 Direction(PlanetaryPageLayout page, int storageX, int storageY) + { + page.Validate(); + if ((uint)storageX >= (uint)page.StorageSize || (uint)storageY >= (uint)page.StorageSize) + throw new ArgumentOutOfRangeException(nameof(storageX)); + var localU = (storageX - page.BorderSize) / (double)(page.InteriorSize - 1); + var localV = (storageY - page.BorderSize) / (double)(page.InteriorSize - 1); + return DirectionAtLocal(page, localU, localV); + } + + public static float3 DirectionAtLocal(PlanetaryPageLayout page, double localU, double localV) + { + page.Validate(); + if (!double.IsFinite(localU) || !double.IsFinite(localV)) throw new ArgumentOutOfRangeException(nameof(localU)); + var count = page.Tile.AxisTileCount; + return PlanetaryTopology.Direction(new( + page.Tile.Face, + -1.0 + 2.0 * (page.Tile.X + localU) / count, + -1.0 + 2.0 * (page.Tile.Y + localV) / count)); + } + + public static float AngularTexelSize(PlanetaryPageLayout page) + { + var center = page.BorderSize + (page.InteriorSize - 1) / 2; + var a = Direction(page, center, center); + var b = Direction(page, Math.Min(center + 1, page.StorageSize - 1), center); + return MathF.Acos(math.clamp(math.dot(a, b), -1.0f, 1.0f)); + } + + public static float NominalAngularTexelSize(PlanetaryPageLayout page) + { + page.Validate(); + return MathF.PI * 0.5f / (page.Tile.AxisTileCount * (page.InteriorSize - 1)); + } + + public static PlanetaryPageSample Bilinear( + ReadOnlySpan samples, + PlanetaryPageLayout page, + float2 local) + { + page.Validate(); + if (samples.Length != page.StorageSize * page.StorageSize) throw new ArgumentException("Page sample count does not match layout.", nameof(samples)); + var texel = local * (page.InteriorSize - 1) + page.BorderSize; + var x0 = Math.Clamp((int)MathF.Floor(texel.x), 0, page.StorageSize - 1); + var y0 = Math.Clamp((int)MathF.Floor(texel.y), 0, page.StorageSize - 1); + var x1 = Math.Min(x0 + 1, page.StorageSize - 1); + var y1 = Math.Min(y0 + 1, page.StorageSize - 1); + var tx = math.frac(texel.x); + var ty = math.frac(texel.y); + var a = Lerp(samples[y0 * page.StorageSize + x0], samples[y0 * page.StorageSize + x1], tx); + var b = Lerp(samples[y1 * page.StorageSize + x0], samples[y1 * page.StorageSize + x1], tx); + return Lerp(a, b, ty); + } + + public static PlanetaryPageSummary Summarize(ReadOnlySpan samples, PlanetaryPageLayout page, float unresolvedHeightBound) + { + page.Validate(); + if (samples.Length != page.StorageSize * page.StorageSize) throw new ArgumentException("Page sample count does not match layout.", nameof(samples)); + var min = float.PositiveInfinity; + var max = float.NegativeInfinity; + var maxSlope = 0.0f; + foreach (var sample in samples) + { + min = MathF.Min(min, sample.Height); + max = MathF.Max(max, sample.Height); + maxSlope = MathF.Max(maxSlope, math.length(sample.TangentGradient)); + } + return new(min, max, maxSlope, MathF.Max(unresolvedHeightBound, 0), AngularTexelSize(page)); + } + + private static PlanetaryPageSample Lerp(PlanetaryPageSample a, PlanetaryPageSample b, float t) => new( + math.lerp(a.Height, b.Height, t), + math.lerp(a.TangentGradient, b.TangentGradient, t), + math.lerp(a.Ridge, b.Ridge, t), + math.lerp(a.Gully, b.Gully, t)); +} + +public readonly record struct PlanetaryBakedPage( + PlanetaryPageLayout Layout, + PlanetaryPageSample[] Samples, + PlanetaryPageSummary Summary); + +public static class PlanetaryPageBaker +{ + public static PlanetaryBakedPage Bake( + in PlanetaryFieldDefinition field, + PlanetaryPageLayout layout, + TSource source) + where TSource : IPlanetaryBaseField + { + layout.Validate(); field.Validate(); + var spacing = PlanetaryPageSampling.NominalAngularTexelSize(layout) * field.Radius; + var parentSpacing = layout.Tile.Level == 0 + ? 0 + : PlanetaryPageSampling.NominalAngularTexelSize(new(layout.Tile.Parent(), layout.InteriorSize, layout.BorderSize)) * field.Radius; + var samples = new PlanetaryPageSample[layout.StorageSize * layout.StorageSize]; + var unresolved = 0.0f; + for (var y = 0; y < layout.StorageSize; y++) + for (var x = 0; x < layout.StorageSize; x++) + { + var direction = PlanetaryPageSampling.Direction(layout, x, y); + var baseSample = source.Sample(direction); + var child = PlanetaryField.Sample(field, direction, baseSample, PlanetaryQueryScale.AtFootprint(spacing)); + var height = child.RadialDisplacement; + var gradient = child.TangentGradient; + if (parentSpacing > 0) + { + var parent = PlanetaryField.Sample(field, direction, baseSample, PlanetaryQueryScale.AtFootprint(parentSpacing)); + height -= parent.RadialDisplacement; + gradient -= parent.TangentGradient; + } + samples[y * layout.StorageSize + x] = new(height, gradient, child.Ridge, child.Gully); + unresolved = MathF.Max(unresolved, child.UnresolvedHeightBound); + } + return new(layout, samples, PlanetaryPageSampling.Summarize(samples, layout, unresolved)); + } +} + +public readonly record struct PlanetaryResidentPage(PlanetaryBakedPage Page, float Blend) +{ + public PlanetaryResidentPage Validate() + { + if (!float.IsFinite(Blend) || Blend is < 0 or > 1) throw new ArgumentOutOfRangeException(nameof(Blend)); + Page.Layout.Validate(); + if (Page.Samples.Length != Page.Layout.StorageSize * Page.Layout.StorageSize) throw new ArgumentException("Page payload does not match layout."); + return this; + } +} + +public static class PlanetaryPageSetSampling +{ + public static bool TrySample(ReadOnlySpan pages, float3 direction, out PlanetaryPageSample sample) + { + PlanetaryTopology.ValidateDirection(direction); + var height = 0.0f; + var gradient = float3.zero; + var ridge = 0.0f; + var gully = 0.0f; + var found = false; + var deepestLevel = -1; + foreach (var resident in pages) + { + resident.Validate(); + if (!PlanetaryTopology.TryLocalCoordinate(direction, resident.Page.Layout.Tile, out var local)) continue; + var contribution = PlanetaryPageSampling.Bilinear(resident.Page.Samples, resident.Page.Layout, local); + var blend = math.saturate(resident.Blend); + height += contribution.Height * blend; + gradient += contribution.TangentGradient * blend; + if (resident.Page.Layout.Tile.Level >= deepestLevel) + { + ridge = found ? math.lerp(ridge, contribution.Ridge, blend) : contribution.Ridge; + gully = found ? math.lerp(gully, contribution.Gully, blend) : contribution.Gully; + deepestLevel = resident.Page.Layout.Tile.Level; + } + found = true; + } + sample = new(height, gradient, ridge, gully); + return found; + } +} diff --git a/src/CultMath/PlanetaryPatch.cs b/src/CultMath/PlanetaryPatch.cs new file mode 100644 index 0000000..0fd3cbb --- /dev/null +++ b/src/CultMath/PlanetaryPatch.cs @@ -0,0 +1,42 @@ +namespace CultMath; + +public readonly record struct PlanetaryPatchVertex(float2 LocalCoordinate, float3 UnitDirection); + +public readonly record struct PlanetaryPatchMesh( + PlanetaryCubeFace Face, + int CellsPerAxis, + PlanetaryPatchVertex[] Vertices, + int[] Indices); + +public static class PlanetaryPatch +{ + public static PlanetaryPatchMesh CreateFace(PlanetaryCubeFace face, int cellsPerAxis) + { + if (!Enum.IsDefined(typeof(PlanetaryCubeFace), face)) throw new ArgumentOutOfRangeException(nameof(face)); + if (cellsPerAxis < 1 || cellsPerAxis > 4096) throw new ArgumentOutOfRangeException(nameof(cellsPerAxis)); + var axisVertices = cellsPerAxis + 1; + var vertices = new PlanetaryPatchVertex[checked(axisVertices * axisVertices)]; + for (var y = 0; y < axisVertices; y++) + for (var x = 0; x < axisVertices; x++) + { + var local = new float2(x / (float)cellsPerAxis, y / (float)cellsPerAxis); + vertices[y * axisVertices + x] = new(local, PlanetaryTopology.Direction(new(face, local.x * 2 - 1, local.y * 2 - 1))); + } + var indices = new int[checked(cellsPerAxis * cellsPerAxis * 6)]; + var write = 0; + for (var y = 0; y < cellsPerAxis; y++) + for (var x = 0; x < cellsPerAxis; x++) + { + var lowerLeft = y * axisVertices + x; + var lowerRight = lowerLeft + 1; + var upperLeft = lowerLeft + axisVertices; + var upperRight = upperLeft + 1; + indices[write++] = lowerLeft; indices[write++] = lowerRight; indices[write++] = upperRight; + indices[write++] = lowerLeft; indices[write++] = upperRight; indices[write++] = upperLeft; + } + return new(face, cellsPerAxis, vertices, indices); + } + + public static PlanetaryPatchMesh[] CreateCubeSphere(int cellsPerAxis) + => ((PlanetaryCubeFace[])Enum.GetValues(typeof(PlanetaryCubeFace))).Select(face => CreateFace(face, cellsPerAxis)).ToArray(); +} diff --git a/src/CultMath/PlanetaryProjection.cs b/src/CultMath/PlanetaryProjection.cs new file mode 100644 index 0000000..f849fe2 --- /dev/null +++ b/src/CultMath/PlanetaryProjection.cs @@ -0,0 +1,320 @@ +namespace CultMath; + +public enum PlanetaryProjectionKind +{ + Equirectangular, + WebMercator, + EqualEarth, + Orthographic, + AzimuthalEquidistant, + AzimuthalEqualArea, + CubeAtlas, + LocalTangent, +} + +public readonly record struct PlanetaryProjectionParameters( + PlanetaryProjectionKind Kind, + double CenterLongitude = 0, + double CenterLatitude = 0, + double Scale = 1) +{ + public PlanetaryProjectionParameters Validate() + { + if (!double.IsFinite(CenterLongitude) || !double.IsFinite(CenterLatitude) || !double.IsFinite(Scale) || Scale <= 0) + throw new ArgumentOutOfRangeException(nameof(Scale)); + if (CenterLatitude is < -Math.PI / 2 or > Math.PI / 2) throw new ArgumentOutOfRangeException(nameof(CenterLatitude)); + return this; + } +} + +public static class PlanetaryProjection +{ + private const double EqualEarthA1 = 1.340264; + private const double EqualEarthA2 = -0.081106; + private const double EqualEarthA3 = 0.000893; + private const double EqualEarthA4 = 0.003796; + private static readonly double EqualEarthM = Math.Sqrt(3) / 2; + private static readonly double EqualEarthXMax = Math.PI / (EqualEarthM * EqualEarthA1); + private static readonly double EqualEarthYMax = EqualEarthY(Math.Asin(EqualEarthM)); + public static readonly double WebMercatorMaximumLatitude = Math.Atan(Math.Sinh(Math.PI)); + + public static bool TryForward(float3 direction, PlanetaryProjectionParameters parameters, out double2 coordinate) + { + parameters.Validate(); + PlanetaryTopology.ValidateDirection(direction); + var d = math.normalize(direction); + if (parameters.Kind == PlanetaryProjectionKind.CubeAtlas) + { + var face = PlanetaryTopology.FaceCoordinate(d); + var index = (int)face.Face; + var column = index % 3; + var row = index / 3; + coordinate = new( + (-1 + (column + (face.U + 1) * 0.5) * (2.0 / 3.0)) / parameters.Scale, + (1 - (row + (face.V + 1) * 0.5)) / parameters.Scale); + return true; + } + + var longitude = Math.Atan2(d.y, d.x); + var latitude = Math.Asin(Math.Clamp(d.z, -1, 1)); + var deltaLongitude = WrapLongitude(longitude - parameters.CenterLongitude); + var centered = CenteredForward(deltaLongitude, latitude, parameters, out coordinate); + coordinate /= parameters.Scale; + return centered; + } + + public static bool TryInverse(double2 coordinate, PlanetaryProjectionParameters parameters, out float3 direction) + { + parameters.Validate(); + if (!double.IsFinite(coordinate.x) || !double.IsFinite(coordinate.y)) + { + direction = float3.zero; + return false; + } + coordinate *= parameters.Scale; + if (parameters.Kind == PlanetaryProjectionKind.CubeAtlas) + return TryInverseCubeAtlas(coordinate, out direction); + + if (!CenteredInverse(coordinate, parameters, out var deltaLongitude, out var latitude)) + { + direction = float3.zero; + return false; + } + var longitude = WrapLongitude(deltaLongitude + parameters.CenterLongitude); + var cosLatitude = Math.Cos(latitude); + direction = math.normalize(new float3( + (float)(cosLatitude * Math.Cos(longitude)), + (float)(cosLatitude * Math.Sin(longitude)), + (float)Math.Sin(latitude))); + return true; + } + + private static bool CenteredForward(double longitude, double latitude, PlanetaryProjectionParameters p, out double2 coordinate) + { + switch (p.Kind) + { + case PlanetaryProjectionKind.Equirectangular: + coordinate = new(longitude / Math.PI, latitude / (Math.PI * 0.5)); + return true; + case PlanetaryProjectionKind.WebMercator: + if (Math.Abs(latitude) > WebMercatorMaximumLatitude) { coordinate = double2.zero; return false; } + coordinate = new(longitude / Math.PI, Math.Asinh(Math.Tan(latitude)) / Math.PI); + return true; + case PlanetaryProjectionKind.EqualEarth: + { + var theta = Math.Asin(EqualEarthM * Math.Sin(latitude)); + var theta2 = theta * theta; + var theta6 = theta2 * theta2 * theta2; + var denominator = EqualEarthM * (EqualEarthA1 + 3 * EqualEarthA2 * theta2 + theta6 * (7 * EqualEarthA3 + 9 * EqualEarthA4 * theta2)); + coordinate = new(longitude * Math.Cos(theta) / denominator / EqualEarthXMax, EqualEarthY(theta) / EqualEarthYMax); + return true; + } + case PlanetaryProjectionKind.Orthographic: + case PlanetaryProjectionKind.AzimuthalEquidistant: + case PlanetaryProjectionKind.AzimuthalEqualArea: + case PlanetaryProjectionKind.LocalTangent: + return ForwardAzimuthal(longitude, latitude, p, out coordinate); + default: + throw new ArgumentOutOfRangeException(nameof(p)); + } + } + + private static bool CenteredInverse(double2 coordinate, PlanetaryProjectionParameters p, out double longitude, out double latitude) + { + longitude = latitude = 0; + switch (p.Kind) + { + case PlanetaryProjectionKind.Equirectangular: + if (Math.Abs(coordinate.y) > 1) return false; + longitude = coordinate.x * Math.PI; latitude = coordinate.y * Math.PI * 0.5; return Math.Abs(coordinate.x) <= 1; + case PlanetaryProjectionKind.WebMercator: + if (Math.Abs(coordinate.x) > 1 || Math.Abs(coordinate.y) > 1) return false; + longitude = coordinate.x * Math.PI; latitude = Math.Atan(Math.Sinh(coordinate.y * Math.PI)); return true; + case PlanetaryProjectionKind.EqualEarth: + { + if (Math.Abs(coordinate.x) > 1 || Math.Abs(coordinate.y) > 1) return false; + var targetY = coordinate.y * EqualEarthYMax; + var theta = targetY / EqualEarthA1; + for (var i = 0; i < 8; i++) + { + var theta2 = theta * theta; + var theta6 = theta2 * theta2 * theta2; + var f = EqualEarthY(theta) - targetY; + var derivative = EqualEarthA1 + 3 * EqualEarthA2 * theta2 + theta6 * (7 * EqualEarthA3 + 9 * EqualEarthA4 * theta2); + theta -= f / derivative; + } + var sinLatitude = Math.Sin(theta) / EqualEarthM; + if (Math.Abs(sinLatitude) > 1.0000000001) return false; + latitude = Math.Asin(Math.Clamp(sinLatitude, -1, 1)); + var thetaSquared = theta * theta; + var thetaSixth = thetaSquared * thetaSquared * thetaSquared; + var denominator = EqualEarthM * (EqualEarthA1 + 3 * EqualEarthA2 * thetaSquared + thetaSixth * (7 * EqualEarthA3 + 9 * EqualEarthA4 * thetaSquared)); + longitude = coordinate.x * EqualEarthXMax * denominator / Math.Cos(theta); + return Math.Abs(longitude) <= Math.PI + 1.0e-10; + } + case PlanetaryProjectionKind.Orthographic: + case PlanetaryProjectionKind.AzimuthalEquidistant: + case PlanetaryProjectionKind.AzimuthalEqualArea: + case PlanetaryProjectionKind.LocalTangent: + return InverseAzimuthal(coordinate, p, out longitude, out latitude); + default: + throw new ArgumentOutOfRangeException(nameof(p)); + } + } + + private static bool ForwardAzimuthal(double longitude, double latitude, PlanetaryProjectionParameters p, out double2 coordinate) + { + var sin0 = Math.Sin(p.CenterLatitude); var cos0 = Math.Cos(p.CenterLatitude); + var sin = Math.Sin(latitude); var cos = Math.Cos(latitude); var cosLon = Math.Cos(longitude); + var cosC = Math.Clamp(sin0 * sin + cos0 * cos * cosLon, -1, 1); + var rawX = cos * Math.Sin(longitude); + var rawY = cos0 * sin - sin0 * cos * cosLon; + double factor; + switch (p.Kind) + { + case PlanetaryProjectionKind.Orthographic: + if (cosC < 0) { coordinate = double2.zero; return false; } + factor = 1; break; + case PlanetaryProjectionKind.AzimuthalEquidistant: + { + var c = Math.Acos(cosC); factor = c < 1.0e-12 ? 1 : c / Math.Sin(c); factor /= Math.PI; break; + } + case PlanetaryProjectionKind.AzimuthalEqualArea: + factor = Math.Sqrt(2 / Math.Max(1 + cosC, 1.0e-20)) * 0.5; break; + case PlanetaryProjectionKind.LocalTangent: + if (cosC <= 0) { coordinate = double2.zero; return false; } + factor = 1 / cosC; break; + default: throw new ArgumentOutOfRangeException(nameof(p)); + } + coordinate = new(rawX * factor, rawY * factor); + return true; + } + + private static bool InverseAzimuthal(double2 coordinate, PlanetaryProjectionParameters p, out double longitude, out double latitude) + { + longitude = latitude = 0; + var rho = Math.Sqrt(coordinate.x * coordinate.x + coordinate.y * coordinate.y); + double c; + switch (p.Kind) + { + case PlanetaryProjectionKind.Orthographic: if (rho > 1) return false; c = Math.Asin(Math.Clamp(rho, 0, 1)); break; + case PlanetaryProjectionKind.AzimuthalEquidistant: if (rho > 1) return false; c = rho * Math.PI; break; + case PlanetaryProjectionKind.AzimuthalEqualArea: if (rho > 1) return false; c = 2 * Math.Asin(Math.Clamp(rho, 0, 1)); break; + case PlanetaryProjectionKind.LocalTangent: c = Math.Atan(rho); break; + default: throw new ArgumentOutOfRangeException(nameof(p)); + } + if (rho < 1.0e-12) { latitude = p.CenterLatitude; return true; } + var sinC = Math.Sin(c); var cosC = Math.Cos(c); + var sin0 = Math.Sin(p.CenterLatitude); var cos0 = Math.Cos(p.CenterLatitude); + latitude = Math.Asin(Math.Clamp(cosC * sin0 + coordinate.y * sinC * cos0 / rho, -1, 1)); + longitude = Math.Atan2(coordinate.x * sinC, rho * cos0 * cosC - coordinate.y * sin0 * sinC); + return true; + } + + private static bool TryInverseCubeAtlas(double2 coordinate, out float3 direction) + { + direction = float3.zero; + if (coordinate.x is < -1 or > 1 || coordinate.y is < -1 or > 1) return false; + var column = Math.Min((int)((coordinate.x + 1) * 1.5), 2); + var row = Math.Min((int)((1 - coordinate.y)), 1); + var u = ((coordinate.x + 1) * 1.5 - column) * 2 - 1; + var v = ((1 - coordinate.y) - row) * 2 - 1; + direction = PlanetaryTopology.Direction(new((PlanetaryCubeFace)(row * 3 + column), u, v)); + return true; + } + + private static double EqualEarthY(double theta) + { + var theta2 = theta * theta; + var theta6 = theta2 * theta2 * theta2; + return theta * (EqualEarthA1 + EqualEarthA2 * theta2 + theta6 * (EqualEarthA3 + EqualEarthA4 * theta2)); + } + + private static double WrapLongitude(double longitude) + { + longitude %= Math.PI * 2; + if (longitude > Math.PI) longitude -= Math.PI * 2; + if (longitude < -Math.PI) longitude += Math.PI * 2; + return longitude; + } +} + +public readonly record struct PlanetaryMapTileLayout( + PlanetaryProjectionParameters Projection, + int Level, + int X, + int Y, + int InteriorSize, + int BorderSize) +{ + public int StorageSize => checked(InteriorSize + BorderSize * 2); + + public PlanetaryMapTileLayout Validate() + { + Projection.Validate(); + if (Level is < 0 or > 30) throw new ArgumentOutOfRangeException(nameof(Level)); + var count = 1 << Level; + if ((uint)X >= (uint)count || (uint)Y >= (uint)count) throw new ArgumentOutOfRangeException(nameof(X)); + if (InteriorSize < 2 || BorderSize < 0 || BorderSize >= InteriorSize) throw new ArgumentOutOfRangeException(nameof(InteriorSize)); + return this; + } +} + +public static class PlanetaryMapTileSampling +{ + public static bool TryDirection(PlanetaryMapTileLayout tile, int storageX, int storageY, out float3 direction) + { + tile.Validate(); + if ((uint)storageX >= (uint)tile.StorageSize || (uint)storageY >= (uint)tile.StorageSize) throw new ArgumentOutOfRangeException(nameof(storageX)); + var count = 1 << tile.Level; + var localU = (storageX - tile.BorderSize) / (double)(tile.InteriorSize - 1); + var localV = (storageY - tile.BorderSize) / (double)(tile.InteriorSize - 1); + var map = new double2(-1 + 2 * (tile.X + localU) / count, 1 - 2 * (tile.Y + localV) / count); + return PlanetaryProjection.TryInverse(map, tile.Projection, out direction); + } +} + +public readonly record struct PlanetaryMapTileKey( + ulong FieldVersion, + uint ProjectionVersion, + uint LayerId, + PlanetaryMapTileLayout Layout, + PlanetaryQueryScale QueryScale) +{ + public PlanetaryMapTileKey Validate() + { + if (FieldVersion == 0 || ProjectionVersion == 0) throw new ArgumentOutOfRangeException(nameof(FieldVersion)); + Layout.Validate(); QueryScale.Validate(); + return this; + } +} + +public readonly record struct PlanetarySurfaceMapTile( + PlanetaryMapTileKey Key, + PlanetarySurfaceSample[] Samples, + bool[] Validity); + +public static class PlanetaryMapTileBaker +{ + public static PlanetarySurfaceMapTile Bake( + in PlanetaryFieldDefinition field, + PlanetaryMapTileKey key, + TSource source) + where TSource : IPlanetaryBaseField + { + field.Validate(); key.Validate(); + if (field.FieldVersion != key.FieldVersion) throw new ArgumentException("Tile field version does not match definition.", nameof(key)); + var count = key.Layout.StorageSize * key.Layout.StorageSize; + var samples = new PlanetarySurfaceSample[count]; + var validity = new bool[count]; + for (var y = 0; y < key.Layout.StorageSize; y++) + for (var x = 0; x < key.Layout.StorageSize; x++) + { + var index = y * key.Layout.StorageSize + x; + if (!PlanetaryMapTileSampling.TryDirection(key.Layout, x, y, out var direction)) continue; + samples[index] = PlanetaryField.Sample(field, direction, source.Sample(direction), key.QueryScale); + validity[index] = true; + } + return new(key, samples, validity); + } +} diff --git a/src/CultMath/PlanetaryQueries.cs b/src/CultMath/PlanetaryQueries.cs new file mode 100644 index 0000000..1f11ca7 --- /dev/null +++ b/src/CultMath/PlanetaryQueries.cs @@ -0,0 +1,191 @@ +namespace CultMath; + +public readonly record struct PlanetaryRadialBounds(float MinimumDisplacement, float MaximumDisplacement) +{ + public PlanetaryRadialBounds Validate() + { + if (!float.IsFinite(MinimumDisplacement) || !float.IsFinite(MaximumDisplacement) || MinimumDisplacement > MaximumDisplacement) + throw new ArgumentOutOfRangeException(nameof(MinimumDisplacement)); + return this; + } +} + +public readonly record struct PlanetaryRayHit( + double Distance, + double3 PlanetLocalPosition, + PlanetarySurfaceSample Surface); + +public readonly record struct PlanetaryPathSample(double Distance, PlanetarySurfaceSample Surface); + +public readonly record struct PlanetaryRegionSummary( + float MinimumDisplacement, + float MaximumDisplacement, + float MaximumSlope, + float MaximumUnresolvedHeight); + +public static class PlanetaryQueries +{ + public static bool TryIntersectRay( + in PlanetaryFieldDefinition field, + double3 rayOrigin, + double3 rayDirection, + double maximumDistance, + in PlanetaryQueryScale scale, + in PlanetaryRadialBounds bounds, + TSource source, + out PlanetaryRayHit hit) + where TSource : IPlanetaryBaseField + { + field.Validate(); scale.Validate(); bounds.Validate(); + ValidateFinite(rayOrigin, nameof(rayOrigin)); ValidateFinite(rayDirection, nameof(rayDirection)); + if (!double.IsFinite(maximumDistance) || maximumDistance <= 0) throw new ArgumentOutOfRangeException(nameof(maximumDistance)); + var directionLength = Length(rayDirection); + if (directionLength < 1.0e-15) throw new ArgumentOutOfRangeException(nameof(rayDirection)); + var direction = rayDirection / directionLength; + var outerRadius = field.Radius + bounds.MaximumDisplacement + scale.MaximumUnresolvedHeight; + if (!double.IsFinite(outerRadius)) outerRadius = field.Radius + bounds.MaximumDisplacement; + if (!TrySphereInterval(rayOrigin, direction, outerRadius, out var enter, out var leave)) { hit = default; return false; } + var start = Math.Max(enter, 0); + var end = Math.Min(leave, maximumDistance); + if (start > end) { hit = default; return false; } + + var previousDistance = start; + var previous = SignedSurfaceDistance(field, rayOrigin + direction * previousDistance, scale, source, out var previousSample); + if (previous <= 0) + { + hit = new(previousDistance, rayOrigin + direction * previousDistance, previousSample); + return true; + } + var minimumStep = Math.Max(scale.FootprintMeters * 0.5, field.Radius * 1.0e-7); + for (var iteration = 0; iteration < 512 && previousDistance < end; iteration++) + { + var nextDistance = Math.Min(previousDistance + Math.Max(previous * 0.5, minimumStep), end); + var next = SignedSurfaceDistance(field, rayOrigin + direction * nextDistance, scale, source, out var nextSample); + if (next <= 0) + { + var low = previousDistance; + var high = nextDistance; + var surface = nextSample; + for (var refine = 0; refine < 32; refine++) + { + var middle = (low + high) * 0.5; + var value = SignedSurfaceDistance(field, rayOrigin + direction * middle, scale, source, out var middleSample); + if (value > 0) low = middle; + else { high = middle; surface = middleSample; } + } + var position = rayOrigin + direction * high; + hit = new(high, position, surface); + return true; + } + if (nextDistance >= end) break; + previousDistance = nextDistance; + previous = next; + previousSample = nextSample; + } + hit = default; + return false; + } + + public static double SegmentClearance( + in PlanetaryFieldDefinition field, + double3 start, + double3 end, + double sampleSpacing, + in PlanetaryQueryScale scale, + TSource source) + where TSource : IPlanetaryBaseField + { + ValidateFinite(start, nameof(start)); ValidateFinite(end, nameof(end)); + if (!double.IsFinite(sampleSpacing) || sampleSpacing <= 0) throw new ArgumentOutOfRangeException(nameof(sampleSpacing)); + var delta = end - start; + var distance = Length(delta); + var steps = Math.Max((int)Math.Ceiling(distance / sampleSpacing), 1); + var minimum = double.PositiveInfinity; + for (var i = 0; i <= steps; i++) + { + var position = start + delta * (i / (double)steps); + minimum = Math.Min(minimum, SignedSurfaceDistance(field, position, scale, source, out _)); + } + return minimum; + } + + public static PlanetaryPathSample[] SampleGreatCircle( + in PlanetaryFieldDefinition field, + float3 startDirection, + float3 endDirection, + int segmentCount, + in PlanetaryQueryScale scale, + TSource source) + where TSource : IPlanetaryBaseField + { + if (segmentCount < 1) throw new ArgumentOutOfRangeException(nameof(segmentCount)); + PlanetaryTopology.ValidateDirection(startDirection); PlanetaryTopology.ValidateDirection(endDirection); + var a = math.normalize(startDirection); var b = math.normalize(endDirection); + var angle = MathF.Acos(math.clamp(math.dot(a, b), -1, 1)); + var sinAngle = MathF.Sin(angle); + var result = new PlanetaryPathSample[segmentCount + 1]; + for (var i = 0; i <= segmentCount; i++) + { + var t = i / (float)segmentCount; + var direction = sinAngle < 1.0e-6f + ? math.normalize(math.lerp(a, b, t)) + : math.normalize(a * (MathF.Sin((1 - t) * angle) / sinAngle) + b * (MathF.Sin(t * angle) / sinAngle)); + var sample = PlanetaryField.Sample(field, direction, source.Sample(direction), scale); + result[i] = new(angle * field.Radius * t, sample); + } + return result; + } + + public static PlanetaryRegionSummary Summarize( + in PlanetaryFieldDefinition field, + ReadOnlySpan directions, + in PlanetaryQueryScale scale, + TSource source) + where TSource : IPlanetaryBaseField + { + if (directions.IsEmpty) throw new ArgumentException("At least one direction is required.", nameof(directions)); + var min = float.PositiveInfinity; var max = float.NegativeInfinity; var slope = 0.0f; var unresolved = 0.0f; + foreach (var direction in directions) + { + var normalized = math.normalize(direction); + var sample = PlanetaryField.Sample(field, normalized, source.Sample(normalized), scale); + min = MathF.Min(min, sample.RadialDisplacement); + max = MathF.Max(max, sample.RadialDisplacement); + slope = MathF.Max(slope, sample.Slope); + unresolved = MathF.Max(unresolved, sample.UnresolvedHeightBound); + } + return new(min, max, slope, unresolved); + } + + private static double SignedSurfaceDistance( + in PlanetaryFieldDefinition field, + double3 position, + in PlanetaryQueryScale scale, + TSource source, + out PlanetarySurfaceSample sample) + where TSource : IPlanetaryBaseField + { + var radius = Length(position); + if (radius < 1.0e-15) { sample = default; return double.NegativeInfinity; } + var direction = new float3((float)(position.x / radius), (float)(position.y / radius), (float)(position.z / radius)); + sample = PlanetaryField.Sample(field, direction, source.Sample(direction), scale); + return radius - sample.SurfaceRadius; + } + + private static bool TrySphereInterval(double3 origin, double3 direction, double radius, out double enter, out double leave) + { + var b = Dot(origin, direction); + var c = Dot(origin, origin) - radius * radius; + var discriminant = b * b - c; + if (discriminant < 0) { enter = leave = 0; return false; } + var root = Math.Sqrt(discriminant); + enter = -b - root; leave = -b + root; return leave >= 0; + } + + private static double Dot(double3 a, double3 b) => a.x * b.x + a.y * b.y + a.z * b.z; + private static double Length(double3 value) => Math.Sqrt(Dot(value, value)); + private static void ValidateFinite(double3 value, string name) + { + if (!double.IsFinite(value.x) || !double.IsFinite(value.y) || !double.IsFinite(value.z)) throw new ArgumentOutOfRangeException(name); + } +} diff --git a/src/CultMath/PlanetaryResidency.cs b/src/CultMath/PlanetaryResidency.cs new file mode 100644 index 0000000..e5dca11 --- /dev/null +++ b/src/CultMath/PlanetaryResidency.cs @@ -0,0 +1,144 @@ +namespace CultMath; + +public readonly record struct PlanetaryLodParameters( + int MaximumLevel, + int PageInteriorSize, + int PageBorderSize, + float ViewportHeightPixels, + float MinimumFootprintMeters) +{ + public PlanetaryLodParameters Validate() + { + if (MaximumLevel is < 0 or > PlanetaryTileAddress.MaxLevel) throw new ArgumentOutOfRangeException(nameof(MaximumLevel)); + _ = new PlanetaryPageLayout(new(PlanetaryCubeFace.PositiveX, 0, 0, 0), PageInteriorSize, PageBorderSize).Validate(); + if (!float.IsFinite(ViewportHeightPixels) || ViewportHeightPixels <= 0) throw new ArgumentOutOfRangeException(nameof(ViewportHeightPixels)); + if (!float.IsFinite(MinimumFootprintMeters) || MinimumFootprintMeters < 0) throw new ArgumentOutOfRangeException(nameof(MinimumFootprintMeters)); + return this; + } +} + +public static class PlanetaryLodSelector +{ + public static int SelectLevel( + in PlanetaryFieldDefinition field, + float cameraDistanceFromCenter, + in PlanetaryLodParameters parameters) + { + field.Validate(); parameters.Validate(); + if (!float.IsFinite(cameraDistanceFromCenter) || cameraDistanceFromCenter < 0) throw new ArgumentOutOfRangeException(nameof(cameraDistanceFromCenter)); + var altitude = MathF.Max(cameraDistanceFromCenter - field.Radius, 0.001f); + var footprint = MathF.Max(altitude / parameters.ViewportHeightPixels * 2, parameters.MinimumFootprintMeters); + var erosion = field.Erosion; + for (var level = 0; level <= parameters.MaximumLevel; level++) + { + var page = new PlanetaryPageLayout(new(PlanetaryCubeFace.PositiveX, level, 0, 0), parameters.PageInteriorSize, parameters.PageBorderSize); + var spacing = PlanetaryPageSampling.NominalAngularTexelSize(page) * field.Radius; + var unresolved = ErosionFrequencyBands.Select( + erosion.Scale * erosion.CellScale, spacing, erosion.Octaves, + erosion.Lacunarity, erosion.Strength * erosion.Scale, erosion.Gain).UnresolvedHeightBound; + if (unresolved <= footprint) return level; + } + return parameters.MaximumLevel; + } + + public static PlanetaryTileAddress[] SelectAncestorChain( + in PlanetaryFieldDefinition field, + float3 cameraDirection, + float cameraDistanceFromCenter, + in PlanetaryLodParameters parameters) + { + PlanetaryTopology.ValidateDirection(cameraDirection); + var level = SelectLevel(field, cameraDistanceFromCenter, parameters); + var result = new PlanetaryTileAddress[6 + level]; + for (var face = 0; face < 6; face++) result[face] = new((PlanetaryCubeFace)face, 0, 0, 0); + var leaf = PlanetaryTopology.TileAt(cameraDirection, level); + for (var currentLevel = level; currentLevel >= 1; currentLevel--) + { + result[6 + currentLevel - 1] = leaf; + leaf = leaf.Parent(); + } + return result; + } +} + +public readonly record struct PlanetaryResidentTile(PlanetaryTileAddress Tile, float Blend, bool Departing); + +public readonly record struct PlanetaryResidencySnapshot( + ulong ContentVersion, + ulong PresentationVersion, + PlanetaryResidentTile[] Tiles); + +/// Owns presentation residency only. It never owns field content. +public sealed class PlanetaryResidualResidency +{ + private readonly Dictionary residents = []; + private readonly Dictionary arrivals = []; + private readonly Dictionary departures = []; + private float lastTime = float.NegativeInfinity; + + public void Reset() + { + residents.Clear(); arrivals.Clear(); departures.Clear(); lastTime = float.NegativeInfinity; + } + + public PlanetaryResidencySnapshot Update(ReadOnlySpan desired, float timeSeconds, float transitionSeconds) + { + if (!float.IsFinite(timeSeconds)) throw new ArgumentOutOfRangeException(nameof(timeSeconds)); + if (!float.IsFinite(transitionSeconds) || transitionSeconds < 0) throw new ArgumentOutOfRangeException(nameof(transitionSeconds)); + if (timeSeconds < lastTime) Reset(); + lastTime = timeSeconds; + + var desiredKeys = new HashSet(); + foreach (var tile in desired) + { + var key = tile.StableKey; + desiredKeys.Add(key); + residents[key] = tile; + if (!arrivals.ContainsKey(key)) arrivals[key] = timeSeconds; + departures.Remove(key); + } + + foreach (var pair in residents.ToArray()) + { + if (pair.Value.Level == 0 || desiredKeys.Contains(pair.Key)) continue; + if (!departures.ContainsKey(pair.Key)) + { + var startingBlend = transitionSeconds == 0 + ? 1 + : Math.Clamp((timeSeconds - arrivals[pair.Key]) / transitionSeconds, 0, 1); + departures[pair.Key] = (timeSeconds, startingBlend); + } + if (transitionSeconds > 0 && timeSeconds - departures[pair.Key].Time < transitionSeconds) continue; + residents.Remove(pair.Key); arrivals.Remove(pair.Key); departures.Remove(pair.Key); + } + + var ordered = residents + .OrderBy(pair => pair.Value.Level) + .ThenBy(pair => pair.Value.Face) + .ThenBy(pair => pair.Value.X) + .ThenBy(pair => pair.Value.Y) + .ToArray(); + var tiles = new PlanetaryResidentTile[ordered.Length]; + for (var i = 0; i < ordered.Length; i++) + { + var pair = ordered[i]; + var departing = departures.TryGetValue(pair.Key, out var departure); + var blend = pair.Value.Level == 0 || transitionSeconds == 0 + ? 1 + : departing + ? departure.StartingBlend * Math.Clamp(1 - (timeSeconds - departure.Time) / transitionSeconds, 0, 1) + : Math.Clamp((timeSeconds - arrivals[pair.Key]) / transitionSeconds, 0, 1); + tiles[i] = new(pair.Value, blend, departing); + } + return new(Hash(ordered.Select(pair => pair.Key)), Hash(tiles.Select(tile => tile.Tile.StableKey ^ (uint)BitConverter.SingleToInt32Bits(tile.Blend))), tiles); + } + + private static ulong Hash(IEnumerable values) + { + const ulong offset = 14695981039346656037UL; + const ulong prime = 1099511628211UL; + var hash = offset; + foreach (var value in values) hash = unchecked((hash ^ value) * prime); + return hash == 0 ? 1 : hash; + } +} diff --git a/src/CultMath/PlanetaryTopology.cs b/src/CultMath/PlanetaryTopology.cs new file mode 100644 index 0000000..c853507 --- /dev/null +++ b/src/CultMath/PlanetaryTopology.cs @@ -0,0 +1,170 @@ +namespace CultMath; + +public enum PlanetaryCubeFace +{ + PositiveX, + NegativeX, + PositiveY, + NegativeY, + PositiveZ, + NegativeZ, +} + +public readonly record struct PlanetaryFaceCoordinate(PlanetaryCubeFace Face, double U, double V) +{ + public PlanetaryFaceCoordinate Validate(bool allowOutsideFace = false) + { + if (!double.IsFinite(U) || !double.IsFinite(V)) throw new ArgumentOutOfRangeException(nameof(U)); + if (!allowOutsideFace && (U is < -1.0 or > 1.0 || V is < -1.0 or > 1.0)) + throw new ArgumentOutOfRangeException(nameof(U), "Face coordinates must be in [-1, 1]."); + if (!Enum.IsDefined(typeof(PlanetaryCubeFace), Face)) throw new ArgumentOutOfRangeException(nameof(Face)); + return this; + } +} + +public readonly record struct PlanetaryTileAddress +{ + public const int MaxLevel = 30; + + public PlanetaryTileAddress(PlanetaryCubeFace face, int level, int x, int y) + { + if (!Enum.IsDefined(typeof(PlanetaryCubeFace), face)) throw new ArgumentOutOfRangeException(nameof(face)); + if (level is < 0 or > MaxLevel) throw new ArgumentOutOfRangeException(nameof(level)); + var count = 1 << level; + if ((uint)x >= (uint)count) throw new ArgumentOutOfRangeException(nameof(x)); + if ((uint)y >= (uint)count) throw new ArgumentOutOfRangeException(nameof(y)); + Face = face; + Level = level; + X = x; + Y = y; + } + + public PlanetaryCubeFace Face { get; } + public int Level { get; } + public int X { get; } + public int Y { get; } + public int AxisTileCount => 1 << Level; + + public PlanetaryTileAddress Parent() => Level == 0 + ? throw new InvalidOperationException("Root cube tile has no parent.") + : new(Face, Level - 1, X >> 1, Y >> 1); + + public PlanetaryTileAddress Child(int childX, int childY) + { + if ((uint)childX > 1) throw new ArgumentOutOfRangeException(nameof(childX)); + if ((uint)childY > 1) throw new ArgumentOutOfRangeException(nameof(childY)); + if (Level == MaxLevel) throw new InvalidOperationException($"Level {MaxLevel} cannot be subdivided."); + return new(Face, Level + 1, (X << 1) | childX, (Y << 1) | childY); + } + + public PlanetaryFaceCoordinate PositionAt(double localU, double localV) + { + if (!double.IsFinite(localU) || localU is < 0.0 or > 1.0) throw new ArgumentOutOfRangeException(nameof(localU)); + if (!double.IsFinite(localV) || localV is < 0.0 or > 1.0) throw new ArgumentOutOfRangeException(nameof(localV)); + var count = AxisTileCount; + return new(Face, -1.0 + 2.0 * (X + localU) / count, -1.0 + 2.0 * (Y + localV) / count); + } + + public ulong StableKey + { + get + { + const ulong offset = 14695981039346656037UL; + const ulong prime = 1099511628211UL; + var hash = offset; + hash = unchecked((hash ^ (uint)Face) * prime); + hash = unchecked((hash ^ (uint)Level) * prime); + hash = unchecked((hash ^ (uint)X) * prime); + hash = unchecked((hash ^ (uint)Y) * prime); + return hash == 0 ? 1 : hash; + } + } + public override string ToString() => $"{Face}/L{Level:D2}/{X}/{Y}"; +} + +public static class PlanetaryTopology +{ + private const float QuarterPi = MathF.PI * 0.25f; + + public static float3 Direction(PlanetaryFaceCoordinate coordinate) + { + coordinate.Validate(allowOutsideFace: true); + var u = MathF.Tan((float)coordinate.U * QuarterPi); + var v = MathF.Tan((float)coordinate.V * QuarterPi); + var cube = coordinate.Face switch + { + PlanetaryCubeFace.PositiveX => new float3(1, v, -u), + PlanetaryCubeFace.NegativeX => new float3(-1, v, u), + PlanetaryCubeFace.PositiveY => new float3(u, 1, -v), + PlanetaryCubeFace.NegativeY => new float3(u, -1, v), + PlanetaryCubeFace.PositiveZ => new float3(u, v, 1), + PlanetaryCubeFace.NegativeZ => new float3(-u, v, -1), + _ => throw new ArgumentOutOfRangeException(nameof(coordinate)), + }; + return math.normalize(cube); + } + + public static PlanetaryFaceCoordinate FaceCoordinate(float3 direction) + { + ValidateDirection(direction); + var d = math.normalize(direction); + var a = math.abs(d); + if (a.x >= a.y && a.x >= a.z) + return d.x >= 0 + ? FromCubeUv(PlanetaryCubeFace.PositiveX, -d.z / a.x, d.y / a.x) + : FromCubeUv(PlanetaryCubeFace.NegativeX, d.z / a.x, d.y / a.x); + if (a.y >= a.z) + return d.y >= 0 + ? FromCubeUv(PlanetaryCubeFace.PositiveY, d.x / a.y, -d.z / a.y) + : FromCubeUv(PlanetaryCubeFace.NegativeY, d.x / a.y, d.z / a.y); + return d.z >= 0 + ? FromCubeUv(PlanetaryCubeFace.PositiveZ, d.x / a.z, d.y / a.z) + : FromCubeUv(PlanetaryCubeFace.NegativeZ, -d.x / a.z, d.y / a.z); + } + + public static PlanetaryTileAddress TileAt(float3 direction, int level) + { + if (level is < 0 or > PlanetaryTileAddress.MaxLevel) throw new ArgumentOutOfRangeException(nameof(level)); + var face = FaceCoordinate(direction); + var count = 1 << level; + var x = Math.Clamp((int)((face.U * 0.5 + 0.5) * count), 0, count - 1); + var y = Math.Clamp((int)((face.V * 0.5 + 0.5) * count), 0, count - 1); + return new(face.Face, level, x, y); + } + + public static bool TryLocalCoordinate(float3 direction, PlanetaryTileAddress tile, out float2 local) + { + var face = FaceCoordinate(direction); + if (face.Face != tile.Face) { local = float2.zero; return false; } + var scaledU = (face.U * 0.5 + 0.5) * tile.AxisTileCount; + var scaledV = (face.V * 0.5 + 0.5) * tile.AxisTileCount; + local = new((float)(scaledU - tile.X), (float)(scaledV - tile.Y)); + return local.x is >= 0 and <= 1 && local.y is >= 0 and <= 1; + } + + public static float3 SurfaceNormal(float3 unitDirection, float3 worldDistanceTangentGradient) + { + ValidateDirection(unitDirection); + ValidateFinite(worldDistanceTangentGradient, nameof(worldDistanceTangentGradient)); + var direction = math.normalize(unitDirection); + var tangent = worldDistanceTangentGradient - direction * math.dot(worldDistanceTangentGradient, direction); + return math.normalize(direction - tangent); + } + + private static PlanetaryFaceCoordinate FromCubeUv(PlanetaryCubeFace face, float cubeU, float cubeV) => new( + face, + MathF.Atan(cubeU) / QuarterPi, + MathF.Atan(cubeV) / QuarterPi); + + internal static void ValidateDirection(float3 direction) + { + ValidateFinite(direction, nameof(direction)); + if (math.length(direction) < 1.0e-12f) throw new ArgumentOutOfRangeException(nameof(direction)); + } + + internal static void ValidateFinite(float3 value, string name) + { + if (!float.IsFinite(value.x) || !float.IsFinite(value.y) || !float.IsFinite(value.z)) + throw new ArgumentOutOfRangeException(name); + } +} diff --git a/tests/CultMath.Tests/HlslMirrorTests.cs b/tests/CultMath.Tests/HlslMirrorTests.cs index 0e904a3..567b038 100644 --- a/tests/CultMath.Tests/HlslMirrorTests.cs +++ b/tests/CultMath.Tests/HlslMirrorTests.cs @@ -65,6 +65,22 @@ public sealed class HlslMirrorTests Assert.Contains("max(length(value), CULTMATH_NORMALIZE_EPSILON)", include); } + [Fact] + public void HlslMirrorPublishesPlanetaryTopologyAndField() + { + var root = GetRepositoryRoot(); + var include = File.ReadAllText(Path.Combine(root, "shaders", "CultMath.hlsl")); + var planetary = File.ReadAllText(Path.Combine(root, "shaders", "Planetary.hlsl")); + Assert.Contains("Planetary.hlsl", include); + Assert.Contains("cultmath_planetary_face_direction", planetary); + Assert.Contains("cultmath_planetary_face_coordinate", planetary); + Assert.Contains("cultmath_planetary_page_local", planetary); + Assert.Contains("cultmath_planetary_field_sample", planetary); + Assert.Contains("cultmath_planetary_equirectangular_forward", planetary); + Assert.Contains("cultmath_planetary_equal_earth_forward", planetary); + Assert.Contains("CultMathPlanetarySurfaceSample", planetary); + } + private static string GetRepositoryRoot() { var directory = new DirectoryInfo(AppContext.BaseDirectory); diff --git a/tests/CultMath.Tests/PlanetaryFieldTests.cs b/tests/CultMath.Tests/PlanetaryFieldTests.cs new file mode 100644 index 0000000..4dc7532 --- /dev/null +++ b/tests/CultMath.Tests/PlanetaryFieldTests.cs @@ -0,0 +1,80 @@ +using Xunit; + +namespace CultMath.Tests; + +public sealed class PlanetaryFieldTests +{ + private static readonly PlanetaryFieldDefinition Definition = new( + 0x7f59_6d11_02ab_9c31UL, + 8.4f, + 42, + AdvancedErosionParameters.Default with { Scale = 0.63f, Strength = 0.12f, Octaves = 7 }); + + [Fact] + public void SampleIsDeterministicScaleAwareAndTangent() + { + var direction = math.normalize(new float3(0.31f, -0.72f, 0.61f)); + var source = new FixtureField(); + var coarse = PlanetaryField.Sample(Definition, direction, source.Sample(direction), PlanetaryQueryScale.AtFootprint(0.1f)); + var fine = PlanetaryField.Sample(Definition, direction, source.Sample(direction), PlanetaryQueryScale.AtFootprint(0.01f)); + var repeated = PlanetaryField.Sample(Definition, direction, source.Sample(direction), PlanetaryQueryScale.AtFootprint(0.01f)); + + Assert.Equal(fine, repeated); + Assert.Equal(Definition.FieldVersion, fine.FieldVersion); + Assert.True(fine.UnresolvedHeightBound <= coarse.UnresolvedHeightBound); + Assert.True(fine.FinestResolvedWavelength <= coarse.FinestResolvedWavelength); + Assert.InRange(MathF.Abs(math.dot(fine.UnitDirection, fine.TangentGradient)), 0, 2.0e-5f); + Assert.Equal(1, math.length(fine.SurfaceNormal), 5); + } + + [Fact] + public void BatchAndDoublePositionUseThePointContract() + { + var directions = new[] + { + math.normalize(new float3(1, 2, 3)), + math.normalize(new float3(-2, 0.5f, 1)), + math.normalize(new float3(0.1f, -0.2f, -1)), + }; + var output = new PlanetarySurfaceSample[directions.Length]; + var source = new FixtureField(); + var scale = PlanetaryQueryScale.AtFootprint(0.05f); + PlanetaryField.SampleBatch(Definition, directions, scale, source, output); + + for (var i = 0; i < directions.Length; i++) + { + var normalized = math.normalize(directions[i]); + var expected = PlanetaryField.Sample(Definition, normalized, source.Sample(normalized), scale); + Assert.InRange(MathF.Abs(expected.RadialDisplacement - output[i].RadialDisplacement), 0, 1.0e-6f); + Assert.InRange(math.distance(expected.TangentGradient, output[i].TangentGradient), 0, 1.0e-6f); + } + + var d = directions[0]; + var position = new double3(d.x * 8.0e11, d.y * 8.0e11, d.z * 8.0e11); + var positioned = PlanetaryField.SamplePosition(Definition, position, scale, source); + Assert.InRange(math.distance(positioned.UnitDirection, d), 0, 2.0e-6f); + } + + [Fact] + public void FieldIdentityChangesWithEveryLoadBearingInput() + { + var erosion = AdvancedErosionParameters.Default; + var baseline = PlanetaryFieldDefinition.Create(3, 8.4f, 7, erosion); + Assert.Equal(baseline.FieldVersion, PlanetaryFieldDefinition.Create(3, 8.4f, 7, erosion).FieldVersion); + Assert.NotEqual(baseline.FieldVersion, PlanetaryFieldDefinition.Create(4, 8.4f, 7, erosion).FieldVersion); + Assert.NotEqual(baseline.FieldVersion, PlanetaryFieldDefinition.Create(3, 8.5f, 7, erosion).FieldVersion); + Assert.NotEqual(baseline.FieldVersion, PlanetaryFieldDefinition.Create(3, 8.4f, 8, erosion).FieldVersion); + Assert.NotEqual(baseline.FieldVersion, PlanetaryFieldDefinition.Create(3, 8.4f, 7, erosion with { Gain = 0.49f }).FieldVersion); + } + + private readonly struct FixtureField : IPlanetaryBaseField + { + public PlanetaryBaseFieldSample Sample(float3 direction) + { + var value = direction.x * 0.6f + direction.y * direction.z * 0.3f; + var fieldGradient = new float3(0.6f, direction.z * 0.3f, direction.y * 0.3f); + fieldGradient -= direction * math.dot(fieldGradient, direction); + return new(value * 0.1f, fieldGradient * 0.1f, value, fieldGradient, math.clamp(value, -1, 1)); + } + } +} diff --git a/tests/CultMath.Tests/PlanetaryPageTests.cs b/tests/CultMath.Tests/PlanetaryPageTests.cs new file mode 100644 index 0000000..c4e4045 --- /dev/null +++ b/tests/CultMath.Tests/PlanetaryPageTests.cs @@ -0,0 +1,57 @@ +using Xunit; + +namespace CultMath.Tests; + +public sealed class PlanetaryPageTests +{ + [Fact] + public void BilinearSamplingAndSummaryUseOnePageContract() + { + var layout = new PlanetaryPageLayout(new(PlanetaryCubeFace.PositiveZ, 0, 0, 0), 3, 1); + var samples = new PlanetaryPageSample[layout.StorageSize * layout.StorageSize]; + for (var y = 0; y < layout.StorageSize; y++) + for (var x = 0; x < layout.StorageSize; x++) + samples[y * layout.StorageSize + x] = new(x + y * 10, new float3(x, y, 0), x / 10.0f, y / 10.0f); + + var sample = PlanetaryPageSampling.Bilinear(samples, layout, new float2(0.25f, 0.75f)); + Assert.Equal(26.5f, sample.Height, 4); + Assert.Equal(new float3(1.5f, 2.5f, 0), sample.TangentGradient); + var summary = PlanetaryPageSampling.Summarize(samples, layout, 0.125f); + Assert.Equal(0, summary.MinimumHeight); + Assert.Equal(44, summary.MaximumHeight); + Assert.Equal(0.125f, summary.UnresolvedHeightBound); + Assert.True(summary.MaximumSlope > 0); + } + + [Fact] + public void ChildPageStoresOnlyResidualAndComposesBackToDirectField() + { + var field = new PlanetaryFieldDefinition(7, 8.4f, 3, AdvancedErosionParameters.Default with + { + Scale = 0.63f, + Strength = 0.12f, + Octaves = 7, + }); + var source = new FixtureField(); + var root = PlanetaryPageBaker.Bake(field, new(new(PlanetaryCubeFace.PositiveZ, 0, 0, 0), 17, 2), source); + var child = PlanetaryPageBaker.Bake(field, new(new(PlanetaryCubeFace.PositiveZ, 1, 1, 0), 17, 2), source); + var direction = PlanetaryPageSampling.DirectionAtLocal(child.Layout, 0.5, 0.5); + var pages = new[] { new PlanetaryResidentPage(root, 1), new PlanetaryResidentPage(child, 1) }; + Assert.True(PlanetaryPageSetSampling.TrySample(pages, direction, out var composed)); + var spacing = PlanetaryPageSampling.NominalAngularTexelSize(child.Layout) * field.Radius; + var direct = PlanetaryField.Sample(field, direction, source.Sample(direction), PlanetaryQueryScale.AtFootprint(spacing)); + Assert.InRange(MathF.Abs(composed.Height - direct.RadialDisplacement), 0, 2.0e-5f); + Assert.InRange(math.distance(composed.TangentGradient, direct.TangentGradient), 0, 3.0e-5f); + } + + private readonly struct FixtureField : IPlanetaryBaseField + { + public PlanetaryBaseFieldSample Sample(float3 direction) + { + var value = direction.x * 0.4f - direction.y * 0.2f + direction.z * 0.5f; + var gradient = new float3(0.4f, -0.2f, 0.5f); + gradient -= direction * math.dot(gradient, direction); + return new(value * 0.08f, gradient * 0.08f, value, gradient, value); + } + } +} diff --git a/tests/CultMath.Tests/PlanetaryPatchTests.cs b/tests/CultMath.Tests/PlanetaryPatchTests.cs new file mode 100644 index 0000000..6564c77 --- /dev/null +++ b/tests/CultMath.Tests/PlanetaryPatchTests.cs @@ -0,0 +1,40 @@ +using Xunit; + +namespace CultMath.Tests; + +public sealed class PlanetaryPatchTests +{ + [Fact] + public void CubeSpherePatchMeshUsesCanonicalTopology() + { + var meshes = PlanetaryPatch.CreateCubeSphere(8); + Assert.Equal(6, meshes.Length); + foreach (var mesh in meshes) + { + Assert.Equal(81, mesh.Vertices.Length); + Assert.Equal(8 * 8 * 6, mesh.Indices.Length); + Assert.All(mesh.Vertices, vertex => Assert.Equal(1, math.length(vertex.UnitDirection), 5)); + foreach (var index in mesh.Indices) Assert.InRange(index, 0, mesh.Vertices.Length - 1); + } + } + + [Fact] + public void NeighborFaceCornersAreSingleValued() + { + var directions = PlanetaryPatch.CreateCubeSphere(1) + .SelectMany(mesh => mesh.Vertices.Select(vertex => vertex.UnitDirection)) + .ToArray(); + var groups = directions.GroupBy(direction => new + { + X = Math.Sign(direction.x), + Y = Math.Sign(direction.y), + Z = Math.Sign(direction.z), + }); + Assert.Equal(8, groups.Count()); + Assert.All(groups, group => + { + var first = group.First(); + Assert.All(group, direction => Assert.InRange(math.distance(first, direction), 0, 1.0e-6f)); + }); + } +} diff --git a/tests/CultMath.Tests/PlanetaryProjectionTests.cs b/tests/CultMath.Tests/PlanetaryProjectionTests.cs new file mode 100644 index 0000000..ad811d6 --- /dev/null +++ b/tests/CultMath.Tests/PlanetaryProjectionTests.cs @@ -0,0 +1,73 @@ +using Xunit; + +namespace CultMath.Tests; + +public sealed class PlanetaryProjectionTests +{ + [Theory] + [InlineData(PlanetaryProjectionKind.Equirectangular)] + [InlineData(PlanetaryProjectionKind.WebMercator)] + [InlineData(PlanetaryProjectionKind.EqualEarth)] + [InlineData(PlanetaryProjectionKind.Orthographic)] + [InlineData(PlanetaryProjectionKind.AzimuthalEquidistant)] + [InlineData(PlanetaryProjectionKind.AzimuthalEqualArea)] + [InlineData(PlanetaryProjectionKind.LocalTangent)] + public void ProjectionRoundTripsVisibleDirections(PlanetaryProjectionKind kind) + { + var parameters = new PlanetaryProjectionParameters(kind, 0.31, -0.22); + foreach (var longitude in new[] { -0.4, -0.1, 0.2, 0.65 }) + foreach (var latitude in new[] { -0.6, -0.15, 0.3, 0.7 }) + { + var direction = Direction(longitude + parameters.CenterLongitude, latitude); + if (!PlanetaryProjection.TryForward(direction, parameters, out var map)) continue; + Assert.True(PlanetaryProjection.TryInverse(map, parameters, out var roundTrip)); + Assert.InRange(math.distance(direction, roundTrip), 0, 3.0e-5f); + } + } + + [Fact] + public void CubeAtlasRoundTripsEveryFaceInterior() + { + var projection = new PlanetaryProjectionParameters(PlanetaryProjectionKind.CubeAtlas); + foreach (var face in Enum.GetValues()) + { + var direction = PlanetaryTopology.Direction(new(face, 0.23, -0.41)); + Assert.True(PlanetaryProjection.TryForward(direction, projection, out var map)); + Assert.True(PlanetaryProjection.TryInverse(map, projection, out var roundTrip)); + Assert.InRange(math.distance(direction, roundTrip), 0, 2.0e-6f); + } + } + + [Fact] + public void MercatorRejectsPolarSingularity() + { + var northPole = new float3(0, 0, 1); + Assert.False(PlanetaryProjection.TryForward(northPole, new(PlanetaryProjectionKind.WebMercator), out _)); + } + + [Fact] + public void MapTileBakerCarriesFieldIdentityAndProjectionValidity() + { + var field = new PlanetaryFieldDefinition(99, 10, 0, AdvancedErosionParameters.Default); + var layout = new PlanetaryMapTileLayout( + new(PlanetaryProjectionKind.Orthographic), 0, 0, 0, 9, 1); + var key = new PlanetaryMapTileKey(field.FieldVersion, 1, 7, layout, PlanetaryQueryScale.AtFootprint(0.1f)); + var tile = PlanetaryMapTileBaker.Bake(field, key, new FlatField()); + Assert.Equal(key, tile.Key); + Assert.Equal(layout.StorageSize * layout.StorageSize, tile.Samples.Length); + Assert.Contains(tile.Validity, valid => valid); + Assert.Contains(tile.Validity, valid => !valid); + Assert.All(tile.Samples.Where((_, index) => tile.Validity[index]), sample => Assert.Equal(field.FieldVersion, sample.FieldVersion)); + } + + private static float3 Direction(double longitude, double latitude) + { + var cos = Math.Cos(latitude); + return new((float)(cos * Math.Cos(longitude)), (float)(cos * Math.Sin(longitude)), (float)Math.Sin(latitude)); + } + + private readonly struct FlatField : IPlanetaryBaseField + { + public PlanetaryBaseFieldSample Sample(float3 unitDirection) => new(0, float3.zero, 0, float3.zero, 0); + } +} diff --git a/tests/CultMath.Tests/PlanetaryQueryTests.cs b/tests/CultMath.Tests/PlanetaryQueryTests.cs new file mode 100644 index 0000000..48d865f --- /dev/null +++ b/tests/CultMath.Tests/PlanetaryQueryTests.cs @@ -0,0 +1,46 @@ +using Xunit; + +namespace CultMath.Tests; + +public sealed class PlanetaryQueryTests +{ + private static readonly PlanetaryFieldDefinition Field = new(12, 10, 0, AdvancedErosionParameters.Default with { Strength = 0, Octaves = 0 }); + private static readonly PlanetaryQueryScale Scale = PlanetaryQueryScale.AtFootprint(0.1f, 0); + + [Fact] + public void RayIntersectionAndClearanceUseTheCpuSurfaceContract() + { + var source = new BulgedSphere(); + Assert.True(PlanetaryQueries.TryIntersectRay( + Field, new double3(0, 0, 20), new double3(0, 0, -1), 30, + Scale, new PlanetaryRadialBounds(-1, 1), source, out var hit)); + Assert.InRange(hit.Distance, 8.99999, 9.00001); + Assert.InRange(hit.PlanetLocalPosition.z, 10.99999, 11.00001); + var clearance = PlanetaryQueries.SegmentClearance( + Field, new double3(-2, 0, 12), new double3(2, 0, 12), 0.1, Scale, source); + Assert.True(clearance > 0.8); + } + + [Fact] + public void GreatCircleAndRegionQueriesReturnVersionedSamples() + { + var source = new BulgedSphere(); + var path = PlanetaryQueries.SampleGreatCircle(Field, new float3(1, 0, 0), new float3(0, 1, 0), 8, Scale, source); + Assert.Equal(9, path.Length); + Assert.Equal(MathF.PI * 5, path[^1].Distance, 4); + Assert.All(path, point => Assert.Equal(Field.FieldVersion, point.Surface.FieldVersion)); + var summary = PlanetaryQueries.Summarize(Field, path.Select(point => point.Surface.UnitDirection).ToArray(), Scale, source); + Assert.True(summary.MaximumDisplacement >= summary.MinimumDisplacement); + Assert.True(summary.MaximumSlope >= 0); + } + + private readonly struct BulgedSphere : IPlanetaryBaseField + { + public PlanetaryBaseFieldSample Sample(float3 direction) + { + var displacement = direction.z; + var gradient = new float3(0, 0, 1) - direction * direction.z; + return new(displacement, gradient, 0, float3.zero, 0); + } + } +} diff --git a/tests/CultMath.Tests/PlanetaryResidencyTests.cs b/tests/CultMath.Tests/PlanetaryResidencyTests.cs new file mode 100644 index 0000000..d2b89b7 --- /dev/null +++ b/tests/CultMath.Tests/PlanetaryResidencyTests.cs @@ -0,0 +1,57 @@ +using Xunit; + +namespace CultMath.Tests; + +public sealed class PlanetaryResidencyTests +{ + private static readonly PlanetaryFieldDefinition Field = new( + 1, 8.4f, 0, AdvancedErosionParameters.Default with { Scale = 0.63f, Strength = 0.12f, Octaves = 7 }); + private static readonly PlanetaryLodParameters Lod = new(5, 17, 2, 720, 8.4f / 65536); + + [Fact] + public void SelectorReturnsSixRootsAndOneAncestorChain() + { + var direction = math.normalize(new float3(-0.2f, 0.91f, 0.37f)); + var selected = PlanetaryLodSelector.SelectAncestorChain(Field, direction, 8.45f, Lod); + Assert.Equal(6, selected.Count(tile => tile.Level == 0)); + var descendants = selected.Where(tile => tile.Level > 0).OrderBy(tile => tile.Level).ToArray(); + for (var i = 1; i < descendants.Length; i++) Assert.Equal(descendants[i - 1], descendants[i].Parent()); + if (descendants.Length > 0) Assert.Equal(PlanetaryTopology.TileAt(direction, descendants[^1].Level), descendants[^1]); + } + + [Fact] + public void ResidencyFadesResidualsWithoutChangingContentDuringPresentation() + { + var root = new PlanetaryTileAddress(PlanetaryCubeFace.PositiveX, 0, 0, 0); + var child = root.Child(1, 0); + var residency = new PlanetaryResidualResidency(); + var arrival = residency.Update(new[] { root, child }, 0, 1); + var midpoint = residency.Update(new[] { root, child }, 0.5f, 1); + var settled = residency.Update(new[] { root, child }, 1, 1); + Assert.Equal(arrival.ContentVersion, midpoint.ContentVersion); + Assert.Equal(midpoint.ContentVersion, settled.ContentVersion); + Assert.Equal(0, arrival.Tiles.Single(tile => tile.Tile == child).Blend); + Assert.Equal(0.5f, midpoint.Tiles.Single(tile => tile.Tile == child).Blend); + Assert.Equal(1, settled.Tiles.Single(tile => tile.Tile == child).Blend); + + var leaving = residency.Update(new[] { root }, 1.25f, 1); + Assert.True(leaving.Tiles.Single(tile => tile.Tile == child).Departing); + Assert.Equal(1, leaving.Tiles.Single(tile => tile.Tile == child).Blend); + var gone = residency.Update(new[] { root }, 2.25f, 1); + Assert.DoesNotContain(gone.Tiles, tile => tile.Tile == child); + } + + [Fact] + public void ReversingAnArrivalDepartsFromItsCurrentWeightWithoutJumping() + { + var root = new PlanetaryTileAddress(PlanetaryCubeFace.PositiveX, 0, 0, 0); + var child = root.Child(0, 1); + var residency = new PlanetaryResidualResidency(); + _ = residency.Update(new[] { root, child }, 0, 1); + var arriving = residency.Update(new[] { root, child }, 0.25f, 1); + var departing = residency.Update(new[] { root }, 0.25f, 1); + Assert.Equal(0.25f, arriving.Tiles.Single(tile => tile.Tile == child).Blend); + Assert.Equal(0.25f, departing.Tiles.Single(tile => tile.Tile == child).Blend); + Assert.True(departing.Tiles.Single(tile => tile.Tile == child).Departing); + } +} diff --git a/tests/CultMath.Tests/PlanetaryTopologyTests.cs b/tests/CultMath.Tests/PlanetaryTopologyTests.cs new file mode 100644 index 0000000..6e6b083 --- /dev/null +++ b/tests/CultMath.Tests/PlanetaryTopologyTests.cs @@ -0,0 +1,55 @@ +using Xunit; + +namespace CultMath.Tests; + +public sealed class PlanetaryTopologyTests +{ + [Fact] + public void FaceDirectionRoundTripsAcrossEveryFace() + { + foreach (var face in Enum.GetValues()) + foreach (var u in new[] { -0.9, -0.25, 0.0, 0.4, 0.9 }) + foreach (var v in new[] { -0.85, -0.1, 0.3, 0.88 }) + { + var direction = PlanetaryTopology.Direction(new(face, u, v)); + var roundTrip = PlanetaryTopology.FaceCoordinate(direction); + Assert.Equal(face, roundTrip.Face); + Assert.Equal(u, roundTrip.U, 5); + Assert.Equal(v, roundTrip.V, 5); + } + } + + [Fact] + public void SiblingPagesShareBitStableBoundaryDirections() + { + var left = new PlanetaryPageLayout(new(PlanetaryCubeFace.PositiveZ, 1, 0, 0), 17, 2); + var right = new PlanetaryPageLayout(new(PlanetaryCubeFace.PositiveZ, 1, 1, 0), 17, 2); + for (var y = 0; y < left.InteriorSize; y++) + { + var a = PlanetaryPageSampling.Direction(left, left.BorderSize + left.InteriorSize - 1, left.BorderSize + y); + var b = PlanetaryPageSampling.Direction(right, right.BorderSize, right.BorderSize + y); + Assert.Equal(a, b); + } + } + + [Fact] + public void DirectionSelectsContainingTileAndLocalCoordinate() + { + var expected = new PlanetaryTileAddress(PlanetaryCubeFace.NegativeY, 5, 19, 7); + var direction = PlanetaryPageSampling.DirectionAtLocal(new(expected, 17, 2), 0.37, 0.61); + Assert.Equal(expected, PlanetaryTopology.TileAt(direction, expected.Level)); + Assert.True(PlanetaryTopology.TryLocalCoordinate(direction, expected, out var local)); + Assert.Equal(0.37f, local.x, 4); + Assert.Equal(0.61f, local.y, 4); + } + + [Fact] + public void BorderSamplesDoNotRequireNeighborResidency() + { + var page = new PlanetaryPageLayout(new(PlanetaryCubeFace.PositiveX, 0, 0, 0), 17, 3); + var edge = PlanetaryPageSampling.Direction(page, page.BorderSize + page.InteriorSize - 1, page.BorderSize + 8); + var border = PlanetaryPageSampling.Direction(page, page.StorageSize - 1, page.BorderSize + 8); + Assert.NotEqual(edge, border); + Assert.Equal(1, math.length(border), 5); + } +} diff --git a/unity/Runtime/CultMath.Unity.asmdef b/unity/Runtime/CultMath.Unity.asmdef new file mode 100644 index 0000000..632c4e4 --- /dev/null +++ b/unity/Runtime/CultMath.Unity.asmdef @@ -0,0 +1,16 @@ +{ + "name": "CultMath.Unity", + "rootNamespace": "CultMath.Unity", + "references": [ + "CultMath" + ], + "includePlatforms": [], + "excludePlatforms": [], + "allowUnsafeCode": false, + "overrideReferences": false, + "precompiledReferences": [], + "autoReferenced": true, + "defineConstraints": [], + "versionDefines": [], + "noEngineReferences": false +} diff --git a/unity/Runtime/PlanetaryPageUpload.cs b/unity/Runtime/PlanetaryPageUpload.cs new file mode 100644 index 0000000..f0b168c --- /dev/null +++ b/unity/Runtime/PlanetaryPageUpload.cs @@ -0,0 +1,81 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Runtime.InteropServices; +using UnityEngine; + +namespace CultMath.Unity; + +[StructLayout(LayoutKind.Sequential)] +public readonly struct PlanetaryPageInputGpu +{ + public PlanetaryPageInputGpu(Vector4 directionRadius, Vector4 sampling) + { + DirectionRadius = directionRadius; + Sampling = sampling; + } + + public readonly Vector4 DirectionRadius; + public readonly Vector4 Sampling; +} + +[StructLayout(LayoutKind.Sequential)] +public readonly struct PlanetaryPageMetadataGpu +{ + public PlanetaryPageMetadataGpu(Vector4 address, Vector4 layout, Vector4 bounds, Vector4 state) + { + Address = address; + Layout = layout; + Bounds = bounds; + State = state; + } + + public readonly Vector4 Address; + public readonly Vector4 Layout; + public readonly Vector4 Bounds; + public readonly Vector4 State; +} + +public readonly record struct PlanetaryPageUpload( + PlanetaryPageInputGpu[] Inputs, + PlanetaryPageMetadataGpu[] Metadata, + int OutputSampleCount, + ulong ContentVersion, + ulong PresentationVersion); + +public static class PlanetaryPageUploadBuilder +{ + public static PlanetaryPageUpload Build( + in PlanetaryResidencySnapshot snapshot, + int interiorSize, + int borderSize, + float radius) + { + if (!float.IsFinite(radius) || radius <= 0) throw new ArgumentOutOfRangeException(nameof(radius)); + var inputs = new List(); + var metadata = new PlanetaryPageMetadataGpu[snapshot.Tiles.Length]; + foreach (var (resident, pageIndex) in snapshot.Tiles.Select((resident, index) => (resident, index))) + { + var layout = new PlanetaryPageLayout(resident.Tile, interiorSize, borderSize).Validate(); + var spacing = PlanetaryPageSampling.NominalAngularTexelSize(layout) * radius; + var parentSpacing = resident.Tile.Level == 0 + ? 0 + : PlanetaryPageSampling.NominalAngularTexelSize(new(resident.Tile.Parent(), interiorSize, borderSize)) * radius; + var outputOffset = inputs.Count; + for (var y = 0; y < layout.StorageSize; y++) + for (var x = 0; x < layout.StorageSize; x++) + { + var direction = PlanetaryPageSampling.Direction(layout, x, y); + inputs.Add(new( + new Vector4(direction.x, direction.y, direction.z, radius), + new Vector4(spacing, parentSpacing, 0, 0))); + } + metadata[pageIndex] = new( + new Vector4((int)resident.Tile.Face, resident.Tile.Level, resident.Tile.X, resident.Tile.Y), + new Vector4(outputOffset, layout.StorageSize, layout.InteriorSize, layout.BorderSize), + Vector4.zero, + new Vector4(1, resident.Blend, PlanetaryPageSampling.NominalAngularTexelSize(layout), spacing)); + } + return new(inputs.ToArray(), metadata, inputs.Count, snapshot.ContentVersion, snapshot.PresentationVersion); + } +} diff --git a/unity/Runtime/PlanetaryPatchMeshAdapter.cs b/unity/Runtime/PlanetaryPatchMeshAdapter.cs new file mode 100644 index 0000000..69a72a2 --- /dev/null +++ b/unity/Runtime/PlanetaryPatchMeshAdapter.cs @@ -0,0 +1,36 @@ +using System; +using UnityEngine; +using UnityEngine.Rendering; + +namespace CultMath.Unity; + +public static class PlanetaryPatchMeshAdapter +{ + public static Mesh CreateFaceMesh(PlanetaryCubeFace face, int cellsPerAxis, float radius = 1) + { + if (!float.IsFinite(radius) || radius <= 0) throw new ArgumentOutOfRangeException(nameof(radius)); + var source = PlanetaryPatch.CreateFace(face, cellsPerAxis); + var vertices = new Vector3[source.Vertices.Length]; + var normals = new Vector3[source.Vertices.Length]; + var uv = new Vector2[source.Vertices.Length]; + for (var i = 0; i < source.Vertices.Length; i++) + { + var direction = source.Vertices[i].UnitDirection; + var normal = new Vector3(direction.x, direction.y, direction.z); + vertices[i] = normal * radius; + normals[i] = normal; + uv[i] = new Vector2(source.Vertices[i].LocalCoordinate.x, source.Vertices[i].LocalCoordinate.y); + } + var mesh = new Mesh + { + name = $"CultMath Planetary {face} {cellsPerAxis}x{cellsPerAxis}", + indexFormat = vertices.Length > ushort.MaxValue ? IndexFormat.UInt32 : IndexFormat.UInt16, + vertices = vertices, + normals = normals, + uv = uv, + triangles = source.Indices, + }; + mesh.RecalculateBounds(); + return mesh; + } +}