Cylinder voxels

Shaders/Private/CesiumCylinder.usf

@@ -0,0 +1,303 @@
+#pragma once
+
+// Copyright 2020-2025 CesiumGS, Inc. and Contributors
+
+/*=============================================================================
+	CesiumCylinder.usf: An implicit cylinder that may be intersected by a ray.
+=============================================================================*/
+
+#include "CesiumRayIntersectionList.usf"
+#include "CesiumLongitudeWedge.usf"
+
+struct Cylinder
+{
+  float3 MinBounds;
+  float3 MaxBounds;
+
+  float RadiusUVScale;
+  float RadiusUVOffset;
+
+  float3 AngleUVExtents; // X = min, Y = max, Z = midpoint
+  float AngleUVScale;
+  float AngleUVOffset;
+
+  bool RadiusHasMinimumBound;
+  bool RadiusIsFlat;
+
+  uint AngleRangeFlag;
+  bool AngleMinAtDiscontinuity;
+  bool AngleMaxAtDiscontinuity;
+  bool AngleIsReversed;
+
+  void Initialize(float3 InMinBounds, float3 InMaxBounds, float4 PackedData0, float4 PackedData1, float4 PackedData2, float4 PackedData3)
+  {
+    MinBounds = InMinBounds; // normalized radius, angle, height
+    MaxBounds = InMaxBounds; // normalized radius, angle, height
+
+    // Flags are packed in CesiumVoxelRendererComponent.cpp
+    RadiusHasMinimumBound = bool(PackedData0.x);
+    RadiusIsFlat = bool(PackedData0.y);
+
+    AngleRangeFlag = round(PackedData1.x);
+    AngleMinAtDiscontinuity = bool(PackedData1.y);
+    AngleMaxAtDiscontinuity = bool(PackedData1.z);
+    AngleIsReversed = bool(PackedData1.w);
+
+    AngleUVExtents = PackedData2.xyz;
+
+    RadiusUVScale = PackedData3.x;
+    RadiusUVOffset = PackedData3.y;
+    AngleUVScale = PackedData3.z;
+    AngleUVOffset = PackedData3.w;
+  }
+
+  /**
+  * Tests whether the input ray intersects the volume defined by two xy-planes at the specified z values.
+  */
+  RayIntersections IntersectHeightBounds(in Ray R, in float2 HeightBounds, in bool IsConvex)
+  {
+    float zPosition = R.Origin.z;
+    float zDirection = R.Direction.z;
+
+    float tmin = (HeightBounds.x - zPosition) / zDirection;
+    float tmax = (HeightBounds.y - zPosition) / zDirection;
+
+    // Normals point outside the volume
+    float signFlip = IsConvex ? 1.0 : -1.0;
+    Intersection min = NewIntersection(tmin, float3(0.0, 0.0, -1.0 * signFlip));
+    Intersection max = NewIntersection(tmax, float3(0.0, 0.0, 1.0 * signFlip));
+
+    if (zDirection < 0.0)
+    {
+      // Ray entered from the top plane down.
+      return NewRayIntersections(max, min);
+    }
+    else
+    {
+      return NewRayIntersections(min, max);
+    }
+  }
+
+  /**
+  * Tests whether the input ray intersects an infinite cylinder with the given radius.
+  */
+  RayIntersections IntersectInfiniteCylinder(in Ray R, in float Radius, in bool IsConvex)
+  {
+    float3 o = R.Origin;
+    float3 d = R.Direction;
+
+    float a = dot(d.xy, d.xy);
+    float b = dot(o.xy, d.xy);
+    float c = dot(o.xy, o.xy) - Radius * Radius;
+    float determinant = b * b - a * c;
+
+    if (determinant < 0.0)
+    {
+      Intersection miss = NewIntersection(NO_HIT, normalize(R.Direction));
+      return NewRayIntersections(miss, miss);
+    }
+
+    determinant = sqrt(determinant);
+    float t1 = (-b - determinant) / a;
+    float t2 = (-b + determinant) / a;
+    float signFlip = IsConvex ? 1.0 : -1.0;
+    Intersection intersect1 = NewIntersection(t1, float3(normalize(o.xy + t1 * d.xy) * signFlip, 0.0));
+    Intersection intersect2 = NewIntersection(t2, float3(normalize(o.xy + t2 * d.xy) * signFlip, 0.0));
+
+    return NewRayIntersections(intersect1, intersect2);
+  }
+
+  /**
+   * Tests whether the given ray intersects a right cylindrical solid of the given radius and height bounds.
+   * The shape is assumed to be convex.
+   */
+  RayIntersections IntersectBoundedCylinder(in Ray R, in float Radius, in float2 MinMaxHeight)
+  {
+    RayIntersections infiniteCylinderIntersection = IntersectInfiniteCylinder(R, Radius, true);
+    RayIntersections heightIntersection = IntersectHeightBounds(R, MinMaxHeight, true);
+    return ResolveIntersections(infiniteCylinderIntersection, heightIntersection);
+  }
+
+  /**
+   * Tests whether the input ray (Unit Space) intersects the cylinder. Outputs the intersections in Unit Space.
+   */
+  void Intersect(in Ray R, out float4 Intersections[INTERSECTIONS_LENGTH], inout IntersectionListState ListState)
+  {
+    ListState.Length = 2;
+
+    RayIntersections OuterResult = IntersectBoundedCylinder(R, MaxBounds.x, float2(MinBounds.z, MaxBounds.z));
+    setShapeIntersections(Intersections, ListState, 0, OuterResult);
+
+    if (OuterResult.Entry.t == NO_HIT)
+    {
+      return;
+    }
+
+    if (RadiusIsFlat)
+    {
+      // When the cylinder is perfectly thin it's necessary to sandwich the
+      // inner cylinder intersection inside the outer cylinder intersection.
+
+      // Without this special case,
+      // [outerMin, outerMax, innerMin, innerMax] will bubble sort to
+      // [outerMin, innerMin, outerMax, innerMax] which will cause the back
+      // side of the cylinder to be invisible because it will think the ray
+      // is still inside the inner (negative) cylinder after exiting the
+      // outer (positive) cylinder.
+
+      // With this special case,
+      // [outerMin, innerMin, innerMax, outerMax] will bubble sort to
+      // [outerMin, innerMin, innerMax, outerMax] which will work correctly.
+
+      // Note: If Sort() changes its sorting function
+      // from bubble sort to something else, this code may need to change.
+      RayIntersections InnerResult = IntersectInfiniteCylinder(R, 1.0, false);
+      setSurfaceIntersection(Intersections, ListState, 0, OuterResult.Entry, true, true); // positive, entering
+      setSurfaceIntersection(Intersections, ListState, 1, InnerResult.Entry, false, true); // negative, entering
+      setSurfaceIntersection(Intersections, ListState, 2, InnerResult.Exit, false, false); // negative, exiting
+      setSurfaceIntersection(Intersections, ListState, 3, OuterResult.Exit, true, false); // positive, exiting
+      ListState.Length += 2; 
+    }
+    else if (RadiusHasMinimumBound)
+    {
+      RayIntersections InnerResult = IntersectInfiniteCylinder(R, MinBounds.x, false);
+      setShapeIntersections(Intersections, ListState, 1, InnerResult);
+      ListState.Length += 2;
+    }
+
+    float2 AngleBounds = float2(MinBounds.y, MaxBounds.y);
+    if (AngleRangeFlag == ANGLE_UNDER_HALF)
+    {
+      // The shape's angle range is under half, so we intersect a NEGATIVE shape that is over half.
+      RayIntersections FirstResult = (RayIntersections) 0;
+      RayIntersections SecondResult = (RayIntersections) 0;
+      IntersectFlippedWedge(R, AngleBounds, FirstResult, SecondResult);
+      switch (ListState.Length)
+      {
+        case 4: // outer cylinder + inner radius
+          setShapeIntersections(Intersections, ListState, 2, FirstResult);
+          setShapeIntersections(Intersections, ListState, 3, SecondResult);
+          break;
+        case 2: // outer cylinder
+        default:
+          setShapeIntersections(Intersections, ListState, 1, FirstResult);
+          setShapeIntersections(Intersections, ListState, 2, SecondResult);
+          break;
+      }
+      ListState.Length += 4;
+    }
+    else if (AngleRangeFlag == ANGLE_OVER_HALF)
+    {
+      // The shape's angle range is over half, so we intersect a NEGATIVE shape that is under half.
+      RayIntersections WedgeResult = IntersectRegularWedge(R, AngleBounds);
+      switch (ListState.Length)
+      {
+        case 4: // outer cylinder + inner radius
+          setShapeIntersections(Intersections, ListState, 2, WedgeResult);
+          break;
+        case 2: // outer cylinder
+        default:
+          setShapeIntersections(Intersections, ListState, 1, WedgeResult);
+          break;
+      }
+      ListState.Length += 2;
+    }
+    else if (AngleRangeFlag == ANGLE_EQUAL_ZERO)
+    {
+      RayIntersections FirstResult = (RayIntersections) 0;
+      RayIntersections SecondResult = (RayIntersections) 0;
+      IntersectHalfPlane(R, AngleBounds.x, FirstResult, SecondResult);
+      switch (ListState.Length)
+      {
+        case 4: // outer cylinder + inner radius
+          setShapeIntersections(Intersections, ListState, 2, FirstResult);
+          setShapeIntersections(Intersections, ListState, 3, SecondResult);
+          break;
+        case 2: // outer cylinder
+        default:
+          setShapeIntersections(Intersections, ListState, 1, FirstResult);
+          setShapeIntersections(Intersections, ListState, 2, SecondResult);
+          break;
+      }
+      ListState.Length += 4;
+    }
+  }
+
+  /**
+  * Scales the input UV coordinates from [0, 1] to their values in UV Shape Space.  
+  */
+  float3 ScaleUVToShapeUVSpace(in float3 UV)
+  {
+    float radius = UV.x;
+    if (RadiusHasMinimumBound)
+    {
+      radius /= RadiusUVScale;
+    }
+
+    // Convert from [0, 1] to radians [-pi, pi]
+    float angle = UV.y * CZM_TWO_PI;
+    if (AngleRangeFlag > 0)
+    {
+      angle /= AngleUVScale;
+    }
+
+    return float3(radius, angle, UV.z);
+  }
+
+ /**
+   * Converts the input position (vanilla UV Space) to its Shape UV Space relative to the
+   * ellipsoid geometry. Also outputs the Jacobian transpose for future use.
+   */
+  float3 ConvertUVToShapeUVSpace(in float3 PositionUV, out float3x3 JacobianT)
+  {
+    // First convert UV to "Unit Space derivative".
+    float3 unitPosition = UVToUnit(PositionUV);
+
+    float radius = length(unitPosition.xy); // [0, 1]
+    float3 radial = normalize(float3(unitPosition.xy, 0.0));
+
+    // Shape space height is defined within [0, 1]
+    float height = PositionUV.z; // [0, 1]
+    float3 z = float3(0.0, 0.0, 1.0);
+
+    float angle = atan2(unitPosition.y, unitPosition.x);
+    float3 east = normalize(float3(-unitPosition.y, unitPosition.x, 0.0));
+
+    JacobianT = float3x3(radial, z, east / length(unitPosition.xy));
+
+    // Then convert Unit Space to Shape UV Space.
+    // Radius: shift and scale to [0, 1]
+    float radiusUV = radius;
+    if (RadiusHasMinimumBound)
+    {
+      radiusUV = radiusUV * RadiusUVScale + RadiusUVOffset;
+    }
+
+    // Angle: shift and scale to [0, 1]
+    float angleUV = (angle + CZM_PI) / CZM_TWO_PI;
+
+    if (AngleRangeFlag > 0)
+    {
+      // Correct the angle when max < min
+      // Technically this should compare against min angle, but it has precision problems so compare against the middle of empty space.
+      if (AngleIsReversed)
+      {
+        angleUV += float(angleUV < AngleUVExtents.z);
+      }
+
+      // Avoid flickering from reading voxels from both sides of the -pi/+pi discontinuity.
+      if (AngleMinAtDiscontinuity)
+      {
+        angleUV = angleUV > AngleUVExtents.z ? AngleUVExtents.x : angleUV;
+      }
+      else if (AngleMaxAtDiscontinuity)
+      {
+        angleUV = angleUV < AngleUVExtents.z ? AngleUVExtents.y : angleUV;
+      }
+
+      angleUV = angleUV * AngleUVScale + AngleUVOffset;
+    }
+
+    return float3(radiusUV, angleUV, height);
+  }
+};