Add existing to tracked

node_modules/@cesium/engine/Source/Shaders/Voxels/IntersectBox.glsl

@@ -0,0 +1,34 @@
+// See IntersectionUtils.glsl for the definitions of Ray, RayShapeIntersection,
+// NO_HIT, Intersections
+
+/* Box defines (set in Scene/VoxelBoxShape.js)
+#define BOX_INTERSECTION_INDEX ### // always 0
+*/
+
+uniform sampler2D u_renderBoundPlanesTexture;
+
+RayShapeIntersection intersectBoundPlanes(in Ray ray) {
+    vec4 lastEntry = vec4(ray.dir, -INF_HIT);
+    vec4 firstExit = vec4(-ray.dir, +INF_HIT);
+    for (int i = 0; i < 6; i++) {
+        vec4 boundPlane = getBoundPlane(u_renderBoundPlanesTexture, i);
+        vec4 intersection = intersectPlane(ray, boundPlane);
+        if (dot(ray.dir, boundPlane.xyz) < 0.0) {
+            lastEntry = intersection.w > lastEntry.w ? intersection : lastEntry;
+        } else {
+            firstExit = intersection.w < firstExit.w ? intersection: firstExit;
+        }
+    }
+
+    if (lastEntry.w < firstExit.w) {
+        return RayShapeIntersection(lastEntry, firstExit);
+    } else {
+        return RayShapeIntersection(vec4(-ray.dir, NO_HIT), vec4(ray.dir, NO_HIT));
+    }
+}
+
+void intersectShape(in Ray rayUV, in  Ray rayEC, inout Intersections ix)
+{
+    RayShapeIntersection intersection = intersectBoundPlanes(rayEC);
+    setShapeIntersection(ix, BOX_INTERSECTION_INDEX, intersection);
+}

node_modules/@cesium/engine/Source/Shaders/Voxels/IntersectBox.js

@@ -0,0 +1,36 @@
+//This file is automatically rebuilt by the Cesium build process.
+export default "// See IntersectionUtils.glsl for the definitions of Ray, RayShapeIntersection,\n\
+// NO_HIT, Intersections\n\
+\n\
+/* Box defines (set in Scene/VoxelBoxShape.js)\n\
+#define BOX_INTERSECTION_INDEX ### // always 0\n\
+*/\n\
+\n\
+uniform sampler2D u_renderBoundPlanesTexture;\n\
+\n\
+RayShapeIntersection intersectBoundPlanes(in Ray ray) {\n\
+    vec4 lastEntry = vec4(ray.dir, -INF_HIT);\n\
+    vec4 firstExit = vec4(-ray.dir, +INF_HIT);\n\
+    for (int i = 0; i < 6; i++) {\n\
+        vec4 boundPlane = getBoundPlane(u_renderBoundPlanesTexture, i);\n\
+        vec4 intersection = intersectPlane(ray, boundPlane);\n\
+        if (dot(ray.dir, boundPlane.xyz) < 0.0) {\n\
+            lastEntry = intersection.w > lastEntry.w ? intersection : lastEntry;\n\
+        } else {\n\
+            firstExit = intersection.w < firstExit.w ? intersection: firstExit;\n\
+        }\n\
+    }\n\
+\n\
+    if (lastEntry.w < firstExit.w) {\n\
+        return RayShapeIntersection(lastEntry, firstExit);\n\
+    } else {\n\
+        return RayShapeIntersection(vec4(-ray.dir, NO_HIT), vec4(ray.dir, NO_HIT));\n\
+    }\n\
+}\n\
+\n\
+void intersectShape(in Ray rayUV, in  Ray rayEC, inout Intersections ix)\n\
+{\n\
+    RayShapeIntersection intersection = intersectBoundPlanes(rayEC);\n\
+    setShapeIntersection(ix, BOX_INTERSECTION_INDEX, intersection);\n\
+}\n\
+";

node_modules/@cesium/engine/Source/Shaders/Voxels/IntersectCylinder.glsl

@@ -0,0 +1,142 @@
+// See IntersectionUtils.glsl for the definitions of Ray, NO_HIT, Intersections,
+// RayShapeIntersection, setSurfaceIntersection, setShapeIntersection,
+// intersectIntersections
+// See IntersectLongitude.glsl for the definitions of intersectHalfPlane,
+// intersectFlippedWedge, intersectRegularWedge
+
+/* Cylinder defines (set in Scene/VoxelCylinderShape.js)
+#define CYLINDER_HAS_RENDER_BOUNDS_RADIUS_MIN
+#define CYLINDER_HAS_RENDER_BOUNDS_RADIUS_FLAT
+#define CYLINDER_HAS_RENDER_BOUNDS_ANGLE
+#define CYLINDER_HAS_RENDER_BOUNDS_ANGLE_RANGE_UNDER_HALF
+#define CYLINDER_HAS_RENDER_BOUNDS_ANGLE_RANGE_OVER_HALF
+#define CYLINDER_HAS_RENDER_BOUNDS_ANGLE_RANGE_EQUAL_ZERO
+
+#define CYLINDER_INTERSECTION_INDEX_RADIUS_MAX
+#define CYLINDER_INTERSECTION_INDEX_RADIUS_MIN
+#define CYLINDER_INTERSECTION_INDEX_ANGLE
+*/
+
+// Cylinder uniforms
+uniform vec2 u_cylinderRenderRadiusMinMax;
+#if defined(CYLINDER_HAS_RENDER_BOUNDS_ANGLE)
+    uniform vec2 u_cylinderRenderAngleMinMax;
+#endif
+
+uniform sampler2D u_renderBoundPlanesTexture;
+
+RayShapeIntersection intersectBoundPlanes(in Ray ray) {
+    vec4 lastEntry = vec4(ray.dir, -INF_HIT);
+    vec4 firstExit = vec4(-ray.dir, +INF_HIT);
+    for (int i = 0; i < 2; i++) {
+        vec4 boundPlane = getBoundPlane(u_renderBoundPlanesTexture, i);
+        vec4 intersection = intersectPlane(ray, boundPlane);
+        if (dot(ray.dir, boundPlane.xyz) < 0.0) {
+            lastEntry = intersection.w > lastEntry.w ? intersection : lastEntry;
+        } else {
+            firstExit = intersection.w < firstExit.w ? intersection: firstExit;
+        }
+    }
+
+    if (lastEntry.w < firstExit.w) {
+        return RayShapeIntersection(lastEntry, firstExit);
+    } else {
+        return RayShapeIntersection(vec4(-ray.dir, NO_HIT), vec4(ray.dir, NO_HIT));
+    }
+}
+
+/**
+ * Find the intersection of a ray with a right cylindrical surface of a given radius
+ * about the z-axis.
+ */
+RayShapeIntersection intersectCylinder(in Ray ray, in float radius, in bool convex)
+{
+    vec2 position = ray.pos.xy;
+    vec2 direction = ray.dir.xy;
+
+    float a = dot(direction, direction);
+    float b = dot(position, direction);
+    float c = dot(position, position) - radius * radius;
+    float determinant = b * b - a * c;
+
+    if (determinant < 0.0) {
+        vec4 miss = vec4(normalize(ray.dir), NO_HIT);
+        return RayShapeIntersection(miss, miss);
+    }
+
+    determinant = sqrt(determinant);
+    float t1 = (-b - determinant) / a;
+    float t2 = (-b + determinant) / a;
+    float signFlip = convex ? 1.0 : -1.0;
+    vec3 normal1 = vec3((position + t1 * direction) * signFlip, 0.0);
+    vec3 normal2 = vec3((position + t2 * direction) * signFlip, 0.0);
+    // Return normals in eye coordinates
+    vec4 intersect1 = vec4(normalize(czm_normal * normal1), t1);
+    vec4 intersect2 = vec4(normalize(czm_normal * normal2), t2);
+
+    return RayShapeIntersection(intersect1, intersect2);
+}
+
+/**
+ * Find the intersection of a ray with a right cylindrical solid of given
+ * radius and height bounds. NOTE: The shape is assumed to be convex.
+ */
+RayShapeIntersection intersectBoundedCylinder(in Ray ray, in Ray rayEC, in float radius)
+{
+    RayShapeIntersection cylinderIntersection = intersectCylinder(ray, radius, true);
+    RayShapeIntersection heightBoundsIntersection = intersectBoundPlanes(rayEC);
+    return intersectIntersections(ray, cylinderIntersection, heightBoundsIntersection);
+}
+
+void intersectShape(in Ray ray, in Ray rayEC, inout Intersections ix)
+{
+    RayShapeIntersection outerIntersect = intersectBoundedCylinder(ray, rayEC, u_cylinderRenderRadiusMinMax.y);
+
+    setShapeIntersection(ix, CYLINDER_INTERSECTION_INDEX_RADIUS_MAX, outerIntersect);
+
+    if (outerIntersect.entry.w == NO_HIT) {
+        return;
+    }
+
+    #if defined(CYLINDER_HAS_RENDER_BOUNDS_RADIUS_FLAT)
+        // When the cylinder is perfectly thin it's necessary to sandwich the
+        // inner cylinder intersection inside the outer cylinder intersection.
+
+        // Without this special case,
+        // [outerMin, outerMax, innerMin, innerMax] will bubble sort to
+        // [outerMin, innerMin, outerMax, innerMax] which will cause the back
+        // side of the cylinder to be invisible because it will think the ray
+        // is still inside the inner (negative) cylinder after exiting the
+        // outer (positive) cylinder.
+
+        // With this special case,
+        // [outerMin, innerMin, innerMax, outerMax] will bubble sort to
+        // [outerMin, innerMin, innerMax, outerMax] which will work correctly.
+
+        // Note: If initializeIntersections() changes its sorting function
+        // from bubble sort to something else, this code may need to change.
+        RayShapeIntersection innerIntersect = intersectCylinder(ray, 1.0, false);
+        setSurfaceIntersection(ix, 0, outerIntersect.entry, true, true);  // positive, enter
+        setSurfaceIntersection(ix, 1, innerIntersect.entry, false, true); // negative, enter
+        setSurfaceIntersection(ix, 2, innerIntersect.exit, false, false); // negative, exit
+        setSurfaceIntersection(ix, 3, outerIntersect.exit, true, false);  // positive, exit
+    #elif defined(CYLINDER_HAS_RENDER_BOUNDS_RADIUS_MIN)
+        RayShapeIntersection innerIntersect = intersectCylinder(ray, u_cylinderRenderRadiusMinMax.x, false);
+        setShapeIntersection(ix, CYLINDER_INTERSECTION_INDEX_RADIUS_MIN, innerIntersect);
+    #endif
+
+    #if defined(CYLINDER_HAS_RENDER_BOUNDS_ANGLE_RANGE_UNDER_HALF)
+        RayShapeIntersection wedgeIntersect = intersectRegularWedge(ray, u_cylinderRenderAngleMinMax);
+        setShapeIntersection(ix, CYLINDER_INTERSECTION_INDEX_ANGLE, wedgeIntersect);
+    #elif defined(CYLINDER_HAS_RENDER_BOUNDS_ANGLE_RANGE_OVER_HALF)
+        RayShapeIntersection wedgeIntersects[2];
+        intersectFlippedWedge(ray, u_cylinderRenderAngleMinMax, wedgeIntersects);
+        setShapeIntersection(ix, CYLINDER_INTERSECTION_INDEX_ANGLE + 0, wedgeIntersects[0]);
+        setShapeIntersection(ix, CYLINDER_INTERSECTION_INDEX_ANGLE + 1, wedgeIntersects[1]);
+    #elif defined(CYLINDER_HAS_RENDER_BOUNDS_ANGLE_RANGE_EQUAL_ZERO)
+        RayShapeIntersection wedgeIntersects[2];
+        intersectHalfPlane(ray, u_cylinderRenderAngleMinMax.x, wedgeIntersects);
+        setShapeIntersection(ix, CYLINDER_INTERSECTION_INDEX_ANGLE + 0, wedgeIntersects[0]);
+        setShapeIntersection(ix, CYLINDER_INTERSECTION_INDEX_ANGLE + 1, wedgeIntersects[1]);
+    #endif
+}

node_modules/@cesium/engine/Source/Shaders/Voxels/IntersectCylinder.js

@@ -0,0 +1,144 @@
+//This file is automatically rebuilt by the Cesium build process.
+export default "// See IntersectionUtils.glsl for the definitions of Ray, NO_HIT, Intersections,\n\
+// RayShapeIntersection, setSurfaceIntersection, setShapeIntersection,\n\
+// intersectIntersections\n\
+// See IntersectLongitude.glsl for the definitions of intersectHalfPlane,\n\
+// intersectFlippedWedge, intersectRegularWedge\n\
+\n\
+/* Cylinder defines (set in Scene/VoxelCylinderShape.js)\n\
+#define CYLINDER_HAS_RENDER_BOUNDS_RADIUS_MIN\n\
+#define CYLINDER_HAS_RENDER_BOUNDS_RADIUS_FLAT\n\
+#define CYLINDER_HAS_RENDER_BOUNDS_ANGLE\n\
+#define CYLINDER_HAS_RENDER_BOUNDS_ANGLE_RANGE_UNDER_HALF\n\
+#define CYLINDER_HAS_RENDER_BOUNDS_ANGLE_RANGE_OVER_HALF\n\
+#define CYLINDER_HAS_RENDER_BOUNDS_ANGLE_RANGE_EQUAL_ZERO\n\
+\n\
+#define CYLINDER_INTERSECTION_INDEX_RADIUS_MAX\n\
+#define CYLINDER_INTERSECTION_INDEX_RADIUS_MIN\n\
+#define CYLINDER_INTERSECTION_INDEX_ANGLE\n\
+*/\n\
+\n\
+// Cylinder uniforms\n\
+uniform vec2 u_cylinderRenderRadiusMinMax;\n\
+#if defined(CYLINDER_HAS_RENDER_BOUNDS_ANGLE)\n\
+    uniform vec2 u_cylinderRenderAngleMinMax;\n\
+#endif\n\
+\n\
+uniform sampler2D u_renderBoundPlanesTexture;\n\
+\n\
+RayShapeIntersection intersectBoundPlanes(in Ray ray) {\n\
+    vec4 lastEntry = vec4(ray.dir, -INF_HIT);\n\
+    vec4 firstExit = vec4(-ray.dir, +INF_HIT);\n\
+    for (int i = 0; i < 2; i++) {\n\
+        vec4 boundPlane = getBoundPlane(u_renderBoundPlanesTexture, i);\n\
+        vec4 intersection = intersectPlane(ray, boundPlane);\n\
+        if (dot(ray.dir, boundPlane.xyz) < 0.0) {\n\
+            lastEntry = intersection.w > lastEntry.w ? intersection : lastEntry;\n\
+        } else {\n\
+            firstExit = intersection.w < firstExit.w ? intersection: firstExit;\n\
+        }\n\
+    }\n\
+\n\
+    if (lastEntry.w < firstExit.w) {\n\
+        return RayShapeIntersection(lastEntry, firstExit);\n\
+    } else {\n\
+        return RayShapeIntersection(vec4(-ray.dir, NO_HIT), vec4(ray.dir, NO_HIT));\n\
+    }\n\
+}\n\
+\n\
+/**\n\
+ * Find the intersection of a ray with a right cylindrical surface of a given radius\n\
+ * about the z-axis.\n\
+ */\n\
+RayShapeIntersection intersectCylinder(in Ray ray, in float radius, in bool convex)\n\
+{\n\
+    vec2 position = ray.pos.xy;\n\
+    vec2 direction = ray.dir.xy;\n\
+\n\
+    float a = dot(direction, direction);\n\
+    float b = dot(position, direction);\n\
+    float c = dot(position, position) - radius * radius;\n\
+    float determinant = b * b - a * c;\n\
+\n\
+    if (determinant < 0.0) {\n\
+        vec4 miss = vec4(normalize(ray.dir), NO_HIT);\n\
+        return RayShapeIntersection(miss, miss);\n\
+    }\n\
+\n\
+    determinant = sqrt(determinant);\n\
+    float t1 = (-b - determinant) / a;\n\
+    float t2 = (-b + determinant) / a;\n\
+    float signFlip = convex ? 1.0 : -1.0;\n\
+    vec3 normal1 = vec3((position + t1 * direction) * signFlip, 0.0);\n\
+    vec3 normal2 = vec3((position + t2 * direction) * signFlip, 0.0);\n\
+    // Return normals in eye coordinates\n\
+    vec4 intersect1 = vec4(normalize(czm_normal * normal1), t1);\n\
+    vec4 intersect2 = vec4(normalize(czm_normal * normal2), t2);\n\
+\n\
+    return RayShapeIntersection(intersect1, intersect2);\n\
+}\n\
+\n\
+/**\n\
+ * Find the intersection of a ray with a right cylindrical solid of given\n\
+ * radius and height bounds. NOTE: The shape is assumed to be convex.\n\
+ */\n\
+RayShapeIntersection intersectBoundedCylinder(in Ray ray, in Ray rayEC, in float radius)\n\
+{\n\
+    RayShapeIntersection cylinderIntersection = intersectCylinder(ray, radius, true);\n\
+    RayShapeIntersection heightBoundsIntersection = intersectBoundPlanes(rayEC);\n\
+    return intersectIntersections(ray, cylinderIntersection, heightBoundsIntersection);\n\
+}\n\
+\n\
+void intersectShape(in Ray ray, in Ray rayEC, inout Intersections ix)\n\
+{\n\
+    RayShapeIntersection outerIntersect = intersectBoundedCylinder(ray, rayEC, u_cylinderRenderRadiusMinMax.y);\n\
+\n\
+    setShapeIntersection(ix, CYLINDER_INTERSECTION_INDEX_RADIUS_MAX, outerIntersect);\n\
+\n\
+    if (outerIntersect.entry.w == NO_HIT) {\n\
+        return;\n\
+    }\n\
+\n\
+    #if defined(CYLINDER_HAS_RENDER_BOUNDS_RADIUS_FLAT)\n\
+        // When the cylinder is perfectly thin it's necessary to sandwich the\n\
+        // inner cylinder intersection inside the outer cylinder intersection.\n\
+\n\
+        // Without this special case,\n\
+        // [outerMin, outerMax, innerMin, innerMax] will bubble sort to\n\
+        // [outerMin, innerMin, outerMax, innerMax] which will cause the back\n\
+        // side of the cylinder to be invisible because it will think the ray\n\
+        // is still inside the inner (negative) cylinder after exiting the\n\
+        // outer (positive) cylinder.\n\
+\n\
+        // With this special case,\n\
+        // [outerMin, innerMin, innerMax, outerMax] will bubble sort to\n\
+        // [outerMin, innerMin, innerMax, outerMax] which will work correctly.\n\
+\n\
+        // Note: If initializeIntersections() changes its sorting function\n\
+        // from bubble sort to something else, this code may need to change.\n\
+        RayShapeIntersection innerIntersect = intersectCylinder(ray, 1.0, false);\n\
+        setSurfaceIntersection(ix, 0, outerIntersect.entry, true, true);  // positive, enter\n\
+        setSurfaceIntersection(ix, 1, innerIntersect.entry, false, true); // negative, enter\n\
+        setSurfaceIntersection(ix, 2, innerIntersect.exit, false, false); // negative, exit\n\
+        setSurfaceIntersection(ix, 3, outerIntersect.exit, true, false);  // positive, exit\n\
+    #elif defined(CYLINDER_HAS_RENDER_BOUNDS_RADIUS_MIN)\n\
+        RayShapeIntersection innerIntersect = intersectCylinder(ray, u_cylinderRenderRadiusMinMax.x, false);\n\
+        setShapeIntersection(ix, CYLINDER_INTERSECTION_INDEX_RADIUS_MIN, innerIntersect);\n\
+    #endif\n\
+\n\
+    #if defined(CYLINDER_HAS_RENDER_BOUNDS_ANGLE_RANGE_UNDER_HALF)\n\
+        RayShapeIntersection wedgeIntersect = intersectRegularWedge(ray, u_cylinderRenderAngleMinMax);\n\
+        setShapeIntersection(ix, CYLINDER_INTERSECTION_INDEX_ANGLE, wedgeIntersect);\n\
+    #elif defined(CYLINDER_HAS_RENDER_BOUNDS_ANGLE_RANGE_OVER_HALF)\n\
+        RayShapeIntersection wedgeIntersects[2];\n\
+        intersectFlippedWedge(ray, u_cylinderRenderAngleMinMax, wedgeIntersects);\n\
+        setShapeIntersection(ix, CYLINDER_INTERSECTION_INDEX_ANGLE + 0, wedgeIntersects[0]);\n\
+        setShapeIntersection(ix, CYLINDER_INTERSECTION_INDEX_ANGLE + 1, wedgeIntersects[1]);\n\
+    #elif defined(CYLINDER_HAS_RENDER_BOUNDS_ANGLE_RANGE_EQUAL_ZERO)\n\
+        RayShapeIntersection wedgeIntersects[2];\n\
+        intersectHalfPlane(ray, u_cylinderRenderAngleMinMax.x, wedgeIntersects);\n\
+        setShapeIntersection(ix, CYLINDER_INTERSECTION_INDEX_ANGLE + 0, wedgeIntersects[0]);\n\
+        setShapeIntersection(ix, CYLINDER_INTERSECTION_INDEX_ANGLE + 1, wedgeIntersects[1]);\n\
+    #endif\n\
+}\n\
+";

node_modules/@cesium/engine/Source/Shaders/Voxels/IntersectDepth.glsl

@@ -0,0 +1,27 @@
+// See IntersectionUtils.glsl for the definitions of Ray, Intersections,
+// setIntersectionPair, INF_HIT, NO_HIT
+
+/* intersectDepth defines (set in Scene/VoxelRenderResources.js)
+#define DEPTH_INTERSECTION_INDEX ###
+*/
+
+void intersectDepth(in vec2 screenCoord, in Ray ray, inout Intersections ix) {
+    float logDepthOrDepth = czm_unpackDepth(texture(czm_globeDepthTexture, screenCoord));
+    float entry;
+    float exit;
+    if (logDepthOrDepth != 0.0) {
+        // Calculate how far the ray must travel before it hits the depth buffer.
+        vec4 eyeCoordinateDepth = czm_screenToEyeCoordinates(screenCoord, logDepthOrDepth);
+        eyeCoordinateDepth /= eyeCoordinateDepth.w;
+        entry = dot(eyeCoordinateDepth.xyz - ray.pos, ray.dir);
+        exit = +INF_HIT;
+    } else {
+        // There's no depth at this location.
+        entry = NO_HIT;
+        exit = NO_HIT;
+    }
+    ix.distanceToDepthBuffer = entry;
+#if defined(DEPTH_TEST)
+    setIntersectionPair(ix, DEPTH_INTERSECTION_INDEX, vec2(entry, exit));
+#endif
+}

node_modules/@cesium/engine/Source/Shaders/Voxels/IntersectDepth.js

@@ -0,0 +1,29 @@
+//This file is automatically rebuilt by the Cesium build process.
+export default "// See IntersectionUtils.glsl for the definitions of Ray, Intersections,\n\
+// setIntersectionPair, INF_HIT, NO_HIT\n\
+\n\
+/* intersectDepth defines (set in Scene/VoxelRenderResources.js)\n\
+#define DEPTH_INTERSECTION_INDEX ###\n\
+*/\n\
+\n\
+void intersectDepth(in vec2 screenCoord, in Ray ray, inout Intersections ix) {\n\
+    float logDepthOrDepth = czm_unpackDepth(texture(czm_globeDepthTexture, screenCoord));\n\
+    float entry;\n\
+    float exit;\n\
+    if (logDepthOrDepth != 0.0) {\n\
+        // Calculate how far the ray must travel before it hits the depth buffer.\n\
+        vec4 eyeCoordinateDepth = czm_screenToEyeCoordinates(screenCoord, logDepthOrDepth);\n\
+        eyeCoordinateDepth /= eyeCoordinateDepth.w;\n\
+        entry = dot(eyeCoordinateDepth.xyz - ray.pos, ray.dir);\n\
+        exit = +INF_HIT;\n\
+    } else {\n\
+        // There's no depth at this location.\n\
+        entry = NO_HIT;\n\
+        exit = NO_HIT;\n\
+    }\n\
+    ix.distanceToDepthBuffer = entry;\n\
+#if defined(DEPTH_TEST)\n\
+    setIntersectionPair(ix, DEPTH_INTERSECTION_INDEX, vec2(entry, exit));\n\
+#endif\n\
+}\n\
+";

node_modules/@cesium/engine/Source/Shaders/Voxels/IntersectEllipsoid.glsl

