320 lines
14 KiB
C#
320 lines
14 KiB
C#
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<TSource>(
|
|
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);
|
|
}
|
|
}
|