diff --git a/Content/Materials/Layers/ML_CesiumVoxel.uasset b/Content/Materials/Layers/ML_CesiumVoxel.uasset
index 811ed7395..28e011b07 100644
Binary files a/Content/Materials/Layers/ML_CesiumVoxel.uasset and b/Content/Materials/Layers/ML_CesiumVoxel.uasset differ
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..3b1bc2ea0
--- /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 LongitudeMinAtDiscontinuity;
+  bool LongitudeMaxAtDiscontinuity;
+  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);
+    LongitudeMinAtDiscontinuity = bool(PackedData0.y);
+    LongitudeMaxAtDiscontinuity = 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 (LongitudeMinAtDiscontinuity)
+    {
+      longitudeUV = longitudeUV > LongitudeUVExtents.z ? LongitudeUVExtents.x : longitudeUV;
+    }
+
+    if (LongitudeMaxAtDiscontinuity)
+    {
+      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..3954f02bc
--- /dev/null
+++ b/Shaders/Private/CesiumRayIntersectionList.usf
@@ -0,0 +1,185 @@
+#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 = NewRayIntersections(NewMissedIntersection(), NewMissedIntersection());
+    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;
+      }
+      
+      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);
+      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 (SurroundingShapeCount == 1 && 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 4af26ff3d..6a413d419 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 2a82091fe..240a3d2a4 100644
--- a/Source/CesiumRuntime/Private/CesiumVoxelRendererComponent.cpp
+++ b/Source/CesiumRuntime/Private/CesiumVoxelRendererComponent.cpp
@@ -83,6 +83,9 @@ EVoxelGridShape getVoxelGridShape(
   if (std::get_if<CesiumGeometry::OrientedBoundingBox>(&boundingVolume)) {
     return EVoxelGridShape::Box;
   }
+  if (std::get_if<CesiumGeospatial::BoundingRegion>(&boundingVolume)) {
+    return EVoxelGridShape::Ellipsoid;
+  }
 
   return EVoxelGridShape::Invalid;
 }
@@ -126,6 +129,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 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;
+
+  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 -pi / pi.
+    const double discontinuityEpsilon =
+        CesiumUtility::Math::Epsilon3; // 0.001 radians = 0.05729578 degrees
+    bool longitudeMinimumAtDiscontinuity = CesiumUtility::Math::equalsEpsilon(
+        minimumLongitude,
+        -defaultMinimumBounds.X,
+        discontinuityEpsilon);
+    bool longitudeMaximumAtDiscontinuity = CesiumUtility::Math::equalsEpsilon(
+        maximumLongitude,
+        defaultMaximumBounds.X,
+        discontinuityEpsilon);
+
+    CartographicAngleDescription longitudeRangeIndicator =
+        interpretLongitudeRange(longitudeRange);
+
+    longitudeFlags.X = longitudeRangeIndicator;
+    longitudeFlags.Y = longitudeMinimumAtDiscontinuity;
+    longitudeFlags.Z = longitudeMaximumAtDiscontinuity;
+    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 the longitude range, describes the proportion of the ellipsoid
+    // that is excluded from that range.
+    double longitudeRangeUVZero = 1.0 - longitudeRange / defaultRange;
+    // Describes the midpoint of the above excluded range.
+    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)
@@ -197,6 +484,7 @@ UCesiumVoxelRendererComponent::CreateVoxelMaterial(
       FName("VoxelMaterial"));
   pVoxelMaterial->SetFlags(
       RF_Transient | RF_DuplicateTransient | RF_TextExportTransient);
+  pVoxelComponent->_pTileset = pTilesetActor;
 
   EVoxelGridShape shape = pVoxelComponent->Options.gridShape;
 
@@ -237,6 +525,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) {
@@ -562,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()) {
@@ -584,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();
@@ -681,6 +978,60 @@ void UCesiumVoxelRendererComponent::UpdateTiles(
   }
 }
 
+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(
@@ -706,8 +1057,24 @@ 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);
+  }
 }
 
-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 34d044bb8..02e2d5ebf 100644
--- a/Source/CesiumRuntime/Private/CesiumVoxelRendererComponent.h
+++ b/Source/CesiumRuntime/Private/CesiumVoxelRendererComponent.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;
 };