Optimum step size and lighting for cylinder- and ellipsoid-shaped voxels

[jjhembd]

Description

This PR implements an optimum step size calculation for cylinder- and ellipsoid-shaped voxels. This generalizes the work done in #11050 for box-shaped voxels, using a new approach based on the partial derivatives of the coordinate transform from world space to shape space.

The partial derivative approach also enables a simpler method to derive the surface normal of the current voxel. This extends the functionality of #11076 to all shape types.

Background

Ray marching is performed in a Cartesian space, which has been scaled to a "UV" space where the range [0, 1] spans the bounds of the shape. At each step of the marching, the voxel shaders convert the Cartesian "UV" ray position into the UV coordinates of the shape, whether BOX (convertUvToBox.glsl), CYLINDER (convertUvToCylinder.glsl), or ELLIPSOID (convertUvToEllipsoid).

We then need to estimate two things:

1. The surface normal of the current voxel cell (at the point where the ray first entered the cell). This is needed to apply lighting effects in the rendering.
2. The optimal distance to step along the ray

The step length could simply be some fraction of the average voxel size. (This is the behavior for CYLINDER and ELLIPSOID prior to this PR.) However, we need to balance two factors:

• For sharp rendering of voxel edges, the step needs to be small enough that the ray position does not skip over any voxels. If the ray is crossing an edge or corner of a voxel, the required step could be arbitrarily small.
• For good performance, the step needs to be reasonably large.

How surface normals and optimum step size are computed

Prior to this PR, both normals and optimum step size for BOX voxels were computed by constructing the bounding box of the current cell, and performing a new ray-shape intersection with that bounding box. This intersection yielded both the normal of the incident face, and the distance to the far side of the voxel. Using the distance to the far side as the next step length ensured that the step would be as long as possible, without skipping across any voxels.

Intersections with individual voxel cells suffer from two problems:

1. The bounding volume for CYLINDER and ELLIPSOID cells is complicated and not convex. A ray can have multiple intersections with the same cell.
2. Intersection tests suffer from precision problems as the cells become smaller.

This PR uses the partial derivatives of the coordinate transform to get a linear estimate of the surface normal and step length.

The derivatives of the coordinate transform describe how fast the shape coordinate is changing for a small change in Cartesian position. As an example, consider the Jacobian matrix of the partial derivatives of the transform from Cartesian coordinates $(x, y, z)$ to ellipsoidal coordinates $(\phi, \theta, h)$:

$$J = \begin{bmatrix} \frac{\partial \phi}{\partial x} & \frac{\partial \phi}{\partial y} & \frac{\partial \phi}{\partial z} & \\\ \frac{\partial \theta}{\partial x} & \frac{\partial \theta}{\partial y} & \frac{\partial \theta}{\partial z} & \\\ \frac{\partial h}{\partial x} & \frac{\partial h}{\partial y} & \frac{\partial h}{\partial z} \end{bmatrix}$$

This matrix is evaluated at the current Cartesian ray position, and then multiplied with the ray direction to obtain the gradient (the rate of change) of the ellipsoidal coordinates along the ray:

$$\nabla (\phi, \theta, h) = J \vec{d}$$

Once we know how fast each ellipsoidal coordinate is changing for a small step along the ray, it is relatively simple to estimate the step length needed to reach the next cell boundary.

The surface normal at the point of entry into the voxel is aligned with the coordinate vector along which the ray last crossed a cell boundary. Conveniently, the rows of the Jacobian matrix are scaled versions of the local ellipsoidal coordinate vectors in Cartesian space. We can therefore simply normalize the appropriate row, and use the sign of the gradient, to obtain the surface normal at the point of intersection. For example, if the most-recently-crossed cell boundary (along the ray direction) is a boundary in longitude ($\phi$), the relevant surface normal is $-sign(J \vec{d}) \frac{J_{123}}{|| J_{123} ||}$

Key code changes

• Removed the computeApproximateStepSize method (and associated specs) from VoxelShape, since all step size calculations are now done in the shader.
• Simplified IntersectBox.glsl to use a simpler logic, more similar to the other shapes. Note that the intersectBox function is no longer used for individual voxel intersections. It is only used for the initial shape intersection to find the shape bounds.
• Moved some shared utility functions into VoxelUtils.glsl.\
• Rewrote the getStepSize function in VoxelFS.glsl to use the same Jacobian-based method for all shape types.
• Reworked convertUvToBox.glsl, convertUvToCylinder.glsl, and convertUvToEllipsoid.glsl to return the Jacobian matrix of partial derivatives, along with the converted shape space coordinate.

For the convertUvTo???.glsl changes: the old convertUvToShapeUvSpace methods converted from Cartesian space, to shape space, to a UV shape space (spanning [0, 1]), all in one method. This has now been broken into two steps:

1. convertUvToShapeSpaceDerivative, which transforms the coordinate to shape space and computes the partial derivatives of the transformation.
2. convertShapeToShapeUvSpace, which converts from the native shape space to a UV space spanning [0, 1].

Testing plan

Load this local Sandcastle and verify that all 3 shape types are rendered with sharp edges and lighting.

Author checklist

• I have submitted a Contributor License Agreement
• I have added my name to CONTRIBUTORS.md
• I have updated CHANGES.md with a short summary of my change
• I have added or updated unit tests to ensure consistent code coverage
• I have update the inline documentation, and included code examples where relevant
• I have performed a self-review of my code

[github-actions]

Thank you for the pull request, @jjhembd!

✅ We can confirm we have a CLA on file for you.

[jjhembd]

This evoluteScale used to be computed in single precision on the shader, in the nearestPointOnEllipse function in convertUvToEllipsoid.glsl. It is the difference of two very similar values, so it is safer to compute it in double precision on the CPU to avoid subtractive cancellation.

[ggetz]

Where does the 0.02 come from? Should this be configurable, or at least declared as a constant?

[ggetz]

I know this was existing code that was just moved to this file in this PR, but would it be possible to better explain what these magic numbers are in a comment?

[jjhembd]

I haven't been able to trace where this hash function came from. It is similar to one of the functions in a popular shadertoy example, but most uses that I have seen use a different value for the hashscale. See for example the hash12 function in this gist.
I added a link to the shadertoy in a comment.

[ggetz]

Are you able to address this TODO?

[ggetz]

Is there a relevant replacement spec for the Jacobian approach? (It's OK if there's not.)

[jjhembd]

No, all step size math has been moved to the shader.

[ggetz]

@jjhembd I'm seeing some artifacts in the center of the screen when zoomed out for both box and cyllinder types, but not ellipsoid:

Box:

Cyllinder:

They seem to flicker, even when the camera is not moving. Maybe some kind of z-fighting?

[jjhembd]

@jjhembd I'm seeing some artifacts in the center of the screen when zoomed out for both box and cyllinder types, but not ellipsoid:

This should be fixed. Let me know if it works on your end now!

[ggetz]

Thanks @jjhembd, the cylinder is looking good.

However, I still seem to see some artifacts with the box shape. They persist despite adjust the step size in the voxel inspector, but do change appearance when the value is changed.

[jjhembd]

Thanks @jjhembd, the cylinder is looking good.

However, I still seem to see some artifacts with the box shape. They persist despite adjust the step size in the voxel inspector, but do change appearance when the value is changed.

As discussed offline: we're not able to reproduce this artifact on a Windows machine. I've opened #11894 to track the bug.

[ggetz]

Thanks @jjhembd! This should be good to go then.