@@ -0,0 +1,330 @@
+// See IntersectionUtils.glsl for the definitions of Ray, NO_HIT, INF_HIT, Intersections,
+// RayShapeIntersection, setSurfaceIntersection, setShapeIntersection
+// See IntersectLongitude.glsl for the definitions of intersectHalfPlane,
+// intersectFlippedWedge, intersectRegularWedge
+
+/* Ellipsoid defines (set in Scene/VoxelEllipsoidShape.js)
+#define ELLIPSOID_HAS_RENDER_BOUNDS_LONGITUDE
+#define ELLIPSOID_HAS_RENDER_BOUNDS_LONGITUDE_RANGE_EQUAL_ZERO
+#define ELLIPSOID_HAS_RENDER_BOUNDS_LONGITUDE_RANGE_UNDER_HALF
+#define ELLIPSOID_HAS_RENDER_BOUNDS_LONGITUDE_RANGE_OVER_HALF
+#define ELLIPSOID_HAS_RENDER_BOUNDS_LATITUDE_MAX_UNDER_HALF
+#define ELLIPSOID_HAS_RENDER_BOUNDS_LATITUDE_MAX_EQUAL_HALF
+#define ELLIPSOID_HAS_RENDER_BOUNDS_LATITUDE_MAX_OVER_HALF
+#define ELLIPSOID_HAS_RENDER_BOUNDS_LATITUDE_MIN_UNDER_HALF
+#define ELLIPSOID_HAS_RENDER_BOUNDS_LATITUDE_MIN_EQUAL_HALF
+#define ELLIPSOID_HAS_RENDER_BOUNDS_LATITUDE_MIN_OVER_HALF
+#define ELLIPSOID_INTERSECTION_INDEX_LONGITUDE
+#define ELLIPSOID_INTERSECTION_INDEX_LATITUDE_MAX
+#define ELLIPSOID_INTERSECTION_INDEX_LATITUDE_MIN
+#define ELLIPSOID_INTERSECTION_INDEX_HEIGHT_MAX
+#define ELLIPSOID_INTERSECTION_INDEX_HEIGHT_MIN
+*/
+
+#if defined(ELLIPSOID_HAS_RENDER_BOUNDS_LONGITUDE)
+    uniform vec2 u_ellipsoidRenderLongitudeMinMax;
+#endif
+uniform float u_eccentricitySquared;
+uniform vec2 u_ellipsoidRenderLatitudeSinMinMax;
+uniform vec2 u_clipMinMaxHeight; // Values are negative: clipHeight - maxShapeHeight
+
+RayShapeIntersection intersectZPlane(in Ray ray, in float z) {
+    float t = -ray.pos.z / ray.dir.z;
+
+    bool startsOutside = sign(ray.pos.z) == sign(z);
+    bool entry = (t >= 0.0) != startsOutside;
+
+    vec4 intersect = vec4(0.0, 0.0, z, t);
+    vec4 farSide = vec4(normalize(ray.dir), INF_HIT);
+
+    if (entry) {
+        return RayShapeIntersection(intersect, farSide);
+    } else {
+        return RayShapeIntersection(-1.0 * farSide, intersect);
+    }
+}
+
+RayShapeIntersection intersectHeight(in Ray ray, in float height, in bool convex)
+{
+    // Scale the ray by the ellipsoid axes to make it a unit sphere
+    // Note: approximating ellipsoid + height as an ellipsoid
+    vec3 radiiCorrection = vec3(1.0) / (u_ellipsoidRadii + height);
+    vec3 position = ray.pos * radiiCorrection;
+    vec3 direction = ray.dir * radiiCorrection;
+
+    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) {
+        vec4 miss = vec4(normalize(direction), NO_HIT);
+        return RayShapeIntersection(miss, miss);
+    }
+
+    determinant = sqrt(determinant);
+
+    // Compute larger root using standard formula
+    float signB = b < 0.0 ? -1.0 : 1.0;
+    // The other root may suffer from subtractive cancellation in the standard formula.
+    // Compute it from the first root instead.
+    float t1 = (-b - signB * determinant) / a;
+    float t2 = c / (a * t1);
+    float tmin = min(t1, t2);
+    float tmax = max(t1, t2);
+
+    float directionScale = convex ? 1.0 : -1.0;
+    vec3 d1 = directionScale * (position + tmin * direction);
+    vec3 d2 = directionScale * (position + tmax * direction);
+
+    // Return normals in eye coordinates. Use spherical approximation for the normal.
+    vec3 normal1 = normalize(czm_normal * d1);
+    vec3 normal2 = normalize(czm_normal * d2);
+
+    return RayShapeIntersection(vec4(normal1, tmin), vec4(normal2, tmax));
+}
+
+/**
+ * Given a circular cone around the z-axis, with apex at the origin,
+ * find the parametric distance(s) along a ray where that ray intersects
+ * the cone.
+ * The cone opening angle is described by the squared cosine of
+ * its half-angle (the angle between the Z-axis and the surface)
+ */
+vec2 intersectDoubleEndedCone(in Ray ray, in float cosSqrHalfAngle)
+{
+    vec3 o = ray.pos;
+    vec3 d = ray.dir;
+    float sinSqrHalfAngle = 1.0 - cosSqrHalfAngle;
+
+    float aSin = d.z * d.z * sinSqrHalfAngle;
+    float aCos = -dot(d.xy, d.xy) * cosSqrHalfAngle;
+    float a = aSin + aCos;
+
+    float bSin = d.z * o.z * sinSqrHalfAngle;
+    float bCos = -dot(o.xy, d.xy) * cosSqrHalfAngle;
+    float b = bSin + bCos;
+
+    float cSin = o.z * o.z * sinSqrHalfAngle;
+    float cCos = -dot(o.xy, o.xy) * cosSqrHalfAngle;
+    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;
+
+    if (determinant < 0.0) {
+        return vec2(NO_HIT);
+    } else if (a == 0.0) {
+        // Ray is parallel to cone surface
+        return (b == 0.0)
+            ? vec2(NO_HIT) // Ray is on cone surface
+            : vec2(-0.5 * c / b, NO_HIT);
+    }
+
+    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);
+    float tmin = min(t1, t2);
+    float tmax = max(t1, t2);
+    return vec2(tmin, tmax);
+}
+
+/**
+ * Given a point on a conical surface, find the surface normal at that point.
+ */
+vec3 getConeNormal(in vec3 p, in bool convex) {
+    // Start with radial component pointing toward z-axis
+    vec2 radial = -abs(p.z) * normalize(p.xy);
+    // Z component points toward opening of cone
+    float zSign = (p.z < 0.0) ? -1.0 : 1.0;
+    float z = length(p.xy) * zSign;
+    // Flip normal if shape is convex
+    float flip = (convex) ? -1.0 : 1.0;
+    return normalize(vec3(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 = u_eccentricitySquared * sinLatitude * sinLatitude;
+    float primeVerticalRadius = u_ellipsoidRadii.x * inversesqrt(1.0 - x2);
+
+    // Compute a shift from the origin to the intersection of the cone with the z-axis
+    return primeVerticalRadius * u_eccentricitySquared * sinLatitude;
+}
+
+void intersectFlippedCone(in Ray ray, in float cosHalfAngle, out RayShapeIntersection intersections[2]) {
+    // Shift the ray to account for the latitude cone not being centered at the Earth center
+    ray.pos.z += getLatitudeConeShift(cosHalfAngle);
+
+    float cosSqrHalfAngle = cosHalfAngle * cosHalfAngle;
+    vec2 intersect = intersectDoubleEndedCone(ray, cosSqrHalfAngle);
+
+    vec4 miss = vec4(normalize(ray.dir), NO_HIT);
+    vec4 farSide = vec4(normalize(ray.dir), INF_HIT);
+
+    // Initialize output with no intersections
+    intersections[0].entry = -1.0 * farSide;
+    intersections[0].exit = farSide;
+    intersections[1].entry = miss;
+    intersections[1].exit = miss;
+
+    if (intersect.x == NO_HIT) {
+        return;
+    }
+
+    // Find the points of intersection
+    float tmin = intersect.x;
+    float tmax = intersect.y;
+    vec3 p0 = ray.pos + tmin * ray.dir;
+    vec3 p1 = ray.pos + tmax * ray.dir;
+
+    vec4 intersect0 = vec4(getConeNormal(p0, true), tmin);
+    vec4 intersect1 = vec4(getConeNormal(p1, true), tmax);
+
+    bool p0InShadowCone = sign(p0.z) != sign(cosHalfAngle);
+    bool p1InShadowCone = sign(p1.z) != sign(cosHalfAngle);
+
+    if (p0InShadowCone && p1InShadowCone) {
+        // no valid intersections
+    } else if (p0InShadowCone) {
+        intersections[0].exit = intersect1;
+    } else if (p1InShadowCone) {
+        intersections[0].entry = intersect0;
+    } else {
+        intersections[0].exit = intersect0;
+        intersections[1].entry = intersect1;
+        intersections[1].exit = farSide;
+    }
+}
+
+RayShapeIntersection intersectRegularCone(in Ray ray, in float cosHalfAngle, in bool convex) {
+    // Shift the ray to account for the latitude cone not being centered at the Earth center
+    ray.pos.z += getLatitudeConeShift(cosHalfAngle);
+
+    float cosSqrHalfAngle = cosHalfAngle * cosHalfAngle;
+    vec2 intersect = intersectDoubleEndedCone(ray, cosSqrHalfAngle);
+
+    vec4 miss = vec4(normalize(ray.dir), NO_HIT);
+    vec4 farSide = vec4(normalize(ray.dir), INF_HIT);
+
+    if (intersect.x == NO_HIT) {
+        return RayShapeIntersection(miss, miss);
+    }
+
+    // Find the points of intersection
+    float tmin = intersect.x;
+    float tmax = intersect.y;
+    vec3 p0 = ray.pos + tmin * ray.dir;
+    vec3 p1 = ray.pos + tmax * ray.dir;
+
+    vec4 intersect0 = vec4(getConeNormal(p0, convex), tmin);
+    vec4 intersect1 = vec4(getConeNormal(p1, convex), tmax);
+
+    bool p0InShadowCone = sign(p0.z) != sign(cosHalfAngle);
+    bool p1InShadowCone = sign(p1.z) != sign(cosHalfAngle);
+
+    if (p0InShadowCone && p1InShadowCone) {
+        return RayShapeIntersection(miss, miss);
+    } else if (p0InShadowCone) {
+        return RayShapeIntersection(intersect1, farSide);
+    } else if (p1InShadowCone) {
+        return RayShapeIntersection(-1.0 * farSide, intersect0);
+    } else {
+        return RayShapeIntersection(intersect0, intersect1);
+    }
+}
+
+void intersectShape(in Ray ray, in Ray rayEC, inout Intersections ix) {    // Outer ellipsoid
+    RayShapeIntersection outerIntersect = intersectHeight(ray, u_clipMinMaxHeight.y, true);
+    setShapeIntersection(ix, ELLIPSOID_INTERSECTION_INDEX_HEIGHT_MAX, outerIntersect);
+
+    // Exit early if the outer ellipsoid was missed.
+    if (outerIntersect.entry.w == NO_HIT) {
+        return;
+    }
+
+    // Inner ellipsoid
+    RayShapeIntersection innerIntersect = intersectHeight(ray, u_clipMinMaxHeight.x, false);
+
+    if (innerIntersect.entry.w == NO_HIT) {
+        setShapeIntersection(ix, ELLIPSOID_INTERSECTION_INDEX_HEIGHT_MIN, innerIntersect);
+    } 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 innerIntersect to outerIntersect 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 initializeIntersections() changes its sorting function
+        // from bubble sort to something else, this code may need to change.
+        innerIntersect.entry.w = max(innerIntersect.entry.w, outerIntersect.entry.w);
+        innerIntersect.exit.w = min(innerIntersect.exit.w, outerIntersect.exit.w);
+        setSurfaceIntersection(ix, 0, outerIntersect.entry, true, true);  // positive, enter
+        setSurfaceIntersection(ix, 1, innerIntersect.entry, false, true); // negative, enter
+        setSurfaceIntersection(ix, 2, innerIntersect.exit, false, false); // negative, exit
+        setSurfaceIntersection(ix, 3, outerIntersect.exit, true, false);  // positive, exit
+    }
+
+    // Bottom cone
+    #if defined(ELLIPSOID_HAS_RENDER_BOUNDS_LATITUDE_MIN_UNDER_HALF)
+        RayShapeIntersection bottomConeIntersection = intersectRegularCone(ray, u_ellipsoidRenderLatitudeSinMinMax.x, false);
+        setShapeIntersection(ix, ELLIPSOID_INTERSECTION_INDEX_LATITUDE_MIN, bottomConeIntersection);
+    #elif defined(ELLIPSOID_HAS_RENDER_BOUNDS_LATITUDE_MIN_EQUAL_HALF)
+        RayShapeIntersection bottomConeIntersection = intersectZPlane(ray, -1.0);
+        setShapeIntersection(ix, ELLIPSOID_INTERSECTION_INDEX_LATITUDE_MIN, bottomConeIntersection);
+    #elif defined(ELLIPSOID_HAS_RENDER_BOUNDS_LATITUDE_MIN_OVER_HALF)
+        RayShapeIntersection bottomConeIntersections[2];
+        intersectFlippedCone(ray, u_ellipsoidRenderLatitudeSinMinMax.x, bottomConeIntersections);
+        setShapeIntersection(ix, ELLIPSOID_INTERSECTION_INDEX_LATITUDE_MIN + 0, bottomConeIntersections[0]);
+        setShapeIntersection(ix, ELLIPSOID_INTERSECTION_INDEX_LATITUDE_MIN + 1, bottomConeIntersections[1]);
+    #endif
+
+    // Top cone
+    #if defined(ELLIPSOID_HAS_RENDER_BOUNDS_LATITUDE_MAX_UNDER_HALF)
+        RayShapeIntersection topConeIntersections[2];
+        intersectFlippedCone(ray, u_ellipsoidRenderLatitudeSinMinMax.y, topConeIntersections);
+        setShapeIntersection(ix, ELLIPSOID_INTERSECTION_INDEX_LATITUDE_MAX + 0, topConeIntersections[0]);
+        setShapeIntersection(ix, ELLIPSOID_INTERSECTION_INDEX_LATITUDE_MAX + 1, topConeIntersections[1]);
+    #elif defined(ELLIPSOID_HAS_RENDER_BOUNDS_LATITUDE_MAX_EQUAL_HALF)
+        RayShapeIntersection topConeIntersection = intersectZPlane(ray, 1.0);
+        setShapeIntersection(ix, ELLIPSOID_INTERSECTION_INDEX_LATITUDE_MAX, topConeIntersection);
+    #elif defined(ELLIPSOID_HAS_RENDER_BOUNDS_LATITUDE_MAX_OVER_HALF)
+        RayShapeIntersection topConeIntersection = intersectRegularCone(ray, u_ellipsoidRenderLatitudeSinMinMax.y, false);
+        setShapeIntersection(ix, ELLIPSOID_INTERSECTION_INDEX_LATITUDE_MAX, topConeIntersection);
+    #endif
+
+    // Wedge
+    #if defined(ELLIPSOID_HAS_RENDER_BOUNDS_LONGITUDE_RANGE_EQUAL_ZERO)
+        RayShapeIntersection wedgeIntersects[2];
+        intersectHalfPlane(ray, u_ellipsoidRenderLongitudeMinMax.x, wedgeIntersects);
+        setShapeIntersection(ix, ELLIPSOID_INTERSECTION_INDEX_LONGITUDE + 0, wedgeIntersects[0]);
+        setShapeIntersection(ix, ELLIPSOID_INTERSECTION_INDEX_LONGITUDE + 1, wedgeIntersects[1]);
+    #elif defined(ELLIPSOID_HAS_RENDER_BOUNDS_LONGITUDE_RANGE_UNDER_HALF)
+        RayShapeIntersection wedgeIntersect = intersectRegularWedge(ray, u_ellipsoidRenderLongitudeMinMax);
+        setShapeIntersection(ix, ELLIPSOID_INTERSECTION_INDEX_LONGITUDE, wedgeIntersect);
+    #elif defined(ELLIPSOID_HAS_RENDER_BOUNDS_LONGITUDE_RANGE_OVER_HALF)
+        RayShapeIntersection wedgeIntersects[2];
+        intersectFlippedWedge(ray, u_ellipsoidRenderLongitudeMinMax, wedgeIntersects);
+        setShapeIntersection(ix, ELLIPSOID_INTERSECTION_INDEX_LONGITUDE + 0, wedgeIntersects[0]);
+        setShapeIntersection(ix, ELLIPSOID_INTERSECTION_INDEX_LONGITUDE + 1, wedgeIntersects[1]);
+    #endif
+}

node_modules/@cesium/engine/Source/Shaders/Voxels/IntersectEllipsoid.js

@@ -0,0 +1,332 @@
+//This file is automatically rebuilt by the Cesium build process.
+export default "// See IntersectionUtils.glsl for the definitions of Ray, NO_HIT, INF_HIT, Intersections,\n\
+// RayShapeIntersection, setSurfaceIntersection, setShapeIntersection\n\
+// See IntersectLongitude.glsl for the definitions of intersectHalfPlane,\n\
+// intersectFlippedWedge, intersectRegularWedge\n\
+\n\
+/* Ellipsoid defines (set in Scene/VoxelEllipsoidShape.js)\n\
+#define ELLIPSOID_HAS_RENDER_BOUNDS_LONGITUDE\n\
+#define ELLIPSOID_HAS_RENDER_BOUNDS_LONGITUDE_RANGE_EQUAL_ZERO\n\
+#define ELLIPSOID_HAS_RENDER_BOUNDS_LONGITUDE_RANGE_UNDER_HALF\n\
+#define ELLIPSOID_HAS_RENDER_BOUNDS_LONGITUDE_RANGE_OVER_HALF\n\
+#define ELLIPSOID_HAS_RENDER_BOUNDS_LATITUDE_MAX_UNDER_HALF\n\
+#define ELLIPSOID_HAS_RENDER_BOUNDS_LATITUDE_MAX_EQUAL_HALF\n\
+#define ELLIPSOID_HAS_RENDER_BOUNDS_LATITUDE_MAX_OVER_HALF\n\
+#define ELLIPSOID_HAS_RENDER_BOUNDS_LATITUDE_MIN_UNDER_HALF\n\
+#define ELLIPSOID_HAS_RENDER_BOUNDS_LATITUDE_MIN_EQUAL_HALF\n\
+#define ELLIPSOID_HAS_RENDER_BOUNDS_LATITUDE_MIN_OVER_HALF\n\
+#define ELLIPSOID_INTERSECTION_INDEX_LONGITUDE\n\
+#define ELLIPSOID_INTERSECTION_INDEX_LATITUDE_MAX\n\
+#define ELLIPSOID_INTERSECTION_INDEX_LATITUDE_MIN\n\
+#define ELLIPSOID_INTERSECTION_INDEX_HEIGHT_MAX\n\
+#define ELLIPSOID_INTERSECTION_INDEX_HEIGHT_MIN\n\
+*/\n\
+\n\
+#if defined(ELLIPSOID_HAS_RENDER_BOUNDS_LONGITUDE)\n\
+    uniform vec2 u_ellipsoidRenderLongitudeMinMax;\n\
+#endif\n\
+uniform float u_eccentricitySquared;\n\
+uniform vec2 u_ellipsoidRenderLatitudeSinMinMax;\n\
+uniform vec2 u_clipMinMaxHeight; // Values are negative: clipHeight - maxShapeHeight\n\
+\n\
+RayShapeIntersection intersectZPlane(in Ray ray, in float z) {\n\
+    float t = -ray.pos.z / ray.dir.z;\n\
+\n\
+    bool startsOutside = sign(ray.pos.z) == sign(z);\n\
+    bool entry = (t >= 0.0) != startsOutside;\n\
+\n\
+    vec4 intersect = vec4(0.0, 0.0, z, t);\n\
+    vec4 farSide = vec4(normalize(ray.dir), INF_HIT);\n\
+\n\
+    if (entry) {\n\
+        return RayShapeIntersection(intersect, farSide);\n\
+    } else {\n\
+        return RayShapeIntersection(-1.0 * farSide, intersect);\n\
+    }\n\
+}\n\
+\n\
+RayShapeIntersection intersectHeight(in Ray ray, in float height, in bool convex)\n\
+{\n\
+    // Scale the ray by the ellipsoid axes to make it a unit sphere\n\
+    // Note: approximating ellipsoid + height as an ellipsoid\n\
+    vec3 radiiCorrection = vec3(1.0) / (u_ellipsoidRadii + height);\n\
+    vec3 position = ray.pos * radiiCorrection;\n\
+    vec3 direction = ray.dir * radiiCorrection;\n\
+\n\
+    float a = dot(direction, direction); // ~ 1.0 (or maybe 4.0 if ray is scaled)\n\
+    float b = dot(direction, position); // roughly inside [-1.0, 1.0] when zoomed in\n\
+    float c = dot(position, position) - 1.0; // ~ 0.0 when zoomed in.\n\
+    float determinant = b * b - a * c; // ~ b * b when zoomed in\n\
+\n\
+    if (determinant < 0.0) {\n\
+        vec4 miss = vec4(normalize(direction), NO_HIT);\n\
+        return RayShapeIntersection(miss, miss);\n\
+    }\n\
+\n\
+    determinant = sqrt(determinant);\n\
+\n\
+    // Compute larger root using standard formula\n\
+    float signB = b < 0.0 ? -1.0 : 1.0;\n\
+    // The other root may suffer from subtractive cancellation in the standard formula.\n\
+    // Compute it from the first root instead.\n\
+    float t1 = (-b - signB * determinant) / a;\n\
+    float t2 = c / (a * t1);\n\
+    float tmin = min(t1, t2);\n\
+    float tmax = max(t1, t2);\n\
+\n\
+    float directionScale = convex ? 1.0 : -1.0;\n\
+    vec3 d1 = directionScale * (position + tmin * direction);\n\
+    vec3 d2 = directionScale * (position + tmax * direction);\n\
+\n\
+    // Return normals in eye coordinates. Use spherical approximation for the normal.\n\
+    vec3 normal1 = normalize(czm_normal * d1);\n\
+    vec3 normal2 = normalize(czm_normal * d2);\n\
+\n\
+    return RayShapeIntersection(vec4(normal1, tmin), vec4(normal2, tmax));\n\
+}\n\
+\n\
+/**\n\
+ * Given a circular cone around the z-axis, with apex at the origin,\n\
+ * find the parametric distance(s) along a ray where that ray intersects\n\
+ * the cone.\n\
+ * The cone opening angle is described by the squared cosine of\n\
+ * its half-angle (the angle between the Z-axis and the surface)\n\
+ */\n\
+vec2 intersectDoubleEndedCone(in Ray ray, in float cosSqrHalfAngle)\n\
+{\n\
+    vec3 o = ray.pos;\n\
+    vec3 d = ray.dir;\n\
+    float sinSqrHalfAngle = 1.0 - cosSqrHalfAngle;\n\
+\n\
+    float aSin = d.z * d.z * sinSqrHalfAngle;\n\
+    float aCos = -dot(d.xy, d.xy) * cosSqrHalfAngle;\n\
+    float a = aSin + aCos;\n\
+\n\
+    float bSin = d.z * o.z * sinSqrHalfAngle;\n\
+    float bCos = -dot(o.xy, d.xy) * cosSqrHalfAngle;\n\
+    float b = bSin + bCos;\n\
+\n\
+    float cSin = o.z * o.z * sinSqrHalfAngle;\n\
+    float cCos = -dot(o.xy, o.xy) * cosSqrHalfAngle;\n\
+    float c = cSin + cCos;\n\
+    // determinant = b * b - a * c. But bSin * bSin = aSin * cSin.\n\
+    // Avoid subtractive cancellation by expanding to eliminate these terms\n\
+    float determinant = 2.0 * bSin * bCos + bCos * bCos - aSin * cCos - aCos * cSin - aCos * cCos;\n\
+\n\
+    if (determinant < 0.0) {\n\
+        return vec2(NO_HIT);\n\
+    } else if (a == 0.0) {\n\
+        // Ray is parallel to cone surface\n\
+        return (b == 0.0)\n\
+            ? vec2(NO_HIT) // Ray is on cone surface\n\
+            : vec2(-0.5 * c / b, NO_HIT);\n\
+    }\n\
+\n\
+    determinant = sqrt(determinant);\n\
+\n\
+    // Compute larger root using standard formula\n\
+    float signB = b < 0.0 ? -1.0 : 1.0;\n\
+    float t1 = (-b - signB * determinant) / a;\n\
+    // The other root may suffer from subtractive cancellation in the standard formula.\n\
+    // Compute it from the first root instead.\n\
+    float t2 = c / (a * t1);\n\
+    float tmin = min(t1, t2);\n\
+    float tmax = max(t1, t2);\n\
+    return vec2(tmin, tmax);\n\
+}\n\
+\n\
+/**\n\
+ * Given a point on a conical surface, find the surface normal at that point.\n\
+ */\n\
+vec3 getConeNormal(in vec3 p, in bool convex) {\n\
+    // Start with radial component pointing toward z-axis\n\
+    vec2 radial = -abs(p.z) * normalize(p.xy);\n\
+    // Z component points toward opening of cone\n\
+    float zSign = (p.z < 0.0) ? -1.0 : 1.0;\n\
+    float z = length(p.xy) * zSign;\n\
+    // Flip normal if shape is convex\n\
+    float flip = (convex) ? -1.0 : 1.0;\n\
+    return normalize(vec3(radial, z) * flip);\n\
+}\n\
+\n\
+/**\n\
+ * Compute the shift between the ellipsoid origin and the apex of a cone of latitude\n\
+ */\n\
+float getLatitudeConeShift(in float sinLatitude) {\n\
+    // Find prime vertical radius of curvature: \n\
+    // the distance along the ellipsoid normal to the intersection with the z-axis\n\
+    float x2 = u_eccentricitySquared * sinLatitude * sinLatitude;\n\
+    float primeVerticalRadius = u_ellipsoidRadii.x * inversesqrt(1.0 - x2);\n\
+\n\
+    // Compute a shift from the origin to the intersection of the cone with the z-axis\n\
+    return primeVerticalRadius * u_eccentricitySquared * sinLatitude;\n\
+}\n\
+\n\
+void intersectFlippedCone(in Ray ray, in float cosHalfAngle, out RayShapeIntersection intersections[2]) {\n\
+    // Shift the ray to account for the latitude cone not being centered at the Earth center\n\
+    ray.pos.z += getLatitudeConeShift(cosHalfAngle);\n\
+\n\
+    float cosSqrHalfAngle = cosHalfAngle * cosHalfAngle;\n\
+    vec2 intersect = intersectDoubleEndedCone(ray, cosSqrHalfAngle);\n\
+\n\
+    vec4 miss = vec4(normalize(ray.dir), NO_HIT);\n\
+    vec4 farSide = vec4(normalize(ray.dir), INF_HIT);\n\
+\n\
+    // Initialize output with no intersections\n\
+    intersections[0].entry = -1.0 * farSide;\n\
+    intersections[0].exit = farSide;\n\
+    intersections[1].entry = miss;\n\
+    intersections[1].exit = miss;\n\
+\n\
+    if (intersect.x == NO_HIT) {\n\
+        return;\n\
+    }\n\
+\n\
+    // Find the points of intersection\n\
+    float tmin = intersect.x;\n\
+    float tmax = intersect.y;\n\
+    vec3 p0 = ray.pos + tmin * ray.dir;\n\
+    vec3 p1 = ray.pos + tmax * ray.dir;\n\
+\n\
+    vec4 intersect0 = vec4(getConeNormal(p0, true), tmin);\n\
+    vec4 intersect1 = vec4(getConeNormal(p1, true), tmax);\n\
+\n\
+    bool p0InShadowCone = sign(p0.z) != sign(cosHalfAngle);\n\
+    bool p1InShadowCone = sign(p1.z) != sign(cosHalfAngle);\n\
+\n\
+    if (p0InShadowCone && p1InShadowCone) {\n\
+        // no valid intersections\n\
+    } else if (p0InShadowCone) {\n\
+        intersections[0].exit = intersect1;\n\
+    } else if (p1InShadowCone) {\n\
+        intersections[0].entry = intersect0;\n\
+    } else {\n\
+        intersections[0].exit = intersect0;\n\
+        intersections[1].entry = intersect1;\n\
+        intersections[1].exit = farSide;\n\
+    }\n\
+}\n\
+\n\
+RayShapeIntersection intersectRegularCone(in Ray ray, in float cosHalfAngle, in bool convex) {\n\
+    // Shift the ray to account for the latitude cone not being centered at the Earth center\n\
+    ray.pos.z += getLatitudeConeShift(cosHalfAngle);\n\
+\n\
+    float cosSqrHalfAngle = cosHalfAngle * cosHalfAngle;\n\
+    vec2 intersect = intersectDoubleEndedCone(ray, cosSqrHalfAngle);\n\
+\n\
+    vec4 miss = vec4(normalize(ray.dir), NO_HIT);\n\
+    vec4 farSide = vec4(normalize(ray.dir), INF_HIT);\n\
+\n\
+    if (intersect.x == NO_HIT) {\n\
+        return RayShapeIntersection(miss, miss);\n\
+    }\n\
+\n\
+    // Find the points of intersection\n\
+    float tmin = intersect.x;\n\
+    float tmax = intersect.y;\n\
+    vec3 p0 = ray.pos + tmin * ray.dir;\n\
+    vec3 p1 = ray.pos + tmax * ray.dir;\n\
+\n\
+    vec4 intersect0 = vec4(getConeNormal(p0, convex), tmin);\n\
+    vec4 intersect1 = vec4(getConeNormal(p1, convex), tmax);\n\
+\n\
+    bool p0InShadowCone = sign(p0.z) != sign(cosHalfAngle);\n\
+    bool p1InShadowCone = sign(p1.z) != sign(cosHalfAngle);\n\
+\n\
+    if (p0InShadowCone && p1InShadowCone) {\n\
+        return RayShapeIntersection(miss, miss);\n\
+    } else if (p0InShadowCone) {\n\
+        return RayShapeIntersection(intersect1, farSide);\n\
+    } else if (p1InShadowCone) {\n\
+        return RayShapeIntersection(-1.0 * farSide, intersect0);\n\
+    } else {\n\
+        return RayShapeIntersection(intersect0, intersect1);\n\
+    }\n\
+}\n\
+\n\
+void intersectShape(in Ray ray, in Ray rayEC, inout Intersections ix) {    // Outer ellipsoid\n\
+    RayShapeIntersection outerIntersect = intersectHeight(ray, u_clipMinMaxHeight.y, true);\n\
+    setShapeIntersection(ix, ELLIPSOID_INTERSECTION_INDEX_HEIGHT_MAX, outerIntersect);\n\
+\n\
+    // Exit early if the outer ellipsoid was missed.\n\
+    if (outerIntersect.entry.w == NO_HIT) {\n\
+        return;\n\
+    }\n\
+\n\
+    // Inner ellipsoid\n\
+    RayShapeIntersection innerIntersect = intersectHeight(ray, u_clipMinMaxHeight.x, false);\n\
+\n\
+    if (innerIntersect.entry.w == NO_HIT) {\n\
+        setShapeIntersection(ix, ELLIPSOID_INTERSECTION_INDEX_HEIGHT_MIN, innerIntersect);\n\
+    } else {\n\
+        // When the ellipsoid is large and thin it's possible for floating point math\n\
+        // to cause the ray to intersect the inner ellipsoid before the outer ellipsoid. \n\
+        // To prevent this from happening, clamp innerIntersect to outerIntersect and\n\
+        // sandwich the inner ellipsoid intersection inside the outer ellipsoid intersection.\n\
+\n\
+        // Without this special case,\n\
+        // [outerMin, outerMax, innerMin, innerMax] will bubble sort to\n\
+        // [outerMin, innerMin, outerMax, innerMax] which will cause the back\n\
+        // side of the ellipsoid to be invisible because it will think the ray\n\
+        // is still inside the inner (negative) ellipsoid after exiting the\n\
+        // outer (positive) ellipsoid.\n\
+\n\
+        // With this special case,\n\
+        // [outerMin, innerMin, innerMax, outerMax] will bubble sort to\n\
+        // [outerMin, innerMin, innerMax, outerMax] which will work correctly.\n\
+\n\
+        // Note: If initializeIntersections() changes its sorting function\n\
+        // from bubble sort to something else, this code may need to change.\n\
+        innerIntersect.entry.w = max(innerIntersect.entry.w, outerIntersect.entry.w);\n\
+        innerIntersect.exit.w = min(innerIntersect.exit.w, outerIntersect.exit.w);\n\
+        setSurfaceIntersection(ix, 0, outerIntersect.entry, true, true);  // positive, enter\n\
+        setSurfaceIntersection(ix, 1, innerIntersect.entry, false, true); // negative, enter\n\
+        setSurfaceIntersection(ix, 2, innerIntersect.exit, false, false); // negative, exit\n\
+        setSurfaceIntersection(ix, 3, outerIntersect.exit, true, false);  // positive, exit\n\
+    }\n\
+\n\
+    // Bottom cone\n\
+    #if defined(ELLIPSOID_HAS_RENDER_BOUNDS_LATITUDE_MIN_UNDER_HALF)\n\
+        RayShapeIntersection bottomConeIntersection = intersectRegularCone(ray, u_ellipsoidRenderLatitudeSinMinMax.x, false);\n\
+        setShapeIntersection(ix, ELLIPSOID_INTERSECTION_INDEX_LATITUDE_MIN, bottomConeIntersection);\n\
+    #elif defined(ELLIPSOID_HAS_RENDER_BOUNDS_LATITUDE_MIN_EQUAL_HALF)\n\
+        RayShapeIntersection bottomConeIntersection = intersectZPlane(ray, -1.0);\n\
+        setShapeIntersection(ix, ELLIPSOID_INTERSECTION_INDEX_LATITUDE_MIN, bottomConeIntersection);\n\
+    #elif defined(ELLIPSOID_HAS_RENDER_BOUNDS_LATITUDE_MIN_OVER_HALF)\n\
+        RayShapeIntersection bottomConeIntersections[2];\n\
+        intersectFlippedCone(ray, u_ellipsoidRenderLatitudeSinMinMax.x, bottomConeIntersections);\n\
+        setShapeIntersection(ix, ELLIPSOID_INTERSECTION_INDEX_LATITUDE_MIN + 0, bottomConeIntersections[0]);\n\
+        setShapeIntersection(ix, ELLIPSOID_INTERSECTION_INDEX_LATITUDE_MIN + 1, bottomConeIntersections[1]);\n\
+    #endif\n\
+\n\
+    // Top cone\n\
+    #if defined(ELLIPSOID_HAS_RENDER_BOUNDS_LATITUDE_MAX_UNDER_HALF)\n\
+        RayShapeIntersection topConeIntersections[2];\n\
+        intersectFlippedCone(ray, u_ellipsoidRenderLatitudeSinMinMax.y, topConeIntersections);\n\
+        setShapeIntersection(ix, ELLIPSOID_INTERSECTION_INDEX_LATITUDE_MAX + 0, topConeIntersections[0]);\n\
+        setShapeIntersection(ix, ELLIPSOID_INTERSECTION_INDEX_LATITUDE_MAX + 1, topConeIntersections[1]);\n\
+    #elif defined(ELLIPSOID_HAS_RENDER_BOUNDS_LATITUDE_MAX_EQUAL_HALF)\n\
+        RayShapeIntersection topConeIntersection = intersectZPlane(ray, 1.0);\n\
+        setShapeIntersection(ix, ELLIPSOID_INTERSECTION_INDEX_LATITUDE_MAX, topConeIntersection);\n\
+    #elif defined(ELLIPSOID_HAS_RENDER_BOUNDS_LATITUDE_MAX_OVER_HALF)\n\
+        RayShapeIntersection topConeIntersection = intersectRegularCone(ray, u_ellipsoidRenderLatitudeSinMinMax.y, false);\n\
+        setShapeIntersection(ix, ELLIPSOID_INTERSECTION_INDEX_LATITUDE_MAX, topConeIntersection);\n\
+    #endif\n\
+\n\
+    // Wedge\n\
+    #if defined(ELLIPSOID_HAS_RENDER_BOUNDS_LONGITUDE_RANGE_EQUAL_ZERO)\n\
+        RayShapeIntersection wedgeIntersects[2];\n\
+        intersectHalfPlane(ray, u_ellipsoidRenderLongitudeMinMax.x, wedgeIntersects);\n\
+        setShapeIntersection(ix, ELLIPSOID_INTERSECTION_INDEX_LONGITUDE + 0, wedgeIntersects[0]);\n\
+        setShapeIntersection(ix, ELLIPSOID_INTERSECTION_INDEX_LONGITUDE + 1, wedgeIntersects[1]);\n\
+    #elif defined(ELLIPSOID_HAS_RENDER_BOUNDS_LONGITUDE_RANGE_UNDER_HALF)\n\
+        RayShapeIntersection wedgeIntersect = intersectRegularWedge(ray, u_ellipsoidRenderLongitudeMinMax);\n\
+        setShapeIntersection(ix, ELLIPSOID_INTERSECTION_INDEX_LONGITUDE, wedgeIntersect);\n\
+    #elif defined(ELLIPSOID_HAS_RENDER_BOUNDS_LONGITUDE_RANGE_OVER_HALF)\n\
+        RayShapeIntersection wedgeIntersects[2];\n\
+        intersectFlippedWedge(ray, u_ellipsoidRenderLongitudeMinMax, wedgeIntersects);\n\
+        setShapeIntersection(ix, ELLIPSOID_INTERSECTION_INDEX_LONGITUDE + 0, wedgeIntersects[0]);\n\
+        setShapeIntersection(ix, ELLIPSOID_INTERSECTION_INDEX_LONGITUDE + 1, wedgeIntersects[1]);\n\
+    #endif\n\
+}\n\
+";

node_modules/@cesium/engine/Source/Shaders/Voxels/IntersectLongitude.glsl

@@ -0,0 +1,114 @@
+// See IntersectionUtils.glsl for the definitions of Ray, NO_HIT, INF_HIT,
+// RayShapeIntersection
+
+vec4 transformNormalToEC(in vec4 intersection) {
+    return vec4(normalize(czm_normal * intersection.xyz), intersection.w);
+}
+
+RayShapeIntersection transformNormalsToEC(in RayShapeIntersection ix) {
+    return RayShapeIntersection(transformNormalToEC(ix.entry), transformNormalToEC(ix.exit));
+}
+
+vec4 intersectLongitude(in Ray ray, in float angle, in bool positiveNormal) {
+    float normalSign = positiveNormal ? 1.0 : -1.0;
+    vec2 planeNormal = vec2(-sin(angle), cos(angle)) * normalSign;
+
+    vec2 position = ray.pos.xy;
+    vec2 direction = ray.dir.xy;
+    float approachRate = dot(direction, planeNormal);
+    float distance = -dot(position, planeNormal);
+
+    float t = (approachRate == 0.0)
+        ? NO_HIT
+        : distance / approachRate;
+
+    return vec4(planeNormal, 0.0, t);
+}
+
+RayShapeIntersection intersectHalfSpace(in Ray ray, in float angle, in bool positiveNormal)
+{
+    vec4 intersection = intersectLongitude(ray, angle, positiveNormal);
+    vec4 farSide = vec4(normalize(ray.dir), INF_HIT);
+
+    bool hitFront = (intersection.w > 0.0) == (dot(ray.pos.xy, intersection.xy) > 0.0);
+    if (!hitFront) {
+        return RayShapeIntersection(intersection, farSide);
+    } else {
+        return RayShapeIntersection(-1.0 * farSide, intersection);
+    }
+}
+
+void intersectFlippedWedge(in Ray ray, in vec2 minMaxAngle, out RayShapeIntersection intersections[2])
+{
+    intersections[0] = transformNormalsToEC(intersectHalfSpace(ray, minMaxAngle.x, false));
+    intersections[1] = transformNormalsToEC(intersectHalfSpace(ray, minMaxAngle.y, true));
+}
+
+bool hitPositiveHalfPlane(in Ray ray, in vec4 intersection, in bool positiveNormal) {
+    float normalSign = positiveNormal ? 1.0 : -1.0;
+    vec2 planeDirection = vec2(intersection.y, -intersection.x) * normalSign;
+    vec2 hit = ray.pos.xy + intersection.w * ray.dir.xy;
+    return dot(hit, planeDirection) > 0.0;
+}
+
+void intersectHalfPlane(in Ray ray, in float angle, out RayShapeIntersection intersections[2]) {
+    vec4 intersection = intersectLongitude(ray, angle, true);
+    vec4 farSide = vec4(normalize(ray.dir), INF_HIT);
+    bool hitPositiveSide = hitPositiveHalfPlane(ray, intersection, true);
+
+    farSide = transformNormalToEC(farSide);
+
+    if (hitPositiveSide) {
+        intersection = transformNormalToEC(intersection);
+        intersections[0].entry = -1.0 * farSide;
+        intersections[0].exit = vec4(-1.0 * intersection.xyz, intersection.w);
+        intersections[1].entry = intersection;
+        intersections[1].exit = farSide;
+    } else {
+        vec4 miss = vec4(normalize(czm_normal * ray.dir), NO_HIT);
+        intersections[0].entry = -1.0 * farSide;
+        intersections[0].exit = farSide;
+        intersections[1].entry = miss;
+        intersections[1].exit = miss;
+    }
+}
+
+RayShapeIntersection intersectRegularWedge(in Ray ray, in vec2 minMaxAngle)
+{
+    // Note: works for maxAngle > minAngle + pi, where the "regular wedge"
+    // is actually a negative volume.
+    // Compute intersections with the two planes.
+    // Normals will point toward the "outside" (negative space)
+    vec4 intersect1 = intersectLongitude(ray, minMaxAngle.x, false);
+    vec4 intersect2 = intersectLongitude(ray, minMaxAngle.y, true);
+
+    // Choose intersection with smallest T as the "first", the other as "last"
+    // Note: first or last could be in the "shadow" wedge, beyond the tip
+    bool inOrder = intersect1.w <= intersect2.w;
+    vec4 first = inOrder ? intersect1 : intersect2;
+    vec4 last = inOrder ? intersect2 : intersect1;
+
+    bool firstIsAhead = first.w >= 0.0;
+    bool startedInsideFirst = dot(ray.pos.xy, first.xy) < 0.0;
+    bool exitFromInside = firstIsAhead == startedInsideFirst;
+    bool lastIsAhead = last.w > 0.0;
+    bool startedOutsideLast = dot(ray.pos.xy, last.xy) >= 0.0;
+    bool enterFromOutside = lastIsAhead == startedOutsideLast;
+
+    vec4 farSide = vec4(normalize(ray.dir), INF_HIT);
+    vec4 miss = vec4(normalize(ray.dir), NO_HIT);
+
+    if (exitFromInside && enterFromOutside) {
+        // Ray crosses both faces of negative wedge, exiting then entering the positive shape
+        return transformNormalsToEC(RayShapeIntersection(first, last));
+    } else if (!exitFromInside && enterFromOutside) {
+        // Ray starts inside wedge. last is in shadow wedge, and first is actually the entry
+        return transformNormalsToEC(RayShapeIntersection(-1.0 * farSide, first));
+    } else if (exitFromInside && !enterFromOutside) {
+        // First intersection was in the shadow wedge, so last is actually the exit
+        return transformNormalsToEC(RayShapeIntersection(last, farSide));
+    } else { // !exitFromInside && !enterFromOutside
+        // Both intersections were in the shadow wedge
+        return transformNormalsToEC(RayShapeIntersection(miss, miss));
+    }
+}

node_modules/@cesium/engine/Source/Shaders/Voxels/IntersectLongitude.js

@@ -0,0 +1,116 @@
+//This file is automatically rebuilt by the Cesium build process.
+export default "// See IntersectionUtils.glsl for the definitions of Ray, NO_HIT, INF_HIT,\n\
+// RayShapeIntersection\n\
+\n\
+vec4 transformNormalToEC(in vec4 intersection) {\n\
+    return vec4(normalize(czm_normal * intersection.xyz), intersection.w);\n\
+}\n\
+\n\
+RayShapeIntersection transformNormalsToEC(in RayShapeIntersection ix) {\n\
+    return RayShapeIntersection(transformNormalToEC(ix.entry), transformNormalToEC(ix.exit));\n\
+}\n\
+\n\
+vec4 intersectLongitude(in Ray ray, in float angle, in bool positiveNormal) {\n\
+    float normalSign = positiveNormal ? 1.0 : -1.0;\n\
+    vec2 planeNormal = vec2(-sin(angle), cos(angle)) * normalSign;\n\
+\n\
+    vec2 position = ray.pos.xy;\n\
+    vec2 direction = ray.dir.xy;\n\
+    float approachRate = dot(direction, planeNormal);\n\
+    float distance = -dot(position, planeNormal);\n\
+\n\
+    float t = (approachRate == 0.0)\n\
+        ? NO_HIT\n\
+        : distance / approachRate;\n\
+\n\
+    return vec4(planeNormal, 0.0, t);\n\
+}\n\
+\n\
+RayShapeIntersection intersectHalfSpace(in Ray ray, in float angle, in bool positiveNormal)\n\
+{\n\
+    vec4 intersection = intersectLongitude(ray, angle, positiveNormal);\n\
+    vec4 farSide = vec4(normalize(ray.dir), INF_HIT);\n\
+\n\
+    bool hitFront = (intersection.w > 0.0) == (dot(ray.pos.xy, intersection.xy) > 0.0);\n\
+    if (!hitFront) {\n\
+        return RayShapeIntersection(intersection, farSide);\n\
+    } else {\n\
+        return RayShapeIntersection(-1.0 * farSide, intersection);\n\
+    }\n\
+}\n\
+\n\
+void intersectFlippedWedge(in Ray ray, in vec2 minMaxAngle, out RayShapeIntersection intersections[2])\n\
+{\n\
+    intersections[0] = transformNormalsToEC(intersectHalfSpace(ray, minMaxAngle.x, false));\n\
+    intersections[1] = transformNormalsToEC(intersectHalfSpace(ray, minMaxAngle.y, true));\n\
+}\n\
+\n\
+bool hitPositiveHalfPlane(in Ray ray, in vec4 intersection, in bool positiveNormal) {\n\
+    float normalSign = positiveNormal ? 1.0 : -1.0;\n\
+    vec2 planeDirection = vec2(intersection.y, -intersection.x) * normalSign;\n\
+    vec2 hit = ray.pos.xy + intersection.w * ray.dir.xy;\n\
+    return dot(hit, planeDirection) > 0.0;\n\
+}\n\
+\n\
+void intersectHalfPlane(in Ray ray, in float angle, out RayShapeIntersection intersections[2]) {\n\
+    vec4 intersection = intersectLongitude(ray, angle, true);\n\
+    vec4 farSide = vec4(normalize(ray.dir), INF_HIT);\n\
+    bool hitPositiveSide = hitPositiveHalfPlane(ray, intersection, true);\n\
+\n\
+    farSide = transformNormalToEC(farSide);\n\
+\n\
+    if (hitPositiveSide) {\n\
+        intersection = transformNormalToEC(intersection);\n\
+        intersections[0].entry = -1.0 * farSide;\n\
+        intersections[0].exit = vec4(-1.0 * intersection.xyz, intersection.w);\n\
+        intersections[1].entry = intersection;\n\
+        intersections[1].exit = farSide;\n\
+    } else {\n\
+        vec4 miss = vec4(normalize(czm_normal * ray.dir), NO_HIT);\n\
+        intersections[0].entry = -1.0 * farSide;\n\
+        intersections[0].exit = farSide;\n\
+        intersections[1].entry = miss;\n\
+        intersections[1].exit = miss;\n\
+    }\n\
+}\n\
+\n\
+RayShapeIntersection intersectRegularWedge(in Ray ray, in vec2 minMaxAngle)\n\
+{\n\
+    // Note: works for maxAngle > minAngle + pi, where the \"regular wedge\"\n\
+    // is actually a negative volume.\n\
+    // Compute intersections with the two planes.\n\
+    // Normals will point toward the \"outside\" (negative space)\n\
+    vec4 intersect1 = intersectLongitude(ray, minMaxAngle.x, false);\n\
+    vec4 intersect2 = intersectLongitude(ray, minMaxAngle.y, true);\n\
+\n\
+    // Choose intersection with smallest T as the \"first\", the other as \"last\"\n\
+    // Note: first or last could be in the \"shadow\" wedge, beyond the tip\n\
+    bool inOrder = intersect1.w <= intersect2.w;\n\
+    vec4 first = inOrder ? intersect1 : intersect2;\n\
+    vec4 last = inOrder ? intersect2 : intersect1;\n\
+\n\
+    bool firstIsAhead = first.w >= 0.0;\n\
+    bool startedInsideFirst = dot(ray.pos.xy, first.xy) < 0.0;\n\
+    bool exitFromInside = firstIsAhead == startedInsideFirst;\n\
+    bool lastIsAhead = last.w > 0.0;\n\
+    bool startedOutsideLast = dot(ray.pos.xy, last.xy) >= 0.0;\n\
+    bool enterFromOutside = lastIsAhead == startedOutsideLast;\n\
+\n\
+    vec4 farSide = vec4(normalize(ray.dir), INF_HIT);\n\
+    vec4 miss = vec4(normalize(ray.dir), NO_HIT);\n\
+\n\
+    if (exitFromInside && enterFromOutside) {\n\
+        // Ray crosses both faces of negative wedge, exiting then entering the positive shape\n\
+        return transformNormalsToEC(RayShapeIntersection(first, last));\n\
+    } else if (!exitFromInside && enterFromOutside) {\n\
+        // Ray starts inside wedge. last is in shadow wedge, and first is actually the entry\n\
+        return transformNormalsToEC(RayShapeIntersection(-1.0 * farSide, first));\n\
+    } else if (exitFromInside && !enterFromOutside) {\n\
+        // First intersection was in the shadow wedge, so last is actually the exit\n\
+        return transformNormalsToEC(RayShapeIntersection(last, farSide));\n\
+    } else { // !exitFromInside && !enterFromOutside\n\
+        // Both intersections were in the shadow wedge\n\
+        return transformNormalsToEC(RayShapeIntersection(miss, miss));\n\
+    }\n\
+}\n\
+";

node_modules/@cesium/engine/Source/Shaders/Voxels/IntersectPlane.glsl

@@ -0,0 +1,80 @@
+// See IntersectionUtils.glsl for the definitions of Ray, Intersections, INF_HIT,
+// NO_HIT, setShapeIntersection
+
+/* Clipping plane defines (set in Scene/VoxelRenderResources.js)
+#define CLIPPING_PLANES_UNION
+#define CLIPPING_PLANES_COUNT
+#define CLIPPING_PLANES_INTERSECTION_INDEX
+*/
+
+uniform sampler2D u_clippingPlanesTexture;
+uniform mat4 u_clippingPlanesMatrix;
+
+// Plane is in Hessian Normal Form
+vec4 intersectPlane(in Ray ray, in vec4 plane) {
+    vec3 n = plane.xyz; // normal
+    float w = plane.w; // -dot(pointOnPlane, normal)
+
+    float a = dot(ray.pos, n);
+    float b = dot(ray.dir, n);
+    float t = -(w + a) / b;
+
+    return vec4(n, t);
+}
+
+#ifdef CLIPPING_PLANES
+void intersectClippingPlanes(in Ray ray, inout Intersections ix) {
+    vec4 backSide = vec4(-ray.dir, -INF_HIT);
+    vec4 farSide = vec4(ray.dir, +INF_HIT);
+    RayShapeIntersection clippingVolume;
+
+    #if (CLIPPING_PLANES_COUNT == 1)
+        // Union and intersection are the same when there's one clipping plane, and the code
+        // is more simplified.
+        vec4 planeUv = getClippingPlane(u_clippingPlanesTexture, 0);
+        vec4 intersection = intersectPlane(ray, planeUv);
+        bool reflects = dot(ray.dir, intersection.xyz) < 0.0;
+        clippingVolume.entry = reflects ? backSide : intersection;
+        clippingVolume.exit = reflects ? intersection : farSide;
+        setShapeIntersection(ix, CLIPPING_PLANES_INTERSECTION_INDEX, clippingVolume);
+    #elif defined(CLIPPING_PLANES_UNION)
+        vec4 firstTransmission = vec4(ray.dir, +INF_HIT);
+        vec4 lastReflection = vec4(-ray.dir, -INF_HIT);
+        for (int i = 0; i < CLIPPING_PLANES_COUNT; i++) {
+            vec4 planeUv = getClippingPlane(u_clippingPlanesTexture, i);
+            vec4 intersection = intersectPlane(ray, planeUv);
+            if (dot(ray.dir, planeUv.xyz) > 0.0) {
+                firstTransmission = intersection.w <= firstTransmission.w ? intersection : firstTransmission;
+            } else {
+                lastReflection = intersection.w >= lastReflection.w ? intersection : lastReflection;
+            }
+        }
+        clippingVolume.entry = backSide;
+        clippingVolume.exit = lastReflection;
+        setShapeIntersection(ix, CLIPPING_PLANES_INTERSECTION_INDEX + 0, clippingVolume);
+        clippingVolume.entry = firstTransmission;
+        clippingVolume.exit = farSide;
+        setShapeIntersection(ix, CLIPPING_PLANES_INTERSECTION_INDEX + 1, clippingVolume);
+    #else // intersection
+        vec4 lastTransmission = vec4(ray.dir, -INF_HIT);
+        vec4 firstReflection = vec4(-ray.dir, +INF_HIT);
+        for (int i = 0; i < CLIPPING_PLANES_COUNT; i++) {
+            vec4 planeUv = getClippingPlane(u_clippingPlanesTexture, i);
+            vec4 intersection = intersectPlane(ray, planeUv);
+            if (dot(ray.dir, planeUv.xyz) > 0.0) {
+                lastTransmission = intersection.w > lastTransmission.w ? intersection : lastTransmission;
+            } else {
+                firstReflection = intersection.w < firstReflection.w ? intersection: firstReflection;
+            }
+        }
+        if (lastTransmission.w < firstReflection.w) {
+            clippingVolume.entry = lastTransmission;
+            clippingVolume.exit = firstReflection;
+        } else {
+            clippingVolume.entry = vec4(-ray.dir, NO_HIT);
+            clippingVolume.exit = vec4(ray.dir, NO_HIT);
+        }
+        setShapeIntersection(ix, CLIPPING_PLANES_INTERSECTION_INDEX, clippingVolume);
+    #endif
+}
+#endif

node_modules/@cesium/engine/Source/Shaders/Voxels/IntersectPlane.js

@@ -0,0 +1,82 @@
+//This file is automatically rebuilt by the Cesium build process.
+export default "// See IntersectionUtils.glsl for the definitions of Ray, Intersections, INF_HIT,\n\
+// NO_HIT, setShapeIntersection\n\
+\n\
+/* Clipping plane defines (set in Scene/VoxelRenderResources.js)\n\
+#define CLIPPING_PLANES_UNION\n\
+#define CLIPPING_PLANES_COUNT\n\
+#define CLIPPING_PLANES_INTERSECTION_INDEX\n\
+*/\n\
+\n\
+uniform sampler2D u_clippingPlanesTexture;\n\
+uniform mat4 u_clippingPlanesMatrix;\n\
+\n\
+// Plane is in Hessian Normal Form\n\
+vec4 intersectPlane(in Ray ray, in vec4 plane) {\n\
+    vec3 n = plane.xyz; // normal\n\
+    float w = plane.w; // -dot(pointOnPlane, normal)\n\
+\n\
+    float a = dot(ray.pos, n);\n\
+    float b = dot(ray.dir, n);\n\
+    float t = -(w + a) / b;\n\
+\n\
+    return vec4(n, t);\n\
+}\n\
+\n\
+#ifdef CLIPPING_PLANES\n\
+void intersectClippingPlanes(in Ray ray, inout Intersections ix) {\n\
+    vec4 backSide = vec4(-ray.dir, -INF_HIT);\n\
+    vec4 farSide = vec4(ray.dir, +INF_HIT);\n\
+    RayShapeIntersection clippingVolume;\n\
+\n\
+    #if (CLIPPING_PLANES_COUNT == 1)\n\
+        // Union and intersection are the same when there's one clipping plane, and the code\n\
+        // is more simplified.\n\
+        vec4 planeUv = getClippingPlane(u_clippingPlanesTexture, 0);\n\
+        vec4 intersection = intersectPlane(ray, planeUv);\n\
+        bool reflects = dot(ray.dir, intersection.xyz) < 0.0;\n\
+        clippingVolume.entry = reflects ? backSide : intersection;\n\
+        clippingVolume.exit = reflects ? intersection : farSide;\n\
+        setShapeIntersection(ix, CLIPPING_PLANES_INTERSECTION_INDEX, clippingVolume);\n\
+    #elif defined(CLIPPING_PLANES_UNION)\n\
+        vec4 firstTransmission = vec4(ray.dir, +INF_HIT);\n\
+        vec4 lastReflection = vec4(-ray.dir, -INF_HIT);\n\
+        for (int i = 0; i < CLIPPING_PLANES_COUNT; i++) {\n\
+            vec4 planeUv = getClippingPlane(u_clippingPlanesTexture, i);\n\
+            vec4 intersection = intersectPlane(ray, planeUv);\n\
+            if (dot(ray.dir, planeUv.xyz) > 0.0) {\n\
+                firstTransmission = intersection.w <= firstTransmission.w ? intersection : firstTransmission;\n\
+            } else {\n\
+                lastReflection = intersection.w >= lastReflection.w ? intersection : lastReflection;\n\
+            }\n\
+        }\n\
+        clippingVolume.entry = backSide;\n\
+        clippingVolume.exit = lastReflection;\n\
+        setShapeIntersection(ix, CLIPPING_PLANES_INTERSECTION_INDEX + 0, clippingVolume);\n\
+        clippingVolume.entry = firstTransmission;\n\
+        clippingVolume.exit = farSide;\n\
+        setShapeIntersection(ix, CLIPPING_PLANES_INTERSECTION_INDEX + 1, clippingVolume);\n\
+    #else // intersection\n\
+        vec4 lastTransmission = vec4(ray.dir, -INF_HIT);\n\
+        vec4 firstReflection = vec4(-ray.dir, +INF_HIT);\n\
+        for (int i = 0; i < CLIPPING_PLANES_COUNT; i++) {\n\
+            vec4 planeUv = getClippingPlane(u_clippingPlanesTexture, i);\n\
+            vec4 intersection = intersectPlane(ray, planeUv);\n\
+            if (dot(ray.dir, planeUv.xyz) > 0.0) {\n\
+                lastTransmission = intersection.w > lastTransmission.w ? intersection : lastTransmission;\n\
+            } else {\n\
+                firstReflection = intersection.w < firstReflection.w ? intersection: firstReflection;\n\
+            }\n\
+        }\n\
+        if (lastTransmission.w < firstReflection.w) {\n\
+            clippingVolume.entry = lastTransmission;\n\
+            clippingVolume.exit = firstReflection;\n\
+        } else {\n\
+            clippingVolume.entry = vec4(-ray.dir, NO_HIT);\n\
+            clippingVolume.exit = vec4(ray.dir, NO_HIT);\n\
+        }\n\
+        setShapeIntersection(ix, CLIPPING_PLANES_INTERSECTION_INDEX, clippingVolume);\n\
+    #endif\n\
+}\n\
+#endif\n\
+";

node_modules/@cesium/engine/Source/Shaders/Voxels/Intersection.glsl

@@ -0,0 +1,50 @@
+// Main intersection function for Voxel scenes.
+// See IntersectBox.glsl, IntersectCylinder.glsl, or IntersectEllipsoid.glsl
+// for the definition of intersectShape. The appropriate function is selected
+// based on the VoxelPrimitive shape type, and added to the shader in
+// Scene/VoxelRenderResources.js.
+// See also IntersectClippingPlane.glsl and IntersectDepth.glsl.
+// See IntersectionUtils.glsl for the definitions of Ray, NO_HIT,
+// getFirstIntersection, initializeIntersections, nextIntersection.
+
+/* Intersection defines (set in Scene/VoxelRenderResources.js)
+#define INTERSECTION_COUNT ###
+*/
+
+RayShapeIntersection intersectScene(in vec2 screenCoord, in Ray ray, in Ray rayEC, out Intersections ix) {
+    // Do a ray-shape intersection to find the exact starting and ending points.
+    intersectShape(ray, rayEC, ix);
+
+    // Exit early if the positive shape was completely missed or behind the ray.
+    RayShapeIntersection intersection = getFirstIntersection(ix);
+    if (intersection.entry.w == NO_HIT) {
+        // Positive shape was completely missed - so exit early.
+        return intersection;
+    }
+
+    // Clipping planes
+    #if defined(CLIPPING_PLANES)
+        intersectClippingPlanes(ray, ix);
+    #endif
+
+    // Depth
+    intersectDepth(screenCoord, rayEC, ix);
+
+    // Find the first intersection that's in front of the ray
+    #if (INTERSECTION_COUNT > 1)
+        initializeIntersections(ix);
+        for (int i = 0; i < INTERSECTION_COUNT; ++i) {
+            intersection = nextIntersection(ix);
+            if (intersection.exit.w > 0.0) {
+                // Set start to 0.0 when ray is inside the shape.
+                intersection.entry.w = max(intersection.entry.w, 0.0);
+                break;
+            }
+        }
+    #else
+        // Set start to 0.0 when ray is inside the shape.
+        intersection.entry.w = max(intersection.entry.w, 0.0);
+    #endif
+
+    return intersection;
+}

node_modules/@cesium/engine/Source/Shaders/Voxels/Intersection.js

@@ -0,0 +1,52 @@
+//This file is automatically rebuilt by the Cesium build process.
+export default "// Main intersection function for Voxel scenes.\n\
+// See IntersectBox.glsl, IntersectCylinder.glsl, or IntersectEllipsoid.glsl\n\
+// for the definition of intersectShape. The appropriate function is selected\n\
+// based on the VoxelPrimitive shape type, and added to the shader in\n\
+// Scene/VoxelRenderResources.js.\n\
+// See also IntersectClippingPlane.glsl and IntersectDepth.glsl.\n\
+// See IntersectionUtils.glsl for the definitions of Ray, NO_HIT,\n\
+// getFirstIntersection, initializeIntersections, nextIntersection.\n\
+\n\
+/* Intersection defines (set in Scene/VoxelRenderResources.js)\n\
+#define INTERSECTION_COUNT ###\n\
+*/\n\
+\n\
+RayShapeIntersection intersectScene(in vec2 screenCoord, in Ray ray, in Ray rayEC, out Intersections ix) {\n\
+    // Do a ray-shape intersection to find the exact starting and ending points.\n\
+    intersectShape(ray, rayEC, ix);\n\
+\n\
+    // Exit early if the positive shape was completely missed or behind the ray.\n\
+    RayShapeIntersection intersection = getFirstIntersection(ix);\n\
+    if (intersection.entry.w == NO_HIT) {\n\
+        // Positive shape was completely missed - so exit early.\n\
+        return intersection;\n\
+    }\n\
+\n\
+    // Clipping planes\n\
+    #if defined(CLIPPING_PLANES)\n\
+        intersectClippingPlanes(ray, ix);\n\
+    #endif\n\
+\n\
+    // Depth\n\
+    intersectDepth(screenCoord, rayEC, ix);\n\
+\n\
+    // Find the first intersection that's in front of the ray\n\
+    #if (INTERSECTION_COUNT > 1)\n\
+        initializeIntersections(ix);\n\
+        for (int i = 0; i < INTERSECTION_COUNT; ++i) {\n\
+            intersection = nextIntersection(ix);\n\
+            if (intersection.exit.w > 0.0) {\n\
+                // Set start to 0.0 when ray is inside the shape.\n\
+                intersection.entry.w = max(intersection.entry.w, 0.0);\n\
+                break;\n\
+            }\n\
+        }\n\
+    #else\n\
+        // Set start to 0.0 when ray is inside the shape.\n\
+        intersection.entry.w = max(intersection.entry.w, 0.0);\n\
+    #endif\n\
+\n\
+    return intersection;\n\
+}\n\
+";

node_modules/@cesium/engine/Source/Shaders/Voxels/IntersectionUtils.glsl

@@ -0,0 +1,167 @@
+/* Intersection defines
+#define INTERSECTION_COUNT ###
+*/
+
+#define NO_HIT (-czm_infinity)
+#define INF_HIT (czm_infinity * 0.5)
+
+struct RayShapeIntersection {
+    vec4 entry;
+    vec4 exit;
+};
+
+vec4 intersectionMin(in vec4 intersect0, in vec4 intersect1)
+{
+    if (intersect0.w == NO_HIT) {
+        return intersect1;
+    } else if (intersect1.w == NO_HIT) {
+        return intersect0;
+    }
+    return (intersect0.w <= intersect1.w) ? intersect0 : intersect1;
+}
+
+vec4 intersectionMax(in vec4 intersect0, in vec4 intersect1)
+{
+    return (intersect0.w >= intersect1.w) ? intersect0 : intersect1;
+}
+
+RayShapeIntersection intersectIntersections(in Ray ray, in RayShapeIntersection intersect0, in RayShapeIntersection intersect1)
+{
+    bool missed = (intersect0.entry.w == NO_HIT) ||
+        (intersect1.entry.w == NO_HIT) ||
+        (intersect0.exit.w < intersect1.entry.w) ||
+        (intersect0.entry.w > intersect1.exit.w);
+    if (missed) {
+        vec4 miss = vec4(normalize(ray.dir), NO_HIT);
+        return RayShapeIntersection(miss, miss);
+    }
+
+    vec4 entry = intersectionMax(intersect0.entry, intersect1.entry);
+    vec4 exit = intersectionMin(intersect0.exit, intersect1.exit);
+
+    return RayShapeIntersection(entry, exit);
+}
+
+struct Intersections {
+    // Don't access these member variables directly - call the functions instead.
+
+    // Store an array of ray-surface intersections. Each intersection is composed of:
+    //  .xyz for the surface normal at the intersection point
+    //  .w for the T value
+    // The scale of the normal 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
+    // INTERSECTION_COUNT is the number of ray-*shape* (volume) intersections,
+    // so we need twice as many to track ray-*surface* intersections
+    vec4 intersections[INTERSECTION_COUNT * 2];
+    float distanceToDepthBuffer;
+
+    #if (INTERSECTION_COUNT > 1)
+        // Maintain state for future nextIntersection calls
+        int index;
+        int surroundCount;
+        bool surroundIsPositive;
+    #endif
+};
+
+RayShapeIntersection getFirstIntersection(in Intersections ix) 
+{
+    return RayShapeIntersection(ix.intersections[0], ix.intersections[1]);
+}
+
+vec4 encodeIntersectionType(vec4 intersection, int index, bool entry)
+{
+    float scale = float(index > 0) * 2.0 + float(!entry) + 1.0;
+    return vec4(intersection.xyz * scale, intersection.w);
+}
+
+// Use defines instead of real functions because WebGL1 cannot access array with non-constant index.
+#define setIntersection(/*inout Intersections*/ ix, /*int*/ index, /*float*/ t, /*bool*/ positive, /*bool*/ enter) (ix).intersections[(index)] = vec4(0.0, float(!positive) * 2.0 + float(!enter) + 1.0, 0.0, (t))
+#define setIntersectionPair(/*inout Intersections*/ ix, /*int*/ index, /*vec2*/ entryExit) (ix).intersections[(index) * 2 + 0] = vec4(0.0, float((index) > 0) * 2.0 + 1.0, 0.0, (entryExit).x); (ix).intersections[(index) * 2 + 1] = vec4(0.0, float((index) > 0) * 2.0 + 2.0, 0.0, (entryExit).y)
+#define setSurfaceIntersection(/*inout Intersections*/ ix, /*int*/ index, /*vec4*/ intersection, /*bool*/ positive, /*bool*/ enter) (ix).intersections[(index)] = encodeIntersectionType((intersection), int(!positive), (enter))
+#define setShapeIntersection(/*inout Intersections*/ ix, /*int*/ index, /*RayShapeIntersection*/ intersection) (ix).intersections[(index) * 2 + 0] = encodeIntersectionType((intersection).entry, (index), true); (ix).intersections[(index) * 2 + 1] = encodeIntersectionType((intersection).exit, (index), false)
+
+#if (INTERSECTION_COUNT > 1)
+void initializeIntersections(inout Intersections ix) {
+    // Sort the intersections from min T to max T with bubble sort.
+    // Note: If this sorting function changes, some of the intersection test may
+    // need to be updated. Search for "bubble sort" to find those areas.
+    const int sortPasses = INTERSECTION_COUNT * 2 - 1;
+    for (int n = sortPasses; n > 0; --n) {
+        for (int i = 0; i < sortPasses; ++i) {
+            // The loop should be: for (i = 0; i < n; ++i) {...} but WebGL1 cannot
+            // loop with non-constant condition, so it has to break early instead
+            if (i >= n) { break; }
+
+            vec4 intersect0 = ix.intersections[i + 0];
+            vec4 intersect1 = ix.intersections[i + 1];
+
+            bool inOrder = intersect0.w <= intersect1.w;
+
+            ix.intersections[i + 0] = inOrder ? intersect0 : intersect1;
+            ix.intersections[i + 1] = inOrder ? intersect1 : intersect0;
+        }
+    }
+
+    // Prepare initial state for nextIntersection
+    ix.index = 0;
+    ix.surroundCount = 0;
+    ix.surroundIsPositive = false;
+}
+#endif
+
+#if (INTERSECTION_COUNT > 1)
+RayShapeIntersection nextIntersection(inout Intersections ix) {
+    vec4 surfaceIntersection = vec4(0.0, 0.0, 0.0, NO_HIT);
+    RayShapeIntersection shapeIntersection = RayShapeIntersection(surfaceIntersection, surfaceIntersection);
+
+    const int passCount = INTERSECTION_COUNT * 2;
+
+    if (ix.index == passCount) {
+        return shapeIntersection;
+    }
+
+    for (int i = 0; i < passCount; ++i) {
+        // The loop should be: for (i = ix.index; i < passCount; ++i) {...} but WebGL1 cannot
+        // loop with non-constant condition, so it has to continue instead.
+        if (i < ix.index) {
+            continue;
+        }
+
+        ix.index = i + 1;
+
+        surfaceIntersection = ix.intersections[i];
+        int intersectionType = int(length(surfaceIntersection.xyz) - 0.5);
+        bool currShapeIsPositive = intersectionType < 2;
+        bool enter = intersectionType % 2 == 0;
+
+        ix.surroundCount += enter ? +1 : -1;
+        ix.surroundIsPositive = currShapeIsPositive ? enter : ix.surroundIsPositive;
+
+        // entering positive or exiting negative
+        if (ix.surroundCount == 1 && ix.surroundIsPositive && enter == currShapeIsPositive) {
+            shapeIntersection.entry = surfaceIntersection;
+        }
+
+        // exiting positive or entering negative after being inside positive
+        bool exitPositive = !enter && currShapeIsPositive && ix.surroundCount == 0;
+        bool enterNegativeFromPositive = enter && !currShapeIsPositive && ix.surroundCount == 2 && ix.surroundIsPositive;
+        if (exitPositive || enterNegativeFromPositive) {
+            shapeIntersection.exit = surfaceIntersection;
+
+            // entry and exit have been found, so the loop can stop
+            if (exitPositive) {
+                // After exiting positive shape there is nothing left to intersect, so jump to the end index.
+                ix.index = passCount;
+            }
+            break;
+        }
+    }
+
+    return shapeIntersection;
+}
+#endif
+
+// NOTE: initializeIntersections, nextIntersection aren't even declared unless INTERSECTION_COUNT > 1

node_modules/@cesium/engine/Source/Shaders/Voxels/IntersectionUtils.js

@@ -0,0 +1,169 @@
+//This file is automatically rebuilt by the Cesium build process.
+export default "/* Intersection defines\n\
+#define INTERSECTION_COUNT ###\n\
+*/\n\
+\n\
+#define NO_HIT (-czm_infinity)\n\
+#define INF_HIT (czm_infinity * 0.5)\n\
+\n\
+struct RayShapeIntersection {\n\
+    vec4 entry;\n\
+    vec4 exit;\n\
+};\n\
+\n\
+vec4 intersectionMin(in vec4 intersect0, in vec4 intersect1)\n\
+{\n\
+    if (intersect0.w == NO_HIT) {\n\
+        return intersect1;\n\
+    } else if (intersect1.w == NO_HIT) {\n\
+        return intersect0;\n\
+    }\n\
+    return (intersect0.w <= intersect1.w) ? intersect0 : intersect1;\n\
+}\n\
+\n\
+vec4 intersectionMax(in vec4 intersect0, in vec4 intersect1)\n\
+{\n\
+    return (intersect0.w >= intersect1.w) ? intersect0 : intersect1;\n\
+}\n\
+\n\
+RayShapeIntersection intersectIntersections(in Ray ray, in RayShapeIntersection intersect0, in RayShapeIntersection intersect1)\n\
+{\n\
+    bool missed = (intersect0.entry.w == NO_HIT) ||\n\
+        (intersect1.entry.w == NO_HIT) ||\n\
+        (intersect0.exit.w < intersect1.entry.w) ||\n\
+        (intersect0.entry.w > intersect1.exit.w);\n\
+    if (missed) {\n\
+        vec4 miss = vec4(normalize(ray.dir), NO_HIT);\n\
+        return RayShapeIntersection(miss, miss);\n\
+    }\n\
+\n\
+    vec4 entry = intersectionMax(intersect0.entry, intersect1.entry);\n\
+    vec4 exit = intersectionMin(intersect0.exit, intersect1.exit);\n\
+\n\
+    return RayShapeIntersection(entry, exit);\n\
+}\n\
+\n\
+struct Intersections {\n\
+    // Don't access these member variables directly - call the functions instead.\n\
+\n\
+    // Store an array of ray-surface intersections. Each intersection is composed of:\n\
+    //  .xyz for the surface normal at the intersection point\n\
+    //  .w for the T value\n\
+    // The scale of the normal encodes the shape intersection type:\n\
+    //  length(intersection.xyz) = 1: positive shape entry\n\
+    //  length(intersection.xyz) = 2: positive shape exit\n\
+    //  length(intersection.xyz) = 3: negative shape entry\n\
+    //  length(intersection.xyz) = 4: negative shape exit\n\
+    // INTERSECTION_COUNT is the number of ray-*shape* (volume) intersections,\n\
+    // so we need twice as many to track ray-*surface* intersections\n\
+    vec4 intersections[INTERSECTION_COUNT * 2];\n\
+    float distanceToDepthBuffer;\n\
+\n\
+    #if (INTERSECTION_COUNT > 1)\n\
+        // Maintain state for future nextIntersection calls\n\
+        int index;\n\
+        int surroundCount;\n\
+        bool surroundIsPositive;\n\
+    #endif\n\
+};\n\
+\n\
+RayShapeIntersection getFirstIntersection(in Intersections ix) \n\
+{\n\
+    return RayShapeIntersection(ix.intersections[0], ix.intersections[1]);\n\
+}\n\
+\n\
+vec4 encodeIntersectionType(vec4 intersection, int index, bool entry)\n\
+{\n\
+    float scale = float(index > 0) * 2.0 + float(!entry) + 1.0;\n\
+    return vec4(intersection.xyz * scale, intersection.w);\n\
+}\n\
+\n\
+// Use defines instead of real functions because WebGL1 cannot access array with non-constant index.\n\
+#define setIntersection(/*inout Intersections*/ ix, /*int*/ index, /*float*/ t, /*bool*/ positive, /*bool*/ enter) (ix).intersections[(index)] = vec4(0.0, float(!positive) * 2.0 + float(!enter) + 1.0, 0.0, (t))\n\
+#define setIntersectionPair(/*inout Intersections*/ ix, /*int*/ index, /*vec2*/ entryExit) (ix).intersections[(index) * 2 + 0] = vec4(0.0, float((index) > 0) * 2.0 + 1.0, 0.0, (entryExit).x); (ix).intersections[(index) * 2 + 1] = vec4(0.0, float((index) > 0) * 2.0 + 2.0, 0.0, (entryExit).y)\n\
+#define setSurfaceIntersection(/*inout Intersections*/ ix, /*int*/ index, /*vec4*/ intersection, /*bool*/ positive, /*bool*/ enter) (ix).intersections[(index)] = encodeIntersectionType((intersection), int(!positive), (enter))\n\
+#define setShapeIntersection(/*inout Intersections*/ ix, /*int*/ index, /*RayShapeIntersection*/ intersection) (ix).intersections[(index) * 2 + 0] = encodeIntersectionType((intersection).entry, (index), true); (ix).intersections[(index) * 2 + 1] = encodeIntersectionType((intersection).exit, (index), false)\n\
+\n\
+#if (INTERSECTION_COUNT > 1)\n\
+void initializeIntersections(inout Intersections ix) {\n\
+    // Sort the intersections from min T to max T with bubble sort.\n\
+    // Note: If this sorting function changes, some of the intersection test may\n\
+    // need to be updated. Search for \"bubble sort\" to find those areas.\n\
+    const int sortPasses = INTERSECTION_COUNT * 2 - 1;\n\
+    for (int n = sortPasses; n > 0; --n) {\n\
+        for (int i = 0; i < sortPasses; ++i) {\n\
+            // The loop should be: for (i = 0; i < n; ++i) {...} but WebGL1 cannot\n\
+            // loop with non-constant condition, so it has to break early instead\n\
+            if (i >= n) { break; }\n\
+\n\
+            vec4 intersect0 = ix.intersections[i + 0];\n\
+            vec4 intersect1 = ix.intersections[i + 1];\n\
+\n\
+            bool inOrder = intersect0.w <= intersect1.w;\n\
+\n\
+            ix.intersections[i + 0] = inOrder ? intersect0 : intersect1;\n\
+            ix.intersections[i + 1] = inOrder ? intersect1 : intersect0;\n\
+        }\n\
+    }\n\
+\n\
+    // Prepare initial state for nextIntersection\n\
+    ix.index = 0;\n\
+    ix.surroundCount = 0;\n\
+    ix.surroundIsPositive = false;\n\
+}\n\
+#endif\n\
+\n\
+#if (INTERSECTION_COUNT > 1)\n\
+RayShapeIntersection nextIntersection(inout Intersections ix) {\n\
+    vec4 surfaceIntersection = vec4(0.0, 0.0, 0.0, NO_HIT);\n\
+    RayShapeIntersection shapeIntersection = RayShapeIntersection(surfaceIntersection, surfaceIntersection);\n\
+\n\
+    const int passCount = INTERSECTION_COUNT * 2;\n\
+\n\
+    if (ix.index == passCount) {\n\
+        return shapeIntersection;\n\
+    }\n\
+\n\
+    for (int i = 0; i < passCount; ++i) {\n\
+        // The loop should be: for (i = ix.index; i < passCount; ++i) {...} but WebGL1 cannot\n\
+        // loop with non-constant condition, so it has to continue instead.\n\
+        if (i < ix.index) {\n\
+            continue;\n\
+        }\n\
+\n\
+        ix.index = i + 1;\n\
+\n\
+        surfaceIntersection = ix.intersections[i];\n\
+        int intersectionType = int(length(surfaceIntersection.xyz) - 0.5);\n\
+        bool currShapeIsPositive = intersectionType < 2;\n\
+        bool enter = intersectionType % 2 == 0;\n\
+\n\
+        ix.surroundCount += enter ? +1 : -1;\n\
+        ix.surroundIsPositive = currShapeIsPositive ? enter : ix.surroundIsPositive;\n\
+\n\
+        // entering positive or exiting negative\n\
+        if (ix.surroundCount == 1 && ix.surroundIsPositive && enter == currShapeIsPositive) {\n\
+            shapeIntersection.entry = surfaceIntersection;\n\
+        }\n\
+\n\
+        // exiting positive or entering negative after being inside positive\n\
+        bool exitPositive = !enter && currShapeIsPositive && ix.surroundCount == 0;\n\
+        bool enterNegativeFromPositive = enter && !currShapeIsPositive && ix.surroundCount == 2 && ix.surroundIsPositive;\n\
+        if (exitPositive || enterNegativeFromPositive) {\n\
+            shapeIntersection.exit = surfaceIntersection;\n\
+\n\
+            // entry and exit have been found, so the loop can stop\n\
+            if (exitPositive) {\n\
+                // After exiting positive shape there is nothing left to intersect, so jump to the end index.\n\
+                ix.index = passCount;\n\
+            }\n\
+            break;\n\
+        }\n\
+    }\n\
+\n\
+    return shapeIntersection;\n\
+}\n\
+#endif\n\
+\n\
+// NOTE: initializeIntersections, nextIntersection aren't even declared unless INTERSECTION_COUNT > 1\n\
+";

node_modules/@cesium/engine/Source/Shaders/Voxels/Megatexture.glsl

@@ -0,0 +1,55 @@
+// See Octree.glsl for the definitions of SampleData
+
+/* Megatexture defines (set in Scene/VoxelRenderResources.js)
+#define SAMPLE_COUNT ###
+#define PADDING
+*/
+
+uniform ivec3 u_megatextureTileCounts; // number of tiles in the megatexture, along each axis
+
+vec3 index1DTo3DTexCoord(int index)
+{
+    int tilesPerZ = u_megatextureTileCounts.x * u_megatextureTileCounts.y;
+    int iz = index / tilesPerZ;
+    int remainder = index - iz * tilesPerZ;
+    int iy = remainder / u_megatextureTileCounts.x;
+    int ix = remainder - iy * u_megatextureTileCounts.x;
+    return vec3(ix, iy, iz) / vec3(u_megatextureTileCounts);
+}
+
+Properties getPropertiesFromMegatexture(in SampleData sampleData) {
+    int tileIndex = sampleData.megatextureIndex;
+
+    vec3 voxelCoord = sampleData.inputCoordinate;
+
+    // UV coordinate of the lower corner of the tile in the megatexture
+    vec3 tileUvOffset = index1DTo3DTexCoord(tileIndex);
+
+    // Voxel location
+    vec3 tileDimensions = vec3(u_inputDimensions);
+    vec3 clampedVoxelCoord = clamp(voxelCoord, vec3(0.5), tileDimensions - vec3(0.5));
+    vec3 voxelUv = clampedVoxelCoord / tileDimensions / vec3(u_megatextureTileCounts);
+
+    return getPropertiesFromMegatextureAtUv(tileUvOffset + voxelUv);
+}
+
+// Convert an array of sample datas to a final weighted properties.
+Properties accumulatePropertiesFromMegatexture(in SampleData sampleDatas[SAMPLE_COUNT]) {
+    #if (SAMPLE_COUNT == 1)
+        return getPropertiesFromMegatexture(sampleDatas[0]);
+    #else
+        // When more than one sample is taken the accumulator needs to start at 0
+        Properties properties = clearProperties();
+        for (int i = 0; i < SAMPLE_COUNT; ++i) {
+            float weight = sampleDatas[i].weight;
+
+            // Avoid reading the megatexture when the weight is 0 as it can be costly.
+            if (weight > 0.0) {
+                Properties tempProperties = getPropertiesFromMegatexture(sampleDatas[i]);
+                tempProperties = scaleProperties(tempProperties, weight);
+                properties = sumProperties(properties, tempProperties);
+            }
+        }
+        return properties;
+    #endif
+}

node_modules/@cesium/engine/Source/Shaders/Voxels/Megatexture.js

@@ -0,0 +1,57 @@
+//This file is automatically rebuilt by the Cesium build process.
+export default "// See Octree.glsl for the definitions of SampleData\n\
+\n\
+/* Megatexture defines (set in Scene/VoxelRenderResources.js)\n\
+#define SAMPLE_COUNT ###\n\
+#define PADDING\n\
+*/\n\
+\n\
+uniform ivec3 u_megatextureTileCounts; // number of tiles in the megatexture, along each axis\n\
+\n\
+vec3 index1DTo3DTexCoord(int index)\n\
+{\n\
+    int tilesPerZ = u_megatextureTileCounts.x * u_megatextureTileCounts.y;\n\
+    int iz = index / tilesPerZ;\n\
+    int remainder = index - iz * tilesPerZ;\n\
+    int iy = remainder / u_megatextureTileCounts.x;\n\
+    int ix = remainder - iy * u_megatextureTileCounts.x;\n\
+    return vec3(ix, iy, iz) / vec3(u_megatextureTileCounts);\n\
+}\n\
+\n\
+Properties getPropertiesFromMegatexture(in SampleData sampleData) {\n\
+    int tileIndex = sampleData.megatextureIndex;\n\
+\n\
+    vec3 voxelCoord = sampleData.inputCoordinate;\n\
+\n\
+    // UV coordinate of the lower corner of the tile in the megatexture\n\
+    vec3 tileUvOffset = index1DTo3DTexCoord(tileIndex);\n\
+\n\
+    // Voxel location\n\
+    vec3 tileDimensions = vec3(u_inputDimensions);\n\
+    vec3 clampedVoxelCoord = clamp(voxelCoord, vec3(0.5), tileDimensions - vec3(0.5));\n\
+    vec3 voxelUv = clampedVoxelCoord / tileDimensions / vec3(u_megatextureTileCounts);\n\
+\n\
+    return getPropertiesFromMegatextureAtUv(tileUvOffset + voxelUv);\n\
+}\n\
+\n\
+// Convert an array of sample datas to a final weighted properties.\n\
+Properties accumulatePropertiesFromMegatexture(in SampleData sampleDatas[SAMPLE_COUNT]) {\n\
+    #if (SAMPLE_COUNT == 1)\n\
+        return getPropertiesFromMegatexture(sampleDatas[0]);\n\
+    #else\n\
+        // When more than one sample is taken the accumulator needs to start at 0\n\
+        Properties properties = clearProperties();\n\
+        for (int i = 0; i < SAMPLE_COUNT; ++i) {\n\
+            float weight = sampleDatas[i].weight;\n\
+\n\
+            // Avoid reading the megatexture when the weight is 0 as it can be costly.\n\
+            if (weight > 0.0) {\n\
+                Properties tempProperties = getPropertiesFromMegatexture(sampleDatas[i]);\n\
+                tempProperties = scaleProperties(tempProperties, weight);\n\
+                properties = sumProperties(properties, tempProperties);\n\
+            }\n\
+        }\n\
+        return properties;\n\
+    #endif\n\
+}\n\
+";

node_modules/@cesium/engine/Source/Shaders/Voxels/Octree.glsl

@@ -0,0 +1,254 @@
+// These octree flags must be in sync with GpuOctreeFlag in VoxelTraversal.js
+#define OCTREE_FLAG_INTERNAL 0
+#define OCTREE_FLAG_LEAF 1
+#define OCTREE_FLAG_PACKED_LEAF_FROM_PARENT 2
+
+#define OCTREE_MAX_LEVELS 32 // Harcoded value because GLSL doesn't like variable length loops
+
+uniform sampler2D u_octreeInternalNodeTexture;
+uniform vec2 u_octreeInternalNodeTexelSizeUv;
+uniform int u_octreeInternalNodeTilesPerRow;
+#if (SAMPLE_COUNT > 1)
+uniform sampler2D u_octreeLeafNodeTexture;
+uniform vec2 u_octreeLeafNodeTexelSizeUv;
+uniform int u_octreeLeafNodeTilesPerRow;
+#endif
+uniform ivec3 u_dimensions; // does not include padding, and is in the z-up orientation
+uniform ivec3 u_inputDimensions; // includes padding, and is in the orientation of the input data
+#if defined(PADDING)
+    uniform ivec3 u_paddingBefore;
+#endif
+
+struct OctreeNodeData {
+    int data;
+    int flag;
+};
+
+struct TraversalData {
+    ivec4 octreeCoords;
+    int parentOctreeIndex;
+};
+
+struct TileAndUvCoordinate {
+    ivec4 tileCoords;
+    vec3 tileUv;
+};
+
+struct SampleData {
+    int megatextureIndex;
+    ivec4 tileCoords;
+    vec3 tileUv;
+    vec3 inputCoordinate;
+    #if (SAMPLE_COUNT > 1)
+        float weight;
+    #endif
+};
+
+int normU8_toInt(in float value) {
+    return int(value * 255.0);
+}
+int normU8x2_toInt(in vec2 value) {
+    return int(value.x * 255.0) + 256 * int(value.y * 255.0);
+}
+float normU8x2_toFloat(in vec2 value) {
+    return float(normU8x2_toInt(value)) / 65535.0;
+}
+
+OctreeNodeData getOctreeNodeData(in vec2 octreeUv) {
+    vec4 texData = texture(u_octreeInternalNodeTexture, octreeUv);
+
+    OctreeNodeData data;
+    data.data = normU8x2_toInt(texData.xy);
+    data.flag = normU8x2_toInt(texData.zw);
+    return data;
+}
+
+OctreeNodeData getOctreeChildData(in int parentOctreeIndex, in ivec3 childCoord) {
+    int childIndex = childCoord.z * 4 + childCoord.y * 2 + childCoord.x;
+    int octreeCoordX = (parentOctreeIndex % u_octreeInternalNodeTilesPerRow) * 9 + 1 + childIndex;
+    int octreeCoordY = parentOctreeIndex / u_octreeInternalNodeTilesPerRow;
+    vec2 octreeUv = u_octreeInternalNodeTexelSizeUv * vec2(float(octreeCoordX) + 0.5, float(octreeCoordY) + 0.5);
+    return getOctreeNodeData(octreeUv);
+}
+
+int getOctreeParentIndex(in int octreeIndex) {
+    int octreeCoordX = (octreeIndex % u_octreeInternalNodeTilesPerRow) * 9;
+    int octreeCoordY = octreeIndex / u_octreeInternalNodeTilesPerRow;
+    vec2 octreeUv = u_octreeInternalNodeTexelSizeUv * vec2(float(octreeCoordX) + 0.5, float(octreeCoordY) + 0.5);
+    vec4 parentData = texture(u_octreeInternalNodeTexture, octreeUv);
+    int parentOctreeIndex = normU8x2_toInt(parentData.xy);
+    return parentOctreeIndex;
+}
+
+vec3 getTileUv(in TileAndUvCoordinate tileAndUv, in ivec4 octreeCoords) {
+    int levelDifference = tileAndUv.tileCoords.w - octreeCoords.w;
+    float scalar = exp2(-1.0 * float(levelDifference));
+    vec3 originShift = vec3(tileAndUv.tileCoords.xyz - (octreeCoords.xyz << levelDifference)) * scalar;
+    return tileAndUv.tileUv * scalar + originShift;
+}
+
+vec3 getClampedTileUv(in TileAndUvCoordinate tileAndUv, in ivec4 octreeCoords) {
+    vec3 tileUv = getTileUv(tileAndUv, octreeCoords);
+    return clamp(tileUv, vec3(0.0), vec3(1.0));
+}
+
+void addSampleCoordinates(in TileAndUvCoordinate tileAndUv, inout SampleData sampleData) {
+    vec3 tileUv = getClampedTileUv(tileAndUv, sampleData.tileCoords);
+
+    vec3 inputCoordinate = tileUv * vec3(u_dimensions);
+#if defined(PADDING)
+    inputCoordinate += vec3(u_paddingBefore);
+#endif
+#if defined(Y_UP_METADATA_ORDER)
+#if defined(SHAPE_BOX)
+    float inputY = inputCoordinate.y;
+    inputCoordinate.y = inputCoordinate.z;
+    // u_inputDimensions.z is the y-up dimension along the 3D Tiles y-axis.
+    inputCoordinate.z = float(u_inputDimensions.z) - inputY;
+#elif defined(SHAPE_CYLINDER)
+    float angle = inputCoordinate.y;
+    float height = inputCoordinate.z;
+    #if (!defined(CYLINDER_HAS_SHAPE_BOUNDS_ANGLE))
+    // Account for the different 0-angle convention in glTF vs 3DTiles
+    if (sampleData.tileCoords.w == 0) {
+        float angleCount = float(u_inputDimensions.z);
+        angle = mod(angle + angleCount / 2.0, angleCount);
+    }
+    #endif
+    inputCoordinate.y = height;
+    inputCoordinate.z = angle;
+#endif
+#endif
+
+    sampleData.tileUv = tileUv;
+    sampleData.inputCoordinate = inputCoordinate;
+}
+
+void getOctreeLeafSampleData(in OctreeNodeData data, in ivec4 octreeCoords, out SampleData sampleData) {
+    sampleData.megatextureIndex = data.data;
+    sampleData.tileCoords = (data.flag == OCTREE_FLAG_PACKED_LEAF_FROM_PARENT)
+        ? ivec4(octreeCoords.xyz / 2, octreeCoords.w - 1)
+        : octreeCoords;
+}
+
+#if (SAMPLE_COUNT > 1)
+void getOctreeLeafSampleDatas(in OctreeNodeData data, in ivec4 octreeCoords, out SampleData sampleDatas[SAMPLE_COUNT]) {
+    int leafIndex = data.data;
+    int leafNodeTexelCount = 2;
+    // Adding 0.5 moves to the center of the texel
+    float leafCoordXStart = float((leafIndex % u_octreeLeafNodeTilesPerRow) * leafNodeTexelCount) + 0.5;
+    float leafCoordY = float(leafIndex / u_octreeLeafNodeTilesPerRow) + 0.5;
+
+    // Get an interpolation weight and a flag to determine whether to read the parent texture
+    vec2 leafUv0 = u_octreeLeafNodeTexelSizeUv * vec2(leafCoordXStart + 0.0, leafCoordY);
+    vec4 leafData0 = texture(u_octreeLeafNodeTexture, leafUv0);
+    float lerp = normU8x2_toFloat(leafData0.xy);
+    sampleDatas[0].weight = 1.0 - lerp;
+    sampleDatas[1].weight = lerp;
+    // TODO: this looks wrong? Should be comparing to OCTREE_FLAG_PACKED_LEAF_FROM_PARENT
+    sampleDatas[0].tileCoords = (normU8_toInt(leafData0.z) == 1)
+        ? ivec4(octreeCoords.xyz / 2, octreeCoords.w - 1)
+        : octreeCoords;
+    sampleDatas[1].tileCoords = (normU8_toInt(leafData0.w) == 1)
+        ? ivec4(octreeCoords.xyz / 2, octreeCoords.w - 1)
+        : octreeCoords;
+
+    // Get megatexture indices for both samples
+    vec2 leafUv1 = u_octreeLeafNodeTexelSizeUv * vec2(leafCoordXStart + 1.0, leafCoordY);
+    vec4 leafData1 = texture(u_octreeLeafNodeTexture, leafUv1);
+    sampleDatas[0].megatextureIndex = normU8x2_toInt(leafData1.xy);
+    sampleDatas[1].megatextureIndex = normU8x2_toInt(leafData1.zw);
+}
+#endif
+
+OctreeNodeData traverseOctreeDownwards(in ivec4 tileCoordinate, inout TraversalData traversalData) {
+    OctreeNodeData childData;
+
+    for (int i = 0; i < OCTREE_MAX_LEVELS; ++i) {
+        // tileCoordinate.xyz is defined at the level of detail tileCoordinate.w.
+        // Find the corresponding coordinate at the level traversalData.octreeCoords.w
+        int level = traversalData.octreeCoords.w + 1;
+        int levelDifference = tileCoordinate.w - level;
+        ivec3 coordinateAtLevel = tileCoordinate.xyz >> levelDifference;
+        traversalData.octreeCoords = ivec4(coordinateAtLevel, level);
+
+        ivec3 childCoordinate = coordinateAtLevel & 1;
+        childData = getOctreeChildData(traversalData.parentOctreeIndex, childCoordinate);
+
+        if (childData.flag != OCTREE_FLAG_INTERNAL) {
+            // leaf tile - stop traversing
+            break;
+        }
+
+        traversalData.parentOctreeIndex = childData.data;
+    }
+
+    return childData;
+}
+
+/**
+* Transform a given position to an octree tile coordinate and a position within that tile,
+* and find the corresponding megatexture index and texture coordinates
+*/
+void traverseOctreeFromBeginning(in TileAndUvCoordinate tileAndUv, out TraversalData traversalData, out SampleData sampleDatas[SAMPLE_COUNT]) {
+    traversalData.octreeCoords = ivec4(0);
+    traversalData.parentOctreeIndex = 0;
+
+    OctreeNodeData nodeData = getOctreeNodeData(vec2(0.0));
+    if (nodeData.flag != OCTREE_FLAG_LEAF) {
+        nodeData = traverseOctreeDownwards(tileAndUv.tileCoords, traversalData);
+    }
+
+    #if (SAMPLE_COUNT == 1)
+        getOctreeLeafSampleData(nodeData, traversalData.octreeCoords, sampleDatas[0]);
+        addSampleCoordinates(tileAndUv, sampleDatas[0]);
+    #else
+        getOctreeLeafSampleDatas(nodeData, traversalData.octreeCoords, sampleDatas);
+        addSampleCoordinates(tileAndUv, sampleDatas[0]);
+        addSampleCoordinates(tileAndUv, sampleDatas[1]);
+    #endif
+}
+
+bool insideTile(in ivec4 tileCoordinate, in ivec4 octreeCoords) {
+    int levelDifference = tileCoordinate.w - octreeCoords.w;
+    if (levelDifference < 0) {
+        return false;
+    }
+    ivec3 coordinateAtLevel = tileCoordinate.xyz >> levelDifference;
+    return coordinateAtLevel == octreeCoords.xyz;
+}
+
+void traverseOctreeFromExisting(in TileAndUvCoordinate tileAndUv, inout TraversalData traversalData, inout SampleData sampleDatas[SAMPLE_COUNT]) {
+    ivec4 tileCoords = tileAndUv.tileCoords;
+    if (insideTile(tileCoords, traversalData.octreeCoords)) {
+        for (int i = 0; i < SAMPLE_COUNT; i++) {
+            addSampleCoordinates(tileAndUv, sampleDatas[i]);
+        }
+        return;
+    }
+
+    // Go up tree until we find a parent tile containing tileCoords.
+    // Assumes all parents are available all they way up to the root.
+    for (int i = 0; i < OCTREE_MAX_LEVELS; ++i) {
+        traversalData.octreeCoords.xyz /= 2;
+        traversalData.octreeCoords.w -= 1;
+
+        if (insideTile(tileCoords, traversalData.octreeCoords)) {
+            break;
+        }
+
+        traversalData.parentOctreeIndex = getOctreeParentIndex(traversalData.parentOctreeIndex);
+    }
+
+    // Go down tree
+    OctreeNodeData nodeData = traverseOctreeDownwards(tileCoords, traversalData);
+
+    #if (SAMPLE_COUNT == 1)
+        getOctreeLeafSampleData(nodeData, traversalData.octreeCoords, sampleDatas[0]);
+        addSampleCoordinates(tileAndUv, sampleDatas[0]);
+    #else
+        getOctreeLeafSampleDatas(nodeData, traversalData.octreeCoords, sampleDatas);
+        addSampleCoordinates(tileAndUv, sampleDatas[0]);
+        addSampleCoordinates(tileAndUv, sampleDatas[1]);
+    #endif
+}

node_modules/@cesium/engine/Source/Shaders/Voxels/Octree.js

@@ -0,0 +1,256 @@
+//This file is automatically rebuilt by the Cesium build process.
+export default "// These octree flags must be in sync with GpuOctreeFlag in VoxelTraversal.js\n\
+#define OCTREE_FLAG_INTERNAL 0\n\
+#define OCTREE_FLAG_LEAF 1\n\
+#define OCTREE_FLAG_PACKED_LEAF_FROM_PARENT 2\n\
+\n\
+#define OCTREE_MAX_LEVELS 32 // Harcoded value because GLSL doesn't like variable length loops\n\
+\n\
+uniform sampler2D u_octreeInternalNodeTexture;\n\
+uniform vec2 u_octreeInternalNodeTexelSizeUv;\n\
+uniform int u_octreeInternalNodeTilesPerRow;\n\
+#if (SAMPLE_COUNT > 1)\n\
+uniform sampler2D u_octreeLeafNodeTexture;\n\
+uniform vec2 u_octreeLeafNodeTexelSizeUv;\n\
+uniform int u_octreeLeafNodeTilesPerRow;\n\
+#endif\n\
+uniform ivec3 u_dimensions; // does not include padding, and is in the z-up orientation\n\
+uniform ivec3 u_inputDimensions; // includes padding, and is in the orientation of the input data\n\
+#if defined(PADDING)\n\
+    uniform ivec3 u_paddingBefore;\n\
+#endif\n\
+\n\
+struct OctreeNodeData {\n\
+    int data;\n\
+    int flag;\n\
+};\n\
+\n\
+struct TraversalData {\n\
+    ivec4 octreeCoords;\n\
+    int parentOctreeIndex;\n\
+};\n\
+\n\
+struct TileAndUvCoordinate {\n\
+    ivec4 tileCoords;\n\
+    vec3 tileUv;\n\
+};\n\
+\n\
+struct SampleData {\n\
+    int megatextureIndex;\n\
+    ivec4 tileCoords;\n\
+    vec3 tileUv;\n\
+    vec3 inputCoordinate;\n\
+    #if (SAMPLE_COUNT > 1)\n\
+        float weight;\n\
+    #endif\n\
+};\n\
+\n\
+int normU8_toInt(in float value) {\n\
+    return int(value * 255.0);\n\
+}\n\
+int normU8x2_toInt(in vec2 value) {\n\
+    return int(value.x * 255.0) + 256 * int(value.y * 255.0);\n\
+}\n\
+float normU8x2_toFloat(in vec2 value) {\n\
+    return float(normU8x2_toInt(value)) / 65535.0;\n\
+}\n\
+\n\
+OctreeNodeData getOctreeNodeData(in vec2 octreeUv) {\n\
+    vec4 texData = texture(u_octreeInternalNodeTexture, octreeUv);\n\
+\n\
+    OctreeNodeData data;\n\
+    data.data = normU8x2_toInt(texData.xy);\n\
+    data.flag = normU8x2_toInt(texData.zw);\n\
+    return data;\n\
+}\n\
+\n\
+OctreeNodeData getOctreeChildData(in int parentOctreeIndex, in ivec3 childCoord) {\n\
+    int childIndex = childCoord.z * 4 + childCoord.y * 2 + childCoord.x;\n\
+    int octreeCoordX = (parentOctreeIndex % u_octreeInternalNodeTilesPerRow) * 9 + 1 + childIndex;\n\
+    int octreeCoordY = parentOctreeIndex / u_octreeInternalNodeTilesPerRow;\n\
+    vec2 octreeUv = u_octreeInternalNodeTexelSizeUv * vec2(float(octreeCoordX) + 0.5, float(octreeCoordY) + 0.5);\n\
+    return getOctreeNodeData(octreeUv);\n\
+}\n\
+\n\
+int getOctreeParentIndex(in int octreeIndex) {\n\
+    int octreeCoordX = (octreeIndex % u_octreeInternalNodeTilesPerRow) * 9;\n\
+    int octreeCoordY = octreeIndex / u_octreeInternalNodeTilesPerRow;\n\
+    vec2 octreeUv = u_octreeInternalNodeTexelSizeUv * vec2(float(octreeCoordX) + 0.5, float(octreeCoordY) + 0.5);\n\
+    vec4 parentData = texture(u_octreeInternalNodeTexture, octreeUv);\n\
+    int parentOctreeIndex = normU8x2_toInt(parentData.xy);\n\
+    return parentOctreeIndex;\n\
+}\n\
+\n\
+vec3 getTileUv(in TileAndUvCoordinate tileAndUv, in ivec4 octreeCoords) {\n\
+    int levelDifference = tileAndUv.tileCoords.w - octreeCoords.w;\n\
+    float scalar = exp2(-1.0 * float(levelDifference));\n\
+    vec3 originShift = vec3(tileAndUv.tileCoords.xyz - (octreeCoords.xyz << levelDifference)) * scalar;\n\
+    return tileAndUv.tileUv * scalar + originShift;\n\
+}\n\
+\n\
+vec3 getClampedTileUv(in TileAndUvCoordinate tileAndUv, in ivec4 octreeCoords) {\n\
+    vec3 tileUv = getTileUv(tileAndUv, octreeCoords);\n\
+    return clamp(tileUv, vec3(0.0), vec3(1.0));\n\
+}\n\
+\n\
+void addSampleCoordinates(in TileAndUvCoordinate tileAndUv, inout SampleData sampleData) {\n\
+    vec3 tileUv = getClampedTileUv(tileAndUv, sampleData.tileCoords);\n\
+\n\
+    vec3 inputCoordinate = tileUv * vec3(u_dimensions);\n\
+#if defined(PADDING)\n\
+    inputCoordinate += vec3(u_paddingBefore);\n\
+#endif\n\
+#if defined(Y_UP_METADATA_ORDER)\n\
+#if defined(SHAPE_BOX)\n\
+    float inputY = inputCoordinate.y;\n\
+    inputCoordinate.y = inputCoordinate.z;\n\
+    // u_inputDimensions.z is the y-up dimension along the 3D Tiles y-axis.\n\
+    inputCoordinate.z = float(u_inputDimensions.z) - inputY;\n\
+#elif defined(SHAPE_CYLINDER)\n\
+    float angle = inputCoordinate.y;\n\
+    float height = inputCoordinate.z;\n\
+    #if (!defined(CYLINDER_HAS_SHAPE_BOUNDS_ANGLE))\n\
+    // Account for the different 0-angle convention in glTF vs 3DTiles\n\
+    if (sampleData.tileCoords.w == 0) {\n\
+        float angleCount = float(u_inputDimensions.z);\n\
+        angle = mod(angle + angleCount / 2.0, angleCount);\n\
+    }\n\
+    #endif\n\
+    inputCoordinate.y = height;\n\
+    inputCoordinate.z = angle;\n\
+#endif\n\
+#endif\n\
+\n\
+    sampleData.tileUv = tileUv;\n\
+    sampleData.inputCoordinate = inputCoordinate;\n\
+}\n\
+\n\
+void getOctreeLeafSampleData(in OctreeNodeData data, in ivec4 octreeCoords, out SampleData sampleData) {\n\
+    sampleData.megatextureIndex = data.data;\n\
+    sampleData.tileCoords = (data.flag == OCTREE_FLAG_PACKED_LEAF_FROM_PARENT)\n\
+        ? ivec4(octreeCoords.xyz / 2, octreeCoords.w - 1)\n\
+        : octreeCoords;\n\
+}\n\
+\n\
+#if (SAMPLE_COUNT > 1)\n\
+void getOctreeLeafSampleDatas(in OctreeNodeData data, in ivec4 octreeCoords, out SampleData sampleDatas[SAMPLE_COUNT]) {\n\
+    int leafIndex = data.data;\n\
+    int leafNodeTexelCount = 2;\n\
+    // Adding 0.5 moves to the center of the texel\n\
+    float leafCoordXStart = float((leafIndex % u_octreeLeafNodeTilesPerRow) * leafNodeTexelCount) + 0.5;\n\
+    float leafCoordY = float(leafIndex / u_octreeLeafNodeTilesPerRow) + 0.5;\n\
+\n\
+    // Get an interpolation weight and a flag to determine whether to read the parent texture\n\
+    vec2 leafUv0 = u_octreeLeafNodeTexelSizeUv * vec2(leafCoordXStart + 0.0, leafCoordY);\n\
+    vec4 leafData0 = texture(u_octreeLeafNodeTexture, leafUv0);\n\
+    float lerp = normU8x2_toFloat(leafData0.xy);\n\
+    sampleDatas[0].weight = 1.0 - lerp;\n\
+    sampleDatas[1].weight = lerp;\n\
+    // TODO: this looks wrong? Should be comparing to OCTREE_FLAG_PACKED_LEAF_FROM_PARENT\n\
+    sampleDatas[0].tileCoords = (normU8_toInt(leafData0.z) == 1)\n\
+        ? ivec4(octreeCoords.xyz / 2, octreeCoords.w - 1)\n\
+        : octreeCoords;\n\
+    sampleDatas[1].tileCoords = (normU8_toInt(leafData0.w) == 1)\n\
+        ? ivec4(octreeCoords.xyz / 2, octreeCoords.w - 1)\n\
+        : octreeCoords;\n\
+\n\
+    // Get megatexture indices for both samples\n\
+    vec2 leafUv1 = u_octreeLeafNodeTexelSizeUv * vec2(leafCoordXStart + 1.0, leafCoordY);\n\
+    vec4 leafData1 = texture(u_octreeLeafNodeTexture, leafUv1);\n\
+    sampleDatas[0].megatextureIndex = normU8x2_toInt(leafData1.xy);\n\
+    sampleDatas[1].megatextureIndex = normU8x2_toInt(leafData1.zw);\n\
+}\n\
+#endif\n\
+\n\
+OctreeNodeData traverseOctreeDownwards(in ivec4 tileCoordinate, inout TraversalData traversalData) {\n\
+    OctreeNodeData childData;\n\
+\n\
+    for (int i = 0; i < OCTREE_MAX_LEVELS; ++i) {\n\
+        // tileCoordinate.xyz is defined at the level of detail tileCoordinate.w.\n\
+        // Find the corresponding coordinate at the level traversalData.octreeCoords.w\n\
+        int level = traversalData.octreeCoords.w + 1;\n\
+        int levelDifference = tileCoordinate.w - level;\n\
+        ivec3 coordinateAtLevel = tileCoordinate.xyz >> levelDifference;\n\
+        traversalData.octreeCoords = ivec4(coordinateAtLevel, level);\n\
+\n\
+        ivec3 childCoordinate = coordinateAtLevel & 1;\n\
+        childData = getOctreeChildData(traversalData.parentOctreeIndex, childCoordinate);\n\
+\n\
+        if (childData.flag != OCTREE_FLAG_INTERNAL) {\n\
+            // leaf tile - stop traversing\n\
+            break;\n\
+        }\n\
+\n\
+        traversalData.parentOctreeIndex = childData.data;\n\
+    }\n\
+\n\
+    return childData;\n\
+}\n\
+\n\
+/**\n\
+* Transform a given position to an octree tile coordinate and a position within that tile,\n\
+* and find the corresponding megatexture index and texture coordinates\n\
+*/\n\
+void traverseOctreeFromBeginning(in TileAndUvCoordinate tileAndUv, out TraversalData traversalData, out SampleData sampleDatas[SAMPLE_COUNT]) {\n\
+    traversalData.octreeCoords = ivec4(0);\n\
+    traversalData.parentOctreeIndex = 0;\n\
+\n\
+    OctreeNodeData nodeData = getOctreeNodeData(vec2(0.0));\n\
+    if (nodeData.flag != OCTREE_FLAG_LEAF) {\n\
+        nodeData = traverseOctreeDownwards(tileAndUv.tileCoords, traversalData);\n\
+    }\n\
+\n\
+    #if (SAMPLE_COUNT == 1)\n\
+        getOctreeLeafSampleData(nodeData, traversalData.octreeCoords, sampleDatas[0]);\n\
+        addSampleCoordinates(tileAndUv, sampleDatas[0]);\n\
+    #else\n\
+        getOctreeLeafSampleDatas(nodeData, traversalData.octreeCoords, sampleDatas);\n\
+        addSampleCoordinates(tileAndUv, sampleDatas[0]);\n\
+        addSampleCoordinates(tileAndUv, sampleDatas[1]);\n\
+    #endif\n\
+}\n\
+\n\
+bool insideTile(in ivec4 tileCoordinate, in ivec4 octreeCoords) {\n\
+    int levelDifference = tileCoordinate.w - octreeCoords.w;\n\
+    if (levelDifference < 0) {\n\
+        return false;\n\
+    }\n\
+    ivec3 coordinateAtLevel = tileCoordinate.xyz >> levelDifference;\n\
+    return coordinateAtLevel == octreeCoords.xyz;\n\
+}\n\
+\n\
+void traverseOctreeFromExisting(in TileAndUvCoordinate tileAndUv, inout TraversalData traversalData, inout SampleData sampleDatas[SAMPLE_COUNT]) {\n\
+    ivec4 tileCoords = tileAndUv.tileCoords;\n\
+    if (insideTile(tileCoords, traversalData.octreeCoords)) {\n\
+        for (int i = 0; i < SAMPLE_COUNT; i++) {\n\
+            addSampleCoordinates(tileAndUv, sampleDatas[i]);\n\
+        }\n\
+        return;\n\
+    }\n\
+\n\
+    // Go up tree until we find a parent tile containing tileCoords.\n\
+    // Assumes all parents are available all they way up to the root.\n\
+    for (int i = 0; i < OCTREE_MAX_LEVELS; ++i) {\n\
+        traversalData.octreeCoords.xyz /= 2;\n\
+        traversalData.octreeCoords.w -= 1;\n\
+\n\
+        if (insideTile(tileCoords, traversalData.octreeCoords)) {\n\
+            break;\n\
+        }\n\
+\n\
+        traversalData.parentOctreeIndex = getOctreeParentIndex(traversalData.parentOctreeIndex);\n\
+    }\n\
+\n\
+    // Go down tree\n\
+    OctreeNodeData nodeData = traverseOctreeDownwards(tileCoords, traversalData);\n\
+\n\
+    #if (SAMPLE_COUNT == 1)\n\
+        getOctreeLeafSampleData(nodeData, traversalData.octreeCoords, sampleDatas[0]);\n\
+        addSampleCoordinates(tileAndUv, sampleDatas[0]);\n\
+    #else\n\
+        getOctreeLeafSampleDatas(nodeData, traversalData.octreeCoords, sampleDatas);\n\
+        addSampleCoordinates(tileAndUv, sampleDatas[0]);\n\
+        addSampleCoordinates(tileAndUv, sampleDatas[1]);\n\
+    #endif\n\
+}\n\
+";

node_modules/@cesium/engine/Source/Shaders/Voxels/VoxelFS.glsl

@@ -0,0 +1,267 @@
+// See Intersection.glsl for the definition of intersectScene
+// See IntersectionUtils.glsl for the definition of nextIntersection
+// See convertLocalToBoxUv.glsl, convertLocalToCylinderUv.glsl, or convertLocalToEllipsoidUv.glsl
+// for the definitions of convertLocalToShapeSpaceDerivative and getTileAndUvCoordinate. 
+// The appropriate functions are selected based on the VoxelPrimitive shape type, 
+// and added to the shader in Scene/VoxelRenderResources.js.
+// See Octree.glsl for the definitions of TraversalData, SampleData,
+// traverseOctreeFromBeginning, and traverseOctreeFromExisting
+// See Megatexture.glsl for the definition of accumulatePropertiesFromMegatexture
+
+#define STEP_COUNT_MAX 1000 // Harcoded value because GLSL doesn't like variable length loops
+#if defined(PICKING_VOXEL)
+    #define ALPHA_ACCUM_MAX 0.1
+#else
+    #define ALPHA_ACCUM_MAX 0.98 // Must be > 0.0 and <= 1.0
+#endif
+
+uniform mat4 u_transformPositionViewToLocal;
+uniform mat3 u_transformDirectionViewToLocal;
+uniform vec3 u_cameraPositionLocal;
+uniform vec3 u_cameraDirectionLocal;
+uniform float u_stepSize;
+
+#if defined(PICKING)
+    uniform vec4 u_pickColor;
+#endif
+
+vec3 getSampleSize(in int level) {
+    vec3 sampleCount = exp2(float(level)) * vec3(u_dimensions);
+    vec3 sampleSizeUv = 1.0 / sampleCount;
+    return scaleShapeUvToShapeSpace(sampleSizeUv);
+}
+
+#define MINIMUM_STEP_SCALAR (0.02)
+#define SHIFT_FRACTION (0.001)
+
+/**
+ * Given a coordinate within a tile, and sample spacings along a ray through
+ * the coordinate, find the distance to the points where the ray entered and
+ * exited the voxel cell, along with the surface normals at those points.
+ * The surface normals are returned in shape space coordinates.
+ */
+RayShapeIntersection getVoxelIntersection(in vec3 tileUv, in vec3 sampleSizeAlongRay) {
+    vec3 voxelCoord = tileUv * vec3(u_dimensions);
+    vec3 directions = sign(sampleSizeAlongRay);
+    vec3 positiveDirections = max(directions, 0.0);
+    vec3 entryCoord = mix(ceil(voxelCoord), floor(voxelCoord), positiveDirections);
+    vec3 exitCoord = entryCoord + directions;
+
+    vec3 distanceFromEntry = -abs((entryCoord - voxelCoord) * sampleSizeAlongRay);
+    float lastEntry = maxComponent(distanceFromEntry);
+    bvec3 isLastEntry = equal(distanceFromEntry, vec3(lastEntry));
+    vec3 entryNormal = -1.0 * vec3(isLastEntry) * directions;
+    vec4 entry = vec4(entryNormal, lastEntry);
+
+    vec3 distanceToExit = abs((exitCoord - voxelCoord) * sampleSizeAlongRay);
+    float firstExit = minComponent(distanceToExit);
+    bvec3 isFirstExit = equal(distanceToExit, vec3(firstExit));
+    vec3 exitNormal = vec3(isFirstExit) * directions;
+    vec4 exit = vec4(exitNormal, firstExit);
+
+    return RayShapeIntersection(entry, exit);
+}
+
+vec4 getStepSize(in SampleData sampleData, in Ray viewRay, in RayShapeIntersection shapeIntersection, in mat3 jacobianT, in float currentT) {
+    vec3 gradient = viewRay.dir * jacobianT;
+    vec3 sampleSizeAlongRay = getSampleSize(sampleData.tileCoords.w) / gradient;
+
+    RayShapeIntersection voxelIntersection = getVoxelIntersection(sampleData.tileUv, sampleSizeAlongRay);
+
+    // Transform normal from shape space to Cartesian space to eye space
+    vec3 voxelNormal = jacobianT * voxelIntersection.entry.xyz;
+    voxelNormal = normalize(czm_normal * voxelNormal);
+
+    // Compare with the shape intersection, to choose the appropriate normal
+    vec4 voxelEntry = vec4(voxelNormal, currentT + voxelIntersection.entry.w);
+    vec4 entry = intersectionMax(shapeIntersection.entry, voxelEntry);
+
+    float fixedStep = minComponent(abs(sampleSizeAlongRay)) * u_stepSize;
+    float shift = fixedStep * SHIFT_FRACTION;
+    float dt = voxelIntersection.exit.w + shift;
+    if ((currentT + dt) > shapeIntersection.exit.w) {
+        // Stop at end of shape
+        dt = shapeIntersection.exit.w - currentT + shift;
+    }
+    float stepSize = clamp(dt, fixedStep * MINIMUM_STEP_SCALAR, fixedStep + shift);
+
+    return vec4(entry.xyz, stepSize);
+}
+
+vec2 packIntToVec2(int value) {
+    float shifted = float(value) / 255.0;
+    float lowBits = fract(shifted);
+    float highBits = floor(shifted) / 255.0;
+    return vec2(highBits, lowBits);
+}
+
+vec2 packFloatToVec2(float value) {
+    float lowBits = fract(value);
+    float highBits = floor(value) / 255.0;
+    return vec2(highBits, lowBits);
+}
+
+int getSampleIndex(in SampleData sampleData) {
+    // tileUv = 1.0 is a valid coordinate but sampleIndex = u_inputDimensions is not.
+    // (tileUv = 1.0 corresponds to the far edge of the last sample, at index = u_inputDimensions - 1).
+    // Clamp to [0, voxelDimensions - 0.5) to avoid numerical error before flooring
+    vec3 maxCoordinate = vec3(u_inputDimensions) - vec3(0.5);
+    vec3 inputCoordinate = clamp(sampleData.inputCoordinate, vec3(0.0), maxCoordinate);
+    ivec3 sampleIndex = ivec3(floor(inputCoordinate));
+    // Convert to a 1D index for lookup in a 1D data array
+    return sampleIndex.x + u_inputDimensions.x * (sampleIndex.y + u_inputDimensions.y * sampleIndex.z);
+}
+
+/**
+ * Compute the view ray at the current fragment, in the local coordinates of the shape.
+ */
+Ray getViewRayLocal() {
+    vec4 eyeCoordinates = czm_windowToEyeCoordinates(gl_FragCoord);
+    vec3 origin;
+    vec3 direction;
+    if (czm_orthographicIn3D == 1.0) {
+        eyeCoordinates.z = 0.0;
+        origin = (u_transformPositionViewToLocal * eyeCoordinates).xyz;
+        direction = u_cameraDirectionLocal;
+    } else {
+        origin = u_cameraPositionLocal;
+        direction = u_transformDirectionViewToLocal * normalize(eyeCoordinates.xyz);
+    }
+    return Ray(origin, direction);
+}
+
+Ray getViewRayEC() {
+    vec4 eyeCoordinates = czm_windowToEyeCoordinates(gl_FragCoord);
+    vec3 viewPosEC = (czm_orthographicIn3D == 1.0)
+        ? vec3(eyeCoordinates.xy, 0.0)
+        : vec3(0.0);
+    vec3 viewDirEC = normalize(eyeCoordinates.xyz);
+    return Ray(viewPosEC, viewDirEC);
+}
+
+void main()
+{
+    Ray viewRayLocal = getViewRayLocal();
+    Ray viewRayEC = getViewRayEC();
+
+    Intersections ix;
+    vec2 screenCoord = (gl_FragCoord.xy - czm_viewport.xy) / czm_viewport.zw; // [0,1]
+    RayShapeIntersection shapeIntersection = intersectScene(screenCoord, viewRayLocal, viewRayEC, ix);
+    // Exit early if the scene was completely missed.
+    if (shapeIntersection.entry.w == NO_HIT) {
+        discard;
+    }
+
+    float currentT = shapeIntersection.entry.w;
+    float endT = shapeIntersection.exit.w;
+
+    vec3 positionEC = viewRayEC.pos + currentT * viewRayEC.dir;
+    TileAndUvCoordinate tileAndUv = getTileAndUvCoordinate(positionEC);
+    vec3 positionLocal = viewRayLocal.pos + currentT * viewRayLocal.dir;
+    mat3 jacobianT = convertLocalToShapeSpaceDerivative(positionLocal);
+
+    // Traverse the tree from the start position
+    TraversalData traversalData;
+    SampleData sampleDatas[SAMPLE_COUNT];
+    traverseOctreeFromBeginning(tileAndUv, traversalData, sampleDatas);
+    vec4 step = getStepSize(sampleDatas[0], viewRayLocal, shapeIntersection, jacobianT, currentT);
+
+    FragmentInput fragmentInput;
+    #if defined(STATISTICS)
+        setStatistics(fragmentInput.metadataStatistics);
+    #endif
+
+    czm_modelMaterial materialOutput;
+    vec4 colorAccum = vec4(0.0);
+
+    for (int stepCount = 0; stepCount < STEP_COUNT_MAX; ++stepCount) {
+        // Read properties from the megatexture based on the traversal state
+        Properties properties = accumulatePropertiesFromMegatexture(sampleDatas);
+
+        // Prepare the custom shader inputs
+        copyPropertiesToMetadata(properties, fragmentInput.metadata);
+
+        fragmentInput.attributes.positionEC = positionEC;
+        // Re-normalize normals: some shape intersections may have been scaled to encode positive/negative shapes
+        fragmentInput.attributes.normalEC = normalize(step.xyz);
+
+        fragmentInput.voxel.viewDirUv = viewRayLocal.dir;
+
+        fragmentInput.voxel.travelDistance = step.w;
+        fragmentInput.voxel.stepCount = stepCount;
+        fragmentInput.voxel.tileIndex = sampleDatas[0].megatextureIndex;
+        fragmentInput.voxel.sampleIndex = getSampleIndex(sampleDatas[0]);
+        fragmentInput.voxel.distanceToDepthBuffer = ix.distanceToDepthBuffer - currentT;
+
+        // Run the custom shader
+        fragmentMain(fragmentInput, materialOutput);
+
+        // Sanitize the custom shader output
+        vec4 color = vec4(materialOutput.diffuse, materialOutput.alpha);
+        color.rgb = max(color.rgb, vec3(0.0));
+        color.a = clamp(color.a, 0.0, 1.0);
+
+        // Pre-multiplied alpha blend
+        colorAccum += (1.0 - colorAccum.a) * vec4(color.rgb * color.a, color.a);
+
+        // Stop traversing if the alpha has been fully saturated
+        if (colorAccum.a > ALPHA_ACCUM_MAX) {
+            colorAccum.a = ALPHA_ACCUM_MAX;
+            break;
+        }
+
+        if (step.w == 0.0) {
+            // Shape is infinitely thin. The ray may have hit the edge of a
+            // foreground voxel. Step ahead slightly to check for more voxels
+            step.w = 0.001;
+        }
+
+        // Keep raymarching
+        currentT += step.w;
+        // Check if there's more intersections.
+        if (currentT > endT) {
+            #if (INTERSECTION_COUNT == 1)
+                break;
+            #else
+                shapeIntersection = nextIntersection(ix);
+                if (shapeIntersection.entry.w == NO_HIT) {
+                    break;
+                } else {
+                    // Found another intersection. Resume raymarching there
+                    currentT = shapeIntersection.entry.w;
+                    endT = shapeIntersection.exit.w;
+                }
+            #endif
+        }
+        positionEC = viewRayEC.pos + currentT * viewRayEC.dir;
+        tileAndUv = getTileAndUvCoordinate(positionEC);
+        positionLocal = viewRayLocal.pos + currentT * viewRayLocal.dir;
+        jacobianT = convertLocalToShapeSpaceDerivative(positionLocal);
+
+        // Traverse the tree from the current ray position.
+        // This is similar to traverseOctreeFromBeginning but is faster when the ray is in the same tile as the previous step.
+        traverseOctreeFromExisting(tileAndUv, traversalData, sampleDatas);
+        step = getStepSize(sampleDatas[0], viewRayLocal, shapeIntersection, jacobianT, currentT);
+    }
+
+    // Convert the alpha from [0,ALPHA_ACCUM_MAX] to [0,1]
+    colorAccum.a /= ALPHA_ACCUM_MAX;
+
+    #if defined(PICKING)
+        // If alpha is 0.0 there is nothing to pick
+        if (colorAccum.a == 0.0) {
+            discard;
+        }
+        out_FragColor = u_pickColor;
+    #elif defined(PICKING_VOXEL)
+        // If alpha is 0.0 there is nothing to pick
+        if (colorAccum.a == 0.0) {
+            discard;
+        }
+        vec2 megatextureId = packIntToVec2(sampleDatas[0].megatextureIndex);
+        vec2 sampleIndex = packIntToVec2(getSampleIndex(sampleDatas[0]));
+        out_FragColor = vec4(megatextureId, sampleIndex);
+    #else
+        out_FragColor = colorAccum;
+    #endif
+}

node_modules/@cesium/engine/Source/Shaders/Voxels/VoxelFS.js

@@ -0,0 +1,269 @@
+//This file is automatically rebuilt by the Cesium build process.
+export default "// See Intersection.glsl for the definition of intersectScene\n\
+// See IntersectionUtils.glsl for the definition of nextIntersection\n\
+// See convertLocalToBoxUv.glsl, convertLocalToCylinderUv.glsl, or convertLocalToEllipsoidUv.glsl\n\
+// for the definitions of convertLocalToShapeSpaceDerivative and getTileAndUvCoordinate. \n\
+// The appropriate functions are selected based on the VoxelPrimitive shape type, \n\
+// and added to the shader in Scene/VoxelRenderResources.js.\n\
+// See Octree.glsl for the definitions of TraversalData, SampleData,\n\
+// traverseOctreeFromBeginning, and traverseOctreeFromExisting\n\
+// See Megatexture.glsl for the definition of accumulatePropertiesFromMegatexture\n\
+\n\
+#define STEP_COUNT_MAX 1000 // Harcoded value because GLSL doesn't like variable length loops\n\
+#if defined(PICKING_VOXEL)\n\
+    #define ALPHA_ACCUM_MAX 0.1\n\
+#else\n\
+    #define ALPHA_ACCUM_MAX 0.98 // Must be > 0.0 and <= 1.0\n\
+#endif\n\
+\n\
+uniform mat4 u_transformPositionViewToLocal;\n\
+uniform mat3 u_transformDirectionViewToLocal;\n\
+uniform vec3 u_cameraPositionLocal;\n\
+uniform vec3 u_cameraDirectionLocal;\n\
+uniform float u_stepSize;\n\
+\n\
+#if defined(PICKING)\n\
+    uniform vec4 u_pickColor;\n\
+#endif\n\
+\n\
+vec3 getSampleSize(in int level) {\n\
+    vec3 sampleCount = exp2(float(level)) * vec3(u_dimensions);\n\
+    vec3 sampleSizeUv = 1.0 / sampleCount;\n\
+    return scaleShapeUvToShapeSpace(sampleSizeUv);\n\
+}\n\
+\n\
+#define MINIMUM_STEP_SCALAR (0.02)\n\
+#define SHIFT_FRACTION (0.001)\n\
+\n\
+/**\n\
+ * Given a coordinate within a tile, and sample spacings along a ray through\n\
+ * the coordinate, find the distance to the points where the ray entered and\n\
+ * exited the voxel cell, along with the surface normals at those points.\n\
+ * The surface normals are returned in shape space coordinates.\n\
+ */\n\
+RayShapeIntersection getVoxelIntersection(in vec3 tileUv, in vec3 sampleSizeAlongRay) {\n\
+    vec3 voxelCoord = tileUv * vec3(u_dimensions);\n\
+    vec3 directions = sign(sampleSizeAlongRay);\n\
+    vec3 positiveDirections = max(directions, 0.0);\n\
+    vec3 entryCoord = mix(ceil(voxelCoord), floor(voxelCoord), positiveDirections);\n\
+    vec3 exitCoord = entryCoord + directions;\n\
+\n\
+    vec3 distanceFromEntry = -abs((entryCoord - voxelCoord) * sampleSizeAlongRay);\n\
+    float lastEntry = maxComponent(distanceFromEntry);\n\
+    bvec3 isLastEntry = equal(distanceFromEntry, vec3(lastEntry));\n\
+    vec3 entryNormal = -1.0 * vec3(isLastEntry) * directions;\n\
+    vec4 entry = vec4(entryNormal, lastEntry);\n\
+\n\
+    vec3 distanceToExit = abs((exitCoord - voxelCoord) * sampleSizeAlongRay);\n\
+    float firstExit = minComponent(distanceToExit);\n\
+    bvec3 isFirstExit = equal(distanceToExit, vec3(firstExit));\n\
+    vec3 exitNormal = vec3(isFirstExit) * directions;\n\
+    vec4 exit = vec4(exitNormal, firstExit);\n\
+\n\
+    return RayShapeIntersection(entry, exit);\n\
+}\n\
+\n\
+vec4 getStepSize(in SampleData sampleData, in Ray viewRay, in RayShapeIntersection shapeIntersection, in mat3 jacobianT, in float currentT) {\n\
+    vec3 gradient = viewRay.dir * jacobianT;\n\
+    vec3 sampleSizeAlongRay = getSampleSize(sampleData.tileCoords.w) / gradient;\n\
+\n\
+    RayShapeIntersection voxelIntersection = getVoxelIntersection(sampleData.tileUv, sampleSizeAlongRay);\n\
+\n\
+    // Transform normal from shape space to Cartesian space to eye space\n\
+    vec3 voxelNormal = jacobianT * voxelIntersection.entry.xyz;\n\
+    voxelNormal = normalize(czm_normal * voxelNormal);\n\
+\n\
+    // Compare with the shape intersection, to choose the appropriate normal\n\
+    vec4 voxelEntry = vec4(voxelNormal, currentT + voxelIntersection.entry.w);\n\
+    vec4 entry = intersectionMax(shapeIntersection.entry, voxelEntry);\n\
+\n\
+    float fixedStep = minComponent(abs(sampleSizeAlongRay)) * u_stepSize;\n\
+    float shift = fixedStep * SHIFT_FRACTION;\n\
+    float dt = voxelIntersection.exit.w + shift;\n\
+    if ((currentT + dt) > shapeIntersection.exit.w) {\n\
+        // Stop at end of shape\n\
+        dt = shapeIntersection.exit.w - currentT + shift;\n\
+    }\n\
+    float stepSize = clamp(dt, fixedStep * MINIMUM_STEP_SCALAR, fixedStep + shift);\n\
+\n\
+    return vec4(entry.xyz, stepSize);\n\
+}\n\
+\n\
+vec2 packIntToVec2(int value) {\n\
+    float shifted = float(value) / 255.0;\n\
+    float lowBits = fract(shifted);\n\
+    float highBits = floor(shifted) / 255.0;\n\
+    return vec2(highBits, lowBits);\n\
+}\n\
+\n\
+vec2 packFloatToVec2(float value) {\n\
+    float lowBits = fract(value);\n\
+    float highBits = floor(value) / 255.0;\n\
+    return vec2(highBits, lowBits);\n\
+}\n\
+\n\
+int getSampleIndex(in SampleData sampleData) {\n\
+    // tileUv = 1.0 is a valid coordinate but sampleIndex = u_inputDimensions is not.\n\
+    // (tileUv = 1.0 corresponds to the far edge of the last sample, at index = u_inputDimensions - 1).\n\
+    // Clamp to [0, voxelDimensions - 0.5) to avoid numerical error before flooring\n\
+    vec3 maxCoordinate = vec3(u_inputDimensions) - vec3(0.5);\n\
+    vec3 inputCoordinate = clamp(sampleData.inputCoordinate, vec3(0.0), maxCoordinate);\n\
+    ivec3 sampleIndex = ivec3(floor(inputCoordinate));\n\
+    // Convert to a 1D index for lookup in a 1D data array\n\
+    return sampleIndex.x + u_inputDimensions.x * (sampleIndex.y + u_inputDimensions.y * sampleIndex.z);\n\
+}\n\
+\n\
+/**\n\
+ * Compute the view ray at the current fragment, in the local coordinates of the shape.\n\
+ */\n\
+Ray getViewRayLocal() {\n\
+    vec4 eyeCoordinates = czm_windowToEyeCoordinates(gl_FragCoord);\n\
+    vec3 origin;\n\
+    vec3 direction;\n\
+    if (czm_orthographicIn3D == 1.0) {\n\
+        eyeCoordinates.z = 0.0;\n\
+        origin = (u_transformPositionViewToLocal * eyeCoordinates).xyz;\n\
+        direction = u_cameraDirectionLocal;\n\
+    } else {\n\
+        origin = u_cameraPositionLocal;\n\
+        direction = u_transformDirectionViewToLocal * normalize(eyeCoordinates.xyz);\n\
+    }\n\
+    return Ray(origin, direction);\n\
+}\n\
+\n\
+Ray getViewRayEC() {\n\
+    vec4 eyeCoordinates = czm_windowToEyeCoordinates(gl_FragCoord);\n\
+    vec3 viewPosEC = (czm_orthographicIn3D == 1.0)\n\
+        ? vec3(eyeCoordinates.xy, 0.0)\n\
+        : vec3(0.0);\n\
+    vec3 viewDirEC = normalize(eyeCoordinates.xyz);\n\
+    return Ray(viewPosEC, viewDirEC);\n\
+}\n\
+\n\
+void main()\n\
+{\n\
+    Ray viewRayLocal = getViewRayLocal();\n\
+    Ray viewRayEC = getViewRayEC();\n\
+\n\
+    Intersections ix;\n\
+    vec2 screenCoord = (gl_FragCoord.xy - czm_viewport.xy) / czm_viewport.zw; // [0,1]\n\
+    RayShapeIntersection shapeIntersection = intersectScene(screenCoord, viewRayLocal, viewRayEC, ix);\n\
+    // Exit early if the scene was completely missed.\n\
+    if (shapeIntersection.entry.w == NO_HIT) {\n\
+        discard;\n\
+    }\n\
+\n\
+    float currentT = shapeIntersection.entry.w;\n\
+    float endT = shapeIntersection.exit.w;\n\
+\n\
+    vec3 positionEC = viewRayEC.pos + currentT * viewRayEC.dir;\n\
+    TileAndUvCoordinate tileAndUv = getTileAndUvCoordinate(positionEC);\n\
+    vec3 positionLocal = viewRayLocal.pos + currentT * viewRayLocal.dir;\n\
+    mat3 jacobianT = convertLocalToShapeSpaceDerivative(positionLocal);\n\
+\n\
+    // Traverse the tree from the start position\n\
+    TraversalData traversalData;\n\
+    SampleData sampleDatas[SAMPLE_COUNT];\n\
+    traverseOctreeFromBeginning(tileAndUv, traversalData, sampleDatas);\n\
+    vec4 step = getStepSize(sampleDatas[0], viewRayLocal, shapeIntersection, jacobianT, currentT);\n\
+\n\
+    FragmentInput fragmentInput;\n\
+    #if defined(STATISTICS)\n\
+        setStatistics(fragmentInput.metadataStatistics);\n\
+    #endif\n\
+\n\
+    czm_modelMaterial materialOutput;\n\
+    vec4 colorAccum = vec4(0.0);\n\
+\n\
+    for (int stepCount = 0; stepCount < STEP_COUNT_MAX; ++stepCount) {\n\
+        // Read properties from the megatexture based on the traversal state\n\
+        Properties properties = accumulatePropertiesFromMegatexture(sampleDatas);\n\
+\n\
+        // Prepare the custom shader inputs\n\
+        copyPropertiesToMetadata(properties, fragmentInput.metadata);\n\
+\n\
+        fragmentInput.attributes.positionEC = positionEC;\n\
+        // Re-normalize normals: some shape intersections may have been scaled to encode positive/negative shapes\n\
+        fragmentInput.attributes.normalEC = normalize(step.xyz);\n\
+\n\
+        fragmentInput.voxel.viewDirUv = viewRayLocal.dir;\n\
+\n\
+        fragmentInput.voxel.travelDistance = step.w;\n\
+        fragmentInput.voxel.stepCount = stepCount;\n\
+        fragmentInput.voxel.tileIndex = sampleDatas[0].megatextureIndex;\n\
+        fragmentInput.voxel.sampleIndex = getSampleIndex(sampleDatas[0]);\n\
+        fragmentInput.voxel.distanceToDepthBuffer = ix.distanceToDepthBuffer - currentT;\n\
+\n\
+        // Run the custom shader\n\
+        fragmentMain(fragmentInput, materialOutput);\n\
+\n\
+        // Sanitize the custom shader output\n\
+        vec4 color = vec4(materialOutput.diffuse, materialOutput.alpha);\n\
+        color.rgb = max(color.rgb, vec3(0.0));\n\
+        color.a = clamp(color.a, 0.0, 1.0);\n\
+\n\
+        // Pre-multiplied alpha blend\n\
+        colorAccum += (1.0 - colorAccum.a) * vec4(color.rgb * color.a, color.a);\n\
+\n\
+        // Stop traversing if the alpha has been fully saturated\n\
+        if (colorAccum.a > ALPHA_ACCUM_MAX) {\n\
+            colorAccum.a = ALPHA_ACCUM_MAX;\n\
+            break;\n\
+        }\n\
+\n\
+        if (step.w == 0.0) {\n\
+            // Shape is infinitely thin. The ray may have hit the edge of a\n\
+            // foreground voxel. Step ahead slightly to check for more voxels\n\
+            step.w = 0.001;\n\
+        }\n\
+\n\
+        // Keep raymarching\n\
+        currentT += step.w;\n\
+        // Check if there's more intersections.\n\
+        if (currentT > endT) {\n\
+            #if (INTERSECTION_COUNT == 1)\n\
+                break;\n\
+            #else\n\
+                shapeIntersection = nextIntersection(ix);\n\
+                if (shapeIntersection.entry.w == NO_HIT) {\n\
+                    break;\n\
+                } else {\n\
+                    // Found another intersection. Resume raymarching there\n\
+                    currentT = shapeIntersection.entry.w;\n\
+                    endT = shapeIntersection.exit.w;\n\
+                }\n\
+            #endif\n\
+        }\n\
+        positionEC = viewRayEC.pos + currentT * viewRayEC.dir;\n\
+        tileAndUv = getTileAndUvCoordinate(positionEC);\n\
+        positionLocal = viewRayLocal.pos + currentT * viewRayLocal.dir;\n\
+        jacobianT = convertLocalToShapeSpaceDerivative(positionLocal);\n\
+\n\
+        // Traverse the tree from the current ray position.\n\
+        // This is similar to traverseOctreeFromBeginning but is faster when the ray is in the same tile as the previous step.\n\
+        traverseOctreeFromExisting(tileAndUv, traversalData, sampleDatas);\n\
+        step = getStepSize(sampleDatas[0], viewRayLocal, shapeIntersection, jacobianT, currentT);\n\
+    }\n\
+\n\
+    // Convert the alpha from [0,ALPHA_ACCUM_MAX] to [0,1]\n\
+    colorAccum.a /= ALPHA_ACCUM_MAX;\n\
+\n\
+    #if defined(PICKING)\n\
+        // If alpha is 0.0 there is nothing to pick\n\
+        if (colorAccum.a == 0.0) {\n\
+            discard;\n\
+        }\n\
+        out_FragColor = u_pickColor;\n\
+    #elif defined(PICKING_VOXEL)\n\
+        // If alpha is 0.0 there is nothing to pick\n\
+        if (colorAccum.a == 0.0) {\n\
+            discard;\n\
+        }\n\
+        vec2 megatextureId = packIntToVec2(sampleDatas[0].megatextureIndex);\n\
+        vec2 sampleIndex = packIntToVec2(getSampleIndex(sampleDatas[0]));\n\
+        out_FragColor = vec4(megatextureId, sampleIndex);\n\
+    #else\n\
+        out_FragColor = colorAccum;\n\
+    #endif\n\
+}\n\
+";

node_modules/@cesium/engine/Source/Shaders/Voxels/VoxelUtils.glsl

@@ -0,0 +1,12 @@
+struct Ray {
+    vec3 pos;
+    vec3 dir;
+};
+
+float minComponent(in vec3 v) {
+    return min(min(v.x, v.y), v.z);
+}
+
+float maxComponent(in vec3 v) {
+    return max(max(v.x, v.y), v.z);
+}

node_modules/@cesium/engine/Source/Shaders/Voxels/VoxelUtils.js

@@ -0,0 +1,14 @@
+//This file is automatically rebuilt by the Cesium build process.
+export default "struct Ray {\n\
+    vec3 pos;\n\
+    vec3 dir;\n\
+};\n\
+\n\
+float minComponent(in vec3 v) {\n\
+    return min(min(v.x, v.y), v.z);\n\
+}\n\
+\n\
+float maxComponent(in vec3 v) {\n\
+    return max(max(v.x, v.y), v.z);\n\
+}\n\
+";

node_modules/@cesium/engine/Source/Shaders/Voxels/VoxelVS.glsl

@@ -0,0 +1,11 @@
+in vec2 position;
+
+uniform vec4 u_ndcSpaceAxisAlignedBoundingBox;
+
+void main() {
+    vec2 aabbMin = u_ndcSpaceAxisAlignedBoundingBox.xy;
+    vec2 aabbMax = u_ndcSpaceAxisAlignedBoundingBox.zw;
+    vec2 translation = 0.5 * (aabbMax + aabbMin);
+    vec2 scale = 0.5 * (aabbMax - aabbMin);
+    gl_Position = vec4(position * scale + translation, 0.0, 1.0);
+}

node_modules/@cesium/engine/Source/Shaders/Voxels/VoxelVS.js

@@ -0,0 +1,13 @@
+//This file is automatically rebuilt by the Cesium build process.
+export default "in vec2 position;\n\
+\n\
+uniform vec4 u_ndcSpaceAxisAlignedBoundingBox;\n\
+\n\
+void main() {\n\
+    vec2 aabbMin = u_ndcSpaceAxisAlignedBoundingBox.xy;\n\
+    vec2 aabbMax = u_ndcSpaceAxisAlignedBoundingBox.zw;\n\
+    vec2 translation = 0.5 * (aabbMax + aabbMin);\n\
+    vec2 scale = 0.5 * (aabbMax - aabbMin);\n\
+    gl_Position = vec4(position * scale + translation, 0.0, 1.0);\n\
+}\n\
+";

node_modules/@cesium/engine/Source/Shaders/Voxels/convertLocalToBoxUv.glsl

@@ -0,0 +1,29 @@
+uniform vec3 u_boxLocalToShapeUvScale;
+
+uniform ivec4 u_cameraTileCoordinates;
+uniform vec3 u_cameraTileUv;
+uniform mat3 u_boxEcToXyz;
+
+mat3 convertLocalToShapeSpaceDerivative(in vec3 positionLocal) {
+    // For BOX, local space = shape space, so the Jacobian is the identity matrix.
+    return mat3(1.0);
+}
+
+vec3 scaleShapeUvToShapeSpace(in vec3 shapeUv) {
+    return shapeUv / u_boxLocalToShapeUvScale;
+}
+
+vec3 convertEcToDeltaTile(in vec3 positionEC) {
+    vec3 dPosition = u_boxEcToXyz * positionEC;
+    return u_boxLocalToShapeUvScale * dPosition * float(1 << u_cameraTileCoordinates.w);
+}
+
+TileAndUvCoordinate getTileAndUvCoordinate(in vec3 positionEC) {
+    vec3 deltaTileCoordinate = convertEcToDeltaTile(positionEC);
+    vec3 tileUvSum = u_cameraTileUv + deltaTileCoordinate;
+    ivec3 tileCoordinate = u_cameraTileCoordinates.xyz + ivec3(floor(tileUvSum));
+    tileCoordinate = min(max(ivec3(0), tileCoordinate), ivec3((1 << u_cameraTileCoordinates.w) - 1));
+    ivec3 tileCoordinateChange = tileCoordinate - u_cameraTileCoordinates.xyz;
+    vec3 tileUv = clamp(tileUvSum - vec3(tileCoordinateChange), 0.0, 1.0);
+    return TileAndUvCoordinate(ivec4(tileCoordinate, u_cameraTileCoordinates.w), tileUv);
+}

node_modules/@cesium/engine/Source/Shaders/Voxels/convertLocalToBoxUv.js

@@ -0,0 +1,31 @@
+//This file is automatically rebuilt by the Cesium build process.
+export default "uniform vec3 u_boxLocalToShapeUvScale;\n\
+\n\
+uniform ivec4 u_cameraTileCoordinates;\n\
+uniform vec3 u_cameraTileUv;\n\
+uniform mat3 u_boxEcToXyz;\n\
+\n\
+mat3 convertLocalToShapeSpaceDerivative(in vec3 positionLocal) {\n\
+    // For BOX, local space = shape space, so the Jacobian is the identity matrix.\n\
+    return mat3(1.0);\n\
+}\n\
+\n\
+vec3 scaleShapeUvToShapeSpace(in vec3 shapeUv) {\n\
+    return shapeUv / u_boxLocalToShapeUvScale;\n\
+}\n\
+\n\
+vec3 convertEcToDeltaTile(in vec3 positionEC) {\n\
+    vec3 dPosition = u_boxEcToXyz * positionEC;\n\
+    return u_boxLocalToShapeUvScale * dPosition * float(1 << u_cameraTileCoordinates.w);\n\
+}\n\
+\n\
+TileAndUvCoordinate getTileAndUvCoordinate(in vec3 positionEC) {\n\
+    vec3 deltaTileCoordinate = convertEcToDeltaTile(positionEC);\n\
+    vec3 tileUvSum = u_cameraTileUv + deltaTileCoordinate;\n\
+    ivec3 tileCoordinate = u_cameraTileCoordinates.xyz + ivec3(floor(tileUvSum));\n\
+    tileCoordinate = min(max(ivec3(0), tileCoordinate), ivec3((1 << u_cameraTileCoordinates.w) - 1));\n\
+    ivec3 tileCoordinateChange = tileCoordinate - u_cameraTileCoordinates.xyz;\n\
+    vec3 tileUv = clamp(tileUvSum - vec3(tileCoordinateChange), 0.0, 1.0);\n\
+    return TileAndUvCoordinate(ivec4(tileCoordinate, u_cameraTileCoordinates.w), tileUv);\n\
+}\n\
+";

node_modules/@cesium/engine/Source/Shaders/Voxels/convertLocalToCylinderUv.glsl

@@ -0,0 +1,95 @@
+uniform vec3 u_cylinderLocalToShapeUvScale; // x = radius scale, y = angle scale, z = height scale
+uniform float u_cylinderShapeUvAngleRangeOrigin;
+uniform mat3 u_cylinderEcToRadialTangentUp;
+uniform ivec4 u_cameraTileCoordinates;
+uniform vec3 u_cameraTileUv;
+uniform vec3 u_cameraShapePosition; // (radial distance, angle, height) of camera in shape space
+
+mat3 convertLocalToShapeSpaceDerivative(in vec3 position) {
+    vec3 radial = normalize(vec3(position.xy, 0.0));
+    vec3 z = vec3(0.0, 0.0, 1.0);
+    vec3 east = normalize(vec3(-position.y, position.x, 0.0));
+    return mat3(radial, east / length(position.xy), z);
+}
+
+vec3 scaleShapeUvToShapeSpace(in vec3 shapeUv) {
+    float radius = shapeUv.x / u_cylinderLocalToShapeUvScale.x;
+    float angle = shapeUv.y * czm_twoPi / u_cylinderLocalToShapeUvScale.y;
+    float height = shapeUv.z / u_cylinderLocalToShapeUvScale.z;
+
+    return vec3(radius, angle, height);
+}
+
+/**
+ * Computes the change in polar coordinates given a change in position.
+ * @param {vec2} dPosition The change in position in Cartesian coordinates.
+ * @param {float} cameraRadialDistance The radial distance of the camera from the origin.
+ * @return {vec2} The change in polar coordinates (radial distance, angle).
+ */
+vec2 computePolarChange(in vec2 dPosition, in float cameraRadialDistance) {
+    float dAngle = atan(dPosition.y, cameraRadialDistance + dPosition.x);
+    // Find the direction of the radial axis at the output angle, in Cartesian coordinates
+    vec2 outputRadialAxis = vec2(cos(dAngle), sin(dAngle));
+    float sinHalfAngle = sin(dAngle / 2.0);
+    float versine = 2.0 * sinHalfAngle * sinHalfAngle;
+    float dRadial = dot(dPosition, outputRadialAxis) - cameraRadialDistance * versine;
+    return vec2(dRadial, dAngle);
+}
+
+vec3 convertEcToDeltaShape(in vec3 positionEC) {
+    // 1. Rotate to radial, tangent, and up coordinates
+    vec3 rtu = u_cylinderEcToRadialTangentUp * positionEC;
+    // 2. Compute change in angular and radial coordinates.
+    vec2 dPolar = computePolarChange(rtu.xy, u_cameraShapePosition.x);
+    return vec3(dPolar.xy, rtu.z);
+}
+
+vec3 convertEcToDeltaTile(in vec3 positionEC) {
+    vec3 deltaShape = convertEcToDeltaShape(positionEC);
+    // Convert to tileset coordinates in [0, 1]
+    float dx = u_cylinderLocalToShapeUvScale.x * deltaShape.x;
+    float dy = deltaShape.y / czm_twoPi;
+#if defined(CYLINDER_HAS_SHAPE_BOUNDS_ANGLE)
+    // Wrap to ensure dy is not crossing through the unoccupied angle range, where
+    // angle to tile coordinate conversions would be more complicated
+    float cameraUvAngle = (u_cameraShapePosition.y + czm_pi) / czm_twoPi;
+    float cameraUvAngleShift = fract(cameraUvAngle - u_cylinderShapeUvAngleRangeOrigin);
+    float rawOutputUvAngle = cameraUvAngleShift + dy;
+    float rotation = floor(rawOutputUvAngle);
+    dy -= rotation;
+#endif
+    dy *= u_cylinderLocalToShapeUvScale.y;
+    float dz = u_cylinderLocalToShapeUvScale.z * deltaShape.z;
+    // Convert to tile coordinate changes
+    return vec3(dx, dy, dz) * float(1 << u_cameraTileCoordinates.w);
+}
+
+TileAndUvCoordinate getTileAndUvCoordinate(in vec3 positionEC) {
+    vec3 deltaTileCoordinate = convertEcToDeltaTile(positionEC);
+    vec3 tileUvSum = u_cameraTileUv + deltaTileCoordinate;
+    ivec3 tileCoordinate = u_cameraTileCoordinates.xyz + ivec3(floor(tileUvSum));
+    int maxTileCoordinate = (1 << u_cameraTileCoordinates.w) - 1;
+    tileCoordinate.x = min(max(0, tileCoordinate.x), maxTileCoordinate);
+    tileCoordinate.z = min(max(0, tileCoordinate.z), maxTileCoordinate);
+#if (!defined(CYLINDER_HAS_SHAPE_BOUNDS_ANGLE))
+    ivec3 tileCoordinateChange = tileCoordinate - u_cameraTileCoordinates.xyz;
+    if (tileCoordinate.y < 0) {
+        tileCoordinate.y += (maxTileCoordinate + 1);
+    } else if (tileCoordinate.y > maxTileCoordinate) {
+        tileCoordinate.y -= (maxTileCoordinate + 1);
+    }
+#else
+    tileCoordinate.y = min(max(0, tileCoordinate.y), maxTileCoordinate);
+    ivec3 tileCoordinateChange = tileCoordinate - u_cameraTileCoordinates.xyz;
+#endif
+    vec3 tileUv = tileUvSum - vec3(tileCoordinateChange);
+    tileUv.x = clamp(tileUv.x, 0.0, 1.0);
+#if (!defined(CYLINDER_HAS_SHAPE_BOUNDS_ANGLE))
+    // If there is only one tile spanning 2*PI angle, the coordinate wraps around
+    tileUv.y = (u_cameraTileCoordinates.w == 0) ? fract(tileUv.y) : clamp(tileUv.y, 0.0, 1.0);
+#else
+    tileUv.y = clamp(tileUv.y, 0.0, 1.0);
+#endif
+    tileUv.z = clamp(tileUv.z, 0.0, 1.0);
+    return TileAndUvCoordinate(ivec4(tileCoordinate, u_cameraTileCoordinates.w), tileUv);
+}

node_modules/@cesium/engine/Source/Shaders/Voxels/convertLocalToCylinderUv.js

@@ -0,0 +1,97 @@
+//This file is automatically rebuilt by the Cesium build process.
+export default "uniform vec3 u_cylinderLocalToShapeUvScale; // x = radius scale, y = angle scale, z = height scale\n\
+uniform float u_cylinderShapeUvAngleRangeOrigin;\n\
+uniform mat3 u_cylinderEcToRadialTangentUp;\n\
+uniform ivec4 u_cameraTileCoordinates;\n\
+uniform vec3 u_cameraTileUv;\n\
+uniform vec3 u_cameraShapePosition; // (radial distance, angle, height) of camera in shape space\n\
+\n\
+mat3 convertLocalToShapeSpaceDerivative(in vec3 position) {\n\
+    vec3 radial = normalize(vec3(position.xy, 0.0));\n\
+    vec3 z = vec3(0.0, 0.0, 1.0);\n\
+    vec3 east = normalize(vec3(-position.y, position.x, 0.0));\n\
+    return mat3(radial, east / length(position.xy), z);\n\
+}\n\
+\n\
+vec3 scaleShapeUvToShapeSpace(in vec3 shapeUv) {\n\
+    float radius = shapeUv.x / u_cylinderLocalToShapeUvScale.x;\n\
+    float angle = shapeUv.y * czm_twoPi / u_cylinderLocalToShapeUvScale.y;\n\
+    float height = shapeUv.z / u_cylinderLocalToShapeUvScale.z;\n\
+\n\
+    return vec3(radius, angle, height);\n\
+}\n\
+\n\
+/**\n\
+ * Computes the change in polar coordinates given a change in position.\n\
+ * @param {vec2} dPosition The change in position in Cartesian coordinates.\n\
+ * @param {float} cameraRadialDistance The radial distance of the camera from the origin.\n\
+ * @return {vec2} The change in polar coordinates (radial distance, angle).\n\
+ */\n\
+vec2 computePolarChange(in vec2 dPosition, in float cameraRadialDistance) {\n\
+    float dAngle = atan(dPosition.y, cameraRadialDistance + dPosition.x);\n\
+    // Find the direction of the radial axis at the output angle, in Cartesian coordinates\n\
+    vec2 outputRadialAxis = vec2(cos(dAngle), sin(dAngle));\n\
+    float sinHalfAngle = sin(dAngle / 2.0);\n\
+    float versine = 2.0 * sinHalfAngle * sinHalfAngle;\n\
+    float dRadial = dot(dPosition, outputRadialAxis) - cameraRadialDistance * versine;\n\
+    return vec2(dRadial, dAngle);\n\
+}\n\
+\n\
+vec3 convertEcToDeltaShape(in vec3 positionEC) {\n\
+    // 1. Rotate to radial, tangent, and up coordinates\n\
+    vec3 rtu = u_cylinderEcToRadialTangentUp * positionEC;\n\
+    // 2. Compute change in angular and radial coordinates.\n\
+    vec2 dPolar = computePolarChange(rtu.xy, u_cameraShapePosition.x);\n\
+    return vec3(dPolar.xy, rtu.z);\n\
+}\n\
+\n\
+vec3 convertEcToDeltaTile(in vec3 positionEC) {\n\
+    vec3 deltaShape = convertEcToDeltaShape(positionEC);\n\
+    // Convert to tileset coordinates in [0, 1]\n\
+    float dx = u_cylinderLocalToShapeUvScale.x * deltaShape.x;\n\
+    float dy = deltaShape.y / czm_twoPi;\n\
+#if defined(CYLINDER_HAS_SHAPE_BOUNDS_ANGLE)\n\
+    // Wrap to ensure dy is not crossing through the unoccupied angle range, where\n\
+    // angle to tile coordinate conversions would be more complicated\n\
+    float cameraUvAngle = (u_cameraShapePosition.y + czm_pi) / czm_twoPi;\n\
+    float cameraUvAngleShift = fract(cameraUvAngle - u_cylinderShapeUvAngleRangeOrigin);\n\
+    float rawOutputUvAngle = cameraUvAngleShift + dy;\n\
+    float rotation = floor(rawOutputUvAngle);\n\
+    dy -= rotation;\n\
+#endif\n\
+    dy *= u_cylinderLocalToShapeUvScale.y;\n\
+    float dz = u_cylinderLocalToShapeUvScale.z * deltaShape.z;\n\
+    // Convert to tile coordinate changes\n\
+    return vec3(dx, dy, dz) * float(1 << u_cameraTileCoordinates.w);\n\
+}\n\
+\n\
+TileAndUvCoordinate getTileAndUvCoordinate(in vec3 positionEC) {\n\
+    vec3 deltaTileCoordinate = convertEcToDeltaTile(positionEC);\n\
+    vec3 tileUvSum = u_cameraTileUv + deltaTileCoordinate;\n\
+    ivec3 tileCoordinate = u_cameraTileCoordinates.xyz + ivec3(floor(tileUvSum));\n\
+    int maxTileCoordinate = (1 << u_cameraTileCoordinates.w) - 1;\n\
+    tileCoordinate.x = min(max(0, tileCoordinate.x), maxTileCoordinate);\n\
+    tileCoordinate.z = min(max(0, tileCoordinate.z), maxTileCoordinate);\n\
+#if (!defined(CYLINDER_HAS_SHAPE_BOUNDS_ANGLE))\n\
+    ivec3 tileCoordinateChange = tileCoordinate - u_cameraTileCoordinates.xyz;\n\
+    if (tileCoordinate.y < 0) {\n\
+        tileCoordinate.y += (maxTileCoordinate + 1);\n\
+    } else if (tileCoordinate.y > maxTileCoordinate) {\n\
+        tileCoordinate.y -= (maxTileCoordinate + 1);\n\
+    }\n\
+#else\n\
+    tileCoordinate.y = min(max(0, tileCoordinate.y), maxTileCoordinate);\n\
+    ivec3 tileCoordinateChange = tileCoordinate - u_cameraTileCoordinates.xyz;\n\
+#endif\n\
+    vec3 tileUv = tileUvSum - vec3(tileCoordinateChange);\n\
+    tileUv.x = clamp(tileUv.x, 0.0, 1.0);\n\
+#if (!defined(CYLINDER_HAS_SHAPE_BOUNDS_ANGLE))\n\
+    // If there is only one tile spanning 2*PI angle, the coordinate wraps around\n\
+    tileUv.y = (u_cameraTileCoordinates.w == 0) ? fract(tileUv.y) : clamp(tileUv.y, 0.0, 1.0);\n\
+#else\n\
+    tileUv.y = clamp(tileUv.y, 0.0, 1.0);\n\
+#endif\n\
+    tileUv.z = clamp(tileUv.z, 0.0, 1.0);\n\
+    return TileAndUvCoordinate(ivec4(tileCoordinate, u_cameraTileCoordinates.w), tileUv);\n\
+}\n\
+";

node_modules/@cesium/engine/Source/Shaders/Voxels/convertLocalToEllipsoidUv.glsl

@@ -0,0 +1,183 @@
+/* Ellipsoid defines (set in Scene/VoxelEllipsoidShape.js)
+#define ELLIPSOID_HAS_SHAPE_BOUNDS_LONGITUDE
+#define ELLIPSOID_HAS_SHAPE_BOUNDS_LATITUDE
+*/
+
+uniform vec3 u_cameraPositionCartographic; // (longitude, latitude, height) in radians and meters
+uniform vec2 u_ellipsoidCurvatureAtLatitude;
+uniform mat3 u_ellipsoidEcToEastNorthUp;
+uniform vec3 u_ellipsoidRadii;
+uniform vec2 u_evoluteScale; // (radii.x ^ 2 - radii.z ^ 2) * vec2(1.0, -1.0) / radii;
+uniform vec3 u_ellipsoidInverseRadiiSquared;
+#if defined(ELLIPSOID_HAS_SHAPE_BOUNDS_LONGITUDE)
+    uniform float u_ellipsoidShapeUvLongitudeRangeOrigin;
+#endif
+uniform vec3 u_ellipsoidLocalToShapeUvScale; // x = longitude scale, y = latitude scale, z = height scale
+
+uniform ivec4 u_cameraTileCoordinates;
+uniform vec3 u_cameraTileUv;
+
+// robust iterative solution without trig functions
+// https://github.com/0xfaded/ellipse_demo/issues/1
+// https://stackoverflow.com/questions/22959698/distance-from-given-point-to-given-ellipse
+// Extended to return radius of curvature along with the point
+vec3 nearestPointAndRadiusOnEllipse(vec2 pos, vec2 radii) {
+    vec2 p = abs(pos);
+    vec2 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
+    vec2 tTrigs = vec2(0.7071067811865476);
+    // Initial guess of point on ellipsoid
+    vec2 v = radii * tTrigs;
+    // Center of curvature of the ellipse at v
+    vec2 evolute = u_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.
+        vec2 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 = u_evoluteScale * tTrigs * tTrigs * tTrigs;
+    }
+
+    return vec3(v * sign(pos), length(v - evolute));
+}
+
+mat3 convertLocalToShapeSpaceDerivative(in vec3 position) {
+    vec3 east = normalize(vec3(-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);
+    vec2 posEllipse = vec2(distanceFromZAxis, position.z);
+    vec3 surfacePointAndRadius = nearestPointAndRadiusOnEllipse(posEllipse, u_ellipsoidRadii.xz);
+    vec2 surfacePoint = surfacePointAndRadius.xy;
+
+    vec2 normal2d = normalize(surfacePoint * u_ellipsoidInverseRadiiSquared.xz);
+    vec3 north = vec3(-normal2d.y * normalize(position.xy), abs(normal2d.x));
+
+    float heightSign = length(posEllipse) < length(surfacePoint) ? -1.0 : 1.0;
+    float height = heightSign * length(posEllipse - surfacePoint);
+    vec3 up = normalize(cross(east, north));
+
+    return mat3(east / distanceFromZAxis, north / (surfacePointAndRadius.z + height), up);
+}
+
+vec3 scaleShapeUvToShapeSpace(in vec3 shapeUv) {
+    // Convert from [0, 1] to radians [-pi, pi]
+    float longitude = shapeUv.x * czm_twoPi;
+    #if defined (ELLIPSOID_HAS_SHAPE_BOUNDS_LONGITUDE)
+        longitude /= u_ellipsoidLocalToShapeUvScale.x;
+    #endif
+
+    // Convert from [0, 1] to radians [-pi/2, pi/2]
+    float latitude = shapeUv.y * czm_pi;
+    #if defined(ELLIPSOID_HAS_SHAPE_BOUNDS_LATITUDE)
+        latitude /= u_ellipsoidLocalToShapeUvScale.y;
+    #endif
+
+    float height = shapeUv.z / u_ellipsoidLocalToShapeUvScale.z;
+
+    return vec3(longitude, latitude, height);
+}
+
+vec3 convertEcToDeltaShape(in vec3 positionEC) {
+    vec3 enu = u_ellipsoidEcToEastNorthUp * positionEC;
+
+    // 1. Compute the change in longitude from the camera to the ENU point
+    // First project the camera and ENU positions to the equatorial XY plane,
+    // positioning the camera on the +x axis, so that enu.x projects along the +y axis
+    float cosLatitude = cos(u_cameraPositionCartographic.y);
+    float sinLatitude = sin(u_cameraPositionCartographic.y);
+    float primeVerticalRadius = 1.0 / u_ellipsoidCurvatureAtLatitude.x;
+    vec2 cameraXY = vec2((primeVerticalRadius + u_cameraPositionCartographic.z) * cosLatitude, 0.0);
+    // Note precision loss in positionXY.x if length(enu) << length(cameraXY)
+    vec2 positionXY = cameraXY + vec2(-enu.y * sinLatitude + enu.z * cosLatitude, enu.x);
+    float dLongitude = atan(positionXY.y, positionXY.x);
+
+    // 2. Find the longitude component of positionXY, by rotating about Z until the y component is zero.
+    // Use the versine  to compute the change in x directly from the change in angle:
+    //   versine(angle) = 2 * sin^2(angle/2)
+    float sinHalfLongitude = sin(dLongitude / 2.0);
+    float dx = length(positionXY) * 2.0 * sinHalfLongitude * sinHalfLongitude;
+    // Rotate longitude component back to ENU North and Up, and remove from enu
+    enu += vec3(-enu.x, -dx * sinLatitude, dx * cosLatitude);
+
+    // 3. Compute the change in latitude from the camera to the ENU point.
+    // First project the camera and ENU positions to the meridional ZX plane,
+    // positioning the camera on the +Z axis, so that enu.y maps to the +X axis.
+    float meridionalRadius = 1.0 / u_ellipsoidCurvatureAtLatitude.y;
+    vec2 cameraZX = vec2(meridionalRadius + u_cameraPositionCartographic.z, 0.0);
+    vec2 positionZX = cameraZX + vec2(enu.z, enu.y);
+    float dLatitude = atan(positionZX.y, positionZX.x);
+
+    // 4. Compute the change in height above the ellipsoid
+    // Find the change in enu.z associated with rotating the point to the latitude of the camera
+    float sinHalfLatitude = sin(dLatitude / 2.0);
+    float dz = length(positionZX) * 2.0 * sinHalfLatitude * sinHalfLatitude;
+    // The remaining change in enu.z is the change in height above the ellipsoid
+    float dHeight = enu.z + dz;
+
+    return vec3(dLongitude, dLatitude, dHeight);
+}
+
+vec3 convertEcToDeltaTile(in vec3 positionEC) {
+    vec3 deltaShape = convertEcToDeltaShape(positionEC);
+    // Convert to tileset coordinates in [0, 1]
+    float dx = deltaShape.x / czm_twoPi;
+
+#if (defined(ELLIPSOID_HAS_SHAPE_BOUNDS_LONGITUDE))
+    // Wrap to ensure dx is not crossing through the unoccupied angle range, where
+    // angle to tile coordinate conversions would be more complicated
+    float cameraUvLongitude = (u_cameraPositionCartographic.x + czm_pi) / czm_twoPi;
+    float cameraUvLongitudeShift = fract(cameraUvLongitude - u_ellipsoidShapeUvLongitudeRangeOrigin);
+    float rawOutputUvLongitude = cameraUvLongitudeShift + dx;
+    float rotation = floor(rawOutputUvLongitude);
+    dx -= rotation;
+    dx *= u_ellipsoidLocalToShapeUvScale.x;
+#endif
+
+    float dy = deltaShape.y / czm_pi;
+#if (defined(ELLIPSOID_HAS_SHAPE_BOUNDS_LATITUDE))
+    dy *= u_ellipsoidLocalToShapeUvScale.y;
+#endif
+
+    float dz = u_ellipsoidLocalToShapeUvScale.z * deltaShape.z;
+    // Convert to tile coordinate changes
+    return vec3(dx, dy, dz) * float(1 << u_cameraTileCoordinates.w);
+}
+
+TileAndUvCoordinate getTileAndUvCoordinate(in vec3 positionEC) {
+    vec3 deltaTileCoordinate = convertEcToDeltaTile(positionEC);
+    vec3 tileUvSum = u_cameraTileUv + deltaTileCoordinate;
+    ivec3 tileCoordinate = u_cameraTileCoordinates.xyz + ivec3(floor(tileUvSum));
+    int maxTileCoordinate = (1 << u_cameraTileCoordinates.w) - 1;
+    tileCoordinate.y = min(max(0, tileCoordinate.y), maxTileCoordinate);
+    tileCoordinate.z = min(max(0, tileCoordinate.z), maxTileCoordinate);
+#if (!defined(ELLIPSOID_HAS_SHAPE_BOUNDS_LONGITUDE))
+    ivec3 tileCoordinateChange = tileCoordinate - u_cameraTileCoordinates.xyz;
+    if (tileCoordinate.x < 0) {
+        tileCoordinate.x += (maxTileCoordinate + 1);
+    } else if (tileCoordinate.x > maxTileCoordinate) {
+        tileCoordinate.x -= (maxTileCoordinate + 1);
+    }
+#else
+    tileCoordinate.x = min(max(0, tileCoordinate.x), maxTileCoordinate);
+    ivec3 tileCoordinateChange = tileCoordinate - u_cameraTileCoordinates.xyz;
+#endif
+    vec3 tileUv = tileUvSum - vec3(tileCoordinateChange);
+#if (!defined(ELLIPSOID_HAS_SHAPE_BOUNDS_LONGITUDE))
+    // If there is only one tile spanning 2*PI angle, the coordinate wraps around
+    tileUv.x = (u_cameraTileCoordinates.w == 0) ? fract(tileUv.x) : clamp(tileUv.x, 0.0, 1.0);
+#else
+    tileUv.x = clamp(tileUv.x, 0.0, 1.0);
+#endif
+    tileUv.y = clamp(tileUv.y, 0.0, 1.0);
+    tileUv.z = clamp(tileUv.z, 0.0, 1.0);
+    return TileAndUvCoordinate(ivec4(tileCoordinate, u_cameraTileCoordinates.w), tileUv);
+}

node_modules/@cesium/engine/Source/Shaders/Voxels/convertLocalToEllipsoidUv.js

@@ -0,0 +1,185 @@
+//This file is automatically rebuilt by the Cesium build process.
+export default "/* Ellipsoid defines (set in Scene/VoxelEllipsoidShape.js)\n\
+#define ELLIPSOID_HAS_SHAPE_BOUNDS_LONGITUDE\n\
+#define ELLIPSOID_HAS_SHAPE_BOUNDS_LATITUDE\n\
+*/\n\
+\n\
+uniform vec3 u_cameraPositionCartographic; // (longitude, latitude, height) in radians and meters\n\
+uniform vec2 u_ellipsoidCurvatureAtLatitude;\n\
+uniform mat3 u_ellipsoidEcToEastNorthUp;\n\
+uniform vec3 u_ellipsoidRadii;\n\
+uniform vec2 u_evoluteScale; // (radii.x ^ 2 - radii.z ^ 2) * vec2(1.0, -1.0) / radii;\n\
+uniform vec3 u_ellipsoidInverseRadiiSquared;\n\
+#if defined(ELLIPSOID_HAS_SHAPE_BOUNDS_LONGITUDE)\n\
+    uniform float u_ellipsoidShapeUvLongitudeRangeOrigin;\n\
+#endif\n\
+uniform vec3 u_ellipsoidLocalToShapeUvScale; // x = longitude scale, y = latitude scale, z = height scale\n\
+\n\
+uniform ivec4 u_cameraTileCoordinates;\n\
+uniform vec3 u_cameraTileUv;\n\
+\n\
+// robust iterative solution without trig functions\n\
+// https://github.com/0xfaded/ellipse_demo/issues/1\n\
+// https://stackoverflow.com/questions/22959698/distance-from-given-point-to-given-ellipse\n\
+// Extended to return radius of curvature along with the point\n\
+vec3 nearestPointAndRadiusOnEllipse(vec2 pos, vec2 radii) {\n\
+    vec2 p = abs(pos);\n\
+    vec2 inverseRadii = 1.0 / radii;\n\
+\n\
+    // We describe the ellipse parametrically: v = radii * vec2(cos(t), sin(t))\n\
+    // but store the cos and sin of t in a vec2 for efficiency.\n\
+    // Initial guess: t = pi/4\n\
+    vec2 tTrigs = vec2(0.7071067811865476);\n\
+    // Initial guess of point on ellipsoid\n\
+    vec2 v = radii * tTrigs;\n\
+    // Center of curvature of the ellipse at v\n\
+    vec2 evolute = u_evoluteScale * tTrigs * tTrigs * tTrigs;\n\
+\n\
+    const int iterations = 3;\n\
+    for (int i = 0; i < iterations; ++i) {\n\
+        // Find the (approximate) intersection of p - evolute with the ellipsoid.\n\
+        vec2 q = normalize(p - evolute) * length(v - evolute);\n\
+        // Update the estimate of t.\n\
+        tTrigs = (q + evolute) * inverseRadii;\n\
+        tTrigs = normalize(clamp(tTrigs, 0.0, 1.0));\n\
+        v = radii * tTrigs;\n\
+        evolute = u_evoluteScale * tTrigs * tTrigs * tTrigs;\n\
+    }\n\
+\n\
+    return vec3(v * sign(pos), length(v - evolute));\n\
+}\n\
+\n\
+mat3 convertLocalToShapeSpaceDerivative(in vec3 position) {\n\
+    vec3 east = normalize(vec3(-position.y, position.x, 0.0));\n\
+\n\
+    // Convert the 3D position to a 2D position relative to the ellipse (radii.x, radii.z)\n\
+    // (assume radii.y == radii.x) and find the nearest point on the ellipse and its normal\n\
+    float distanceFromZAxis = length(position.xy);\n\
+    vec2 posEllipse = vec2(distanceFromZAxis, position.z);\n\
+    vec3 surfacePointAndRadius = nearestPointAndRadiusOnEllipse(posEllipse, u_ellipsoidRadii.xz);\n\
+    vec2 surfacePoint = surfacePointAndRadius.xy;\n\
+\n\
+    vec2 normal2d = normalize(surfacePoint * u_ellipsoidInverseRadiiSquared.xz);\n\
+    vec3 north = vec3(-normal2d.y * normalize(position.xy), abs(normal2d.x));\n\
+\n\
+    float heightSign = length(posEllipse) < length(surfacePoint) ? -1.0 : 1.0;\n\
+    float height = heightSign * length(posEllipse - surfacePoint);\n\
+    vec3 up = normalize(cross(east, north));\n\
+\n\
+    return mat3(east / distanceFromZAxis, north / (surfacePointAndRadius.z + height), up);\n\
+}\n\
+\n\
+vec3 scaleShapeUvToShapeSpace(in vec3 shapeUv) {\n\
+    // Convert from [0, 1] to radians [-pi, pi]\n\
+    float longitude = shapeUv.x * czm_twoPi;\n\
+    #if defined (ELLIPSOID_HAS_SHAPE_BOUNDS_LONGITUDE)\n\
+        longitude /= u_ellipsoidLocalToShapeUvScale.x;\n\
+    #endif\n\
+\n\
+    // Convert from [0, 1] to radians [-pi/2, pi/2]\n\
+    float latitude = shapeUv.y * czm_pi;\n\
+    #if defined(ELLIPSOID_HAS_SHAPE_BOUNDS_LATITUDE)\n\
+        latitude /= u_ellipsoidLocalToShapeUvScale.y;\n\
+    #endif\n\
+\n\
+    float height = shapeUv.z / u_ellipsoidLocalToShapeUvScale.z;\n\
+\n\
+    return vec3(longitude, latitude, height);\n\
+}\n\
+\n\
+vec3 convertEcToDeltaShape(in vec3 positionEC) {\n\
+    vec3 enu = u_ellipsoidEcToEastNorthUp * positionEC;\n\
+\n\
+    // 1. Compute the change in longitude from the camera to the ENU point\n\
+    // First project the camera and ENU positions to the equatorial XY plane,\n\
+    // positioning the camera on the +x axis, so that enu.x projects along the +y axis\n\
+    float cosLatitude = cos(u_cameraPositionCartographic.y);\n\
+    float sinLatitude = sin(u_cameraPositionCartographic.y);\n\
+    float primeVerticalRadius = 1.0 / u_ellipsoidCurvatureAtLatitude.x;\n\
+    vec2 cameraXY = vec2((primeVerticalRadius + u_cameraPositionCartographic.z) * cosLatitude, 0.0);\n\
+    // Note precision loss in positionXY.x if length(enu) << length(cameraXY)\n\
+    vec2 positionXY = cameraXY + vec2(-enu.y * sinLatitude + enu.z * cosLatitude, enu.x);\n\
+    float dLongitude = atan(positionXY.y, positionXY.x);\n\
+\n\
+    // 2. Find the longitude component of positionXY, by rotating about Z until the y component is zero.\n\
+    // Use the versine  to compute the change in x directly from the change in angle:\n\
+    //   versine(angle) = 2 * sin^2(angle/2)\n\
+    float sinHalfLongitude = sin(dLongitude / 2.0);\n\
+    float dx = length(positionXY) * 2.0 * sinHalfLongitude * sinHalfLongitude;\n\
+    // Rotate longitude component back to ENU North and Up, and remove from enu\n\
+    enu += vec3(-enu.x, -dx * sinLatitude, dx * cosLatitude);\n\
+\n\
+    // 3. Compute the change in latitude from the camera to the ENU point.\n\
+    // First project the camera and ENU positions to the meridional ZX plane,\n\
+    // positioning the camera on the +Z axis, so that enu.y maps to the +X axis.\n\
+    float meridionalRadius = 1.0 / u_ellipsoidCurvatureAtLatitude.y;\n\
+    vec2 cameraZX = vec2(meridionalRadius + u_cameraPositionCartographic.z, 0.0);\n\
+    vec2 positionZX = cameraZX + vec2(enu.z, enu.y);\n\
+    float dLatitude = atan(positionZX.y, positionZX.x);\n\
+\n\
+    // 4. Compute the change in height above the ellipsoid\n\
+    // Find the change in enu.z associated with rotating the point to the latitude of the camera\n\
+    float sinHalfLatitude = sin(dLatitude / 2.0);\n\
+    float dz = length(positionZX) * 2.0 * sinHalfLatitude * sinHalfLatitude;\n\
+    // The remaining change in enu.z is the change in height above the ellipsoid\n\
+    float dHeight = enu.z + dz;\n\
+\n\
+    return vec3(dLongitude, dLatitude, dHeight);\n\
+}\n\
+\n\
+vec3 convertEcToDeltaTile(in vec3 positionEC) {\n\
+    vec3 deltaShape = convertEcToDeltaShape(positionEC);\n\
+    // Convert to tileset coordinates in [0, 1]\n\
+    float dx = deltaShape.x / czm_twoPi;\n\
+\n\
+#if (defined(ELLIPSOID_HAS_SHAPE_BOUNDS_LONGITUDE))\n\
+    // Wrap to ensure dx is not crossing through the unoccupied angle range, where\n\
+    // angle to tile coordinate conversions would be more complicated\n\
+    float cameraUvLongitude = (u_cameraPositionCartographic.x + czm_pi) / czm_twoPi;\n\
+    float cameraUvLongitudeShift = fract(cameraUvLongitude - u_ellipsoidShapeUvLongitudeRangeOrigin);\n\
+    float rawOutputUvLongitude = cameraUvLongitudeShift + dx;\n\
+    float rotation = floor(rawOutputUvLongitude);\n\
+    dx -= rotation;\n\
+    dx *= u_ellipsoidLocalToShapeUvScale.x;\n\
+#endif\n\
+\n\
+    float dy = deltaShape.y / czm_pi;\n\
+#if (defined(ELLIPSOID_HAS_SHAPE_BOUNDS_LATITUDE))\n\
+    dy *= u_ellipsoidLocalToShapeUvScale.y;\n\
+#endif\n\
+\n\
+    float dz = u_ellipsoidLocalToShapeUvScale.z * deltaShape.z;\n\
+    // Convert to tile coordinate changes\n\
+    return vec3(dx, dy, dz) * float(1 << u_cameraTileCoordinates.w);\n\
+}\n\
+\n\
+TileAndUvCoordinate getTileAndUvCoordinate(in vec3 positionEC) {\n\
+    vec3 deltaTileCoordinate = convertEcToDeltaTile(positionEC);\n\
+    vec3 tileUvSum = u_cameraTileUv + deltaTileCoordinate;\n\
+    ivec3 tileCoordinate = u_cameraTileCoordinates.xyz + ivec3(floor(tileUvSum));\n\
+    int maxTileCoordinate = (1 << u_cameraTileCoordinates.w) - 1;\n\
+    tileCoordinate.y = min(max(0, tileCoordinate.y), maxTileCoordinate);\n\
+    tileCoordinate.z = min(max(0, tileCoordinate.z), maxTileCoordinate);\n\
+#if (!defined(ELLIPSOID_HAS_SHAPE_BOUNDS_LONGITUDE))\n\
+    ivec3 tileCoordinateChange = tileCoordinate - u_cameraTileCoordinates.xyz;\n\
+    if (tileCoordinate.x < 0) {\n\
+        tileCoordinate.x += (maxTileCoordinate + 1);\n\
+    } else if (tileCoordinate.x > maxTileCoordinate) {\n\
+        tileCoordinate.x -= (maxTileCoordinate + 1);\n\
+    }\n\
+#else\n\
+    tileCoordinate.x = min(max(0, tileCoordinate.x), maxTileCoordinate);\n\
+    ivec3 tileCoordinateChange = tileCoordinate - u_cameraTileCoordinates.xyz;\n\
+#endif\n\
+    vec3 tileUv = tileUvSum - vec3(tileCoordinateChange);\n\
+#if (!defined(ELLIPSOID_HAS_SHAPE_BOUNDS_LONGITUDE))\n\
+    // If there is only one tile spanning 2*PI angle, the coordinate wraps around\n\
+    tileUv.x = (u_cameraTileCoordinates.w == 0) ? fract(tileUv.x) : clamp(tileUv.x, 0.0, 1.0);\n\
+#else\n\
+    tileUv.x = clamp(tileUv.x, 0.0, 1.0);\n\
+#endif\n\
+    tileUv.y = clamp(tileUv.y, 0.0, 1.0);\n\
+    tileUv.z = clamp(tileUv.z, 0.0, 1.0);\n\
+    return TileAndUvCoordinate(ivec4(tileCoordinate, u_cameraTileCoordinates.w), tileUv);\n\
+}\n\
+";