From 2fbf4287699f3a8bb9f1e3ac1d5dc0a67b72cb73 Mon Sep 17 00:00:00 2001
From: Janine Liu <janine@cesium.com>
Date: Mon, 10 Mar 2025 16:55:18 -0400
Subject: [PATCH 1/8] Add support for ellipsoid voxel rendering

---
 Content/Materials/CesiumVoxelTemplate.hlsl    | 877 +++++++++++++++++-
 .../Private/CesiumVoxelRendererComponent.cpp  | 364 +++++++-
 2 files changed, 1236 insertions(+), 5 deletions(-)

diff --git a/Content/Materials/CesiumVoxelTemplate.hlsl b/Content/Materials/CesiumVoxelTemplate.hlsl
index fc1f12eb2..a880c2f75 100644
--- a/Content/Materials/CesiumVoxelTemplate.hlsl
+++ b/Content/Materials/CesiumVoxelTemplate.hlsl
@@ -363,6 +363,26 @@ struct ShapeUtility
   float3 MinBounds;
   float3 MaxBounds;
 
+  // Ellipsoid variables
+  float3 RadiiUV;
+  float3 InverseRadiiUVSquared;
+  float InverseHeightDifferenceUV;
+
+  float EccentricitySquared;
+  float2 EvoluteScale;
+  float3 LongitudeUVExtents; // X = min, Y = max, Z = midpoint
+  float LongitudeUVScale;
+  float LongitudeUVOffset;
+  float LatitudeUVScale;
+  float LatitudeUVOffset;
+  
+  uint LongitudeRangeFlag;
+  bool LongitudeMinHasDiscontinuity;
+  bool LongitudeMaxHasDiscontinuity;
+  bool LongitudeIsReversed;
+  uint LatitudeMinFlag;
+  uint LatitudeMaxFlag;
+
   /**
   * Converts a position from Unit Shape Space coordinates to UV Space.
   * [-1, -1] => [0, 1]
@@ -395,6 +415,31 @@ struct ShapeUtility
       MinBounds = float3(-1, -1, -1);
       MaxBounds = float3(1, 1, 1);
     }
+    else if (ShapeConstant == ELLIPSOID)
+    {
+      MinBounds = InMinBounds; // longitude, sine(latitude), height
+      MaxBounds = InMaxBounds; // longitude, sine(latitude), height
+
+      // Flags are packed in CesiumVoxelRendererComponent.cpp
+      LongitudeRangeFlag = round(PackedData0.x);
+      LongitudeMinHasDiscontinuity = bool(PackedData0.y);
+      LongitudeMaxHasDiscontinuity = bool(PackedData0.z);
+      LongitudeIsReversed = bool(PackedData0.w);
+      
+      LatitudeMinFlag = round(PackedData1.x);
+      LatitudeMaxFlag = round(PackedData1.y);
+
+      EvoluteScale = PackedData1.zw;
+      RadiiUV = PackedData2.xyz;
+      InverseRadiiUVSquared = PackedData3.xyz;
+      InverseHeightDifferenceUV = PackedData3.w;
+      EccentricitySquared = PackedData4.w;
+      LongitudeUVExtents = PackedData4.xyz;
+      LongitudeUVScale = PackedData5.x;
+      LongitudeUVOffset = PackedData5.y;
+      LatitudeUVScale = PackedData5.z;
+      LatitudeUVOffset = PackedData5.w;
+    }
   }
   
   /**
@@ -440,6 +485,646 @@ struct ShapeUtility
     
     setShapeIntersections(Intersections, ListState, 0, result);
   }
+
+  /**
+  * Tests whether the input ray intersects the ellipsoid at the given height, relative to original radii.
+  */
+  RayIntersections IntersectEllipsoidAtHeight(in Ray R, in float RelativeHeight, in bool IsConvex)
+  {
+    // Scale the ray by the ellipsoid axes to make it a unit sphere.
+    // Note: This approximates ellipsoid + height as an ellipsoid.
+    float3 radiiCorrection = RadiiUV / (RadiiUV + RelativeHeight);
+    float3 position = R.Origin * radiiCorrection;
+    float3 direction = R.Direction * radiiCorrection;
+
+    // Intersect the ellipsoid defined by x^2/a^2 + y^2/b^2 + z^2/c^2 = 1.
+    // Solve for the values of t where the ray intersects the ellipsoid, by substituting
+    // the ray equation o + t * d into the ellipsoid x/y values. Results in a quadratic
+    // equation on t where we can find the determinant.
+    float a = dot(direction, direction); // ~ 1.0 (or maybe 4.0 if ray is scaled)
+    float b = dot(direction, position); // roughly inside [-1.0, 1.0] when zoomed in
+    float c = dot(position, position) - 1.0; // ~ 0.0 when zoomed in.
+    float determinant = b * b - a * c; // ~ b * b when zoomed in
+
+    if (determinant < 0.0)
+    {
+      // Output a "normal" so that the shape entry information can be encoded.
+      // (See EncodeIntersection)
+      Intersection miss = Utils.NewIntersection(NO_HIT, normalize(direction));
+      return Utils.NewRayIntersections(miss, miss);
+    }
+
+    determinant = sqrt(determinant);
+    RayIntersections result = (RayIntersections) 0;
+
+    // Compute the larger t using standard formula.
+    // The other t may suffer from subtractive cancellation in the standard formula,
+    // so it is computed from the first t instead.
+    float t1 = (-b - sign(b) * determinant) / a;
+    float t2 = c / (a * t1);
+    result.Entry.t = min(t1, t2);
+    result.Exit.t = max(t1, t2);
+  
+    float directionScale = IsConvex ? 1.0 : -1.0;
+    result.Entry.Normal = directionScale * normalize(position + result.Entry.t * direction);
+    result.Exit.Normal = directionScale * normalize(position + result.Exit.t * direction);
+    return result;
+  }
+  
+  /**
+  * Intersects the plane at the specified longitude angle.
+  */
+  Intersection IntersectLongitudePlane(in Ray R, in float Angle, in bool PositiveNormal)
+  {
+    float normalSign = PositiveNormal ? 1.0 : -1.0;
+    float2 planeNormal = normalSign * float2(-sin(Angle), cos(Angle));
+
+    float approachRate = dot(R.Direction.xy, planeNormal);
+    float distance = -dot(R.Origin.xy, planeNormal);
+
+    float t = (approachRate == 0.0)
+        ? NO_HIT
+        : distance / approachRate;
+
+    return Utils.NewIntersection(t, float3(planeNormal, 0));
+  }
+  
+  /**
+  * Intersects the space on one side of the specified longitude angle.
+  */
+  RayIntersections IntersectHalfSpace(in Ray R, in float Angle, in bool PositiveNormal)
+  {
+    Intersection intersection = IntersectLongitudePlane(R, Angle, PositiveNormal);
+    Intersection farSide = Utils.NewIntersection(INF_HIT, normalize(R.Direction));
+
+    bool hitFront = (intersection.t > 0.0) == (dot(R.Origin.xy, intersection.Normal.xy) > 0.0);
+    if (!hitFront)
+    {
+      return Utils.NewRayIntersections(intersection, farSide);
+    }
+    else
+    {
+      return Utils.NewRayIntersections(Utils.Multiply(farSide, -1), intersection);
+    }
+  }
+
+  /**
+  * Intersects a "flipped" wedge formed by the negative space of the specified angle
+  * bounds and returns up to four intersections. The "flipped" wedge is the union of
+  * two half-spaces defined at the given angles and represents a *negative* volume
+  * of over > 180 degrees.
+  */
+  void IntersectFlippedWedge(in Ray R, in float2 AngleBounds, out RayIntersections FirstIntersection, out RayIntersections SecondIntersection)
+  {
+    FirstIntersection = IntersectHalfSpace(R, AngleBounds.x, false);
+    SecondIntersection = IntersectHalfSpace(R, AngleBounds.y, true);
+  }
+
+  /**
+  * Intersects the wedge formed by the negative space of the min/max longitude, where
+  * maxAngle > minAngle + pi. The wedge is represented by two planes at such angles.
+  * There is an opposite "shadow wedge", i.e. the wedge formed at the *opposite* side
+  * of the planes' intersection, that must be specially handled.
+  */
+  RayIntersections IntersectRegularWedge(in Ray R, in float2 AngleBounds)
+  {
+    // Normals will point toward the "outside" (into the negative space)
+    Intersection intersect1 = IntersectLongitudePlane(R, AngleBounds.x, false);
+    Intersection intersect2 = IntersectLongitudePlane(R, AngleBounds.y, true);
+
+    // Note: the intersections could be in the "shadow" wedge, beyond the tip of
+    // the actual wedge.
+    Intersection first = intersect1;
+    Intersection last = intersect2;
+    if (first.t > last.t)
+    {
+      first = intersect2;
+      last = intersect1;
+    }
+
+    bool firstIntersectionAheadOfRay = first.t >= 0.0;
+    bool startedInsideFirst = dot(R.Origin.xy, first.Normal.xy) < 0.0;
+    bool isExitingFromInside = firstIntersectionAheadOfRay == startedInsideFirst;
+    bool lastIntersectionAheadOfRay = last.t > 0.0;
+    bool startedOutsideLast = dot(R.Origin.xy, last.Normal.xy) >= 0.0;
+    bool isEnteringFromOutside = lastIntersectionAheadOfRay == startedOutsideLast;
+
+    Intersection farSide = Utils.NewIntersection(INF_HIT, normalize(R.Direction));
+    Intersection miss = Utils.NewIntersection(NO_HIT, normalize(R.Direction));
+
+    if (isExitingFromInside && isEnteringFromOutside)
+    {
+      // Ray crosses both faces of negative wedge, exiting then entering the positive shape
+      return Utils.NewRayIntersections(first, last);
+    }
+    else if (!isExitingFromInside && isEnteringFromOutside)
+    {
+      // Ray starts inside wedge. last is in shadow wedge, and first is actually the entry
+      return Utils.NewRayIntersections(Utils.Multiply(farSide, -1), first);
+    }
+    else if (isExitingFromInside && !isEnteringFromOutside)
+    {
+      // First intersection was in the shadow wedge, so last is actually the exit
+      return Utils.NewRayIntersections(last, farSide);
+    }
+    else
+    { // !exitFromInside && !enterFromOutside
+      // Both intersections were in the shadow wedge
+      return Utils.NewRayIntersections(miss, miss);
+    }
+  }
+  
+  bool HitsPositiveHalfPlane(in Ray R, in Intersection InputIntersection, in bool positiveNormal)
+  {
+    float normalSign = positiveNormal ? 1.0 : -1.0;
+    float2 planeDirection = float2(InputIntersection.Normal.y, -InputIntersection.Normal.x) * normalSign;
+    float2 hit = R.Origin.xy + InputIntersection.t * R.Direction.xy;
+    return dot(hit, planeDirection) > 0.0;
+  }
+
+  void IntersectHalfPlane(in Ray R, in float Angle, out RayIntersections FirstIntersection, out RayIntersections SecondIntersection)
+  {
+    Intersection intersection = IntersectLongitudePlane(R, Angle, true);
+    Intersection farSide = Utils.NewIntersection(INF_HIT, normalize(R.Direction));
+
+    if (HitsPositiveHalfPlane(R, intersection, true))
+    {
+      FirstIntersection.Entry = Utils.Multiply(farSide, -1);
+      FirstIntersection.Exit = Utils.NewIntersection(intersection.t, float3(-1.0 * intersection.Normal.xy, 0.0));
+      SecondIntersection.Entry = intersection;
+      SecondIntersection.Exit = farSide;
+    }
+    else
+    {
+      Intersection miss = Utils.NewIntersection(NO_HIT, normalize(R.Direction));
+      FirstIntersection.Entry = Utils.Multiply(farSide, -1);
+      FirstIntersection.Exit = farSide;
+      SecondIntersection.Entry = miss;
+      SecondIntersection.Exit = miss;
+    }
+  }
+
+  /**
+  * Given a circular quadric cone around the z-axis, with apex at the origin,
+  * find the parametric distance(s) along a ray where that ray intersects
+  * the cone. Computation of the normal is deferred to GetFlippedCone.
+  * 
+  * The cone opening angle is described by the squared cosine of its half-angle
+  * (the angle between the Z-axis and the surface)
+  */
+  RayIntersections IntersectQuadricCone(in Ray R, in float cosSquaredHalfAngle)
+  {
+    float3 o = R.Origin;
+    float3 d = R.Direction;
+    
+    float sinSquaredHalfAngle = 1.0 - cosSquaredHalfAngle;
+
+    float aSin = d.z * d.z * sinSquaredHalfAngle;
+    float aCos = -dot(d.xy, d.xy) * cosSquaredHalfAngle;
+    float a = aSin + aCos;
+
+    float bSin = d.z * o.z * sinSquaredHalfAngle;
+    float bCos = -dot(o.xy, d.xy) * cosSquaredHalfAngle;
+    float b = bSin + bCos;
+
+    float cSin = o.z * o.z * sinSquaredHalfAngle;
+    float cCos = -dot(o.xy, o.xy) * cosSquaredHalfAngle;
+    float c = cSin + cCos;
+    // determinant = b * b - a * c. But bSin * bSin = aSin * cSin.
+    // Avoid subtractive cancellation by expanding to eliminate these terms
+    float determinant = 2.0 * bSin * bCos + bCos * bCos - aSin * cCos - aCos * cSin - aCos * cCos;
+
+    RayIntersections result = (RayIntersections) 0;
+    result.Entry.t = NO_HIT;
+    result.Exit.t = NO_HIT;
+
+    if (determinant < 0.0)
+    {
+      return result;
+    }
+
+    if (a == 0.0)
+    {
+      // Ray is parallel to cone surface if b == 0.0.
+      // Otherwise, ray has a single intersection on cone surface
+      if (b != 0.0)
+      {
+        result.Entry.t = -0.5 * c / b;
+      }
+    }
+
+    determinant = sqrt(determinant);
+
+    // Compute larger root using standard formula
+    float signB = b < 0.0 ? -1.0 : 1.0;
+    float t1 = (-b - signB * determinant) / a;
+    // The other root may suffer from subtractive cancellation in the standard formula.
+    // Compute it from the first root instead.
+    float t2 = c / (a * t1);
+    result.Entry.t = min(t1, t2);
+    result.Exit.t = max(t1, t2);
+    return result;
+  }
+
+  /**
+   * Given a point on a conical surface, find the surface normal at that point.
+   */
+  float3 GetConeNormal(in float3 Point, in bool IsConvex)
+  {
+    // Start with radial component pointing toward z-axis
+    float2 radial = -abs(Point.z) * normalize(Point.xy);
+    // Z component points toward opening of cone
+    float zSign = (Point.z < 0.0) ? -1.0 : 1.0;
+    float z = length(Point.xy) * zSign;
+    // Flip normal if the shape is convex
+    float flip = (IsConvex) ? -1.0 : 1.0;
+    return normalize(float3(radial, z) * flip);
+  }
+
+  /**
+   * Compute the shift between the ellipsoid origin and the apex of a cone of latitude
+   */
+  float GetLatitudeConeShift(in float sinLatitude)
+  {
+    // Find prime vertical radius of curvature: 
+    // the distance along the ellipsoid normal to the intersection with the z-axis
+    float x2 = EccentricitySquared * sinLatitude * sinLatitude;
+    float primeVerticalRadius = rsqrt(1.0 - x2);
+
+    // Compute a shift from the origin to the intersection of the cone with the z-axis
+    return primeVerticalRadius * EccentricitySquared * sinLatitude;
+  }
+
+  /**
+  * Tests if the input ray intersects the upper half of the quadric cone, i.e., the cone where its apex
+  * points downward.
+  */ 
+  void IntersectFlippedCone(in Ray R, in float CosHalfAngle, out RayIntersections FirstIntersections, out RayIntersections SecondIntersections)
+  {
+    // Undo the scaling from ellipsoid to sphere
+    float3 position = R.Origin * RadiiUV;
+    float3 direction = R.Direction * RadiiUV;
+    // Shift the ray to account for the latitude cone not being centered at the Earth center
+    position.z += GetLatitudeConeShift(CosHalfAngle);
+
+    RayIntersections quadricIntersections = IntersectQuadricCone(R, CosHalfAngle * CosHalfAngle);
+    
+    Intersection miss = Utils.NewIntersection(NO_HIT, normalize(direction));
+    Intersection farSide = Utils.NewIntersection(INF_HIT, normalize(direction));
+
+    // Initialize output with no intersections
+    FirstIntersections = (RayIntersections) 0;
+    FirstIntersections.Entry = Utils.Multiply(farSide, -1);
+    FirstIntersections.Exit = farSide;
+    
+    SecondIntersections = Utils.NewRayIntersections(miss, miss);
+    
+    if (quadricIntersections.Entry.t == NO_HIT)
+    {
+      return;
+    }
+    
+    // Find the points of intersection
+    float3 p0 = position + quadricIntersections.Entry.t * direction;
+    float3 p1 = position + quadricIntersections.Exit.t * direction;
+
+    quadricIntersections.Entry.Normal = GetConeNormal(p0, true);
+    quadricIntersections.Exit.Normal = GetConeNormal(p1, true);
+
+    // shadow cone = the half of the quadric cone that is being discarded.
+    bool p0InShadowCone = sign(p0.z) != sign(CosHalfAngle);
+    bool p1InShadowCone = sign(p1.z) != sign(CosHalfAngle);
+
+    if (p0InShadowCone && p1InShadowCone)
+    {
+      // both points in shadow cone; no valid intersections
+      return;
+    }
+    else if (p0InShadowCone)
+    {
+      FirstIntersections.Exit = quadricIntersections.Exit;
+    }
+    else if (p1InShadowCone)
+    {
+      FirstIntersections.Entry = quadricIntersections.Entry;
+    }
+    else
+    {
+      FirstIntersections.Exit = quadricIntersections.Entry;
+      SecondIntersections.Entry = quadricIntersections.Exit;
+      SecondIntersections.Exit = farSide;
+    }
+  }
+  
+  /**
+  * Tests if the input ray intersects the upper half of the quadric cone, i.e., the cone where its apex
+  * points downward.
+  */ 
+  RayIntersections IntersectUprightCone(in Ray R, in float CosHalfAngle)
+  {
+    // Undo the scaling from ellipsoid to sphere
+    float3 position = R.Origin * RadiiUV;
+    float3 direction = R.Direction * RadiiUV;
+    // Shift the ray to account for the latitude cone not being centered at the Earth center
+    position.z += GetLatitudeConeShift(CosHalfAngle);
+
+    RayIntersections quadricIntersections = IntersectQuadricCone(R, CosHalfAngle * CosHalfAngle);
+    
+    Intersection miss = Utils.NewIntersection(NO_HIT, normalize(direction));
+    Intersection farSide = Utils.NewIntersection(INF_HIT, normalize(direction));
+    
+    if (quadricIntersections.Entry.t == NO_HIT)
+    {
+      return Utils.NewRayIntersections(miss, miss);
+    }
+    
+    // Find the points of intersection
+    float3 p0 = position + quadricIntersections.Entry.t * direction;
+    float3 p1 = position + quadricIntersections.Exit.t * direction;
+
+    quadricIntersections.Entry.Normal = GetConeNormal(p0, false);
+    quadricIntersections.Exit.Normal = GetConeNormal(p1, false);
+
+    bool p0InShadowCone = sign(p0.z) != sign(CosHalfAngle);
+    bool p1InShadowCone = sign(p1.z) != sign(CosHalfAngle);
+
+    if (p0InShadowCone && p1InShadowCone)
+    {
+      return Utils.NewRayIntersections(miss, miss);
+    }
+    else if (p0InShadowCone)
+    {
+      return Utils.NewRayIntersections(quadricIntersections.Exit, farSide);
+    }
+    else if (p1InShadowCone)
+    {
+      return Utils.NewRayIntersections(Utils.Multiply(farSide, -1), quadricIntersections.Entry);
+    }
+    else
+    {
+      return Utils.NewRayIntersections(quadricIntersections.Entry, quadricIntersections.Exit);
+    }
+  }
+
+  /**
+  * Intersects an XY-plane at Z = 0. The given Z describes the direction of plane's normal.
+  */
+  RayIntersections IntersectZPlane(in Ray R, in float Z)
+  {
+    float t = -R.Origin.z / R.Direction.z;
+
+    bool startsOutside = sign(R.Origin.z) == sign(Z);
+    bool isEntering = (t >= 0.0) != startsOutside;
+
+    Intersection intersect = Utils.NewIntersection(t, float3(0.0, 0.0, Z));
+    Intersection farSide = Utils.NewIntersection(INF_HIT, normalize(R.Direction));
+
+    if (isEntering)
+    {
+      return Utils.NewRayIntersections(intersect, farSide);
+    }
+    else
+    {
+      return Utils.NewRayIntersections(Utils.Multiply(farSide, -1), intersect);
+    }
+  }
+  
+#define ANGLE_EQUAL_ZERO 1
+#define ANGLE_UNDER_HALF 2
+#define ANGLE_HALF 3
+#define ANGLE_OVER_HALF 4
+
+  /**
+  * Tests whether the input ray (Unit Space) intersects the ellipsoid region. Outputs the intersections in Unit Space.
+  */
+  void IntersectEllipsoid(in Ray R, inout float4 Intersections[INTERSECTIONS_LENGTH])
+  {
+    ListState.Length = 2;
+
+    // The region can be thought of as the space between two ellipsoids.
+    RayIntersections OuterResult = IntersectEllipsoidAtHeight(R, MaxBounds.z, true);
+    setShapeIntersections(Intersections, ListState, 0, OuterResult);
+    if (OuterResult.Entry.t == NO_HIT)
+    {
+      return;
+    }
+
+    ListState.Length += 2;
+    RayIntersections InnerResult = IntersectEllipsoidAtHeight(R, MinBounds.z, false);
+    if (InnerResult.Entry.t == NO_HIT)
+    {
+      setShapeIntersections(Intersections, ListState, 1, InnerResult);
+    }
+    else
+    {
+      // When the ellipsoid is large and thin it's possible for floating point math
+      // to cause the ray to intersect the inner ellipsoid before the outer ellipsoid. 
+      // To prevent this from happening, clamp the inner intersections to the outer ones
+      // and sandwich the inner ellipsoid intersection inside the outer ellipsoid 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 ellipsoid to be invisible because it will think the ray
+      // is still inside the inner (negative) ellipsoid after exiting the
+      // outer (positive) ellipsoid.
+      //
+      // 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.
+      InnerResult.Entry.t = max(InnerResult.Entry.t, OuterResult.Entry.t);
+      InnerResult.Exit.t = min(InnerResult.Exit.t, OuterResult.Exit.t);
+      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
+    }
+
+    // For min latitude, intersect a NEGATIVE cone at the bottom, based on the latitude value.
+    if (LatitudeMinFlag == ANGLE_UNDER_HALF)
+    {
+      ListState.Length += 2;
+      RayIntersections BottomConeResult = IntersectUprightCone(R, MinBounds.y);
+      setShapeIntersections(Intersections, ListState, 2, BottomConeResult);
+    }
+    else if (LatitudeMinFlag == ANGLE_HALF)
+    {
+      ListState.Length += 2;
+      RayIntersections PlaneResult = IntersectZPlane(R, -1.0);
+      setShapeIntersections(Intersections, ListState, 2, PlaneResult);
+    }
+    else if (LatitudeMinFlag == ANGLE_OVER_HALF)
+    {
+      ListState.Length += 4;
+      RayIntersections FirstIntersection;
+      RayIntersections SecondIntersection;
+      IntersectFlippedCone(R, MinBounds.y, FirstIntersection, SecondIntersection);
+      setShapeIntersections(Intersections, ListState, 2, FirstIntersection);
+      setShapeIntersections(Intersections, ListState, 3, SecondIntersection);
+    }
+
+    // For max latitude, intersect a NEGATIVE cone at the top, based on the latitude value.
+    if (LatitudeMaxFlag == ANGLE_UNDER_HALF)
+    {
+      RayIntersections FirstIntersection;
+      RayIntersections SecondIntersection;
+      IntersectFlippedCone(R, MaxBounds.y, FirstIntersection, SecondIntersection);
+      // The array index depends on how many intersections were previously found.
+      [branch]
+      switch (ListState.Length)
+      {
+        case 6: // maxHeight + minHeight + minLatitude <= half
+          setShapeIntersections(Intersections, ListState, 3, FirstIntersection);
+          setShapeIntersections(Intersections, ListState, 4, SecondIntersection);
+          break;
+        case 8: // maxHeight + minHeight + minLatitude > half
+          // It makes no sense for min latitude to be in the top half, AND max latitude in the bottom half.
+          // The region is invalid, so make things easier and mark the shape as not hit.
+          Intersections[0].w = NO_HIT;
+          Intersections[1].w = NO_HIT;
+          ListState.Length = 0;
+          return;
+        case 4: // maxHeight + minHeight
+        default:
+          setShapeIntersections(Intersections, ListState, 2, FirstIntersection);
+          setShapeIntersections(Intersections, ListState, 3, SecondIntersection);
+          break;
+      }
+      ListState.Length += 4;
+    }
+    else if (LatitudeMaxFlag == ANGLE_HALF)
+    {
+      RayIntersections PlaneResult = IntersectZPlane(R, 1.0);
+      // The array index depends on how many intersections were previously found.
+      [branch]
+      switch (ListState.Length)
+      {
+        case 6: // maxHeight + minHeight + minLatitude <= half
+          // Not worth checking if minLatitude == half, just let the sort algorithm figure it out
+          setShapeIntersections(Intersections, ListState, 3, PlaneResult);
+          break;
+        case 8: // maxHeight + minHeight + minLatitude > half
+          // It makes no sense for min latitude to be over half, AND max latitude to be half.
+          // The region is invalid, so make things easier and mark the shape as not hit.
+          Intersections[0].w = NO_HIT;
+          Intersections[1].w = NO_HIT;
+          ListState.Length = 0;
+          return;
+        case 4: // maxHeight + minHeight
+        default:
+          setShapeIntersections(Intersections, ListState, 2, PlaneResult);
+          break;
+      }
+      ListState.Length += 2;
+    }
+    else if (LatitudeMaxFlag == ANGLE_OVER_HALF)
+    {
+      RayIntersections TopConeResult = IntersectUprightCone(R, MaxBounds.y);
+      [branch]
+      switch (ListState.Length)
+      {
+        case 6: // maxHeight + minHeight + minLatitude <= half
+          setShapeIntersections(Intersections, ListState, 3, TopConeResult);
+          break;
+        case 8: // maxHeight + minHeight + minLatitude > half
+          setShapeIntersections(Intersections, ListState, 4, TopConeResult);
+          break;
+        case 4: // maxHeight + minHeight
+        default:
+          setShapeIntersections(Intersections, ListState, 2, TopConeResult);
+          break;
+      }
+      ListState.Length += 2;
+    }
+    
+    float2 LongitudeBounds = float2(MinBounds.x, MaxBounds.x);
+    if (LongitudeRangeFlag == ANGLE_OVER_HALF)
+    {
+      // The shape's longitude range is over half, so we intersect a NEGATIVE shape that is under half.
+      RayIntersections WedgeResult = IntersectRegularWedge(R, LongitudeBounds);
+      // The array index depends on how many intersections were previously found.
+      [branch]
+      switch (ListState.Length)
+      {
+        case 6: // maxHeight + minHeight + (minLatitude <= half || maxLatitude >= half)
+          setShapeIntersections(Intersections, ListState, 3, WedgeResult);
+          break;
+        case 8: // maxHeight + minHeight + (minLatitude > half || maxLatitude < half || minLatitude <= half && maxLatitude >= half)
+          setShapeIntersections(Intersections, ListState, 4, WedgeResult);
+          break;
+        case 10: // maxHeight + minHeight + (minLatitude > half && maxLatitude > half || minLatitude < half + maxLatitude < half)
+          setShapeIntersections(Intersections, ListState, 5, WedgeResult);
+          break;
+        case 4: // maxHeight + minHeight
+        default:
+          setShapeIntersections(Intersections, ListState, 2, WedgeResult);
+          break;
+      }
+      ListState.Length += 2;
+    }
+    else if (LongitudeRangeFlag == ANGLE_UNDER_HALF)
+    {
+      // The shape's longitude range is under half, so we intersect a NEGATIVE shape that is over half.
+      RayIntersections FirstResult = (RayIntersections) 0;
+      RayIntersections SecondResult = (RayIntersections) 0;
+      IntersectFlippedWedge(R, LongitudeBounds, FirstResult, SecondResult);
+      // The array index depends on how many intersections were previously found.
+      [branch]
+      switch (ListState.Length)
+      {
+        case 6: // maxHeight + minHeight + (minLatitude <= half || maxLatitude >= half)
+          setShapeIntersections(Intersections, ListState, 3, FirstResult);
+          setShapeIntersections(Intersections, ListState, 4, SecondResult);
+          break;
+        case 8: // maxHeight + minHeight + (minLatitude > half || maxLatitude < half || minLatitude <= half && maxLatitude >= half)
+          setShapeIntersections(Intersections, ListState, 4, FirstResult);
+          setShapeIntersections(Intersections, ListState, 5, SecondResult);
+          break;
+        case 10: // maxHeight + minHeight + (minLatitude > half && maxLatitude > half || minLatitude < half + maxLatitude < half)
+          setShapeIntersections(Intersections, ListState, 5, FirstResult);
+          setShapeIntersections(Intersections, ListState, 6, SecondResult);
+          break;
+        case 4: // maxHeight + minHeight
+        default:
+          setShapeIntersections(Intersections, ListState, 2, FirstResult);
+          setShapeIntersections(Intersections, ListState, 3, SecondResult);
+          break;
+      }
+      ListState.Length += 4;
+    }
+    else if (LongitudeRangeFlag == ANGLE_EQUAL_ZERO)
+    {
+      RayIntersections FirstResult = (RayIntersections) 0;
+      RayIntersections SecondResult = (RayIntersections) 0;
+      IntersectHalfPlane(R, LongitudeBounds.x, FirstResult, SecondResult);
+      // The array index depends on how many intersections were previously found.
+      [branch]
+      switch (ListState.Length)
+      {
+        case 6: // maxHeight + minHeight + (minLatitude <= half || maxLatitude >= half)
+          setShapeIntersections(Intersections, ListState, 3, FirstResult);
+          setShapeIntersections(Intersections, ListState, 4, SecondResult);
+          break;
+        case 8: // maxHeight + minHeight + (minLatitude > half || maxLatitude < half || minLatitude <= half && maxLatitude >= half)
+          setShapeIntersections(Intersections, ListState, 4, FirstResult);
+          setShapeIntersections(Intersections, ListState, 5, SecondResult);
+          break;
+        case 10: // maxHeight + minHeight + (minLatitude > half && maxLatitude > half || minLatitude < half + maxLatitude < half)
+          setShapeIntersections(Intersections, ListState, 5, FirstResult);
+          setShapeIntersections(Intersections, ListState, 6, SecondResult);
+          break;
+        case 4: // maxHeight + minHeight
+        default:
+          setShapeIntersections(Intersections, ListState, 2, FirstResult);
+          setShapeIntersections(Intersections, ListState, 3, SecondResult);
+          break;
+      }
+      ListState.Length += 4;
+    }
+  }
   
   /**
   * Tests whether the input ray (Unit Space) intersects the shape.
@@ -453,11 +1138,43 @@ struct ShapeUtility
       case BOX:
         IntersectBox(R, Intersections);
         break;
+      case ELLIPSOID:
+        IntersectEllipsoid(R, Intersections);
+        break;
       default:
         return Utils.NewRayIntersections(Utils.NewIntersection(NO_HIT, 0), Utils.NewIntersection(NO_HIT, 0));
     }
     
     RayIntersections result = ListState.GetFirstIntersections(Intersections);
+    if (result.Entry.t == NO_HIT)
+    {
+      return result;
+    }
+    
+    // Don't bother with sorting if the positive shape was missed.
+    // Box intersection is straightforward and doesn't require sorting.
+    // Currently, cylinders do not require sorting, but they will once clipping / exaggeration
+    // is supported.
+    if (ShapeConstant == ELLIPSOID)
+    {
+      ListState.Sort(Intersections);
+      for (int i = 0; i < SHAPE_INTERSECTIONS; ++i)
+      {
+        result = ListState.GetNextIntersections(Intersections);
+        if (result.Exit.t > 0.0)
+        {
+          // Set start to 0.0 when ray is inside the shape.
+          result.Entry.t = max(result.Entry.t, 0.0);
+          break;
+        }
+      }
+    }
+    else
+    {
+      // Set start to 0.0 when ray is inside the shape.
+      result.Entry.t = max(result.Entry.t, 0.0);
+    }
+
     return result;
   }
   
@@ -466,7 +1183,32 @@ struct ShapeUtility
   */
   float3 ScaleUVToShapeUVSpace(in float3 UV)
   {
-    return UV;
+    if (ShapeConstant == ELLIPSOID)
+    {
+      // Convert from [0, 1] to radians [-pi, pi]
+      float longitude = UV.x * CZM_TWO_PI;
+      if (LongitudeRangeFlag > 0)
+      {
+        longitude /= LongitudeUVScale;
+      }
+
+      // Convert from [0, 1] to radians [-pi/2, pi/2]
+      float latitude = UV.y * CZM_PI;
+      if (LatitudeMinFlag > 0 || LatitudeMaxFlag > 0)
+      {
+        latitude /= LatitudeUVScale;
+      }
+    
+      float height = UV.z / InverseHeightDifferenceUV;
+
+      return float3(longitude, latitude, height);
+    }
+    else
+    {
+      // For CYLINDER: Once clipping / exaggeration is supported, an offset + scale
+      // should be applied to its radius / height / angle. (See CesiumJS)
+      return UV;
+    }
   }
 
   /**
@@ -481,6 +1223,115 @@ struct ShapeUtility
     JacobianT = float3x3(0.5f, 0, 0, 0, 0.5f, 0, 0, 0, 0.5f);
     return PositionUV;
   }
+
+  /**
+    * Given a specified point, gets the nearest point (XY) on the ellipse using a robust
+    * iterative solution without trig functions, as well as the radius of curvature
+    * at that point (Z).
+    * https://github.com/0xfaded/ellipse_demo/issues/1
+    * https://stackoverflow.com/questions/22959698/distance-from-given-point-to-given-ellipse
+    */
+  float3 GetNearestPointAndRadiusOnEllipse(float2 pos, float2 radii)
+  {
+    float2 p = abs(pos);
+    float2 inverseRadii = 1.0 / radii;
+
+    // We describe the ellipse parametrically: v = radii * vec2(cos(t), sin(t))
+    // but store the cos and sin of t in a vec2 for efficiency.
+    // Initial guess: t = pi/4
+    float2 tTrigs = 0.7071067811865476;
+    // Initial guess of point on ellipsoid
+    float2 v = radii * tTrigs;
+    // Center of curvature of the ellipse at v
+    float2 evolute = EvoluteScale * tTrigs * tTrigs * tTrigs;
+
+    const int iterations = 3;
+    for (int i = 0; i < iterations; ++i)
+    {
+      // Find the (approximate) intersection of p - evolute with the ellipsoid.
+      float2 q = normalize(p - evolute) * length(v - evolute);
+      // Update the estimate of t.
+      tTrigs = (q + evolute) * inverseRadii;
+      tTrigs = normalize(clamp(tTrigs, 0.0, 1.0));
+      v = radii * tTrigs;
+      evolute = EvoluteScale * tTrigs * tTrigs * tTrigs;
+    }
+
+    return float3(v * sign(pos), length(v - evolute));
+  }
+  
+  /**
+   * 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 ConvertUVToShapeUVSpaceEllipsoid(in float3 PositionUV, out float3x3 JacobianT)
+  {
+    // First convert UV to "Unit Space derivative".
+    // Get the position in unit space, undo the scaling from ellipsoid to sphere
+    float3 position = UVToUnit(PositionUV);
+    position *= RadiiUV;
+
+    float longitude = atan2(position.y, position.x);
+    float3 east = normalize(float3(-position.y, position.x, 0.0));
+
+    // Convert the 3D position to a 2D position relative to the ellipse (radii.x, radii.z)
+    // (assume radii.y == radii.x) and find the nearest point on the ellipse and its normal
+    float distanceFromZAxis = length(position.xy);
+    float2 positionEllipse = float2(distanceFromZAxis, position.z);
+    float3 surfacePointAndRadius = GetNearestPointAndRadiusOnEllipse(positionEllipse, RadiiUV.xz);
+    float2 surfacePoint = surfacePointAndRadius.xy;
+
+    float2 normal2d = normalize(surfacePoint * InverseRadiiUVSquared.xz);
+    float latitude = atan2(normal2d.y, normal2d.x);
+    float3 north = float3(-normal2d.y * normalize(position.xy), abs(normal2d.x));
+
+    float heightSign = length(positionEllipse) < length(surfacePoint) ? -1.0 : 1.0;
+    float height = heightSign * length(positionEllipse - surfacePoint);
+    float3 up = normalize(cross(east, north));
+
+    JacobianT = float3x3(east / distanceFromZAxis, north / (surfacePointAndRadius.z + height), up);
+    // Transpose because HLSL matrices are constructed in row-major order.
+    JacobianT = transpose(JacobianT);
+
+    // Then convert Unit Space to Shape UV Space
+    // Longitude: shift & scale to [0, 1]
+    float longitudeUV = (longitude + CZM_PI) / CZM_TWO_PI;
+
+    // Correct the angle when max < min
+    // Technically this should compare against min longitude - but it has precision problems so compare against the middle of empty space.
+    if (LongitudeIsReversed)
+    {
+      longitudeUV += float(longitudeUV < LongitudeUVExtents.z);
+    }
+
+    // Avoid flickering from reading voxels from both sides of the -pi/+pi discontinuity.
+    if (LongitudeMinHasDiscontinuity)
+    {
+      longitudeUV = longitudeUV > LongitudeUVExtents.z ? LongitudeUVExtents.x : longitudeUV;
+    }
+
+    if (LongitudeMaxHasDiscontinuity)
+    {
+      longitudeUV = longitudeUV < LongitudeUVExtents.z ? LongitudeUVExtents.y : longitudeUV;
+    }
+
+    if (LongitudeRangeFlag > 0)
+    {
+      longitudeUV = longitudeUV * LongitudeUVScale + LongitudeUVOffset;
+    }
+
+    // Latitude: shift and scale to [0, 1]
+    float latitudeUV = (latitude + CZM_PI_OVER_TWO) / CZM_PI;
+    if (LatitudeMinFlag > 0 || LatitudeMaxFlag > 0)
+    {
+      latitudeUV = latitudeUV * LatitudeUVScale + LatitudeUVOffset;
+    }
+
+    // Height: scale to the range [0, 1]
+    float heightUV = 1.0 + height * InverseHeightDifferenceUV;
+
+    return float3(longitudeUV, latitudeUV, heightUV);
+  }
   
   /**
   * Converts the input position (vanilla UV Space) to its Shape UV Space relative to the
@@ -492,6 +1343,8 @@ struct ShapeUtility
     {
       case BOX:
         return ConvertUVToShapeUVSpaceBox(UVPosition, JacobianT);
+      case ELLIPSOID:
+        return ConvertUVToShapeUVSpaceEllipsoid(UVPosition, JacobianT);
       default:
         // Default return
         JacobianT = float3x3(1, 0, 0, 0, 1, 0, 0, 0, 1);
@@ -937,6 +1790,7 @@ switch (ShapeConstant) {
     DataTextures.PaddingBefore = round(PaddingBefore.xzy);
     DataTextures.PaddingAfter = round(PaddingAfter.xzy);
     break;
+  case ELLIPSOID:
   default:
     DataTextures.GridDimensions = round(GridDimensions);
     DataTextures.PaddingBefore = round(PaddingBefore);
@@ -954,7 +1808,13 @@ float3 PositionUV = R.Origin + CurrentT * R.Direction;
 float3x3 JacobianT;
 float3 PositionShapeUVSpace = VoxelOctree.ShapeUtils.ConvertUVToShapeUVSpace(PositionUV, JacobianT);
 
+// For ellipsoids, the UV direction has been scaled to a space where the ellipsoid is a sphere.
+// Undo this scaling to get the raw direction.
 float3 RawDirection = R.Direction;
+if (VoxelOctree.ShapeUtils.ShapeConstant == ELLIPSOID)
+{
+  RawDirection *= VoxelOctree.ShapeUtils.RadiiUV;
+}
 
 OctreeTraversal Traversal;
 TileSample Sample;
@@ -989,7 +1849,20 @@ for (int step = 0; step < STEP_COUNT_MAX; step++) {
   // Keep raymarching
   CurrentT += NextIntersection.t;
   if (CurrentT > EndT) {
-    break;
+    if (ShapeConstant == ELLIPSOID) {
+      // For CYLINDER: once clipping / exaggeration is supported, this must be done for cylinders too.
+      Intersections = VoxelOctree.ShapeUtils.ListState.GetNextIntersections(IntersectionList);
+      if (Intersections.Entry.t == NO_HIT) {
+        break;
+      } else {
+        // Found another intersection. Resume raymarching there
+        CurrentT = Intersections.Entry.t;
+        EndT = Intersections.Exit.t;
+      }
+    } else {
+        break;
+    }
+
   }
 
   PositionUV = R.Origin + CurrentT * R.Direction;
diff --git a/Source/CesiumRuntime/Private/CesiumVoxelRendererComponent.cpp b/Source/CesiumRuntime/Private/CesiumVoxelRendererComponent.cpp
index ca70b9bf8..6d73f1668 100644
--- a/Source/CesiumRuntime/Private/CesiumVoxelRendererComponent.cpp
+++ b/Source/CesiumRuntime/Private/CesiumVoxelRendererComponent.cpp
@@ -67,11 +67,12 @@ void UCesiumVoxelRendererComponent::BeginDestroy() {
 namespace {
 EVoxelGridShape getVoxelGridShape(
     const Cesium3DTilesSelection::BoundingVolume& boundingVolume) {
-  const CesiumGeometry::OrientedBoundingBox* pBox =
-      std::get_if<CesiumGeometry::OrientedBoundingBox>(&boundingVolume);
-  if (pBox) {
+  if (std::get_if<CesiumGeometry::OrientedBoundingBox>(&boundingVolume)) {
     return EVoxelGridShape::Box;
   }
+  if (std::get_if<CesiumGeospatial::BoundingRegion>(&boundingVolume)) {
+    return EVoxelGridShape::Ellipsoid;
+  }
 
   return EVoxelGridShape::Invalid;
 }
@@ -111,6 +112,290 @@ void setVoxelBoxProperties(
       FVector4(0, 0, 0.02, 0));
 }
 
+/**
+ * @brief Describes the quality of a radian value relative to the axis it is
+ * defined in. This determines the math for the ray-intersection tested against
+ * that value in the voxel shader.
+ */
+enum CartographicAngleDescription : int8 {
+  None = 0,
+  Zero = 1,
+  UnderHalf = 2,
+  Half = 3,
+  OverHalf = 4
+};
+
+CartographicAngleDescription interpretLongitudeRange(double value) {
+  const double longitudeEpsilon = CesiumUtility::Math::Epsilon10;
+
+  if (value >= CesiumUtility::Math::OnePi - longitudeEpsilon &&
+      value < CesiumUtility::Math::TwoPi - longitudeEpsilon) {
+    // longitude range > PI
+    return CartographicAngleDescription::OverHalf;
+  }
+  if (value > longitudeEpsilon &&
+      value < CesiumUtility::Math::OnePi - longitudeEpsilon) {
+    // longitude range < PI
+    return CartographicAngleDescription::UnderHalf;
+  }
+  if (value < longitudeEpsilon) {
+    // longitude range ~= 0
+    return CartographicAngleDescription::Zero;
+  }
+
+  return CartographicAngleDescription::None;
+}
+
+CartographicAngleDescription interpretLatitudeValue(double value) {
+  const double latitudeEpsilon = CesiumUtility::Math::Epsilon10;
+  const double zeroLatitudeEpsilon =
+      CesiumUtility::Math::Epsilon3; // 0.001 radians = 0.05729578 degrees
+
+  if (value > -CesiumUtility::Math::OnePi + latitudeEpsilon &&
+      value < -zeroLatitudeEpsilon) {
+    // latitude between (-PI, 0)
+    return CartographicAngleDescription::UnderHalf;
+  }
+  if (value >= -zeroLatitudeEpsilon && value <= +zeroLatitudeEpsilon) {
+    // latitude ~= 0
+    return CartographicAngleDescription::Half;
+  }
+  if (value > zeroLatitudeEpsilon) {
+    // latitude between (0, PI)
+    return CartographicAngleDescription::OverHalf;
+  }
+
+  return CartographicAngleDescription::None;
+}
+
+FVector getEllipsoidRadii(const ACesiumGeoreference* pGeoreference) {
+  FVector radii =
+      VecMath::createVector(CesiumGeospatial::Ellipsoid::WGS84.getRadii());
+  if (pGeoreference) {
+    UCesiumEllipsoid* pEllipsoid = pGeoreference->GetEllipsoid();
+    radii = pEllipsoid ? pEllipsoid->GetRadii() : radii;
+  }
+
+  return radii;
+}
+
+void setVoxelEllipsoidProperties(
+    UCesiumVoxelRendererComponent* pVoxelComponent,
+    UMaterialInstanceDynamic* pVoxelMaterial,
+    const CesiumGeospatial::BoundingRegion& region,
+    ACesium3DTileset* pTileset) {
+  // Although the ellipsoid corresponds to the size & location of the Earth, the
+  // cube is scaled to fit the region, which may be much smaller. This prevents
+  // unnecessary pixels from running the voxel raymarching shader.
+  const CesiumGeometry::OrientedBoundingBox box = region.getBoundingBox();
+  glm::dmat3 halfAxes = box.getHalfAxes();
+  pVoxelComponent->HighPrecisionTransform = glm::dmat4(
+      glm::dvec4(halfAxes[0], 0) * 0.02,
+      glm::dvec4(halfAxes[1], 0) * 0.02,
+      glm::dvec4(halfAxes[2], 0) * 0.02,
+      glm::dvec4(box.getCenter(), 1));
+
+  FVector radii = getEllipsoidRadii(pTileset->ResolveGeoreference());
+  // The default bounds define the minimum and maximum extents for the shape's
+  // actual bounds, in the order of (longitude, latitude, height). The longitude
+  // and latitude bounds describe the angular range of the full ellipsoid.
+  const FVector defaultMinimumBounds(
+      -CesiumUtility::Math::OnePi,
+      -CesiumUtility::Math::PiOverTwo,
+      -radii.GetMin());
+  const FVector defaultMaximumBounds(
+      CesiumUtility::Math::OnePi,
+      CesiumUtility::Math::PiOverTwo,
+      10 * radii.GetMin());
+
+  const CesiumGeospatial::GlobeRectangle& rectangle = region.getRectangle();
+  double minimumLongitude = rectangle.getWest();
+  double maximumLongitude = rectangle.getEast();
+  double minimumLatitude = rectangle.getSouth();
+  double maximumLatitude = rectangle.getNorth();
+
+  // Don't let the minimum height extend past the center of the Earth.
+  double minimumHeight =
+      glm::max(region.getMinimumHeight(), defaultMinimumBounds.Z);
+  double maximumHeight = region.getMaximumHeight();
+
+  // Defines the extents of the longitude in UV space. In other words, this
+  // expresses the minimum, maximum, and midpoint values of the longitude range
+  // in UV space.
+  FVector longitudeUVExtents = FVector::Zero();
+  double longitudeUVScale = 1;
+  double longitudeUVOffset = 0;
+
+  FIntVector4 longitudeFlags(0);
+
+  // Longitude
+  {
+    const double defaultRange = CesiumUtility::Math::TwoPi;
+    bool isLongitudeReversed = maximumLongitude < minimumLongitude;
+    double longitudeRange = maximumLongitude - minimumLongitude +
+                            isLongitudeReversed * defaultRange;
+
+    // Refers to the discontinuity at longitude 0 / 2pi.
+    const double discontinuityEpsilon =
+        CesiumUtility::Math::Epsilon3; // 0.001 radians = 0.05729578 degrees
+    bool longitudeMinimumHasDiscontinuity = CesiumUtility::Math::equalsEpsilon(
+        minimumLongitude,
+        0,
+        discontinuityEpsilon);
+    bool longitudeMaximumHasDiscontinuity = CesiumUtility::Math::equalsEpsilon(
+        maximumLongitude,
+        CesiumUtility::Math::TwoPi,
+        discontinuityEpsilon);
+
+    CartographicAngleDescription longitudeRangeIndicator =
+        interpretLongitudeRange(longitudeRange);
+
+    longitudeFlags.X = longitudeRangeIndicator;
+    longitudeFlags.Y = longitudeMinimumHasDiscontinuity;
+    longitudeFlags.Z = longitudeMaximumHasDiscontinuity;
+    longitudeFlags.W = isLongitudeReversed;
+
+    // Compute the extents of the longitude range in UV Shape Space.
+    double minimumLongitudeUV =
+        (minimumLongitude - defaultMinimumBounds.X) / defaultRange;
+    double maximumLongitudeUV =
+        (maximumLongitude - defaultMinimumBounds.X) / defaultRange;
+    // Given a longitude range, represents the actual value where "0" would be
+    // in UV coordinates.
+    double longitudeRangeUVZero = 1.0 - longitudeRange / defaultRange;
+    // TODO: document this
+    double longitudeRangeUVZeroMid =
+        glm::fract(maximumLongitudeUV + 0.5 * longitudeRangeUVZero);
+
+    longitudeUVExtents = FVector(
+        minimumLongitudeUV,
+        maximumLongitudeUV,
+        longitudeRangeUVZeroMid);
+
+    const double longitudeEpsilon = CesiumUtility::Math::Epsilon10;
+    if (longitudeRange > longitudeEpsilon) {
+      longitudeUVScale = defaultRange / longitudeRange;
+      longitudeUVOffset =
+          -(minimumLongitude - defaultMinimumBounds.X) / longitudeRange;
+    }
+  }
+
+  // Latitude
+  CartographicAngleDescription latitudeMinValueFlag =
+      interpretLatitudeValue(minimumLatitude);
+  CartographicAngleDescription latitudeMaxValueFlag =
+      interpretLatitudeValue(maximumLatitude);
+  double latitudeUVScale = 1;
+  double latitudeUVOffset = 0;
+
+  {
+    const double latitudeEpsilon = CesiumUtility::Math::Epsilon10;
+    const double defaultRange = CesiumUtility::Math::OnePi;
+    double latitudeRange = maximumLatitude - minimumLatitude;
+    if (latitudeRange >= latitudeEpsilon) {
+      latitudeUVScale = defaultRange / latitudeRange;
+      latitudeUVOffset =
+          (defaultMinimumBounds.Y - minimumLatitude) / latitudeRange;
+    }
+  }
+
+  // Compute the farthest a point can be from the center of the ellipsoid.
+  const FVector outerExtent = radii + maximumHeight;
+  const double maximumExtent = outerExtent.GetMax();
+
+  const FVector radiiUV = outerExtent / maximumExtent;
+  const double axisRatio = radiiUV.Z / radiiUV.X;
+  const double eccentricitySquared = 1.0 - axisRatio * axisRatio;
+  const FVector2D evoluteScale(
+      (radiiUV.X * radiiUV.X - radiiUV.Z * radiiUV.Z) / radiiUV.X,
+      (radiiUV.Z * radiiUV.Z - radiiUV.X * radiiUV.X) / radiiUV.Z);
+
+  // Used to compute geodetic surface normal.
+  const FVector inverseRadiiSquaredUV = FVector::One() / (radiiUV * radiiUV);
+  // The percent of space that is between the inner and outer ellipsoid.
+  const double thickness = (maximumHeight - minimumHeight) / maximumExtent;
+  const double inverseHeightDifferenceUV =
+      maximumHeight != minimumHeight ? 1.0 / thickness : 0;
+
+  // Ray-intersection math for latitude requires sin(latitude).
+  // The actual latitude values aren't used by other parts of the shader, so
+  // passing sin(latitude) here saves space.
+  // Shape Min Bounds = Region Min (xyz)
+  // X = longitude, Y = sin(latitude), Z = height relative to the maximum extent
+  pVoxelMaterial->SetVectorParameterValueByInfo(
+      FMaterialParameterInfo(
+          UTF8_TO_TCHAR("Shape Min Bounds"),
+          EMaterialParameterAssociation::LayerParameter,
+          0),
+      FVector(
+          minimumLongitude,
+          FMath::Sin(minimumLatitude),
+          (minimumHeight - maximumHeight) / maximumExtent));
+
+  // Shape Max Bounds = Region Max (xyz)
+  // X = longitude, Y = sin(latitude), Z = height relative to the maximum extent
+  // Since clipping isn't supported, Z resolves to 0.
+  pVoxelMaterial->SetVectorParameterValueByInfo(
+      FMaterialParameterInfo(
+          UTF8_TO_TCHAR("Shape Max Bounds"),
+          EMaterialParameterAssociation::LayerParameter,
+          0),
+      FVector(maximumLongitude, FMath::Sin(maximumLatitude), 0));
+
+  // Data is packed across multiple vec4s to conserve space.
+  // 0 = Longitude Range Flags (xyzw)
+  pVoxelMaterial->SetVectorParameterValueByInfo(
+      FMaterialParameterInfo(
+          UTF8_TO_TCHAR("Shape Packed Data 0"),
+          EMaterialParameterAssociation::LayerParameter,
+          0),
+      FVector4(longitudeFlags));
+  // 1 = Min Latitude Flag (x), Max Latitude Flag (y), Evolute Scale (zw)
+  pVoxelMaterial->SetVectorParameterValueByInfo(
+      FMaterialParameterInfo(
+          UTF8_TO_TCHAR("Shape Packed Data 1"),
+          EMaterialParameterAssociation::LayerParameter,
+          0),
+      FVector4(
+          latitudeMinValueFlag,
+          latitudeMaxValueFlag,
+          evoluteScale.X,
+          evoluteScale.Y));
+  // 2 = Radii UV (xyz)
+  pVoxelMaterial->SetVectorParameterValueByInfo(
+      FMaterialParameterInfo(
+          UTF8_TO_TCHAR("Shape Packed Data 2"),
+          EMaterialParameterAssociation::LayerParameter,
+          0),
+      FVector4(radiiUV, 0));
+  // 3 = Inverse Radii UV Squared (xyz), Inverse Height Difference UV (w)
+  pVoxelMaterial->SetVectorParameterValueByInfo(
+      FMaterialParameterInfo(
+          UTF8_TO_TCHAR("Shape Packed Data 3"),
+          EMaterialParameterAssociation::LayerParameter,
+          0),
+      FVector4(inverseRadiiSquaredUV, inverseHeightDifferenceUV));
+  // 4 = Longitude UV extents (xyz), Eccentricity Squared (w)
+  pVoxelMaterial->SetVectorParameterValueByInfo(
+      FMaterialParameterInfo(
+          UTF8_TO_TCHAR("Shape Packed Data 4"),
+          EMaterialParameterAssociation::LayerParameter,
+          0),
+      FVector4(longitudeUVExtents, eccentricitySquared));
+  // 5 = UV -> Shape UV Transforms (scale / offset)
+  // Longitude (xy), Latitude (zw)
+  pVoxelMaterial->SetVectorParameterValueByInfo(
+      FMaterialParameterInfo(
+          UTF8_TO_TCHAR("Shape Packed Data 5"),
+          EMaterialParameterAssociation::LayerParameter,
+          0),
+      FVector4(
+          longitudeUVScale,
+          longitudeUVOffset,
+          latitudeUVScale,
+          latitudeUVOffset));
+}
+
 FCesiumMetadataValue
 getMetadataValue(const std::optional<CesiumUtility::JsonValue>& jsonValue) {
   if (!jsonValue)
@@ -244,6 +529,15 @@ UCesiumVoxelRendererComponent::CreateVoxelMaterial(
         std::get_if<CesiumGeometry::OrientedBoundingBox>(&boundingVolume);
     assert(pBox != nullptr);
     setVoxelBoxProperties(pVoxelComponent, pVoxelMaterial, *pBox);
+  } else if (shape == EVoxelGridShape::Ellipsoid) {
+    const CesiumGeospatial::BoundingRegion* pRegion =
+        std::get_if<CesiumGeospatial::BoundingRegion>(&boundingVolume);
+    assert(pRegion != nullptr);
+    setVoxelEllipsoidProperties(
+        pVoxelComponent,
+        pVoxelMaterial,
+        *pRegion,
+        pTilesetActor);
   }
 
   if (pDescription && pVoxelClass) {
@@ -516,6 +810,60 @@ void UCesiumVoxelRendererComponent::UpdateTiles(
   this->_pResources->Update(VisibleTiles, VisibleTileScreenSpaceErrors);
 }
 
+namespace {
+/**
+ * Updates the input voxel material to account for origin shifting or ellipsoid
+ * changes from the tileset's georeference.
+ */
+void updateEllipsoidVoxelParameters(
+    UMaterialInstanceDynamic* pMaterial,
+    const ACesiumGeoreference* pGeoreference) {
+  if (!pMaterial || !pGeoreference) {
+    return;
+  }
+  FVector radii = getEllipsoidRadii(pGeoreference);
+  FMatrix unrealToEcef =
+      pGeoreference->ComputeUnrealToEarthCenteredEarthFixedTransformation();
+  glm::dmat4 transformToUnit = glm::dmat4(
+                                   glm::dvec4(1.0 / (radii.X), 0, 0, 0),
+                                   glm::dvec4(0, 1.0 / (radii.Y), 0, 0),
+                                   glm::dvec4(0, 0, 1.0 / (radii.Z), 0),
+                                   glm::dvec4(0, 0, 0, 1)) *
+                               VecMath::createMatrix4D(unrealToEcef);
+
+  pMaterial->SetVectorParameterValueByInfo(
+      FMaterialParameterInfo(
+          UTF8_TO_TCHAR("Shape TransformToUnit Row 0"),
+          EMaterialParameterAssociation::LayerParameter,
+          0),
+      FVector4(
+          transformToUnit[0][0],
+          transformToUnit[1][0],
+          transformToUnit[2][0],
+          transformToUnit[3][0]));
+  pMaterial->SetVectorParameterValueByInfo(
+      FMaterialParameterInfo(
+          UTF8_TO_TCHAR("Shape TransformToUnit Row 1"),
+          EMaterialParameterAssociation::LayerParameter,
+          0),
+      FVector4(
+          transformToUnit[0][1],
+          transformToUnit[1][1],
+          transformToUnit[2][1],
+          transformToUnit[3][1]));
+  pMaterial->SetVectorParameterValueByInfo(
+      FMaterialParameterInfo(
+          UTF8_TO_TCHAR("Shape TransformToUnit Row 2"),
+          EMaterialParameterAssociation::LayerParameter,
+          0),
+      FVector4(
+          transformToUnit[0][2],
+          transformToUnit[1][2],
+          transformToUnit[2][2],
+          transformToUnit[3][2]));
+}
+} // namespace
+
 void UCesiumVoxelRendererComponent::UpdateTransformFromCesium(
     const glm::dmat4& CesiumToUnrealTransform) {
   FTransform transform = FTransform(VecMath::createMatrix(
@@ -541,4 +889,14 @@ void UCesiumVoxelRendererComponent::UpdateTransformFromCesium(
     this->MeshComponent->UpdateComponentToWorld();
     this->MeshComponent->MarkRenderTransformDirty();
   }
+
+  if (this->Options.gridShape == EVoxelGridShape::Ellipsoid) {
+    // Ellipsoid voxels are rendered specially due to the ellipsoid radii and
+    // georeference, so the material must be updated here.
+    UMaterialInstanceDynamic* pMaterial =
+        Cast<UMaterialInstanceDynamic>(this->MeshComponent->GetMaterial(0));
+    ACesiumGeoreference* pGeoreference = this->_pTileset->ResolveGeoreference();
+
+    updateEllipsoidVoxelParameters(pMaterial, pGeoreference);
+  }
 }

From b984e008d3768507a86a0248bfe0d058abde4a16 Mon Sep 17 00:00:00 2001
From: Janine Liu <janine@cesium.com>
Date: Wed, 12 Mar 2025 15:08:00 -0400
Subject: [PATCH 2/8] Bugfixes and shader comments

---
 Content/Materials/CesiumVoxelTemplate.hlsl        | 14 +++++++-------
 .../CesiumRuntime/Private/CesiumGltfComponent.cpp |  1 +
 Source/CesiumRuntime/Private/VoxelOctree.cpp      | 15 +++++++++------
 3 files changed, 17 insertions(+), 13 deletions(-)

diff --git a/Content/Materials/CesiumVoxelTemplate.hlsl b/Content/Materials/CesiumVoxelTemplate.hlsl
index a880c2f75..81a64d494 100644
--- a/Content/Materials/CesiumVoxelTemplate.hlsl
+++ b/Content/Materials/CesiumVoxelTemplate.hlsl
@@ -756,8 +756,8 @@ struct ShapeUtility
   }
 
   /**
-  * Tests if the input ray intersects the upper half of the quadric cone, i.e., the cone where its apex
-  * points downward.
+  * Tests if the input ray intersects the negative space outside of the cone. In other words,
+  * this captures the volume outside the cone, excluding the part of the ray that is inside the cone.
   */ 
   void IntersectFlippedCone(in Ray R, in float CosHalfAngle, out RayIntersections FirstIntersections, out RayIntersections SecondIntersections)
   {
@@ -817,10 +817,10 @@ struct ShapeUtility
   }
   
   /**
-  * Tests if the input ray intersects the upper half of the quadric cone, i.e., the cone where its apex
-  * points downward.
+  * Tests if the input ray intersects the volume inside the quadric cone. The cone's orientation is
+  * dependent on the sign of the input angle.
   */ 
-  RayIntersections IntersectUprightCone(in Ray R, in float CosHalfAngle)
+  RayIntersections IntersectCone(in Ray R, in float CosHalfAngle)
   {
     // Undo the scaling from ellipsoid to sphere
     float3 position = R.Origin * RadiiUV;
@@ -947,7 +947,7 @@ struct ShapeUtility
     if (LatitudeMinFlag == ANGLE_UNDER_HALF)
     {
       ListState.Length += 2;
-      RayIntersections BottomConeResult = IntersectUprightCone(R, MinBounds.y);
+      RayIntersections BottomConeResult = IntersectCone(R, MinBounds.y);
       setShapeIntersections(Intersections, ListState, 2, BottomConeResult);
     }
     else if (LatitudeMinFlag == ANGLE_HALF)
@@ -1022,7 +1022,7 @@ struct ShapeUtility
     }
     else if (LatitudeMaxFlag == ANGLE_OVER_HALF)
     {
-      RayIntersections TopConeResult = IntersectUprightCone(R, MaxBounds.y);
+      RayIntersections TopConeResult = IntersectCone(R, MaxBounds.y);
       [branch]
       switch (ListState.Length)
       {
diff --git a/Source/CesiumRuntime/Private/CesiumGltfComponent.cpp b/Source/CesiumRuntime/Private/CesiumGltfComponent.cpp
index 9d7514ed6..b96a821db 100644
--- a/Source/CesiumRuntime/Private/CesiumGltfComponent.cpp
+++ b/Source/CesiumRuntime/Private/CesiumGltfComponent.cpp
@@ -2026,6 +2026,7 @@ static void loadVoxels(
         ValidatedVoxelBuffer{pBuffer, pBufferView});
   }
 
+  primitiveResult.meshIndex = options.pMeshOptions->meshIndex;
   primitiveResult.primitiveIndex = options.primitiveIndex;
   primitiveResult.transform = transform * yInvertMatrix;
   primitiveResult.Metadata =
diff --git a/Source/CesiumRuntime/Private/VoxelOctree.cpp b/Source/CesiumRuntime/Private/VoxelOctree.cpp
index b10a052ae..7d0c98188 100644
--- a/Source/CesiumRuntime/Private/VoxelOctree.cpp
+++ b/Source/CesiumRuntime/Private/VoxelOctree.cpp
@@ -113,12 +113,15 @@ void UCesiumVoxelOctreeTexture::Update(const std::vector<std::byte>& data) {
   ENQUEUE_RENDER_COMMAND(Cesium_UpdateResource)
   ([pResource = this->_pResource, &data, region, sourcePitch](
        FRHICommandListImmediate& RHICmdList) {
-    RHIUpdateTexture2D(
-        pResource->TextureRHI,
-        0,
-        region,
-        sourcePitch,
-        reinterpret_cast<const uint8*>(data.data()));
+    // The data could be emptied if the tileset refreshes mid-update.
+    if (!data.empty()) {
+      RHIUpdateTexture2D(
+          pResource->TextureRHI,
+          0,
+          region,
+          sourcePitch,
+          reinterpret_cast<const uint8*>(data.data()));
+    }
   });
 }
 

From 55e1f1f8723cc116a257788203f4dd44c8265356 Mon Sep 17 00:00:00 2001
From: Janine Liu <janine@cesium.com>
Date: Wed, 9 Jul 2025 15:25:39 -0400
Subject: [PATCH 3/8] Update for ellipsoid voxels

---
 .../Materials/Layers/ML_CesiumVoxel.uasset    | Bin 37165 -> 47932 bytes
 Shaders/Private/CesiumBox.usf                 |  15 +-
 Shaders/Private/CesiumEllipsoidRegion.usf     | 835 ++++++++++++++++++
 Shaders/Private/CesiumRayIntersectionList.usf | 182 ++++
 Shaders/Private/CesiumShape.usf               |  94 +-
 Shaders/Private/CesiumVectorUtility.usf       |  18 +
 Shaders/Private/CesiumVoxelTemplate.usf       |  34 +-
 .../CesiumRuntime/Private/Cesium3DTileset.cpp |   2 +-
 .../Private/CesiumVoxelRendererComponent.cpp  |  17 +-
 .../Private/CesiumVoxelRendererComponent.h    |  10 +-
 10 files changed, 1157 insertions(+), 50 deletions(-)
 create mode 100644 Shaders/Private/CesiumEllipsoidRegion.usf
 create mode 100644 Shaders/Private/CesiumRayIntersectionList.usf

diff --git a/Content/Materials/Layers/ML_CesiumVoxel.uasset b/Content/Materials/Layers/ML_CesiumVoxel.uasset
index 811ed7395b8086693c3f5df093ab4183cd2e228c..28e011b07474b0ae709dd285423b8cd819d0ee4f 100644
GIT binary patch
delta 14400
zcmbtb3p|wB`+wh&OES4Mx)8mg7$k*A)F|VY`_P4In;4A73}y(eW=gt<(94D*m3t||
zb}zN7+h!BeZl%<=m2O*Gy0`t$nVDBJLw@`F{ol{$9p^pI_dL&eF6W%*ocC3`hjii%
zDJO5(Q5u5aTQux+5O9@(`vwY22r^D5_l`xtI8Z+wJk{$W$QW@gg5+92I)W5iN$#t&
z5rjMYaHr;c%>4-O{Tx=t{1;B1xz)B899;{t>c_c1O=*WrA)it_6lEBVoP!{0$03ie
z*B5ZX&ZL7XjpGe%XuGOV1C;%lVJlhTu1o>b8hKN^;I4e#G{i+@vq8RYI#N*4O_i^+
zL8cx*;wxWgi<mz_6Xomdkm5Jt+9Sns!}03jr-a&|%!p#fvM?WJ5;m8YzzydkK65`O
z%MjRHg79ISay_NLg{1Eo>^x?T<n6y_JUDF{VqCw`bx>L;HXZq}o&97`B4&g5WN2g#
zPP0WkZZtIxPP0S4oNw#E52_pD>FVKDnnaI*ic$?-QKh3_tJz#>wZ+bhfepRtN0iB%
zKpS{H27cRtPh{((n{cmDrv3H3pgsfE%hhZdP?I^JCTl>=)&Vuy18Q;x)Z`AR**2gC
zA5fDwpeBDn%?|@=whyQ&7*JC<pr#1ajKXtAQP7QIT8YG@;@=Ih1}(v-vGs9lwgy?2
zhQfAJ{s9&y8cYL)&ESUPOby(A(oqaYl=#B0I*!7BbJW07Db^DSGzpMo(Q0#tLMtht
z7tR|BJxl>z;W`vrSpnTg9}2CafOd5s3azSu#yy5YLwm^D<eukH=-~?JY4a7)OLZmo
z&Hx0$JLOW4$XgNJDuX7xL!c-5D59N5OL%XQ;T^3mL0uYvu9ZO(ykn<G(6);dd8bmw
z+i|4Sj|!2NDI;&OG}6k_6***$k@PN027MD63|yZ#Nd6-Iars$_-U<kWor0vif6Gxs
z$BdQm&Xz%ckb*3?DWbcjXvl`(O`s>^L!k*f5oj!LC^Rub2((`QP-tR=5NOH|L!p%c
z5ymBf9=UxeG%-R5v|7PXXkvsAXj0)&XkvsA=uZRCsrnM{<N`K?y;n$m`KCyb_cJLP
zx{W|LNYRUSDx#N;m+*#eBhb^O=n=#u*&kJ{QZ#fMfmW5Gop&knb~cdkhHfL!94Xqh
zR1sY(MMJj{Xd5ZIZnq-Zz)->)x{W}4NYQ`nQAB6SpothH&=zGwp+^8B^w0q*Z$2@7
z^xFiE>>4Tcg~G^&y@MGM<ZwxfnOF%P;3dc5+z=ed01P?2lwzJ#!RFwS<M3<<4lo0Q
z91cq{w`yQ>aLI8XqGT{FVDbVvyq9vo_rvDklH<@fL?4LJOxQj1fP}$U4zgY%Ml%sz
zod*@s-4nh>05O^g^v$}V(8Opa(D&+xLKCBzK=&RR3Jq-pZvy?IVJI}Q-Vx}xjYFZ~
zQUrP5m0J0qCPg%i#cnC}FNIiyQx+sa4nz>&IjqQ`c%ro)%x8+F1lU6hDaZ9k@ZBzY
zc#r)s{Fb8{NPWgznMz5|WQ4(%#5@Cs&loA@;L-l(8Z^+{$aWktWxz2x7q$b}xntl5
zM`>g5oS0#F4n+-rW@IMBKqc%KL{5|`o63}}WXeV|Wn-DLiA;HtOnI`j4EqY@5PzZm
z!0tk10w`#M+fYh`9f2XO18&$r3*2z1!{@Q!9tUn1#4t2rz!6~t8IA%s+)034ATk=<
zW56vd6!LuvED8rctS15ah<^VAf{d3b>&cYOWXe>TvbjvzK&EUVQ?`^T!$d{4L&Jc4
z1WqEbM{p{DMA*xz;D$by1w+edOY2}6hOI1}FztWySNQ*B0NcaPLN`H0pfzE~2nC=;
zVW**+p^Kp_31u<pf7siJ;D#E(po7{&(?fqlvqPUklS98kb3@-kQ$znkGeaLk6T{em
z=7sSA0~W@J5x8NT7=xRL6&R>6K;dwm45banUrwYnOrErkBn&~k%=U;3i>d5U<x>z!
zR(r7zHxIcKVpgiWtntH&?4sIvk#P=}tv&v@7F2|+%V=M!8o<nr{(@wAoYKFhUtawM
zId%P0l+x&ouGxCmR6iZtb>{jfcI8huiq@^U@CDgGMpzx!pLcdFS{j<!%Cii+mXopB
z{TG@t<w#knihYPMX6A=>e-dMR?Tno#qtMrRjd@5`d!Zli7&{LMpi|IEl-BcIS?9Y%
z6TEpRqHd+Kn@Q87&^*SQY_Ay^$f$D1f;Eg5{p?y~cJ8s@9aCmMc=u0ltmyA&&AOe7
zqPn_7)4rSG^2ysVRoP<ec^g*Q(mSsxNq&EwCLOMIGS#_zKv)qur`u)iKSJ|vmmB-B
z^4NN>n;wULJfFXEEcVdSUi)3qsVevASS)bE>(Uif(YdMTs5K*kYD)blu1$JSLduBV
zDd-J6FxH>+ecAGB7FlC_k6+oba<t#1yqGx!`~F$9_ZQ>aqM0kND4+d(2HqI?@!%`X
zD^;IY7I)fOblXp{4c#mJf|P8&^abfVt@paAr*jAT%(o}$=2mu7@3CF!PZP9$9d)tb
zwP5VUw+2r(U9;6cNj=?ed^x<wqtaY>cFU%V<-U!d_qG?gU2U;8pxeat?TzazRfe|f
ziH}VQd`3YoQ9rtu>(i0x8D^!VlGBNWe%~FZjhIh2uU<jQ#G6S%kuehw4N2U5!tyY!
z=XLg{=yN+N_A9mKDN!f9ZRNO~2>Q`><4<~f(*40iq3U&LU7_BC{YLk4T3q|86HD&o
z8)f82q3g=%SK8A=GoN@J+}ioG_g+emxm^z9r-QYKn<L}Af!U&8sb_7T(R7TtPr6E!
zr-hwQF&VbXp;S4-KvihkmA}OBL{QKkk!i(69k0x7X5q|rW#QW13yjOX_Uk{PjgL6m
z^J3Fm$K?*=3mkX78~#uJxeIDNnOZA<uU~R~vy;QIZMJWkr`&k+$i#5to!iCDuWG1&
zyIq}Sk`dh5Q#}^zpl%xTS5RKd$IcH3l^Rv5@{|%oY4j(>PzrB0S2Be9)uqZBK0J08
zNz0I?vGEIHP*b$Dbk6CV$d73sy6(>p{H@)T{@$;B|G5)EYN_UzBeSFDnFa*bx;2@%
zDRq_os>0sl5jt69u=AnqU2n==?|Cf_2^AT!PiUjJmpd>PmOJ>9G7_g-vy~nEhiM%5
zHPf}Rx|^|am|;pAI^WaFOh}@yX?@=j`h#u#f#-f+t0yNLL@5=%oUm~eV;-MF>aKd&
z?(!)0^OGx=lT601v-mA;$*L`lbXvlQa3u0=^Xk{UmldKv;tuC}9_hIGyl87d`mLy2
zY^oW_i)ycKpMivIyP?^*WZ0vVo^PrSt6x8~ex9Mn`&;4d+@I)=@-z@@RUr~w_OQ~c
zrlOYazhRMAF`1U{*g-mL=UuUHiud}58A$$)cS6Si8f|qBd#-(b#YDOTgJxW1+Y^PJ
zEjdZIW^Au29kH^F(K2n#qsj~<ww$sq>yl$GZ6xL5!mMV&6PJAVleU*VdRB<^RQz<!
z>D|>=a~7`kD;?HCY0PU5J0MIph*qWNjx+JU^wTiyGrG}x6P{DiKa1Y{AlUa-q!VT=
z;H*yjCvLG_@2Sb{D%Ir;Mi1^GAx^9MbRVRX1GR2YZ~9FBiALU&VMbaRb>jDHe!dky
z2EA(W@95xg9v*G3ReM*(Df80}@v2E<8y45x5}-_TR@LZA22H1SQRa$Lp|8)PWqGK1
z=azt4Cx4>=U&GgD{3y+_!LKLeI|uOF7<A{9wqYc<^C7W`eRGQ4W-&H4A#>uKTalM7
z_D#UK*(nLOf^rIKIi;!ow`)aZ=AAJ|(Z6V#dork^ybV7TJOi<YSksfaCzsnk3rV7A
z8?Kt1syXK~@9atEP;#H}K<N3=q0g2d54aptUpLEeSJmezT5M-vJ44&Bqt~sY*RS&2
z@v`7Phf=b>;CJ<@r5E=N@7inF`1;A>yil6jnO+8Fas5@z%c8^MuC$9L&1>;)(y4qa
zG|$y^dh^VC#I($yC;pYB@(?87$E!)bD96xD+pIDpQ!mN3IY%$0uS%HiaWaZ@tVG!%
zX!%P&3i7hAT$OUQWn8IhMu1Owsp^{tWozqxPONgNR7!hpM*5S|)_fp(MuC9RoV(Mk
zY6N$t?fQqcB$0Y`ZRqC)`ms$U&&(Xv_?Z(f)waz_Yin1!FhbquY{lvIp?7Z0J!zbq
zw?<`Ouv)4Q$0l}H$cBohFG3_c{S({THhrdJXEgh9a+O(=!Q#wo4y_@XFKMIOI?L;j
zDy@!-?@#Sykj^^$W(SWhXV6A>X5{Eqjo;d5bb(HqlIEwf&$6*Ek3qU%XH*hhmlw2U
z`U8=%0rmF>xhpscl$Rq?G7uNvO10sSfBL&})0pFWNTwmxrPa3{d$vnhT}?sT(ATAE
zqbGBJPK^#03AL}Zm$gz*jb9Fr?d|o}bxJc!JzDemSR3i8rx#n)sMI>iyss5W&fb(a
z`FJ97?z=g~tgiDJDtEM&Y9@R^+UQ{&TKRUBV@AB|w3|}@b8p}7wG|_&kE6|l9rj<E
zaxS`EuUtR4>W0Y5E$(`0xiIr3-r2eMz2<33pKWtpMOjF?k4MGYQA#aNf~vKZX*8X!
zg_UUmO1*w61v+WIj92QQ?~m3nX=AiL?I}|!&4j<ANh&+ZZe#~N1h+j}b*3^xG@)oo
zhm9Ls^t>r_xABX8<El*Vwj!8qU_`!)DH4df`BDB!BZ)xp3g@!dJ~W&>tLfw3=c@lW
zT1;+Cc&J&sH23qTHNnN+uC|YU{hjhbi=7Bo+iWkS!#3}?n*y)Y>T89r^*@?-`UC~_
zT-$V$f+pFIcplBN3z@Dz?x}ggu+^WbYaNU7eP|VBj26;l)!St%yH&*fEYZ}}vIUyv
z>7|Uf22mxfyj?nC)0%)lCI8*;=}vk`o?fioSTMbiXlxgtiV{sHA`__c7+DiN)6e*J
z&Lw4|Q-56H<~0Wgb|+C@ima(NZc(Ubh0Suk;rpLOKY4o0r?WI?u9jM#ljcJYP3MyZ
zOO{<R>ov)|Q&{u5{d^Z5*H@$U7^FN?LTC6NtQ6YWF*biXiaxu2KyQOB(s|*U@=Wr?
zB5gyV`k}0#iEEmoXH+O(<Mcf!3s$Y5FIZ6B);6qWW4%g=lIGnx>C>RE)JHmYtSu+|
zjdiAdsuqb}22OlmjePe!BUd-xp>(+M9J6Y)X63~5AxSwO+LY_qXExRbX{54Or4%8D
zAE5TP@m~#Ueq{4o5MP#awf(5hPC<moq}9=;m0lEY7)3hAw()7TZL%I8a+Ejr_Ol;1
zT=Sbg<IO{S{xl=un7F-h6<28^Z$&@VNzTpeGl)+28Sy^cu)Eh*`G#*^>4S`Qo9}Ub
z7H8FD?9{8*Z&Q99U94Aad~H3Ky6MN-z<t5N(fJ`F%j(KBe~X}uC{ir@QL-*#XEdRx
zfV^xSaFzOmA4<rGXOYAuXivx5s#ji`6clp`d|q`%Sl$X~K~6>4wmlyA+lJk&nIZJb
zJhdk8qRFDGg}*5WZ1ZE%#@%>RX2~hkt+#Ei2|HYUU_#WsqKeW4{Wsfl(*C5aR<90P
zZy)t`Z&YOfO*>$I-POi6665;wj>IDW#1~Pd%Q|6y&Pn$Rvd>OwpqzO(XZ8Ee#_3^x
z?k(;#TaobCw&1)q7IudF&piaIri)(1f*S1$`~{=EkR!E$7Qw+eDV;U1n-7t_vM$ur
zTr+aGWrha5p!Rf9i%*h0Y<g)rO1`RgV^wydeWj`oc4|`E_I7q`P>n_cu6V_M7s7qK
zv=e!`wmtV|M#SSVdt2n*waf3a7mS>JUg^H@m%uss!X;gHYa=#{RkO+cf_%UF^U7Qd
z-`O!&drjZ$!mJZTOYaqoEGxTh2Bb#q$v~)CogKYdooyt4hnbWB+UUF6lYUD5>tSH9
zT`$Ar<c?R>QQ<xJ%C~xL?4>!BDc|3=-y);c)BLO_O;4@QyeVL-SK+dPF-mcPZg-r6
z7yeN@U;mTvtjO}*9kXhs;6-=-)WuSeLwb333pIT;X`_{Zf0VAC5}9_pJx)uV_8=(m
zax;3-v<X=>b?suk19d-^wxHyh2EyZCkn8(DyDy8pykMicn?~^$MEz&RHoos!!^^+6
z6m0wXK7M-}Sjf5_vGd-1tSAI&djL)ImUGlCYM#H6@U3PDd)Uscd3k*w8<FnH_N>bG
zU3A3OU~kdVpSp^=*=D5R?Ywym%S+uaqRVGGw%3iZ@mtsuw9umC@5;@obPJ6W+dg(a
zF+CZvFRiU}T_6>$IuyrnH@$2@wK7k7u5#C#Jof&3cRLKf-C(>kM%Q|$zlvyE@q-0P
zN0d%Z`|F<2`*}{P4NY%Nj3HMw<m~nkB_k7)4<#0}?|oQZV({X0_jvV0?Y1S4U(hO7
zCTBmY{LrP_-PdyW^B+uix`Fwj*Od%(hmX?M#t$n6rAo1TWBV$Aa;C<l9gQFSPqJSl
ztw5s+ppZMA;*Wy7*}L*!cd0T+HPEcCZTZLAnf|D+Z{v;lgi=+oQ$e!c`q=6?>QGzX
z`b-xYPKtGG-@H>0y)^@=n-urSE!)xWuOb!&DKF?-A=Jw{R+Q>?C#enf1Px#J`c_nW
zB_NZG1?6RM^95!n-IR{rEnT>W-CUL3U~5-PYsj%XgrW)f-I-=6hKrJo@oamw{^v9H
zEKrV%F~j{FsNzGH`CzZtf5D9=yd=#hDm0&npRr$zvuA2ZS(MEDk1R|TSeS}g9C0LA
zJahbyEUXk*SOFGJC=IV)J)Rs3PBz>nT?TyyYH1lLXL9UxBHcx%qD=EKy}NOgGwwg>
zX`~>*NJc^(*s#k>_&sK7V(k6ru9fqX8OQ~VyZ=Ulv4R9+841ncEJI$xo%UP#cPEg1
zsz*+$3s%!F{5KLz6eO6)NPve5@)F|IqK^7cx=P<(p5n6oZ1~cDBVm$)gh~Ap3^c&$
zhP(j7MG1e{yz`o07wWNjhF9gse<NVBf`G~W0>+b<g99LW3DckcVZgBK^x5^mHEXrz
zF;yr5ora@Z-SDNei16Ariz=xj$AZ%gIR?1B6QnVn#-f6T6sK+!72u0zdx6kCz@lrf
zmj>00{@{Y+Ngcejb;Tp*XaN*>#uL|Sffu`W+G&HK@e;Nr?jJQy+i6g0WU8d?_XDUe
zX))TV18SkX?JjifHG@)g=hC%f2c@JZnc)5rYWV576mo(T0ne2v;^WS@bB)loxZb=a
zK>JnmaOuFdaix;s;T1eTf|meiyZ#JlR&+WVhUdhU;&T$lfp$4@RX~4beR#kHFM;Hb
z`1c|_SZXnDc;hw$6;8#1qS$;aoE6Ig-(FxmE;f(FXD7rg;3cs*eqn-m7RxGuA7P2D
zNDyF5jKd3Kaxh~yH;j`I&cec&TrN+5tzcpCEIyBu$O^}zSn(_?vij7in5O`XN#KjA
z95xqxXM%yxJ6Q3Uk*Bn^8=GU!!u*+nD83N}KJ#E=^8lU~gK>F)H>eOqv9K7XfECYX
za;$JGcLUtjLqH;{HyrW!RZR#)cz83sv2Y%X51I%t4tpgFW5&lblQF&^J|PUqXT_{w
z#q%vOK2OXn9DXbk#pGiG9_%y_77x0`6o`oeUhEVOE0M*KbUz3P=d)8lNB9_%0~kax
z<HN&1&%!~|FlNF`J_`#^<}zd0z$6?X5?}!@C>}J3XNkK6D7;myaKJHw!($3;u^2X=
zkE6>+34Q&(^Y9F!&g!Rla)F|JR+xaz<9f6Cf>2N2AbKE!J}=1A&zJF?H{I7g$YUwy
zOjch_R>#DQEdR~#;9z~&kwMy%)l+BT>jKUQEiv(jIX=YKDa14J5sM{oT>r;IvJcZ_
zAEwJb*vLNE%0AdZJCM}_f#(7-XUr^+*$@=CNL9xYQ?yFM{>5ranqf>}e7d){r$58b
z)77z|ZixvA#)BI#9%FGh>{vdJ9S+ok0TjFd3unhm>?-#C3Kolt@qyETyFeT;F@7vF
zjD@X&kt=S|&q3^ZK>ixj!Q5~dNgyO312#7ja)@NX0MY`3fVt|cEGxWy+2m#H2+T}u
z$ay?2U%=!FFlXnlJ*P^Y8OU7a+TTZWXG}Wq+~e8d;>K2itVj?R(lRrg%?@4wESB2v
zW?9HE;eU{UiQ@$Wm()V@7JGa8y3zw3!EgtPhXH%VLRF?>tR!Y+BrBc?V-PCOPl%1>
z#S2*B(3Ier7|-E6E(?nQ#tuv7u(jc+hC_-6!W@dhaN%O3nJ|VK0x){5+*tx&R+2#0
zm3%YVCkTZRAa0D|b&Jg-B*hi+Eau97LMrer@O!^HQiZ^*f-%AP2VMe`#@9h)g$1(s
z2{F(u$uZ3Muqd&+pf3qY82<fI6IF>Q3@=|gVTz=jib)Uy1<H0PwP8Wn1QM8@6Iu4+
z#mQLZzABhnweWz*@%Ud6#9Cw)N!Mmd=U=&1OBudaK=El2v>{iBiQ4glwZ@yln>U`l
z23+UxDQ;@R<#^-A9K@>_K0Z^CqMr`eor&5D2Icz8S}_%8M-lzEjG*IBBBf8LDqLx`
z7Z27JRTI-f2Bi%buh20BAKX4s(u-w-q46!NB-EvYQe(cC^zmVTA9v8`+N%b|?wVtX
zUt&{$o9?h9N&f11TcW;99#RlZRfh|zB-zFTA~5_(G}IMl;wCH0q@|ZDVM?OUaRJXu
zIkwW;W1O*|1oICv{)saOl$~NEG#V)T#<}39H>%;=<96e1JecW2o~>Ni2GUG<_0l>t
zmZ&Qc?8dF*)nxE?ae&s~8>fVlv@{etzc(KCsx^KL;A2HM#?|;nQi+8y4*(K&!@mtB
zf>`Os0~aR4u%O_fi4&#5%M+<2nkL|EDy_g)4U+k<lNvVx4@hwkb2zwauC#$(670t8
zWQj2d>Nxx?sFx^BNtQ*Jp=KbVbVA$2JB1#QqsZrfBU9wjR5AL67~KqB%gLkFR?DN!
zWzePIeVja+^u0XV0EDG6esI+|5=Mg7CFg)wck)d>DKxQw0!aQ6xCp%AlZU@am6u=%
z;lLIK`0vTm>{sic`N<W+)IPsL65sLgZT3@XicFZ$44pzWBGbStKRE$fhYVM|tXtww
z%h7N;-unX$hmCc<#<0{;-{I)IvD)F%Cz$fVOX7}oYXYT?`MgFkpR^J30qyQgV|AFd
zswsh=k<doN;rmZ`sp8Wj^-)cHya>*EaWTgDD$!iD6u*)<o{WL~Rj#x6h_!<Rll*3O
zD$rGCFN_~dGr|8_J4xC@=XKO?D2L~)g9CN<x&X8pcbNsYA_H9KnkMH9kocRRtMSe>
z>)|O9We_1wUT=*i;c@BC`0Dk<?xb$LTs!gh!WguqrCX!P0P893+#l%>SG);EK6j%T
z?z=$?6&oI3=At%S-sp832%}qUlx^i8ZE<Tri+3A|7W_@$G=lSb1wmTZotGrLwS@T9
zo2yv|T`ydMYj2%FE*-F2zIlMQ?o`U~`E_;f`57-Pl3*B#U#$1<Eu<!bf5iKfI8gB8
zf>f6o8_NM%?+Ok`jw7%{xK-d|31DgdnmI*+)CHz=i$FT(4e|*A+%AA*5_97)BgyKE
z<*LEA-;!!d$U;D%C-OK60v5)AIU1sdqx7+O=_XyVjm<Sha}b2M5kz9`*mOlB-1_ne
zy+7Areu;ypm7s9vV_*b(jRE`l@Z<S=XW#$DeV@vP;`4j>zd;*t;ZPdP2L8|x-v5gY
z7KcEsh`baFgb5tn-vlohLc2SD!8FMrj)sE4B^OUx$7kKs4tll1tLV?w)e#*J=ZoF^
z0;Cjj15j&&?l_hEF?^4-^n+esHbw^x>1Oc~$_JwU8#jZj5BOQEqbD~J?qQ@Fg3Ox{
zmjJdxAQGTo3983v7CSOZfVr|GB0x3<b5H=4m|~>dsh33mm16A`^<OYy&EKnRDD!Kc
z|B$`tsM|QP{#pahJFa+I;=?($-bGi(pR=*wCO#DUW|)h&@NjDQZ}k7FzBi8>DT#ET
z56E~~Tmc_5i(zxYMvIvQcU$4HJdmOqfV^KE=wV?163-4}3z7*_ebWqh4^1)1_lOqq
zFJ-5Cp9(tTq_VKrL~H{80S7XTU-TD-Uw$@!Z)H%0`LDF-Aw!@U8vDfmV1i)2WIrW{
zXL9*qZxDkeGP!IHhY6h<#^c3@v$<e5z=t_H5BrV}|LtQ+Jc|P+heUAz$Y4N7EfP-=
z{RS29*Dwo=7ZJf{3CQZ=iK2gEkx-dnBOswG3@c0V+NFNoSbX3k2`V1c6I^nG3a4!*
zrrhJ3zd7ZO%pej(_YA6Z?Wr4(#hlIf76a`i6445>eYX(Ia^#Yn-CF{_u>xv`4$3Gv
zlVB7(ApdF0{ARBb2W7S-i(s~VK&o>+>zm9X24kj$+8b_hKX8HGw3)o?39HTs1}vVH
zL8L{WYmY>ITKv6N-Tj-U>A{DYFi^jxI%`k?W~#|EzI<Bmx8BS}%a7YqD47%nq)<Kd
zN6~N5pXm)DQLnpJS9T30C~{CiHP(F-y3?-Aub<?rRG=Qt6AM~C;L!WH@KxMKhev*!
z)?GU{q0A&?$e~vSC6$5uznOg`voBmu`^QQklN-sBU1<9$Z~=C8XY<2&TmhS#0Hz%<
zlH}UL|HE>X2f&$1hB-JP15<psbcmgvHIG~e@;SLh2mgNSAX%T0)Q}1I>8-{p&)}08
zel>m+{&A~@uo*s11veN+|H~(N|IaV-=KX&^$rFE(7dr|by5hUDsmfpyM)dG$tHzS8
tK?yFY@DiurTBb`_Pj4p9_Sqa=@=IJj>W^FI_@E*@C1*VN-JbKq{{h*yKcxTw

delta 8609
zcmc&(cUV(N*S|?9p@gDArB^jUQPfaGN>Bokr3VmM3yK9$5ev#9novZI1PKF5C$xk_
zP!zDU3N~~tNFwekxPrT@;G(PI%ZiHc9mr)V?)yC7U*9}T?%dz}&Y3gkoSA#?aOpqz
z>Qa10mcl|i1OaLqHj@yvWT5pi4z>^^R||Arm>D$OL=ZEyF^`NOQ-`_`MB5zpAc*1g
z;p%ULAZkAh{$g5pz>Ckky3VlZ@B5BR`C?kLgIf;ol&;60@vTrZ)RO_e$Fa=PQHhUm
z#8KrdTA)R-DE1|zwguO!4D?j51Dv`Mc%-(jHd6hZA}LNT(1Vf>1J`5vBS<%j{4fED
zpMdO{fFw*n_D(?dO+eTa5Y7Z7aRQPw0ZE>Kq)b5ePe4*9AO|KO2O)xs1D=y+j`0!~
zn1)OS+B#(1?jf1<Q9;->B$E!dO_>ae{fNNtm^x{6KMu-T4WhorgiO)K4`_t}CZg3>
zo<lH{Lk0M|O$H%u>fnoxwJFBNqg|uf2krkCn=p?3XyJJFUfp5c@n{ETNrjv=W6^ka
z?<h8A3C8|u0(+(2u<kvhbVsQU>x^+6ulw#OHm18@zF;1^k7s|;F|c_r@BK8+s|lmV
zTQ2wNgQ)Qus`ZDh+dGQ=L}i#|xp_R>&|sLoZxlOFu3K@}cs5VYW{+YI$Vtn0k7s|8
zvpL^jv#^ExPd1EVV}8PPd&7i+<(q6Q8!)y93_2{|WGkQ@Fw-!$D-0tn-()MI9i!Qc
zU=U#WCR+*Z7|phUL4ak1z1L*e--#%Q<>gwrzvsdrz%s(_m9vvj5M$TM*?brTSVq_@
zO^0<SqaenfFK4?Y4sl1Z@5<RJqu9!F_AWReEF-#|%!YOEM?p;YW;uH=91xZfcD0<H
zih>w>ft=k92ZUvWZD>BM`v3}JY^I!T09C*;!sd-)V|jwH*FhDqe3Omk3C6BEI86RF
z8_N^SCw+&;vq>}mu^Y=1jQxrGzt~uwU~GKa|6(hnoiM8~whC`N+sVwr>dm;_0!|L|
zCbpMpVO#fze+(_H^-NJFrV?%xujK0M(?O!Y9td+)04v;x=yq)|=Y_nhchq1I!t7T=
z+f(JC6C4>4H(;RRmNeZXa&Z{9#q#!9_}B;9jHy7dQ~_|N5W$vN6!tU}gmyuu?x?n&
zybZ^I{n!=S4I>6$NuUu^&<gvYJdC?KS|K0$5qclqlhCSxR!y{Op;a5La-Hm<1{C}g
zohb|e?10K9k7}c@KqEMOXTmu8quK_e+J>XrI-}Y~@-|L>RY1_SJrr>ifHQ-Np)tm2
zg_#ELF)}m?&KP2_4ZSeBALIR>u5fRGgTbd7Gz88XD#i?eGargIdWEirE`?r##=_9S
z!C=^+K`>;{c<3SMcIX}GdgvMGe&`h#0_YJK2Ivi(HrQrKt`*Kn!n0X)6$aUlAg!=6
z+gTw{6rS~rwYlI&phlE(boTuvl9Cf1M(NKwbYo4DhmoYbMS&cg`O@LAhSln~q9ddp
zUaRaY^6RES8ji=<WK4=G8(iu8?A0Jm+S%)c-`D&VNhKlcrH@|Aq#-N6Ke6%b(Rn@1
z47al11$b_#ClVejWJ~j_rz6uB#!Dk_wQ-CY<Sc~O?>)ch*PT)dQp%;^J*oLT{Ug<g
zb54ghdbL+i%R|<azm)RIS~&{cK0@QHXa46}2;KS>%TG&RZ$0w*7}0)CJ8mT-mA7S6
z>a0JST@D@a6;y2d`H#(Yhc6DUJ)-iZ|5@+{QA*be8`)h`_F#{?=KGJqJsJ<!D3ec2
zp?)7ZCHU8#>Tm^v-&4C>CNFzP?DbJl^?RC?YMI(6p4z@PPf+hEQJkl)_R1H>u<?0X
zXc=@>74gkV3!257>T~ga{-!Cu)3dg?<QyLi>AkK`>=$T+_BY0y@h<4q9enlq>igb-
zM~f;WmOYLRjX3MU{)#wF&i;zL+Go?fsZnf8cwBhp?O6lqCTXmT+Rr_v1vA0CmoMjm
z*NYx$i(lHsv)JVxIk&FxH3GjP`|<N5S=BBPeL{ZSIu=gGvNI0#zS+hx=pfB6P?@Q8
zETT`K<nQ~~iGu#&viunc`9rst=wydLv7$6P@D$E{=}WdEGB4Fvr&bh^-}F;#ff4^J
z!V8h^e8~CXG-ti};K%BNkGT0f)4y(LFRny1IJ=HhsQeRAO%!sf<;?2u>WVymXO@<S
zA`xQk1$}~N8t&c;R}=VNZjyOEI(p62&PQz=b3^L$bcUL)h<F!a^3Ls%UY{3!Gu}d3
z*2U2&*u{SFG<(JFuF?uJ&dHXU8S;yTU-MN(W;5;jT8BVMV*?M_miHpL$-KRr^GrR*
z0Zq1<1IzSH%#$0H8j5x9{0*)i7MuHS#&=Y9w=kVFFJJyhvUcp@Ch+TgZkc7r+FQC4
zX}8zv#7nymw-wVBr(2PcjwD{vQ7+<L3KqqY**)QVpVL*R<vo>FDVR8NbTtPT1<@{A
zv)Ly1`*DM@=O@X|&<;J=#|jQ#rAR-}T6KeMxn!n)K(%|lWeY(Tn0w`x(8_en4US>K
zZ+S67)fk~^SAm0+VIM0{o7GxGJ>FVC!8IMzokmM4px~8Oij1>$7R>479++nD6oFKf
zv(bs0UEc7yyZR@+z!mS}S-&i3enmHUbnvr%UY>q_*pyxm=S<BdZ{Dm*2yZlbVb+ks
zeST?WsU=lrm4Pe%6*pk_x0Jr^+>_?=_Q>9TZx64<uW8Y=tSsDga{hk({fO1lEy|DQ
z;Mtl1dnHy)^?#OV-7H9WbE~`~_3#7AG$f3zfgE3Y-zPiBJCs|ZUX$&hmuML)R7{D=
z^@%zyYxar22XF2C+|9t@UNUJVnx_PdxiLb{o@5`T91ZeY8B1|tJQq3Bpy<e&{F)ml
zP*(M4QrgSJWTew8JFeYDPhH~FEnMqCAR{gF<HP}TdVSK8E+N!DB&YOfCTH3!c9BPD
zlG}3aOon=uwX1=4u*7av7ss)7RqMU%eSf9pG<WD5^|0^Fn$ci)?E?KF?aaCr?Q}Jd
z^uf|gEhWMfTkbCD(i$nQ4oz<FEmQLy#ND#C+eJ~SMLF&}+MfldBm~aqI|fBsN=Tw7
z_1~v+X3zA@jZzfLS`>~4n@_8p%m%?zxkz>tqq41tqD>KL#7oh+;2M{GDtjlA%!yZx
z5S<J~ILx14`DolcUwPBj$4G{Jxz<l3+$2Aa5h#*-BD$|r1S)H`RHxsldT1J4w)a+1
zA2L8G%~4@SMyCHvMmXwW%Z|QkBX~a)CAw<A*uJsl{a2*z^DL>X^ZXn3inHmn&pvM4
zRhIK4JZzVfi~m3ir?ax6TxjMbYUK%weR;c?$9w1klK#$LRObqA9am{9TOAr7TXow?
zNPJP@mAmS4;j*-cmG|O$>61;ngrSF6pDt%J4-^$HV&U(E>o~|%u85CvX0QUR7T255
zh*_-0%0^MLmkJJ@SehO?{)0fZCyw0dFYYR$;GVyt+m{>o(={uU@_ht~;;_y03?<q=
z$_mDw$ax0MdRAVzmxF4UQvsz?_pCwQ;vH>?bhW?_$Ju8h&iaRChFjp;(>;1DJ18$q
z`e>ZD@0vi&2l2_yDgFIg>HOx~E`%USinH*H?m6G;C9Ke|$W3-rYgnY}_Sj09^ilN&
z|0suAUCFk*dfl#^DJ`;!8l+N7`u@FfG&OOJk;q-bjS8aSw)j?f3dIM^8<?lH>4(Fl
zjgp_{>X_Z8CWh(77e$@R3OsDtq7YtTZxlIY#!Y`;t+SN^1#(Xer>M%sGA?niH{z~I
z7W^U7Ht8)Ei`isc3$9N<oMIb&ZBJfynDMn8ec}b~5qRpc=y}rCG9AZc%HH#m=b}&~
zG;|A^^a7(u>ar6NggXgti8fDal;(f@{A$h_$AEcV`X#R~4-l*3*G>}Fb#=1!QY+5b
z7dqEj(;6E^;!1H~S1skN?r~;R)|TM1=m`988x4w*ve}wDWU}Q=#<z=vIaA0oZWP(Z
zEKk9U`n<x8<fw;a?Y!%#chJ}~j&>zWZ+mfGnKC!NILA-%nhQmYtJL~jTCT*Iu2evB
zC)w*EfF_WINfQOie;lLTd5{;js5<y1v#+aag8&>NAttnd4Ow^}M4+ha$Qn>rMO2--
zRRnBJ`?vEtPIef_8vp&Wirf0AeQ9+E|I6onL4{s!w3kmdbDyT@tt@Nhm2zlzr`y_g
zOe&6w-QYia3-|2gAVqtdwYJxhRivg*dAYvT!D*@H&e;yh9|l<XraZnkU0BY#j-R1a
zu2jMbAkZfz`zmt}8n9$iykF$YeQW%ZeSI6)^r=I$5G;Ae%|Pd?7h^Mx_ACNVk90-T
z6cyT#CdHM-A}%%tBK{lZdOUyij=fZl<<VDH#7amPpydjbCwpzUuF$rtKK4H=XZ*1u
zW3d+TgQMmPrlwOvj_)4Fva;qoE|Y`0D~?ppc>XHTREBO9q4k?fRxjbRg6!F=f66fG
zAg!=bkHCrf=PIeq+n5_tJc^sD%7TciL2y}A*oxr$tD=ICPeJRW3Q35kWRF(lGwZ?0
zf3whVMQ8~Iq=^+IUj-w<zQ8N=II-5}rVPcwmtyT@=3FFm0i2caw1Gxs*An`aMv;<H
zz2CwQ6&*bh_+N9B<9+RIcS+8RHgb`I+j|0{XNoJY>vnN0?m8^E%gEbqch^fzxBRYM
zEUg}~&-?vS_?|P5^@ALnpFB|pcO1;vq#MV!2_Mik+N1uQ^aKCBELNgJn_H9oq+tf<
z;?mYE1Fog&%tPmGHrPrmW1_BYrB8hz*%@t_=@jDLm2)`ip=EtgusEm^e>*SE24Bw<
za?nV*{}zaoxVhL^>!Dk^wx9T8<g*rqToJ(^-S9a^Xp$YF3ErYX$uhFEW=HNoTiKd6
zMv<;JXK-e8cx~*f{0a04ms;LPyLZJhEpv}5a_>>SY52P4>O(?BmSKZ;%qNeA6%TC)
zmA}?qN-?%1{2b(2ZhDcS&v4)(T<%ICPso{6ba_{PU2Ae(Sc_tg#o(lU+0Id^oxy7h
zmNs@?F8prE3|e18;7(2}p8d{gIq&Z(S$wTN|0##2Em<E`Vo>mKe<o+ntU&I=pmj#|
zmJQaklLO4RTb^!UrrW6^L5}Lu#!A+rc)g5`tsD!b0*xqZl0Bilh49`-nHxcnNmG@h
z5dZt0F-CPAvWu0t=Y$IL`USD`@@H@iLeCd#_mni&CcaznsnsQW1-s2AaYO1f;)0s5
z$d68g_08E9?;=Zxm*b1Q2wQ0v_gfn>>|)B$X9W5hiwNzW(gwm~X-F0!ncs%K*S61p
zGtFhsM_H~UvQ@ZrkKK!amRpV2zca|zY~^fYs@<9%vP`Kqz_TWE8(ojsRS_(qPWL+L
zrIsjRKKVl;`5Y?_XS30bpz_J@<i8#U+kIp(jX1MYD@Z|CrO7Rd9N9$)sc)-5rQiR8
z;63%o{L-TvRFwlF@Qn$<?#HXFBsl_8k@?Af_P~c?@uutVbxNA}$urqar^vq_60OnP
zuyO63r53+|r<q(T2Yny+@QqH*r7M!r2SYohotouG=+_KoM=BI>J>$PTMxv$N?$U1G
zQjVVUD$<<JSKiTeSvbp|&X$Oeu8;oP`9nZKVPuF}uIPA|WL9>!@1y~4u<j4j#cjUv
zcb6IRn5X$n5(A{y)3zrxozCnMrL%Kb1owr`{sFU#KI9raj}WIkZ{yGuiu&Hnuf6tc
zM_%}^*NvZMur)WgC>7^><R`Y(RJ3urEUaF~Red_#{W@-U(^q6?4E1LZ*B!?jkwxF%
zVKZFXYEHa4F<0}AzN!PIJLp)2Lyr362;;dk3Yi@gq}?f^-6;~C)~)6%!d~aBxlctc
z1lwXy^lvX2XNb;Q(P>}r{GtX5E~1o+ez%lYMAxm&LRg^MwxBQq2lo{7owSx~J1toG
zB|<;KwG&Eq0f!29Mx~=m+w*yu-fD-%^Dap#0>#-k#WUYlJ`sy#TVySGe+C)pL%(7~
zp0UVu^PLst(OJlC3jFl?;GA<g9yB{IRx+KAAXav@S}vz?xb<MoK6~JtM+Ta%#G$*{
zz|D0n+DcqX0ry?0I8*Ssa2o7cLIKKbt<l^i6r3(_bzKg+7pu!9jFya*2m#|I^oAtt
z-7twQZct*ISAbl@pI*cF?qG#C#^HNYhmJkIL%%gg3jFwlWnHbzy90+?2TuEXqB{5q
z=D@_yOnc{WFJddd4$ojen*px)Y02laEoC++4$#%!@l9@Kkp+5)ZVE2>WBm8Y3{YJ~
zR2i#tM*ssjmyJi?1~9avzo~aCivijaFrAuvX9L|pGHTJBKnA$Cj|lbz2BEWg5D4TR
z;RR7i@Q45&mcX(OTn)O0_L-)cgX6(wz%iK!Y!Byx=w-S?^0tRdfOQCI^15OA&@UZ+
zraKzd`i0oyfX?>@z%MkWwkh-;5qw<NP;0RvQn|Kx3z1NJA;ts`Htn#eeX}D(L22#w
z^=mdQs$Cdwj7Ke6o?wTLiureiT>DsGByt$w^Qm#Zz$t<By>4UX%mDNx%meoh5CQX;
z!H}<SaO5aX1N~Qq`4Gk&!8Av(3&|Ky3&XTGe$%4e{<A@03K`}{%0|4F3gBnygODIB
zwKnaiszT=urxp<SnSdYnL)|d{`}dd2Th~%yR=!SkMenT-(BvrZ0Jbu|a)_YqfIW(9
z9>jV#AB4T8X|A|XF!ykcTt;w6j)H+hr{tK7i(zvPV;C_F#cEB{2sk`lZJg1&c+eUF
zZ|abKbMS;`HZ+EywkUlG9tHiAyucj3sa$;=pNhKz5(^E$k9=(2@A=<YQrqPo99}E^
zNfsah=z{!As$5Yo+A;;&S<qUzAM~;yW3jyV1sL0_%7<&;3LQEl%K)V3VncUkjpeP%
z9OgY6;vLMxcu%s&^61A#@kl1<;g!dTy)ar2kB;TXCJ)(<*TO^hk+jg$v9SX8kBt@h
z;6F3~o<1x9kIeoND9Q1k^iQJ52VM>Y0n@;bNkkR%5mQX_<K+&ymOp;s4&?6`a-F$q
z49X|Mb2cK(Uz{?yTwtyO&*LEG3)TRApi~HVg~X-iz`1ZS-b8&U2?`4__p}y5>xVbz
z2YzEVaWpgpjs-YVG+l0qLNWE<j9*raWqU?36$r`<NzW#DR!uP0QRvQy+XV6j)+*Tl
zv$g`6z#8x>h#(|a7lf1=lKw&u6aVoZjL+_pv0C=Cp%&#*Ymgwq<c2ei3XYdj6ri^O
z=Y$SAhu%_q+>o!eSEC>-|M+`Ibs6^fcvnUpD@Xw=%e97u!57}>KtaLF<(QyD#RS10
zPsjze*9;FjZd8+#m~j5dfA#O{lbC;J)(rU<9IC@YSX%RMV%ovq9C6O6VHFwmLn==F
zYf5^t$IsN}G!!*bN6&}HhNb58zZn5k>oD6LYvs8yVG*6Jg;Cw+<xZ-Z5b^KMePfd`
zlh15DNIz$;R1NPa;M93zHU~EU|BkWIg{KV~Gp2!>^HfE+645u6%4sAMv=xq4ScaYo
lS|{2D(@!UNGm;#mRk9<FxPV0$ym2w$&;<kZS9sy*e*t-%hp7Mn

diff --git a/Shaders/Private/CesiumBox.usf b/Shaders/Private/CesiumBox.usf
index 238ba13bc..7c09732b2 100644
--- a/Shaders/Private/CesiumBox.usf
+++ b/Shaders/Private/CesiumBox.usf
@@ -6,7 +6,7 @@
 	CesiumBox.usf: An implicit box shape that may be intersected by a ray.
 =============================================================================*/
 
-#include "CesiumRayIntersection.usf"
+#include "CesiumRayIntersectionList.usf"
 
 struct Box
 {
@@ -16,10 +16,11 @@ struct Box
   /**
    * Tests whether the input ray (Unit Space) intersects the box. Outputs the intersections in Unit Space.
    */
-  RayIntersections Intersect(in Ray R)
+  void Intersect(in Ray R, inout float4 Intersections[INTERSECTIONS_LENGTH], inout IntersectionListState ListState)
   {
+    ListState.Length = 2;
+    
     // Consider the box as the intersection of the space between 3 pairs of parallel planes.
-
     // Compute the distance along the ray to each plane.
     float3 t0 = (MinBounds - R.Origin) / R.Direction;
     float3 t1 = (MaxBounds - R.Origin) / R.Direction;
@@ -36,8 +37,10 @@ struct Box
 
     if (entryT > exitT)
     {
-      Intersection miss = NewMissedIntersection();
-      return NewRayIntersections(miss, miss);
+      result.Entry.t = NO_HIT;
+      result.Exit.t = NO_HIT;
+      setShapeIntersections(Intersections, ListState, 0, result);
+      return;
     }
 
     // Compute normals
@@ -50,7 +53,7 @@ struct Box
     result.Exit.Normal = float3(isFirstExit) * directions;
     result.Exit.t = exitT;
 
-    return result;
+    setShapeIntersections(Intersections, ListState, 0, result);
   }
 
   /**
diff --git a/Shaders/Private/CesiumEllipsoidRegion.usf b/Shaders/Private/CesiumEllipsoidRegion.usf
new file mode 100644
index 000000000..56b08b6c1
--- /dev/null
+++ b/Shaders/Private/CesiumEllipsoidRegion.usf
@@ -0,0 +1,835 @@
+#pragma once
+
+// Copyright 2020-2025 CesiumGS, Inc. and Contributors
+
+/*=============================================================================
+	CesiumEllipsoidRegion.usf: An implicit ellipsoid region that may be intersected by a ray.
+=============================================================================*/
+
+#include "CesiumRayIntersectionList.usf"
+
+struct EllipsoidRegion
+{
+  float3 MinBounds;
+  float3 MaxBounds;
+  
+  float3 RadiiUV;
+  float3 InverseRadiiUVSquared;
+  float InverseHeightDifferenceUV;
+
+  float EccentricitySquared;
+  float2 EvoluteScale;
+  float3 LongitudeUVExtents; // X = min, Y = max, Z = midpoint
+  float LongitudeUVScale;
+  float LongitudeUVOffset;
+  float LatitudeUVScale;
+  float LatitudeUVOffset;
+
+  uint LongitudeRangeFlag;
+  bool LongitudeMinHasDiscontinuity;
+  bool LongitudeMaxHasDiscontinuity;
+  bool LongitudeIsReversed;
+  uint LatitudeMinFlag;
+  uint LatitudeMaxFlag;
+
+  void Initialize(float3 InMinBounds, float3 InMaxBounds, float4 PackedData0, float4 PackedData1, float4 PackedData2, float4 PackedData3, float4 PackedData4, float4 PackedData5)
+  {
+    MinBounds = InMinBounds; // longitude, sine(latitude), height
+    MaxBounds = InMaxBounds; // longitude, sine(latitude), height
+
+    // Flags are packed in CesiumVoxelRendererComponent.cpp
+    LongitudeRangeFlag = round(PackedData0.x);
+    LongitudeMinHasDiscontinuity = bool(PackedData0.y);
+    LongitudeMaxHasDiscontinuity = bool(PackedData0.z);
+    LongitudeIsReversed = bool(PackedData0.w);
+
+    LatitudeMinFlag = round(PackedData1.x);
+    LatitudeMaxFlag = round(PackedData1.y);
+
+    EvoluteScale = PackedData1.zw;
+    RadiiUV = PackedData2.xyz;
+    InverseRadiiUVSquared = PackedData3.xyz;
+    InverseHeightDifferenceUV = PackedData3.w;
+    EccentricitySquared = PackedData4.w;
+    LongitudeUVExtents = PackedData4.xyz;
+    LongitudeUVScale = PackedData5.x;
+    LongitudeUVOffset = PackedData5.y;
+    LatitudeUVScale = PackedData5.z;
+    LatitudeUVOffset = PackedData5.w;
+  }
+
+  /**
+  * Tests whether the input ray intersects the ellipsoid at the given height, relative to original radii.
+  */
+  RayIntersections IntersectEllipsoidAtHeight(in Ray R, in float RelativeHeight, in bool IsConvex)
+  {
+    // Scale the ray by the ellipsoid axes to make it a unit sphere.
+    // Note: This approximates ellipsoid + height as an ellipsoid.
+    float3 radiiCorrection = RadiiUV / (RadiiUV + RelativeHeight);
+    float3 position = R.Origin * radiiCorrection;
+    float3 direction = R.Direction * radiiCorrection;
+
+    // Intersect the ellipsoid defined by x^2/a^2 + y^2/b^2 + z^2/c^2 = 1.
+    // Solve for the values of t where the ray intersects the ellipsoid, by substituting
+    // the ray equation o + t * d into the ellipsoid x/y values. Results in a quadratic
+    // equation on t where we can find the determinant.
+    float a = dot(direction, direction); // ~ 1.0 (or maybe 4.0 if ray is scaled)
+    float b = dot(direction, position); // roughly inside [-1.0, 1.0] when zoomed in
+    float c = dot(position, position) - 1.0; // ~ 0.0 when zoomed in.
+    float determinant = b * b - a * c; // ~ b * b when zoomed in
+
+    if (determinant < 0.0)
+    {
+      // Output a "normal" so that the shape entry information can be encoded.
+      // (See EncodeIntersection)
+      Intersection miss = NewIntersection(NO_HIT, normalize(direction));
+      return NewRayIntersections(miss, miss);
+    }
+
+    determinant = sqrt(determinant);
+    RayIntersections result = (RayIntersections) 0;
+
+    // Compute the larger t using standard formula.
+    // The other t may suffer from subtractive cancellation in the standard formula,
+    // so it is computed from the first t instead.
+    float t1 = (-b - sign(b) * determinant) / a;
+    float t2 = c / (a * t1);
+    result.Entry.t = min(t1, t2);
+    result.Exit.t = max(t1, t2);
+
+    float directionScale = IsConvex ? 1.0 : -1.0;
+    result.Entry.Normal = directionScale * normalize(position + result.Entry.t * direction);
+    result.Exit.Normal = directionScale * normalize(position + result.Exit.t * direction);
+    return result;
+  }
+
+  /**
+  * Intersects the plane at the specified longitude angle.
+  */
+  Intersection IntersectLongitudePlane(in Ray R, in float Angle, in bool PositiveNormal)
+  {
+    float normalSign = PositiveNormal ? 1.0 : -1.0;
+    float2 planeNormal = normalSign * float2(-sin(Angle), cos(Angle));
+
+    float approachRate = dot(R.Direction.xy, planeNormal);
+    float distance = -dot(R.Origin.xy, planeNormal);
+
+    float t = (approachRate == 0.0)
+        ? NO_HIT
+        : distance / approachRate;
+
+    return NewIntersection(t, float3(planeNormal, 0));
+  }
+
+  /**
+  * Intersects the space on one side of the specified longitude angle.
+  */
+  RayIntersections IntersectHalfSpace(in Ray R, in float Angle, in bool PositiveNormal)
+  {
+    Intersection intersection = IntersectLongitudePlane(R, Angle, PositiveNormal);
+    Intersection farSide = NewIntersection(INF_HIT, normalize(R.Direction));
+
+    bool hitFront = (intersection.t > 0.0) == (dot(R.Origin.xy, intersection.Normal.xy) > 0.0);
+    if (!hitFront)
+    {
+      return NewRayIntersections(intersection, farSide);
+    }
+    else
+    {
+      return NewRayIntersections(Multiply(farSide, -1), intersection);
+    }
+  }
+
+  /**
+  * Intersects a "flipped" wedge formed by the negative space of the specified angle
+  * bounds and returns up to four intersections. The "flipped" wedge is the union of
+  * two half-spaces defined at the given angles and represents a *negative* volume
+  * of over > 180 degrees.
+  */
+  void IntersectFlippedWedge(in Ray R, in float2 AngleBounds, out RayIntersections FirstIntersection, out RayIntersections SecondIntersection)
+  {
+    FirstIntersection = IntersectHalfSpace(R, AngleBounds.x, false);
+    SecondIntersection = IntersectHalfSpace(R, AngleBounds.y, true);
+  }
+
+  /**
+  * Intersects the wedge formed by the negative space of the min/max longitude, where
+  * maxAngle > minAngle + pi. The wedge is represented by two planes at such angles.
+  * There is an opposite "shadow wedge", i.e. the wedge formed at the *opposite* side
+  * of the planes' intersection, that must be specially handled.
+  */
+  RayIntersections IntersectRegularWedge(in Ray R, in float2 AngleBounds)
+  {
+    // Normals will point toward the "outside" (into the negative space)
+    Intersection intersect1 = IntersectLongitudePlane(R, AngleBounds.x, false);
+    Intersection intersect2 = IntersectLongitudePlane(R, AngleBounds.y, true);
+
+    // Note: the intersections could be in the "shadow" wedge, beyond the tip of
+    // the actual wedge.
+    Intersection first = intersect1;
+    Intersection last = intersect2;
+    if (first.t > last.t)
+    {
+      first = intersect2;
+      last = intersect1;
+    }
+
+    bool firstIntersectionAheadOfRay = first.t >= 0.0;
+    bool startedInsideFirst = dot(R.Origin.xy, first.Normal.xy) < 0.0;
+    bool isExitingFromInside = firstIntersectionAheadOfRay == startedInsideFirst;
+    bool lastIntersectionAheadOfRay = last.t > 0.0;
+    bool startedOutsideLast = dot(R.Origin.xy, last.Normal.xy) >= 0.0;
+    bool isEnteringFromOutside = lastIntersectionAheadOfRay == startedOutsideLast;
+
+    Intersection farSide = NewIntersection(INF_HIT, normalize(R.Direction));
+    Intersection miss = NewIntersection(NO_HIT, normalize(R.Direction));
+
+    if (isExitingFromInside && isEnteringFromOutside)
+    {
+      // Ray crosses both faces of negative wedge, exiting then entering the positive shape
+      return NewRayIntersections(first, last);
+    }
+    else if (!isExitingFromInside && isEnteringFromOutside)
+    {
+      // Ray starts inside wedge. last is in shadow wedge, and first is actually the entry
+      return NewRayIntersections(Multiply(farSide, -1), first);
+    }
+    else if (isExitingFromInside && !isEnteringFromOutside)
+    {
+      // First intersection was in the shadow wedge, so last is actually the exit
+      return NewRayIntersections(last, farSide);
+    }
+    else
+    { // !exitFromInside && !enterFromOutside
+      // Both intersections were in the shadow wedge
+      return NewRayIntersections(miss, miss);
+    }
+  }
+
+  bool HitsPositiveHalfPlane(in Ray R, in Intersection InputIntersection, in bool positiveNormal)
+  {
+    float normalSign = positiveNormal ? 1.0 : -1.0;
+    float2 planeDirection = float2(InputIntersection.Normal.y, -InputIntersection.Normal.x) * normalSign;
+    float2 hit = R.Origin.xy + InputIntersection.t * R.Direction.xy;
+    return dot(hit, planeDirection) > 0.0;
+  }
+
+  void IntersectHalfPlane(in Ray R, in float Angle, out RayIntersections FirstIntersection, out RayIntersections SecondIntersection)
+  {
+    Intersection intersection = IntersectLongitudePlane(R, Angle, true);
+    Intersection farSide = NewIntersection(INF_HIT, normalize(R.Direction));
+
+    if (HitsPositiveHalfPlane(R, intersection, true))
+    {
+      FirstIntersection.Entry = Multiply(farSide, -1);
+      FirstIntersection.Exit = NewIntersection(intersection.t, float3(-1.0 * intersection.Normal.xy, 0.0));
+      SecondIntersection.Entry = intersection;
+      SecondIntersection.Exit = farSide;
+    }
+    else
+    {
+      Intersection miss = NewIntersection(NO_HIT, normalize(R.Direction));
+      FirstIntersection.Entry = Multiply(farSide, -1);
+      FirstIntersection.Exit = farSide;
+      SecondIntersection.Entry = miss;
+      SecondIntersection.Exit = miss;
+    }
+  }
+
+  /**
+  * Given a circular quadric cone around the z-axis, with apex at the origin,
+  * find the parametric distance(s) along a ray where that ray intersects
+  * the cone. Computation of the normal is deferred to GetFlippedCone.
+  * 
+  * The cone opening angle is described by the squared cosine of its half-angle
+  * (the angle between the Z-axis and the surface)
+  */
+  RayIntersections IntersectQuadricCone(in Ray R, in float cosSquaredHalfAngle)
+  {
+    float3 o = R.Origin;
+    float3 d = R.Direction;
+
+    float sinSquaredHalfAngle = 1.0 - cosSquaredHalfAngle;
+
+    float aSin = d.z * d.z * sinSquaredHalfAngle;
+    float aCos = -dot(d.xy, d.xy) * cosSquaredHalfAngle;
+    float a = aSin + aCos;
+
+    float bSin = d.z * o.z * sinSquaredHalfAngle;
+    float bCos = -dot(o.xy, d.xy) * cosSquaredHalfAngle;
+    float b = bSin + bCos;
+
+    float cSin = o.z * o.z * sinSquaredHalfAngle;
+    float cCos = -dot(o.xy, o.xy) * cosSquaredHalfAngle;
+    float c = cSin + cCos;
+    // determinant = b * b - a * c. But bSin * bSin = aSin * cSin.
+    // Avoid subtractive cancellation by expanding to eliminate these terms
+    float determinant = 2.0 * bSin * bCos + bCos * bCos - aSin * cCos - aCos * cSin - aCos * cCos;
+
+    RayIntersections result = (RayIntersections) 0;
+    result.Entry.t = NO_HIT;
+    result.Exit.t = NO_HIT;
+
+    if (determinant < 0.0)
+    {
+      return result;
+    }
+
+    if (a == 0.0)
+    {
+      // Ray is parallel to cone surface if b == 0.0.
+      // Otherwise, ray has a single intersection on cone surface
+      if (b != 0.0)
+      {
+        result.Entry.t = -0.5 * c / b;
+      }
+    }
+
+    determinant = sqrt(determinant);
+
+    // Compute larger root using standard formula
+    float signB = b < 0.0 ? -1.0 : 1.0;
+    float t1 = (-b - signB * determinant) / a;
+    // The other root may suffer from subtractive cancellation in the standard formula.
+    // Compute it from the first root instead.
+    float t2 = c / (a * t1);
+    result.Entry.t = min(t1, t2);
+    result.Exit.t = max(t1, t2);
+    return result;
+  }
+
+  /**
+   * Given a point on a conical surface, find the surface normal at that point.
+   */
+  float3 GetConeNormal(in float3 Point, in bool IsConvex)
+  {
+    // Start with radial component pointing toward z-axis
+    float2 radial = -abs(Point.z) * normalize(Point.xy);
+    // Z component points toward opening of cone
+    float zSign = (Point.z < 0.0) ? -1.0 : 1.0;
+    float z = length(Point.xy) * zSign;
+    // Flip normal if the shape is convex
+    float flip = (IsConvex) ? -1.0 : 1.0;
+    return normalize(float3(radial, z) * flip);
+  }
+
+  /**
+   * Compute the shift between the ellipsoid origin and the apex of a cone of latitude
+   */
+  float GetLatitudeConeShift(in float sinLatitude)
+  {
+    // Find prime vertical radius of curvature: 
+    // the distance along the ellipsoid normal to the intersection with the z-axis
+    float x2 = EccentricitySquared * sinLatitude * sinLatitude;
+    float primeVerticalRadius = rsqrt(1.0 - x2);
+
+    // Compute a shift from the origin to the intersection of the cone with the z-axis
+    return primeVerticalRadius * EccentricitySquared * sinLatitude;
+  }
+
+  /**
+  * Tests if the input ray intersects the negative space outside of the cone. In other words,
+  * this captures the volume outside the cone, excluding the part of the ray that is inside the cone.
+  */ 
+  void IntersectFlippedCone(in Ray R, in float CosHalfAngle, out RayIntersections FirstIntersections, out RayIntersections SecondIntersections)
+  {
+    // Undo the scaling from ellipsoid to sphere
+    float3 position = R.Origin * RadiiUV;
+    float3 direction = R.Direction * RadiiUV;
+    // Shift the ray to account for the latitude cone not being centered at the Earth center
+    position.z += GetLatitudeConeShift(CosHalfAngle);
+
+    RayIntersections quadricIntersections = IntersectQuadricCone(R, CosHalfAngle * CosHalfAngle);
+
+    Intersection miss = NewIntersection(NO_HIT, normalize(direction));
+    Intersection farSide = NewIntersection(INF_HIT, normalize(direction));
+
+    // Initialize output with no intersections
+    FirstIntersections = (RayIntersections) 0;
+    FirstIntersections.Entry = Multiply(farSide, -1);
+    FirstIntersections.Exit = farSide;
+
+    SecondIntersections = NewRayIntersections(miss, miss);
+
+    if (quadricIntersections.Entry.t == NO_HIT)
+    {
+      return;
+    }
+
+    // Find the points of intersection
+    float3 p0 = position + quadricIntersections.Entry.t * direction;
+    float3 p1 = position + quadricIntersections.Exit.t * direction;
+
+    quadricIntersections.Entry.Normal = GetConeNormal(p0, true);
+    quadricIntersections.Exit.Normal = GetConeNormal(p1, true);
+
+    // shadow cone = the half of the quadric cone that is being discarded.
+    bool p0InShadowCone = sign(p0.z) != sign(CosHalfAngle);
+    bool p1InShadowCone = sign(p1.z) != sign(CosHalfAngle);
+
+    if (p0InShadowCone && p1InShadowCone)
+    {
+      // both points in shadow cone; no valid intersections
+      return;
+    }
+    else if (p0InShadowCone)
+    {
+      FirstIntersections.Exit = quadricIntersections.Exit;
+    }
+    else if (p1InShadowCone)
+    {
+      FirstIntersections.Entry = quadricIntersections.Entry;
+    }
+    else
+    {
+      FirstIntersections.Exit = quadricIntersections.Entry;
+      SecondIntersections.Entry = quadricIntersections.Exit;
+      SecondIntersections.Exit = farSide;
+    }
+  }
+
+  /**
+  * Tests if the input ray intersects the volume inside the quadric cone. The cone's orientation is
+  * dependent on the sign of the input angle.
+  */ 
+  RayIntersections IntersectCone(in Ray R, in float CosHalfAngle)
+  {
+    // Undo the scaling from ellipsoid to sphere
+    float3 position = R.Origin * RadiiUV;
+    float3 direction = R.Direction * RadiiUV;
+    // Shift the ray to account for the latitude cone not being centered at the Earth center
+    position.z += GetLatitudeConeShift(CosHalfAngle);
+
+    RayIntersections quadricIntersections = IntersectQuadricCone(R, CosHalfAngle * CosHalfAngle);
+
+    Intersection miss = NewIntersection(NO_HIT, normalize(direction));
+    Intersection farSide = NewIntersection(INF_HIT, normalize(direction));
+
+    if (quadricIntersections.Entry.t == NO_HIT)
+    {
+      return NewRayIntersections(miss, miss);
+    }
+
+    // Find the points of intersection
+    float3 p0 = position + quadricIntersections.Entry.t * direction;
+    float3 p1 = position + quadricIntersections.Exit.t * direction;
+
+    quadricIntersections.Entry.Normal = GetConeNormal(p0, false);
+    quadricIntersections.Exit.Normal = GetConeNormal(p1, false);
+
+    bool p0InShadowCone = sign(p0.z) != sign(CosHalfAngle);
+    bool p1InShadowCone = sign(p1.z) != sign(CosHalfAngle);
+
+    if (p0InShadowCone && p1InShadowCone)
+    {
+      return NewRayIntersections(miss, miss);
+    }
+    else if (p0InShadowCone)
+    {
+      return NewRayIntersections(quadricIntersections.Exit, farSide);
+    }
+    else if (p1InShadowCone)
+    {
+      return NewRayIntersections(Multiply(farSide, -1), quadricIntersections.Entry);
+    }
+    else
+    {
+      return NewRayIntersections(quadricIntersections.Entry, quadricIntersections.Exit);
+    }
+  }
+
+  /**
+  * Intersects an XY-plane at Z = 0. The given Z describes the direction of plane's normal.
+  */
+  RayIntersections IntersectZPlane(in Ray R, in float Z)
+  {
+    float t = -R.Origin.z / R.Direction.z;
+
+    bool startsOutside = sign(R.Origin.z) == sign(Z);
+    bool isEntering = (t >= 0.0) != startsOutside;
+
+    Intersection intersect = NewIntersection(t, float3(0.0, 0.0, Z));
+    Intersection farSide = NewIntersection(INF_HIT, normalize(R.Direction));
+
+    if (isEntering)
+    {
+      return NewRayIntersections(intersect, farSide);
+    }
+    else
+    {
+      return NewRayIntersections(Multiply(farSide, -1), intersect);
+    }
+  }
+  
+  
+#define ANGLE_EQUAL_ZERO 1
+#define ANGLE_UNDER_HALF 2
+#define ANGLE_HALF 3
+#define ANGLE_OVER_HALF 4
+
+  /**
+  * Tests whether the input ray (Unit Space) intersects the ellipsoid region. Outputs the intersections in Unit Space.
+  */
+  void Intersect(in Ray R, inout float4 Intersections[INTERSECTIONS_LENGTH], inout IntersectionListState ListState)
+  {
+    ListState.Length = 2;
+
+    // The region can be thought of as the space between two ellipsoids.
+    RayIntersections OuterResult = IntersectEllipsoidAtHeight(R, MaxBounds.z, true);
+    setShapeIntersections(Intersections, ListState, 0, OuterResult);
+    if (OuterResult.Entry.t == NO_HIT)
+    {
+      return;
+    }
+
+    ListState.Length += 2;
+    RayIntersections InnerResult = IntersectEllipsoidAtHeight(R, MinBounds.z, false);
+    if (InnerResult.Entry.t == NO_HIT)
+    {
+      setShapeIntersections(Intersections, ListState, 1, InnerResult);
+    }
+    else
+    {
+      // When the ellipsoid is large and thin it's possible for floating point math
+      // to cause the ray to intersect the inner ellipsoid before the outer ellipsoid. 
+      // To prevent this from happening, clamp the inner intersections to the outer ones
+      // and sandwich the inner ellipsoid intersection inside the outer ellipsoid 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 ellipsoid to be invisible because it will think the ray
+      // is still inside the inner (negative) ellipsoid after exiting the
+      // outer (positive) ellipsoid.
+      //
+      // 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.
+      InnerResult.Entry.t = max(InnerResult.Entry.t, OuterResult.Entry.t);
+      InnerResult.Exit.t = min(InnerResult.Exit.t, OuterResult.Exit.t);
+      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
+    }
+
+    // For min latitude, intersect a NEGATIVE cone at the bottom, based on the latitude value.
+    if (LatitudeMinFlag == ANGLE_UNDER_HALF)
+    {
+      ListState.Length += 2;
+      RayIntersections BottomConeResult = IntersectCone(R, MinBounds.y);
+      setShapeIntersections(Intersections, ListState, 2, BottomConeResult);
+    }
+    else if (LatitudeMinFlag == ANGLE_HALF)
+    {
+      ListState.Length += 2;
+      RayIntersections PlaneResult = IntersectZPlane(R, -1.0);
+      setShapeIntersections(Intersections, ListState, 2, PlaneResult);
+    }
+    else if (LatitudeMinFlag == ANGLE_OVER_HALF)
+    {
+      ListState.Length += 4;
+      RayIntersections FirstIntersection;
+      RayIntersections SecondIntersection;
+      IntersectFlippedCone(R, MinBounds.y, FirstIntersection, SecondIntersection);
+      setShapeIntersections(Intersections, ListState, 2, FirstIntersection);
+      setShapeIntersections(Intersections, ListState, 3, SecondIntersection);
+    }
+
+    // For max latitude, intersect a NEGATIVE cone at the top, based on the latitude value.
+    if (LatitudeMaxFlag == ANGLE_UNDER_HALF)
+    {
+      RayIntersections FirstIntersection;
+      RayIntersections SecondIntersection;
+      IntersectFlippedCone(R, MaxBounds.y, FirstIntersection, SecondIntersection);
+      // The array index depends on how many intersections were previously found.
+      [branch]
+      switch (ListState.Length)
+      {
+        case 6: // maxHeight + minHeight + minLatitude <= half
+          setShapeIntersections(Intersections, ListState, 3, FirstIntersection);
+          setShapeIntersections(Intersections, ListState, 4, SecondIntersection);
+          break;
+        case 8: // maxHeight + minHeight + minLatitude > half
+          // It makes no sense for min latitude to be in the top half, AND max latitude in the bottom half.
+          // The region is invalid, so make things easier and mark the shape as not hit.
+          Intersections[0].w = NO_HIT;
+          Intersections[1].w = NO_HIT;
+          ListState.Length = 0;
+          return;
+        case 4: // maxHeight + minHeight
+        default:
+          setShapeIntersections(Intersections, ListState, 2, FirstIntersection);
+          setShapeIntersections(Intersections, ListState, 3, SecondIntersection);
+          break;
+      }
+      ListState.Length += 4;
+    }
+    else if (LatitudeMaxFlag == ANGLE_HALF)
+    {
+      RayIntersections PlaneResult = IntersectZPlane(R, 1.0);
+      // The array index depends on how many intersections were previously found.
+      [branch]
+      switch (ListState.Length)
+      {
+        case 6: // maxHeight + minHeight + minLatitude <= half
+          // Not worth checking if minLatitude == half, just let the sort algorithm figure it out
+          setShapeIntersections(Intersections, ListState, 3, PlaneResult);
+          break;
+        case 8: // maxHeight + minHeight + minLatitude > half
+          // It makes no sense for min latitude to be over half, AND max latitude to be half.
+          // The region is invalid, so make things easier and mark the shape as not hit.
+          Intersections[0].w = NO_HIT;
+          Intersections[1].w = NO_HIT;
+          ListState.Length = 0;
+          return;
+        case 4: // maxHeight + minHeight
+        default:
+          setShapeIntersections(Intersections, ListState, 2, PlaneResult);
+          break;
+      }
+      ListState.Length += 2;
+    }
+    else if (LatitudeMaxFlag == ANGLE_OVER_HALF)
+    {
+      RayIntersections TopConeResult = IntersectCone(R, MaxBounds.y);
+      [branch]
+      switch (ListState.Length)
+      {
+        case 6: // maxHeight + minHeight + minLatitude <= half
+          setShapeIntersections(Intersections, ListState, 3, TopConeResult);
+          break;
+        case 8: // maxHeight + minHeight + minLatitude > half
+          setShapeIntersections(Intersections, ListState, 4, TopConeResult);
+          break;
+        case 4: // maxHeight + minHeight
+        default:
+          setShapeIntersections(Intersections, ListState, 2, TopConeResult);
+          break;
+      }
+      ListState.Length += 2;
+    }
+
+    float2 LongitudeBounds = float2(MinBounds.x, MaxBounds.x);
+    if (LongitudeRangeFlag == ANGLE_OVER_HALF)
+    {
+      // The shape's longitude range is over half, so we intersect a NEGATIVE shape that is under half.
+      RayIntersections WedgeResult = IntersectRegularWedge(R, LongitudeBounds);
+      // The array index depends on how many intersections were previously found.
+      [branch]
+      switch (ListState.Length)
+      {
+        case 6: // maxHeight + minHeight + (minLatitude <= half || maxLatitude >= half)
+          setShapeIntersections(Intersections, ListState, 3, WedgeResult);
+          break;
+        case 8: // maxHeight + minHeight + (minLatitude > half || maxLatitude < half || minLatitude <= half && maxLatitude >= half)
+          setShapeIntersections(Intersections, ListState, 4, WedgeResult);
+          break;
+        case 10: // maxHeight + minHeight + (minLatitude > half && maxLatitude > half || minLatitude < half + maxLatitude < half)
+          setShapeIntersections(Intersections, ListState, 5, WedgeResult);
+          break;
+        case 4: // maxHeight + minHeight
+        default:
+          setShapeIntersections(Intersections, ListState, 2, WedgeResult);
+          break;
+      }
+      ListState.Length += 2;
+    }
+    else if (LongitudeRangeFlag == ANGLE_UNDER_HALF)
+    {
+      // The shape's longitude range is under half, so we intersect a NEGATIVE shape that is over half.
+      RayIntersections FirstResult = (RayIntersections) 0;
+      RayIntersections SecondResult = (RayIntersections) 0;
+      IntersectFlippedWedge(R, LongitudeBounds, FirstResult, SecondResult);
+      // The array index depends on how many intersections were previously found.
+      [branch]
+      switch (ListState.Length)
+      {
+        case 6: // maxHeight + minHeight + (minLatitude <= half || maxLatitude >= half)
+          setShapeIntersections(Intersections, ListState, 3, FirstResult);
+          setShapeIntersections(Intersections, ListState, 4, SecondResult);
+          break;
+        case 8: // maxHeight + minHeight + (minLatitude > half || maxLatitude < half || minLatitude <= half && maxLatitude >= half)
+          setShapeIntersections(Intersections, ListState, 4, FirstResult);
+          setShapeIntersections(Intersections, ListState, 5, SecondResult);
+          break;
+        case 10: // maxHeight + minHeight + (minLatitude > half && maxLatitude > half || minLatitude < half + maxLatitude < half)
+          setShapeIntersections(Intersections, ListState, 5, FirstResult);
+          setShapeIntersections(Intersections, ListState, 6, SecondResult);
+          break;
+        case 4: // maxHeight + minHeight
+        default:
+          setShapeIntersections(Intersections, ListState, 2, FirstResult);
+          setShapeIntersections(Intersections, ListState, 3, SecondResult);
+          break;
+      }
+      ListState.Length += 4;
+    }
+    else if (LongitudeRangeFlag == ANGLE_EQUAL_ZERO)
+    {
+      RayIntersections FirstResult = (RayIntersections) 0;
+      RayIntersections SecondResult = (RayIntersections) 0;
+      IntersectHalfPlane(R, LongitudeBounds.x, FirstResult, SecondResult);
+      // The array index depends on how many intersections were previously found.
+      [branch]
+      switch (ListState.Length)
+      {
+        case 6: // maxHeight + minHeight + (minLatitude <= half || maxLatitude >= half)
+          setShapeIntersections(Intersections, ListState, 3, FirstResult);
+          setShapeIntersections(Intersections, ListState, 4, SecondResult);
+          break;
+        case 8: // maxHeight + minHeight + (minLatitude > half || maxLatitude < half || minLatitude <= half && maxLatitude >= half)
+          setShapeIntersections(Intersections, ListState, 4, FirstResult);
+          setShapeIntersections(Intersections, ListState, 5, SecondResult);
+          break;
+        case 10: // maxHeight + minHeight + (minLatitude > half && maxLatitude > half || minLatitude < half + maxLatitude < half)
+          setShapeIntersections(Intersections, ListState, 5, FirstResult);
+          setShapeIntersections(Intersections, ListState, 6, SecondResult);
+          break;
+        case 4: // maxHeight + minHeight
+        default:
+          setShapeIntersections(Intersections, ListState, 2, FirstResult);
+          setShapeIntersections(Intersections, ListState, 3, 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)
+  {
+    // Convert from [0, 1] to radians [-pi, pi]
+    float longitude = UV.x * CZM_TWO_PI;
+    if (LongitudeRangeFlag > 0)
+    {
+      longitude /= LongitudeUVScale;
+    }
+
+      // Convert from [0, 1] to radians [-pi/2, pi/2]
+    float latitude = UV.y * CZM_PI;
+    if (LatitudeMinFlag > 0 || LatitudeMaxFlag > 0)
+    {
+      latitude /= LatitudeUVScale;
+    }
+    
+    float height = UV.z / InverseHeightDifferenceUV;
+
+    return float3(longitude, latitude, height);
+  }
+
+  /**
+    * Given a specified point, gets the nearest point (XY) on the ellipse using a robust
+    * iterative solution without trig functions, as well as the radius of curvature
+    * at that point (Z).
+    * https://github.com/0xfaded/ellipse_demo/issues/1
+    * https://stackoverflow.com/questions/22959698/distance-from-given-point-to-given-ellipse
+    */
+  float3 GetNearestPointAndRadiusOnEllipse(float2 pos, float2 radii)
+  {
+    float2 p = abs(pos);
+    float2 inverseRadii = 1.0 / radii;
+
+    // We describe the ellipse parametrically: v = radii * vec2(cos(t), sin(t))
+    // but store the cos and sin of t in a vec2 for efficiency.
+    // Initial guess: t = pi/4
+    float2 tTrigs = 0.7071067811865476;
+    // Initial guess of point on ellipsoid
+    float2 v = radii * tTrigs;
+    // Center of curvature of the ellipse at v
+    float2 evolute = EvoluteScale * tTrigs * tTrigs * tTrigs;
+
+    const int iterations = 3;
+    for (int i = 0; i < iterations; ++i)
+    {
+      // Find the (approximate) intersection of p - evolute with the ellipsoid.
+      float2 q = normalize(p - evolute) * length(v - evolute);
+      // Update the estimate of t.
+      tTrigs = (q + evolute) * inverseRadii;
+      tTrigs = normalize(clamp(tTrigs, 0.0, 1.0));
+      v = radii * tTrigs;
+      evolute = EvoluteScale * tTrigs * tTrigs * tTrigs;
+    }
+
+    return float3(v * sign(pos), length(v - evolute));
+  }
+  
+  
+ /**
+   * 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".
+    // Get the position in unit space, undo the scaling from ellipsoid to sphere
+    float3 position = UVToUnit(PositionUV);
+    position *= RadiiUV;
+
+    float longitude = atan2(position.y, position.x);
+    float3 east = normalize(float3(-position.y, position.x, 0.0));
+
+    // Convert the 3D position to a 2D position relative to the ellipse (radii.x, radii.z)
+    // (assume radii.y == radii.x) and find the nearest point on the ellipse and its normal
+    float distanceFromZAxis = length(position.xy);
+    float2 positionEllipse = float2(distanceFromZAxis, position.z);
+    float3 surfacePointAndRadius = GetNearestPointAndRadiusOnEllipse(positionEllipse, RadiiUV.xz);
+    float2 surfacePoint = surfacePointAndRadius.xy;
+
+    float2 normal2d = normalize(surfacePoint * InverseRadiiUVSquared.xz);
+    float latitude = atan2(normal2d.y, normal2d.x);
+    float3 north = float3(-normal2d.y * normalize(position.xy), abs(normal2d.x));
+
+    float heightSign = length(positionEllipse) < length(surfacePoint) ? -1.0 : 1.0;
+    float height = heightSign * length(positionEllipse - surfacePoint);
+    float3 up = normalize(cross(east, north));
+
+    JacobianT = float3x3(east / distanceFromZAxis, north / (surfacePointAndRadius.z + height), up);
+    // Transpose because HLSL matrices are constructed in row-major order.
+    JacobianT = transpose(JacobianT);
+
+    // Then convert Unit Space to Shape UV Space
+    // Longitude: shift & scale to [0, 1]
+    float longitudeUV = (longitude + CZM_PI) / CZM_TWO_PI;
+
+    // Correct the angle when max < min
+    // Technically this should compare against min longitude - but it has precision problems so compare against the middle of empty space.
+    if (LongitudeIsReversed)
+    {
+      longitudeUV += float(longitudeUV < LongitudeUVExtents.z);
+    }
+
+    // Avoid flickering from reading voxels from both sides of the -pi/+pi discontinuity.
+    if (LongitudeMinHasDiscontinuity)
+    {
+      longitudeUV = longitudeUV > LongitudeUVExtents.z ? LongitudeUVExtents.x : longitudeUV;
+    }
+
+    if (LongitudeMaxHasDiscontinuity)
+    {
+      longitudeUV = longitudeUV < LongitudeUVExtents.z ? LongitudeUVExtents.y : longitudeUV;
+    }
+
+    if (LongitudeRangeFlag > 0)
+    {
+      longitudeUV = longitudeUV * LongitudeUVScale + LongitudeUVOffset;
+    }
+
+    // Latitude: shift and scale to [0, 1]
+    float latitudeUV = (latitude + CZM_PI_OVER_TWO) / CZM_PI;
+    if (LatitudeMinFlag > 0 || LatitudeMaxFlag > 0)
+    {
+      latitudeUV = latitudeUV * LatitudeUVScale + LatitudeUVOffset;
+    }
+
+    // Height: scale to the range [0, 1]
+    float heightUV = 1.0 + height * InverseHeightDifferenceUV;
+
+    return float3(longitudeUV, latitudeUV, heightUV);
+  }
+};
diff --git a/Shaders/Private/CesiumRayIntersectionList.usf b/Shaders/Private/CesiumRayIntersectionList.usf
new file mode 100644
index 000000000..318484600
--- /dev/null
+++ b/Shaders/Private/CesiumRayIntersectionList.usf
@@ -0,0 +1,182 @@
+#pragma once
+
+// Copyright 2020-2024 CesiumGS, Inc. and Contributors
+
+#include "CesiumRayIntersection.usf"
+
+/*===========================
+  CesiumRayIntersectionList.usf: Utility for tracking ray intersections across a complex shape.
+=============================*/
+
+// SHAPE_INTERSECTIONS is the number of ray-*shape* intersections (i.e., the volume intersection pairs),
+// INTERSECTIONS_LENGTH is the number of ray-*surface* intersections.
+#define SHAPE_INTERSECTIONS 7
+#define INTERSECTIONS_LENGTH SHAPE_INTERSECTIONS * 2
+
+/**
+* Holds state while iterating through the top-level IntersectionList. Used to interact with an array of encoded ray-surface intersections.
+* (See EncodeIntersection).
+*/
+struct IntersectionListState
+{
+  int Length; // Set based on ShapeConstant
+  
+  // Don't access the other member variables directly - call the functions on this struct instead.
+  int Index; // Emulates dynamic indexing.
+  
+  // The variables below relate to the shapes that the intersection is inside (counting when it's on the surface itself)
+  // For example, given a hollow ellipsoid volume, count = 1 on the outer ellipsoid, 2 on the inner ellipsoid.
+  int SurroundingShapeCount;
+  bool IsInsidePositiveShape; // will be true as long as it is inside any positive shape.
+  
+  /**
+  * Intersections are encoded as float4s:
+  * - .xyz for the surface normal at the intersection point
+  * - .w for the T value
+  * The normal's scale encodes the shape intersection type: 
+  * length(intersection.xyz) = 1: positive shape entry  
+  * length(intersection.xyz) = 2: positive shape exit
+  * length(intersection.xyz) = 3: negative shape entry  
+  * length(intersection.xyz) = 4: negative shape exit
+  *
+  * When the voxel volume is hollow, the "positive" shape is the original volume.
+  * The "negative" shape is subtracted from the positive shape.
+  */
+  float4 EncodeIntersection(in Intersection Input, bool IsPositive, bool IsEntering)
+  {
+    float scale = float(!IsPositive) * 2.0 + float(!IsEntering) + 1.0;
+    return float4(Input.Normal * scale, Input.t);
+  }
+
+  /**
+  * Sort the intersections from min T to max T with bubble sort. Also prepares for iteration
+  * over the intersections.
+  *
+  * Note: If this sorting function changes, some of the intersection tests may need to be updated.
+  * Search for "Sort()" to find those areas.
+  */
+  void Sort(inout float4 Data[INTERSECTIONS_LENGTH])
+  {
+    const int sortPasses = INTERSECTIONS_LENGTH - 1;
+    for (int n = sortPasses; n > 0; --n)
+    {
+      // Skip to n = Length - 1
+      if (n >= Length)
+      {
+        continue;
+      }
+      
+      for (int i = 0; i < sortPasses; ++i)
+      {
+      // The loop should be: for (i = 0; i < n; ++i) {...} but since loops with
+      // non-constant conditions are not allowed, this breaks early instead.
+        if (i >= n)
+        {
+          break;
+        }
+
+        float4 first = Data[i + 0];
+        float4 second = Data[i + 1];
+
+        bool inOrder = first.w <= second.w;
+        Data[i + 0] = inOrder ? first : second;
+        Data[i + 1] = inOrder ? second : first;
+      }
+    }
+
+    // Prepare initial state for GetNextIntersections()
+    Index = 0;
+    SurroundingShapeCount = 0;
+    IsInsidePositiveShape = false;
+  }
+
+  RayIntersections GetFirstIntersections(in float4 Data[INTERSECTIONS_LENGTH])
+  {
+    RayIntersections result = (RayIntersections) 0;
+    result.Entry.t = Data[0].w;
+    result.Entry.Normal = normalize(Data[0].xyz);
+    result.Exit.t = Data[1].w;
+    result.Exit.Normal = normalize(Data[1].xyz);
+
+    return result;
+  }
+
+  /**
+  * Gets the intersection at the current value of Index, while managing the state of the ray's
+  * trajectory with respect to the intersected shapes.
+  */
+  RayIntersections GetNextIntersections(in float4 Data[INTERSECTIONS_LENGTH])
+  {
+    RayIntersections result = (RayIntersections) 0;
+    result.Entry.t = NO_HIT;
+    result.Exit.t = NO_HIT;
+
+    if (Index >= Length)
+    {
+      return result;
+    }
+
+    float4 surfaceIntersection = float4(0, 0, 0, NO_HIT);
+
+    for (int i = 0; i < INTERSECTIONS_LENGTH; ++i)
+    {
+      // The loop should be: for (i = index; i < loopCount; ++i) {...} but it's not possible
+      // to loop with non-constant condition. Instead, continue until i = index.
+      if (i < Index)
+      {
+        continue;
+      }
+
+      Index = i + 1;
+
+      surfaceIntersection = Data[i];
+      // Maps from [1-4] -> [0-3] (see EncodeIntersection for the types)
+      int intersectionType = int(length(surfaceIntersection.xyz) - 0.5);
+      bool isCurrentShapePositive = intersectionType < 2;
+      bool isEnteringShape = (intersectionType % 2) == 0;
+
+      SurroundingShapeCount += isEnteringShape ? +1 : -1;
+      IsInsidePositiveShape = isCurrentShapePositive ? isEnteringShape : IsInsidePositiveShape;
+
+      // True if entering positive shape or exiting negative shape
+      if (IsInsidePositiveShape && isEnteringShape == isCurrentShapePositive)
+      {
+        result.Entry.t = surfaceIntersection.w;
+        result.Entry.Normal = normalize(surfaceIntersection.xyz);
+      }
+
+      // True if exiting the outermost positive shape
+      bool isExitingOutermostShape = !isEnteringShape && isCurrentShapePositive && SurroundingShapeCount == 0;
+      // True if entering negative shape while being inside a positive one
+      bool isEnteringNegativeFromPositive = isEnteringShape && !isCurrentShapePositive && SurroundingShapeCount == 2 && IsInsidePositiveShape;
+      
+      if (isExitingOutermostShape || isEnteringNegativeFromPositive)
+      {
+        result.Exit.t = surfaceIntersection.w;
+        result.Exit.Normal = normalize(surfaceIntersection.xyz);
+        // Entry and exit have been found, so the loop can stop
+        if (isExitingOutermostShape)
+        {
+          // After exiting the outermost positive shape, there is nothing left to intersect. Jump to the end.
+          Index = INTERSECTIONS_LENGTH;
+        }
+        break;
+      }
+      // Otherwise, keep searching for the correct exit.
+    }
+    
+    return result;
+  }
+};
+
+// Use defines instead of real functions to get array access with a non-constant index.
+
+/**
+* Encodes and stores a single intersection.
+*/
+#define setSurfaceIntersection(/*inout float4[]*/ list, /*in IntersectionListState*/ state, /*int*/ index, /*Intersection*/ intersection, /*bool*/ isPositive, /*bool*/ isEntering) (list)[(index)] = (state).EncodeIntersection((intersection), (isPositive), (isEntering))
+
+/**
+* Encodes and stores the given shape intersections, i.e., the intersections where a ray enters and exits a volume.
+*/
+#define setShapeIntersections(/*inout float4[]*/ list, /*in IntersectionListState*/ state, /*int*/ pairIndex, /*RayIntersection*/ intersections) (list)[(pairIndex) * 2 + 0] = (state).EncodeIntersection((intersections).Entry, (pairIndex) == 0, true); (list)[(pairIndex) * 2 + 1] = (state).EncodeIntersection((intersections).Exit, (pairIndex) == 0, false)
diff --git a/Shaders/Private/CesiumShape.usf b/Shaders/Private/CesiumShape.usf
index 02d60d887..49998e2ab 100644
--- a/Shaders/Private/CesiumShape.usf
+++ b/Shaders/Private/CesiumShape.usf
@@ -8,64 +8,83 @@
 
 #include "CesiumShapeConstants.usf"
 #include "CesiumBox.usf"
-
-/**
-* Converts a position from Unit Shape Space coordinates to UV Space.
-* [-1, -1] => [0, 1]
-*/
-float3 UnitToUV(float3 UnitPosition)
-{
-  return 0.5 * UnitPosition + 0.5;
-}
-
-/**
-* Converts a position from UV Space coordinates to Unit Shape Space.
-* [0, 1] => [-1, -1]
-*/
-float3 UVToUnit(float3 UVPosition)
-{
-  return 2.0 * UVPosition - 1.0;
-}
+#include "CesiumEllipsoidRegion.usf"
 
 struct Shape
 {
   int ShapeConstant;
   Box BoxShape;
+  EllipsoidRegion RegionShape;
+  IntersectionListState ListState;
   
   /**
   * Interpret the input parameters according to the voxel grid shape.
   */
-  void Initialize(in int InShapeConstant)
+  void Initialize(in int InShapeConstant, float3 InMinBounds, float3 InMaxBounds, float4 PackedData0, float4 PackedData1, float4 PackedData2, float4 PackedData3, float4 PackedData4, float4 PackedData5)
   {
     ShapeConstant = InShapeConstant;
-    
-    if (ShapeConstant == BOX)
+
+    [branch]
+    switch (ShapeConstant)
     {
+      case BOX:
       // Initialize with default unit box bounds.
-      BoxShape.MinBounds = -1;
-      BoxShape.MaxBounds = 1;
+        BoxShape.MinBounds = -1;
+        BoxShape.MaxBounds = 1;
+        break;
+      case ELLIPSOID:
+        RegionShape.Initialize(InMinBounds, InMaxBounds, PackedData0, PackedData1, PackedData2, PackedData3, PackedData4, PackedData5);
+        break;
     }
   }
     
   /**
   * Tests whether the input ray (Unit Space) intersects the shape.
   */
-  RayIntersections Intersect(in Ray R)
-  {
-    RayIntersections result;
-    
+  RayIntersections Intersect(in Ray R, out float4 Intersections[INTERSECTIONS_LENGTH])
+  { 
     [branch]
     switch (ShapeConstant)
     {
       case BOX:
-        result = BoxShape.Intersect(R);
+        BoxShape.Intersect(R, Intersections, ListState);
+        break;
+      case ELLIPSOID:
+        RegionShape.Intersect(R, Intersections, ListState);
         break;
       default:
         return NewRayIntersections(NewMissedIntersection(), NewMissedIntersection());
     }
-    
-    // Set start to 0.0 when ray is inside the shape.
-    result.Entry.t = max(result.Entry.t, 0.0);
+
+    RayIntersections result = ListState.GetFirstIntersections(Intersections);
+    if (result.Entry.t == NO_HIT)
+    {
+      return result;
+    }
+
+    // Don't bother with sorting if the positive shape was missed.
+    // Box intersection is straightforward and doesn't require sorting.
+    // Currently, cylinders do not require sorting, but they will once clipping / exaggeration
+    // is supported.
+    if (ShapeConstant == ELLIPSOID)
+    {
+      ListState.Sort(Intersections);
+      for (int i = 0; i < SHAPE_INTERSECTIONS; ++i)
+      {
+        result = ListState.GetNextIntersections(Intersections);
+        if (result.Exit.t > 0.0)
+        {
+          // Set start to 0.0 when ray is inside the shape.
+          result.Entry.t = max(result.Entry.t, 0.0);
+          break;
+        }
+      }
+    }
+    else
+    {
+      // Set start to 0.0 when ray is inside the shape.
+      result.Entry.t = max(result.Entry.t, 0.0);
+    }
 
     return result;
   }
@@ -75,8 +94,15 @@ struct Shape
   */
   float3 ScaleUVToShapeUVSpace(in float3 UV)
   {
-    // This is trivial for boxes, but will become relevant for more complex shapes.
-    return UV;
+    [branch]
+    switch (ShapeConstant)
+    {
+      case ELLIPSOID:
+        return RegionShape.ScaleUVToShapeUVSpace(UV);
+      case BOX:
+      default:
+        return UV;
+    }
   }
   
   /**
@@ -89,6 +115,8 @@ struct Shape
     {
       case BOX:
         return BoxShape.ConvertUVToShapeUVSpace(UVPosition, JacobianT);
+      case ELLIPSOID:
+        return RegionShape.ConvertUVToShapeUVSpace(UVPosition, JacobianT);
       default:
         // Default return
         JacobianT = float3x3(1, 0, 0, 0, 1, 0, 0, 0, 1);
diff --git a/Shaders/Private/CesiumVectorUtility.usf b/Shaders/Private/CesiumVectorUtility.usf
index 5fb160b4f..3dbeddb0b 100644
--- a/Shaders/Private/CesiumVectorUtility.usf
+++ b/Shaders/Private/CesiumVectorUtility.usf
@@ -34,3 +34,21 @@ int ConstructInt(in float2 value)
 {
   return int(value.x * 255.0) + (256 * int(value.y * 255.0));
 }
+
+/**
+* Converts a position from Unit Shape Space coordinates to UV Space.
+* [-1, -1] => [0, 1]
+*/
+float3 UnitToUV(float3 UnitPosition)
+{
+  return 0.5 * UnitPosition + 0.5;
+}
+
+/**
+* Converts a position from UV Space coordinates to Unit Shape Space.
+* [0, 1] => [-1, -1]
+*/
+float3 UVToUnit(float3 UVPosition)
+{
+  return 2.0 * UVPosition - 1.0;
+}
diff --git a/Shaders/Private/CesiumVoxelTemplate.usf b/Shaders/Private/CesiumVoxelTemplate.usf
index d28fe9613..1353202d3 100644
--- a/Shaders/Private/CesiumVoxelTemplate.usf
+++ b/Shaders/Private/CesiumVoxelTemplate.usf
@@ -96,12 +96,21 @@ struct VoxelMegatextures
   MAIN FUNCTION BODY
 =============================*/
 
+// HLSL does not like array struct members, so the data has to be stored at the top-level.
+// The size is also hardcoded because dynamically sized arrays are not allowed.
+float4 miss = float4(0, 0, 0, NO_HIT);
+float4 IntersectionList[INTERSECTIONS_LENGTH] =
+{
+  miss, miss, miss, miss, miss, miss, miss,
+  miss, miss, miss, miss, miss, miss, miss,
+};
+
 #define STEP_COUNT_MAX 1000
 #define ALPHA_ACCUMULATION_MAX 0.98 // Must be > 0.0 and <= 1.0
 
 VoxelOctree Octree;
 Octree.GridShape = (Shape) 0;
-Octree.GridShape.Initialize(ShapeConstant);
+Octree.GridShape.Initialize(ShapeConstant, ShapeMinBounds, ShapeMaxBounds, PackedData0, PackedData1, PackedData2, PackedData3, PackedData4, PackedData5);
 
 Ray R = (Ray) 0;
 R.Origin = RayOrigin;
@@ -132,7 +141,7 @@ R.Direction = RayDirection;
 // Shape UV Space: Voxel space from [0, 1] conforming to the actual voxel volume. The volume could be a box, a part of a cylinder, or a region on an ellipsoid.
 // (Vanilla) UV Space: Voxel space within an untransformed voxel octree. This can be envisioned as a simple cube spanning [0, 1] in three dimensions.
 
-RayIntersections Intersections = Octree.GridShape.Intersect(R);
+RayIntersections Intersections = Octree.GridShape.Intersect(R, IntersectionList);
 if (Intersections.Entry.t == NO_HIT) {
   return 0;
 }
@@ -157,6 +166,7 @@ switch (ShapeConstant) {
     DataTextures.PaddingBefore = round(PaddingBefore.xzy);
     DataTextures.PaddingAfter = round(PaddingAfter.xzy);
     break;
+  case ELLIPSOID:
   default:
     DataTextures.GridDimensions = round(GridDimensions);
     DataTextures.PaddingBefore = round(PaddingBefore);
@@ -174,7 +184,13 @@ float3 PositionUV = R.Origin + CurrentT * R.Direction;
 float3x3 JacobianT;
 float3 PositionShapeUVSpace = Octree.GridShape.ConvertUVToShapeUVSpace(PositionUV, JacobianT);
 
+// For ellipsoids, the UV direction has been scaled to a space where the ellipsoid is a sphere.
+// Undo this scaling to get the raw direction.
 float3 RawDirection = R.Direction;
+if (ShapeConstant == ELLIPSOID)
+{
+  RawDirection *= Octree.GridShape.RegionShape.RadiiUV;
+}
 
 OctreeTraversal Traversal;
 TileSample Sample;
@@ -209,7 +225,19 @@ for (int step = 0; step < STEP_COUNT_MAX; step++) {
   // Keep raymarching
   CurrentT += NextIntersection.t;
   if (CurrentT > EndT) {
-    break;
+    if (ShapeConstant == ELLIPSOID) {
+      // For CYLINDER: once clipping / exaggeration is supported, this must be done for cylinders too.
+      Intersections = Octree.GridShape.ListState.GetNextIntersections(IntersectionList);
+      if (Intersections.Entry.t == NO_HIT) {
+        break;
+      } else {
+        // Found another intersection. Resume raymarching there
+        CurrentT = Intersections.Entry.t;
+        EndT = Intersections.Exit.t;
+      } 
+    } else {
+        break;
+    }
   }
 
   PositionUV = R.Origin + CurrentT * R.Direction;
diff --git a/Source/CesiumRuntime/Private/Cesium3DTileset.cpp b/Source/CesiumRuntime/Private/Cesium3DTileset.cpp
index 2f76ce756..a77b2a16d 100644
--- a/Source/CesiumRuntime/Private/Cesium3DTileset.cpp
+++ b/Source/CesiumRuntime/Private/Cesium3DTileset.cpp
@@ -674,7 +674,7 @@ void ACesium3DTileset::OnFocusEditorViewportOnThis() {
   };
 
   const Cesium3DTilesSelection::Tile* pRootTile =
-      this->_pTileset->getRootTile();
+      this->_pTileset ? this->_pTileset->getRootTile() : nullptr;
   if (!pRootTile) {
     return;
   }
diff --git a/Source/CesiumRuntime/Private/CesiumVoxelRendererComponent.cpp b/Source/CesiumRuntime/Private/CesiumVoxelRendererComponent.cpp
index a9874a7e8..823945c64 100644
--- a/Source/CesiumRuntime/Private/CesiumVoxelRendererComponent.cpp
+++ b/Source/CesiumRuntime/Private/CesiumVoxelRendererComponent.cpp
@@ -484,6 +484,7 @@ UCesiumVoxelRendererComponent::CreateVoxelMaterial(
       FName("VoxelMaterial"));
   pVoxelMaterial->SetFlags(
       RF_Transient | RF_DuplicateTransient | RF_TextExportTransient);
+  pVoxelComponent->_pTileset = pTilesetActor;
 
   EVoxelGridShape shape = pVoxelComponent->Options.gridShape;
 
@@ -858,7 +859,7 @@ void UCesiumVoxelRendererComponent::UpdateTiles(
 
         // Don't create the missing node just yet? It may not be added to the
         // tree depending on the priority of other nodes.
-        priorityQueue.push({pVoxel, sse, computePriority(sse)});
+        priorityQueue.push({pVoxel, sse, computePriority(pVoxel->TileId, sse)});
       });
 
   if (this->_visibleTileQueue.empty()) {
@@ -880,8 +881,8 @@ void UCesiumVoxelRendererComponent::UpdateTiles(
         if (!pRight) {
           return true;
         }
-        return computePriority(pLeft->lastKnownScreenSpaceError) >
-               computePriority(pRight->lastKnownScreenSpaceError);
+        return computePriority(lhs, pLeft->lastKnownScreenSpaceError) >
+               computePriority(rhs, pRight->lastKnownScreenSpaceError);
       });
 
   size_t existingNodeCount = this->_loadedNodeIds.size();
@@ -1068,6 +1069,12 @@ void UCesiumVoxelRendererComponent::UpdateTransformFromCesium(
   }
 }
 
-double UCesiumVoxelRendererComponent::computePriority(double sse) {
-  return 10.0 * sse / (sse + 1.0);
+double UCesiumVoxelRendererComponent::computePriority(
+    const CesiumGeometry::OctreeTileID& tileId,
+    double sse) {
+  // This heuristic is intentionally biased towards tiles with lower levels.
+  // Without this, tilesets with many leaf tiles will kick all of the lower
+  // level detail tiles from the megatexture, resulting in holes or other
+  // artifacts.
+  return sse / (sse + 1.0 + tileId.level);
 }
diff --git a/Source/CesiumRuntime/Private/CesiumVoxelRendererComponent.h b/Source/CesiumRuntime/Private/CesiumVoxelRendererComponent.h
index 411797df1..1dfa57828 100644
--- a/Source/CesiumRuntime/Private/CesiumVoxelRendererComponent.h
+++ b/Source/CesiumRuntime/Private/CesiumVoxelRendererComponent.h
@@ -9,8 +9,8 @@
 #include "CustomDepthParameters.h"
 #include "Materials/MaterialInstanceDynamic.h"
 #include "Templates/UniquePtr.h"
-#include "VoxelMegatextures.h"
 #include "VoxelGridShape.h"
+#include "VoxelMegatextures.h"
 #include "VoxelOctree.h"
 
 #include <Cesium3DTilesSelection/Tile.h>
@@ -101,7 +101,8 @@ class UCesiumVoxelRendererComponent : public USceneComponent {
       const FCesiumVoxelClassDescription* pDescription,
       const Cesium3DTilesSelection::BoundingVolume& boundingVolume);
 
-  static double computePriority(double sse);
+  static double
+  computePriority(const CesiumGeometry::OctreeTileID& tileId, double sse);
 
   struct VoxelTileUpdateInfo {
     const UCesiumGltfVoxelComponent* pComponent;
@@ -126,4 +127,9 @@ class UCesiumVoxelRendererComponent : public USceneComponent {
   std::vector<CesiumGeometry::OctreeTileID> _loadedNodeIds;
   MaxPriorityQueue _visibleTileQueue;
   bool _needsOctreeUpdate;
+
+  /**
+   * The tileset that owns this voxel renderer.
+   */
+  ACesium3DTileset* _pTileset = nullptr;
 };

From 1e0ca35b8bae49b883c0529953b50d9ee0e505e6 Mon Sep 17 00:00:00 2001
From: Janine Liu <janine@cesium.com>
Date: Mon, 4 Aug 2025 14:35:24 -0400
Subject: [PATCH 4/8] Update comments for correctness

---
 Shaders/Private/CesiumEllipsoidRegion.usf      |  8 ++++----
 .../Private/CesiumVoxelRendererComponent.cpp   | 18 +++++++++---------
 2 files changed, 13 insertions(+), 13 deletions(-)

diff --git a/Shaders/Private/CesiumEllipsoidRegion.usf b/Shaders/Private/CesiumEllipsoidRegion.usf
index 56b08b6c1..9520e9ae4 100644
--- a/Shaders/Private/CesiumEllipsoidRegion.usf
+++ b/Shaders/Private/CesiumEllipsoidRegion.usf
@@ -26,8 +26,8 @@ struct EllipsoidRegion
   float LatitudeUVOffset;
 
   uint LongitudeRangeFlag;
-  bool LongitudeMinHasDiscontinuity;
-  bool LongitudeMaxHasDiscontinuity;
+  bool LongitudeMinAtDiscontinuity;
+  bool LongitudeMaxAtDiscontinuity;
   bool LongitudeIsReversed;
   uint LatitudeMinFlag;
   uint LatitudeMaxFlag;
@@ -39,8 +39,8 @@ struct EllipsoidRegion
 
     // Flags are packed in CesiumVoxelRendererComponent.cpp
     LongitudeRangeFlag = round(PackedData0.x);
-    LongitudeMinHasDiscontinuity = bool(PackedData0.y);
-    LongitudeMaxHasDiscontinuity = bool(PackedData0.z);
+    LongitudeMinAtDiscontinuity = bool(PackedData0.y);
+    LongitudeMaxAtDiscontinuity = bool(PackedData0.z);
     LongitudeIsReversed = bool(PackedData0.w);
 
     LatitudeMinFlag = round(PackedData1.x);
diff --git a/Source/CesiumRuntime/Private/CesiumVoxelRendererComponent.cpp b/Source/CesiumRuntime/Private/CesiumVoxelRendererComponent.cpp
index f14f6c482..10cea421c 100644
--- a/Source/CesiumRuntime/Private/CesiumVoxelRendererComponent.cpp
+++ b/Source/CesiumRuntime/Private/CesiumVoxelRendererComponent.cpp
@@ -237,8 +237,8 @@ void setVoxelEllipsoidProperties(
   double maximumHeight = region.getMaximumHeight();
 
   // Defines the extents of the longitude in UV space. In other words, this
-  // expresses the minimum, maximum, and midpoint values of the longitude range
-  // in UV space.
+  // expresses the minimum and maximum values of the longitude range, as well as
+  // the midpoint of the negative space.
   FVector longitudeUVExtents = FVector::Zero();
   double longitudeUVScale = 1;
   double longitudeUVOffset = 0;
@@ -255,11 +255,11 @@ void setVoxelEllipsoidProperties(
     // Refers to the discontinuity at longitude 0 / 2pi.
     const double discontinuityEpsilon =
         CesiumUtility::Math::Epsilon3; // 0.001 radians = 0.05729578 degrees
-    bool longitudeMinimumHasDiscontinuity = CesiumUtility::Math::equalsEpsilon(
+    bool longitudeMinimumAtDiscontinuity = CesiumUtility::Math::equalsEpsilon(
         minimumLongitude,
         0,
         discontinuityEpsilon);
-    bool longitudeMaximumHasDiscontinuity = CesiumUtility::Math::equalsEpsilon(
+    bool longitudeMaximumAtDiscontinuity = CesiumUtility::Math::equalsEpsilon(
         maximumLongitude,
         CesiumUtility::Math::TwoPi,
         discontinuityEpsilon);
@@ -268,8 +268,8 @@ void setVoxelEllipsoidProperties(
         interpretLongitudeRange(longitudeRange);
 
     longitudeFlags.X = longitudeRangeIndicator;
-    longitudeFlags.Y = longitudeMinimumHasDiscontinuity;
-    longitudeFlags.Z = longitudeMaximumHasDiscontinuity;
+    longitudeFlags.Y = longitudeMinimumAtDiscontinuity;
+    longitudeFlags.Z = longitudeMaximumAtDiscontinuity;
     longitudeFlags.W = isLongitudeReversed;
 
     // Compute the extents of the longitude range in UV Shape Space.
@@ -277,10 +277,10 @@ void setVoxelEllipsoidProperties(
         (minimumLongitude - defaultMinimumBounds.X) / defaultRange;
     double maximumLongitudeUV =
         (maximumLongitude - defaultMinimumBounds.X) / defaultRange;
-    // Given a longitude range, represents the actual value where "0" would be
-    // in UV coordinates.
+    // Given the longitude range, describes the proportion of the ellipsoid
+    // that is excluded from that range.
     double longitudeRangeUVZero = 1.0 - longitudeRange / defaultRange;
-    // TODO: document this
+    // Describes the midpoint of the above excluded range.
     double longitudeRangeUVZeroMid =
         glm::fract(maximumLongitudeUV + 0.5 * longitudeRangeUVZero);
 

From 359fda73f040bff4a2710cf5c12ac66a09bafbd4 Mon Sep 17 00:00:00 2001
From: Janine Liu <janine@cesium.com>
Date: Mon, 4 Aug 2025 15:01:26 -0400
Subject: [PATCH 5/8] Fix old names

---
 Shaders/Private/CesiumEllipsoidRegion.usf | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Shaders/Private/CesiumEllipsoidRegion.usf b/Shaders/Private/CesiumEllipsoidRegion.usf
index 9520e9ae4..3b1bc2ea0 100644
--- a/Shaders/Private/CesiumEllipsoidRegion.usf
+++ b/Shaders/Private/CesiumEllipsoidRegion.usf
@@ -805,12 +805,12 @@ struct EllipsoidRegion
     }
 
     // Avoid flickering from reading voxels from both sides of the -pi/+pi discontinuity.
-    if (LongitudeMinHasDiscontinuity)
+    if (LongitudeMinAtDiscontinuity)
     {
       longitudeUV = longitudeUV > LongitudeUVExtents.z ? LongitudeUVExtents.x : longitudeUV;
     }
 
-    if (LongitudeMaxHasDiscontinuity)
+    if (LongitudeMaxAtDiscontinuity)
     {
       longitudeUV = longitudeUV < LongitudeUVExtents.z ? LongitudeUVExtents.y : longitudeUV;
     }

From 7bb09efa5f1526fe480a2f2ff70a073d67ebe988 Mon Sep 17 00:00:00 2001
From: Janine Liu <janine@cesium.com>
Date: Tue, 5 Aug 2025 14:52:14 -0400
Subject: [PATCH 6/8] Fix intersection list error

---
 Shaders/Private/CesiumRayIntersectionList.usf | 7 ++-----
 1 file changed, 2 insertions(+), 5 deletions(-)

diff --git a/Shaders/Private/CesiumRayIntersectionList.usf b/Shaders/Private/CesiumRayIntersectionList.usf
index 318484600..4bc8d6fa2 100644
--- a/Shaders/Private/CesiumRayIntersectionList.usf
+++ b/Shaders/Private/CesiumRayIntersectionList.usf
@@ -107,10 +107,7 @@ struct IntersectionListState
   */
   RayIntersections GetNextIntersections(in float4 Data[INTERSECTIONS_LENGTH])
   {
-    RayIntersections result = (RayIntersections) 0;
-    result.Entry.t = NO_HIT;
-    result.Exit.t = NO_HIT;
-
+    RayIntersections result = NewRayIntersections(NewMissedIntersection(), NewMissedIntersection());
     if (Index >= Length)
     {
       return result;
@@ -139,7 +136,7 @@ struct IntersectionListState
       IsInsidePositiveShape = isCurrentShapePositive ? isEnteringShape : IsInsidePositiveShape;
 
       // True if entering positive shape or exiting negative shape
-      if (IsInsidePositiveShape && isEnteringShape == isCurrentShapePositive)
+      if (SurroundingShapeCount == 1 && IsInsidePositiveShape && isEnteringShape == isCurrentShapePositive)
       {
         result.Entry.t = surfaceIntersection.w;
         result.Entry.Normal = normalize(surfaceIntersection.xyz);

From e1ef2106d5c7280ad3b992cf7f65345ac18d66f4 Mon Sep 17 00:00:00 2001
From: Janine Liu <janine@cesium.com>
Date: Tue, 19 Aug 2025 14:34:48 -0400
Subject: [PATCH 7/8] Fix region discontinuity values

---
 .../CesiumRuntime/Private/CesiumVoxelRendererComponent.cpp  | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/Source/CesiumRuntime/Private/CesiumVoxelRendererComponent.cpp b/Source/CesiumRuntime/Private/CesiumVoxelRendererComponent.cpp
index 10cea421c..02a6294f2 100644
--- a/Source/CesiumRuntime/Private/CesiumVoxelRendererComponent.cpp
+++ b/Source/CesiumRuntime/Private/CesiumVoxelRendererComponent.cpp
@@ -252,16 +252,16 @@ void setVoxelEllipsoidProperties(
     double longitudeRange = maximumLongitude - minimumLongitude +
                             isLongitudeReversed * defaultRange;
 
-    // Refers to the discontinuity at longitude 0 / 2pi.
+    // Refers to the discontinuity at longitude -pi / pi.
     const double discontinuityEpsilon =
         CesiumUtility::Math::Epsilon3; // 0.001 radians = 0.05729578 degrees
     bool longitudeMinimumAtDiscontinuity = CesiumUtility::Math::equalsEpsilon(
         minimumLongitude,
-        0,
+        -defaultMinimumBounds.X,
         discontinuityEpsilon);
     bool longitudeMaximumAtDiscontinuity = CesiumUtility::Math::equalsEpsilon(
         maximumLongitude,
-        CesiumUtility::Math::TwoPi,
+        defaultMaximumBounds.X,
         discontinuityEpsilon);
 
     CartographicAngleDescription longitudeRangeIndicator =

From eae4464a254585790c1a71b6425634434a9e6ce0 Mon Sep 17 00:00:00 2001
From: Janine Liu <janine@cesium.com>
Date: Tue, 19 Aug 2025 15:17:50 -0400
Subject: [PATCH 8/8] Avoid flickering artifacts

---
 Shaders/Private/CesiumRayIntersectionList.usf | 10 ++++++++--
 1 file changed, 8 insertions(+), 2 deletions(-)

diff --git a/Shaders/Private/CesiumRayIntersectionList.usf b/Shaders/Private/CesiumRayIntersectionList.usf
index 4bc8d6fa2..3954f02bc 100644
--- a/Shaders/Private/CesiumRayIntersectionList.usf
+++ b/Shaders/Private/CesiumRayIntersectionList.usf
@@ -123,9 +123,15 @@ struct IntersectionListState
       {
         continue;
       }
-
+      
+      if (i >= Length)
+      {
+        Index = INTERSECTIONS_LENGTH;
+        break;
+      }
+      
       Index = i + 1;
-
+      
       surfaceIntersection = Data[i];
       // Maps from [1-4] -> [0-3] (see EncodeIntersection for the types)
       int intersectionType = int(length(surfaceIntersection.xyz) - 0.5);