Math: line intersections

Lately a lot of my shaders have used a bit of raymarching, where I primarily rely on finding the intersections between a line (ray) and a shape. In this post I will go over the general approach for calculating these intersections and provide a list of HLSL implementations. The shapes discussed here will include planes, quadric surfaces (e.g. a sphere) and some useful combinations of them.

This tutorial was made with Unity 2020.3.1f1 and the Universal Render Pipeline (URP). Note that the general concept can be applied to any shader or coding environment. You can download the finished HLSL file, including an example shader, as Unity package through my git repository. If you’ve never imported a Unity package via a git repository this documentation will be helpful.

Acknowledgements

As this is mostly a tutorial on mathematics and its implementation my main acknowledgement would be my university level classes in Calculus and Linear Algebra. For those interested I shall list the Calculus book I’ve used, but know that I’ve put plenty of links to the relevant wiki articles throughout this tutorial as well.

  • Calculus: A Complete Course (8th edition), by Robert Adams and Christopher Essex.

Contents

  1. Finding intersections
    1. Space transformation
    2. HLSL include file
    3. Plane intersection
    4. Visualise shader
  2. Quadric surfaces
    1. Elliptic cylinder
    2. Ellipsoid
    3. Hyperboloid
    4. Paraboloid
  3. Conditions and combinations
    1. Ellipse
    2. Capped cylinder
    3. Capped ellipsoid
    4. Capped half hyperboloid
    5. Rectangle
    6. Box
  4. Conclusion
  5. References

1 Finding intersections

The method for finding line intersections that I will discuss here relies on knowing a general equation for a shape and substituting it’s coordinates with a parameterised line.

Starting with what a parameterised line is. It is defined solely by an origin and a direction that scales with a parameter (a number that can be anywhere between negative and positive infinity). The points on this line are given by the equation \overrightarrow{P} = \overrightarrow{O_L} + \overrightarrow{D_L} t. Here \overrightarrow{O_L} is the origin, \overrightarrow{D_L} the direction and t the parameter. Our parameter t thus denotes the distance from our line’s origin along the direction and can be anywhere between negative and positive infinity.

For those unfamiliar with the notation, the arrows denote vectors with an x, y and z component, like \overrightarrow{P} = \begin{pmatrix} P_x \\ P_y \\ P_z \end{pmatrix} = \begin{pmatrix} O_x + D_x t \\ O_y + D_y t \\ O_z + D_z t \end{pmatrix}.

Now to find the intersections between this line and a shape we need to know the shape’s equation, substitute it’s x, y and z component with \overrightarrow{P}, and solve for t. For example, the equation of a sphere would be x^2 + y^2 + z^2 = r^2, with r the radius. To find the intersections we would substitute x \rightarrow P_x, y \rightarrow P_y and z \rightarrow P_z. This would result in an equation that has only one variable, our parameter t. Solving this equation for t would then give us the distance between the intersection and the line’s origin \overrightarrow{O_L} (if there is an intersection).

So to summarise, the method consists of the following 3 steps:

  1. Know/find the shape equation
  2. Substitute coordinates with the parameterised line equation
  3. Solve the equation for parameter t, the distance to the intersection

Note that this method will also return negative distances, i.e. when the intersection is behind the line’s origin. In most use cases you would either discard these negative distances or use it to say something about your location, e.g. if the shape is convex and there is one intersection behind you and one in front you are within the shape’s volume.

1.1 Space transformation

Most of our shapes will point in some direction, like a plane that has its normal pointing upward along the y-axis. For simplicity’s sake we will define each shape equation in its own local space where the shape is centred on the origin and y is always the up axis, z is the forward axis and x is the right axis.

In order to translate and rotate our shapes to any place or direction we will be transforming our input line (that is likely defined in world space) to our shape’s local space, somewhat similar to transforming from world space to object space except we don’t adjust the scale. In order to do this we will provide two orthonormal vectors that form a basis of our shape space, the forward and up direction, in addition to the origin of the shape. These vectors must be defined within the same space as our input line. With this transformation we can easily define our shapes in their local x, y and z coordinates and use a transformation matrix to rotate and the origin to translate the shape to any location/rotation we want.

In short, our equations will be defined in their own space, where it is centred on the origin and always pointing upward. However, our input is defined in world space (or any other space), which consists of the line’s origin, the line’s direction, the shape’s origin and the shape’s direction (up and forward vectors). Finally, we use the shape’s origin and direction to transform our line to shape space.

I know this might be a bit much and the next sections will show it implemented in code, but let me know if you want a more in-depth tutorial on matrix transformations. In the meanwhile you can also check out this tutorial by Catlike Coding on matrix transformations.

1.2 HLSL include file

The actual HLSL implemenation of our line intersections will be in the form of a HLSL include file that we can add to any shader. Create a new text file and change its name (including format) to “LineIntersections.hlsl” and add the following lines.

#ifndef LINE_INTERSECTIONS_INCLUDED
#define LINE_INTERSECTIONS_INCLUDED

float3x3 constructTransitionMatrix(float3 forwardDir, float3 upDir)
{
    float3 rightDir = cross(forwardDir, upDir);
    float3x3 result = {rightDir, upDir, forwardDir};
    return result;
}

#endif

The constructTransitionMatrix function uses the shape’s forward and up direction to construct a rotation matrix to transform from our input space to the shape’s space. Notice the definition lines starting with a #. These are used to ensure that we don’t accidently include our functions twice in a shader. Basically, it checks if LINE_INTERSECTIONS_INCLUDED is defined, if not, define it and add our functions.

1.3 Plane intersection

Now onto the actual stuff, finding intersections. We’ll start with the simplest shape, as it always has exactly one intersection, the flat plane. We can define a plane with the following equation [1]:

\displaystyle s_x x + s_y y + s_z z + s_t = 0

Here the shape parameters \begin{pmatrix} s_x \\ s_y \\ s_z \end{pmatrix} determine the normal of the plane and s_t an offset along the normal. Now in order to find the intersection with the line we need to substitute our \begin{pmatrix} x \\ y \\ z \end{pmatrix} coordinates with the line parameterisation \overrightarrow{P} = \begin{pmatrix} O_x + D_x t \\ O_y + D_y t \\ O_z + D_z t \end{pmatrix}. Because this is a simple linear equation we can solve it for t to find the distance between the line’s origin and the intersection. This process is written out below.

\displaystyle s_x (O_x + D_x t) + s_y (O_y + D_y t) + s_z (O_z + D_z t) + s_t = 0

\displaystyle (s_x D_x + s_y D_y + s_z D_z) t + s_x O_x + s_y O_y + s_z O_z + s_t = 0

\displaystyle t = -\frac{s_x O_x + s_y O_y + s_z O_z + s_t}{s_x D_x + s_y D_y + s_z D_z}

Because we’re transforming our line to the shape’s local space, which already defines the shape’s direction and origin/offset, we don’t have to provide the shape parameters as they overlap. In our shape space the plane always points up (the y-direction), this allows us to simplify the equation to t = -\frac{O_y}{D_y}.

The code below shows how this is implemented in HLSL, where we transform the input line to shape space and calculate t. As input we provide the line origin, line direction, shape origin, and shape up direction. Add the following lines to our HLSL file below the constructTransitionMatrix function but before #endif.

// Based on plane equation from https://en.wikipedia.org/wiki/Plane_(geometry)#Point%E2%80%93normal_form_and_general_form_of_the_equation_of_a_plane
float intersectPlane(float3 lineOrigin, float3 lineDir, float3 shapeOrigin, float3 shapeUpDir)
{
    // Transform line origin and direction from world space to the shape space
    float3x3 transitionMatrix = construcTransitionMatrix(float3(0,0,0), shapeUpDir);
    float3 lO = mul(transitionMatrix, lineOrigin - shapeOrigin);
    float3 lD = mul(transitionMatrix, lineDir);

    float denominator = lD.y;
    float numerator = lO.y;
    
    return - numerator / denominator;
}

In the previous section I noted we must always provide two orthonormal vectors, but because our plane is infinite and thus symmetric within the plane we can get away with only providing one vector, the normal or up direction.

1.4 Visualise shader

Having this LineIntersections.hlsl is great and all, but we can’t see anything yet as we do not have a shader or material. Let’s change that by making a new shader called VisualiseIntersection.shader and fill it with the code block below.

Shader "KelvinvanHoorn/VisualiseIntersection"
{
    Properties
    {
        _ShapeParams ("Shape parameters", vector) = (1,1,1,1)
        _DistanceScale ("Distance scale", float) = 100
        _CapHeights ("Cap heights", vector) = (0,1,0,0)
    }
    SubShader
    {
        Tags { "RenderType" = "Opaque" "RenderPipeline" = "UniversalRenderPipeline" "Queue" = "Geometry" }
        Cull Front
 
        Pass
        {
            HLSLPROGRAM
            #pragma vertex vert
            #pragma fragment frag
 
            #include "Packages/com.unity.render-pipelines.universal/ShaderLibrary/Core.hlsl"
            #include "LineIntersections.hlsl"          
 
            struct Attributes
            {
                float4 posOS    : POSITION;
            };
 
            struct v2f
            {
                float4 posCS        : SV_POSITION;
                float3 posWS        : TEXCOORD0;
 
                float3 objectOrigin : TEXCOORD1;
            };

            float4 _ShapeParams;
            float _DistanceScale;
            float2 _CapHeights;
 
            v2f vert(Attributes IN)
            {
                v2f OUT = (v2f)0;
 
                VertexPositionInputs vertexInput = GetVertexPositionInputs(IN.posOS.xyz);
 
                OUT.posCS = vertexInput.positionCS;
                OUT.posWS = vertexInput.positionWS;
 
                // Object information
                OUT.objectOrigin = UNITY_MATRIX_M._m03_m13_m23;

                return OUT;
            }
 
            float4 frag (v2f IN) : SV_Target
            {
                // Line information
                float3 lineOrigin = _WorldSpaceCameraPos;
                float3 lineDir = normalize(IN.posWS - _WorldSpaceCameraPos);

                // Shape information
                float3 shapeUpDir = normalize(mul(unity_ObjectToWorld, float4(0,1,0,0)).xyz);
                float3 shapeForwardDir = normalize(mul(unity_ObjectToWorld, float4(0,0,1,0)).xyz);
                float3 shapeOrigin = IN.objectOrigin;

                float r = 0;
                float g = 0;
                float b = 0;

                // Intersect information
                float intersect = intersectPlane(lineOrigin, lineDir, shapeOrigin, shapeUpDir);
                if(intersect < maxLineDst && intersect > 0)
                {
                    float3 pos = lineOrigin - shapeOrigin + lineDir * intersect;
                    r = saturate(abs(pos.x)/_DistanceScale);
                    g = saturate(abs(pos.y)/_DistanceScale);
                    b = saturate(abs(pos.z)/_DistanceScale);
                }
     
                return float4(r, g, b, 1);
            }
            ENDHLSL
        }
    }
}

This tutorial is primarily focussed on the mathematics and implementation behind finding line intersections, I will thus only shortly gloss over what the shader does.

We use 3 properties to influence our shapes, the shape parameters which usually say something about the scale, the distance scaling which is used to colour the shape and the cap heights which will be used in section 3 for the cap positions.

We include our LineIntersections.hlsl file at line 21. The #include requires the path to our HLSL file, if this is in the same folder we can simply put the name there. Otherwise it would look something like the following.

#include "Assets/folder/subfolder/LineIntersections.hlsl"
// or if you're using the package 
#include "Packages/com.kelvinvanhoorn.line-intersections/ShaderLibrary/LineIntersections.hlsl"

We start our line at the camera’s world space location, with its direction being towards the fragment world position of the gameObject to which we apply the material. The shape’s origin is set to the gameObject’s origin and the shape’s forward and up direction are set to the gameObject’s local forward and up direction.

Finally we get the intersection information to determine the output colour. If there is an intersection and it is in front of the camera we colour it depending on its distance with respect to the shape’s origin, otherwise it is black to show the bounds of our gameObject.

I’ve applied the shader to a standard Unity sphere with all scales set to 50. The image below shows a visualisation of the plane intersection. The red and blue line correspond with the world x-axis and z-axis respectively.

Plane visualisation on a 50 scale sphere object with distance scale = 50, the red and blue line correspond to the x and z axis of the scene

2 Quadric surfaces

In the following sections I’ll go through the surface equations of quadrics and use the quadratic formula to find their intersections. Quadrics describe surfaces that are expressed as a polynomial of second degree, i.e. quadratic functions like the following.

\displaystyle a x^2 + b x + c = 0

When we substitute the line equation into a quadric equation it will also be quadratic in t, and thus like any quadratic function it will have at most 2 intersections. Moreover, it means we can refactor it and use the quadratic formula to solve it.

You’ve probably already encountered the quadratic formula in a mathematics class. Basically we refactor a quadratic equation to a general form, calculate the discriminant (d) and use that to calculate the intersections. Written out this method consists of the following equations.

\displaystyle a t^2 + b t + c = 0

\displaystyle d = b^2 - 4 a c

\displaystyle t = \frac{-b \pm \sqrt{d}}{2 a}

Note that depending on the discriminant there are either 2, 1 or no intersections, which respectively corresponds to d > 0, d = 0 or d < 0.

Implementing this in HLSL is a straightforward copy of the mathematical equations, except that we need to decide on what to return if there are no intersections. In my case I’ve decided to return the maximum float value as we should never deal with an intersection that far away. Add the following lines to our HLSL file just above the constructTransitionMatrix function. We will use this function to solve all our quadric surface intersections.

static const float maxLineDst = 3.402823466e+38;

float2 quadraticFormula(float a, float b, float c)
{
    float discriminant = b * b - 4 * a * c;
    if (discriminant >= 0) 
    {
        float s = sqrt(discriminant);

        float minDst = (-b - s) / (2 * a);
        float maxDst = (-b + s) / (2 * a);
 
        return float2(minDst, maxDst);
    }
    // returns maxLineDst if no intersection
    return float2(maxLineDst, maxLineDst);
}

2.1 Elliptic cylinder

Elliptic cylinder [2]

The first quadric to discuss is the elliptic cylinder, an infinite cylinder with two focal points, i.e. the x-z plane cross-section is an ellipse. The equation for an elliptic cylinder is as follows [3]:

\displaystyle \frac{x^2}{s_x^2} + \frac{z^2}{s_z^2} = 1.

Here the scale parameters s_x and s_z act similar to the radius of a circular cylinder. Substituting the line equation results in:

\displaystyle \frac{(O_x + D_x t)^2}{s_x^2} + \frac{(O_z + D_z t)^2}{s_z^2} - 1 = 0.

Refactoring this equation results in a quadratic function of form a t^2 + b t + c = 0, with the following coefficients:

\displaystyle a = \frac{D_x^2}{s_x^2} + \frac{D_z^2}{s_z^2},

\displaystyle b = 2(\frac{O_x D_x}{s_x^2} + \frac{O_z D_z}{s_z^2}),

\displaystyle c = \frac{O_x^2}{s_x^2} + \frac{O_z^2}{s_z^2} - 1.

We can then use the quadratic formula from the previous section to find up to two intersections. The code below shows how this is implemented in HLSL. As input we provide the line origin, line direction, shape origin, shape forward direction, shape up direction, and the scale parameters s_x and s_z. Add the following lines to our HLSL file below the intersectPlane function.

// Based on quadrics equations from https://en.wikipedia.org/wiki/Quadric
float2 intersectEllipticCylinder(float3 lineOrigin, float3 lineDir, float3 shapeOrigin, float3 shapeForwardDir, float3 shapeUpDir, float2 scaleParams)
{
    // Transform line origin and direction from world space to the shape space
    float3x3 transitionMatrix = construcTransitionMatrix(shapeForwardDir, shapeUpDir);
    float3 lO = mul(transitionMatrix, lineOrigin - shapeOrigin);
    float3 lD = mul(transitionMatrix, lineDir);

    // Squared reciprocals of the shape parameters s_x and s_z
    float rsx = 1.0 / (scaleParams.x*scaleParams.x);
    float rsz = 1.0 / (scaleParams.y*scaleParams.y);

    float a = lD.x*lD.x * rsx + lD.z*lD.z * rsz;
    float b = 2 * (lO.x * lD.x * rsx + lO.z * lD.z * rsz);
    float c = lO.x*lO.x * rsx + lO.z*lO.z * rsz - 1;
    
    return quadraticFormula(a, b, c);
}

To visualise this intersection change the intersection information in our visualise shader to the following. We now determine the colour on the first intersection in front of the camera, as we get 2 values now instead of 1.

// Intersect information
float2 intersect = intersectEllipticCylinder(lineOrigin, lineDir, shapeOrigin, shapeForwardDir, shapeUpDir, _ShapeParams.xz);
float frontIntersection = intersect.x >= 0 ? intersect.x : intersect.y;
if(frontIntersection < maxLineDst && frontIntersection >= 0)
{
    float3 pos = lineOrigin - shapeOrigin + lineDir * frontIntersection;
    r = saturate(abs(pos.x)/_DistanceScale);
    g = saturate(abs(pos.y)/_DistanceScale);
    b = saturate(abs(pos.z)/_DistanceScale);
}
Elliptic cylinder on a 50 scale sphere object with parameters sx = 10, sz = 15 and distance scale = 15

2.2 Ellipsoid

Ellipsoid [2]

The ellipsoid is my most used quadric surface, as it includes the sphere where all axes have the same radius. The equation for an ellipsoid is given by the following [3]:

\displaystyle \frac{x^2}{s_x^2} + \frac{y^2}{s_y^2} + \frac{z^2}{s_z^2} = 1.

Here the scale parameters s_x, s_y and s_z act similar to the radius of a sphere. Substituting the line equation results in:

\displaystyle \frac{(O_x + D_x t)^2}{s_x^2} + \frac{(O_y + D_y t)^2}{s_y^2} + \frac{(O_z + D_z t)^2}{s_z^2} - 1 = 0.

Refactoring this equation results in a quadratic function of form a t^2 + b t + c = 0, with the following coefficients:

\displaystyle a = \frac{D_x^2}{s_x^2} + \frac{D_y^2}{s_y^2} + \frac{D_z^2}{s_z^2},

\displaystyle b = 2(\frac{O_x D_x}{s_x^2} + \frac{O_y D_y}{s_y^2} + \frac{O_z D_z}{s_z^2}),

\displaystyle c = \frac{O_x^2}{s_x^2} + \frac{O_y^2}{s_y^2} + \frac{O_z^2}{s_z^2} - 1.

Using the quadratic formula we can find up to two intersections. The code below shows how this is implemented in HLSL. As input we provide the line origin, line direction, shape origin, shape forward direction, shape up direction, and the scale parameters s_x, s_y and s_z. Add the following lines to our HLSL file below the intersectEllipticCylinder function.

// Based on quadrics equations from https://en.wikipedia.org/wiki/Quadric
float2 intersectEllipsoid(float3 lineOrigin, float3 lineDir, float3 shapeOrigin, float3 shapeForwardDir, float3 shapeUpDir, float3 scaleParams)
{
    // Transform line origin and direction from world space to the shape space
    float3x3 transitionMatrix = constructTransitionMatrix(shapeForwardDir, shapeUpDir);
    float3 lO = mul(transitionMatrix, lineOrigin - shapeOrigin);
    float3 lD = mul(transitionMatrix, lineDir);

    // Squared reciprocals of the shape parameters
    float rsx = 1.0 / (scaleParams.x*scaleParams.x);
    float rsy = 1.0 / (scaleParams.y*scaleParams.y);
    float rsz = 1.0 / (scaleParams.z*scaleParams.z);

    float a = lD.x*lD.x * rsx + lD.y*lD.y * rsy + lD.z*lD.z * rsz;
    float b = 2 * (lO.x * lD.x * rsx + lO.y * lD.y * rsy + lO.z * lD.z * rsz);
    float c = lO.x*lO.x * rsx + lO.y*lO.y * rsy + lO.z*lO.z * rsz - 1;

    return quadraticFormula(a, b, c);
}

To visualise this intersection change the intersection information in our visualise shader to the following.

// Intersect information
float2 intersect = intersectEllipsoid(lineOrigin, lineDir, shapeOrigin, shapeForwardDir, shapeUpDir, _ShapeParams.xyz);
float frontIntersection = intersect.x >= 0 ? intersect.x : intersect.y;
if(frontIntersection < maxLineDst && frontIntersection >= 0)
{
    float3 pos = lineOrigin - shapeOrigin + lineDir * frontIntersection;
    r = saturate(abs(pos.x)/_DistanceScale);
    g = saturate(abs(pos.y)/_DistanceScale);
    b = saturate(abs(pos.z)/_DistanceScale);
}
Ellipsoid on a 50 scale sphere object with parameters sx = 10, sy = 15, sz = 20 and distance scale = 20

2.3 Hyperboloid

Hyperbolic hyperboloid [2]
st > 0
Elliptic hyperboloid [2]
st < 0
Elliptic cone [2]
st = 0

The hyperboloid equation is very similar to the ellipsoid, except that the y component is negative and that we have an additional parameter for type (st). Depending on this type parameter we can have a hyperbolic (st > 0) or elliptic (st < 0) hyperboloid or an elliptic cone (st = 0). The equation for hyperboloids is given by the following [3]:

\displaystyle \frac{x^2}{s_x^2} - \frac{y^2}{s_y^2} + \frac{z^2}{s_z^2} = s_t.

By substituting the line equation we get:

\displaystyle \frac{(O_x + D_x t)^2}{s_x^2} - \frac{(O_y + D_y t)^2}{s_y^2} + \frac{(O_z + D_z t)^2}{s_z^2} - s_t = 0.

Refactoring this equation results in a quadratic function of form a t^2 + b t + c = 0, with the following coefficients:

\displaystyle a = \frac{D_x^2}{s_x^2} - \frac{D_y^2}{s_y^2} + \frac{D_z^2}{s_z^2},

\displaystyle b = 2(\frac{O_x D_x}{s_x^2} - \frac{O_y D_y}{s_y^2} + \frac{O_z D_z}{s_z^2}),

\displaystyle c = \frac{O_x^2}{s_x^2} - \frac{O_y^2}{s_y^2} + \frac{O_z^2}{s_z^2} - s_t.

Using the quadratic formula we can find up to two intersections. The code below shows how this is implemented in HLSL. As input we provide the line origin, line direction, shape origin, shape forward direction, shape up direction, the scale parameters s_x, s_y and s_z, and the type parameter s_t. Add the following lines to our HLSL file below the intersectEllipsoid function.

// Based on quadrics equations from https://en.wikipedia.org/wiki/Quadric
float2 intersectHyperboloid(float3 lineOrigin, float3 lineDir, float3 shapeOrigin, float3 shapeForwardDir, float3 shapeUpDir, float3 scaleParams, float typeParam)
{
    // Transform line origin and direction from world space to the shape space
    float3x3 transitionMatrix = constructTransitionMatrix(shapeForwardDir, shapeUpDir);
    float3 lO = mul(transitionMatrix, lineOrigin - shapeOrigin);
    float3 lD = mul(transitionMatrix, lineDir);

    // Squared reciprocals of the shape parameters
    float rsx = 1.0 / (scaleParams.x*scaleParams.x);
    float rsy = 1.0 / (scaleParams.y*scaleParams.y);
    float rsz = 1.0 / (scaleParams.z*scaleParams.z);

    float a = lD.x*lD.x * rsx - lD.y*lD.y * rsy + lD.z*lD.z * rsz;
    float b = 2 * (lO.x * lD.x * rsx - lO.y * lD.y * rsy + lO.z * lD.z * rsz);
    float c = lO.x*lO.x * rsx - lO.y*lO.y * rsy + lO.z*lO.z * rsz - typeParam;

    return quadraticFormula(a, b, c);
}

To visualise this intersection change the intersection information in our visualise shader to the following. I’ve used _ShapeParams.w for the type parameter.

// Intersect information
float2 intersect = intersectHyperboloid(lineOrigin, lineDir, shapeOrigin, shapeForwardDir, shapeUpDir, _ShapeParams.xyz, _ShapeParams.w);
float frontIntersection = intersect.x >= 0 ? intersect.x : intersect.y;
if(frontIntersection < maxLineDst && frontIntersection >= 0)
{
    float3 pos = lineOrigin - shapeOrigin + lineDir * frontIntersection;
    r = saturate(abs(pos.x)/_DistanceScale);
    g = saturate(abs(pos.y)/_DistanceScale);
    b = saturate(abs(pos.z)/_DistanceScale);
}
Hyperboloid on a 50 scale sphere object with parameters sx = 10, sy = 15, sz = 15, st = 0.5 and distance scale = 15

2.4 Paraboloid

Elliptic paraboloid [2]
st = 1

Hyperbolic paraboloid [2]
st = -1

I’m not sure how useful the paraboloids are for raytracing, but they are a quadric so I’ve added them here for completeness. Paraboloids also have a type parameter where st = 1 describes an elliptic paraboloid and st = -1 describes a hyperbolic paraboloid. The equation for hyperboloids is given by the following [3]:

\displaystyle \frac{x^2}{s_x^2} + sign(s_t) \frac{z^2}{s_z^2} - y = 0

The sign function returns -1 if negative, 1 if positive and 0 if 0. By substituting the line equation we get:

\displaystyle \frac{(O_x + D_x t)^2}{s_x^2} + sign(s_t) \frac{(O_z + D_z t)^2}{s_z^2} - (O_y + D_y t) = 0

Refactoring this equation results in a quadratic function of form a t^2 + b t + c = 0, with the following coefficients:

\displaystyle a = \frac{D_x^2}{s_x^2} + sign(s_t) \frac{D_z^2}{s_z^2},

\displaystyle b = 2(\frac{O_x D_x}{s_x^2} + sign(s_t) \frac{O_z D_z}{s_z^2}) - D_y,

\displaystyle c = \frac{O_x^2}{s_x^2} + sign(s_t) \frac{O_z^2}{s_z^2} - O_y.

Using the quadratic formula we can find up to two intersections. The code below shows how this is implemented in HLSL. As input we provide the line origin, line direction, shape origin, shape forward direction, shape up direction, the scale parameters s_x and s_z, and the type parameter s_t. Add the following lines to our HLSL file below the intersectHyperboloid function.

// Based on quadrics equations from https://en.wikipedia.org/wiki/Quadric
float2 intersectParaboloid(float3 lineOrigin, float3 lineDir, float3 shapeOrigin, float3 shapeForwardDir, float3 shapeUpDir, float2 scaleParams, float typeParam)
{
    // Transform line origin and direction from world space to the shape space
    float3x3 transitionMatrix = constructTransitionMatrix(shapeForwardDir, shapeUpDir);
    float3 lO = mul(transitionMatrix, lineOrigin - shapeOrigin);
    float3 lD = mul(transitionMatrix, lineDir);

    // Squared reciprocals of the shape parameters
    float rsx = 1.0 / (scaleParams.x*scaleParams.x);
    float rsz = 1.0 / (scaleParams.y*scaleParams.y);

    // Sign of type parameter
    float sst = sign(typeParam);

    float a = lD.x*lD.x * rsx + sst * lD.z*lD.z * rsz;
    float b = 2 * (lO.x * lD.x * rsx + sst * lO.z * lD.z * rsz) - lD.y ;
    float c = lO.x*lO.x * rsx + sst * lO.z*lO.z * rsz - lO.y;

    return quadraticFormula(a, b, c);
}

To visualise this intersection change the intersection information in our visualise shader to the following.

// Intersect information
float2 intersect = intersectParaboloid(lineOrigin, lineDir, shapeOrigin, shapeForwardDir, shapeUpDir, _ShapeParams.xz, _ShapeParams.w);
float frontIntersection = intersect.x >= 0 ? intersect.x : intersect.y;
if(frontIntersection < maxLineDst && frontIntersection >= 0)
{
    float3 pos = lineOrigin - shapeOrigin + lineDir * frontIntersection;
    r = saturate(abs(pos.x)/_DistanceScale);
    g = saturate(abs(pos.y)/_DistanceScale);
    b = saturate(abs(pos.z)/_DistanceScale);
}
Paraboloid on a 50 scale sphere object with parameters sx = 10, sz = 5, st = 1 and distance scale = 10

3 Conditions and combinations

We’ve seen quite a few useful shapes, but we can make them even more versatile. One thing we can do is add additional conditions, like limiting the plane to be within a certain radius to form a circle or ellipse. Another thing would be to combine shapes, like using planes/ellipses to add caps to the cylinder, ellipsoid and hyperboloid. The following sections will shortly discuss some conditions/combinations and show how it is implemented in HLSL.

3.1 Ellipse

The ellipse is the simplest of combinations as it is a plane where we check if the intersection is within the bounds of an ellipse. The surface of an ellipse is described by following equation:

\displaystyle \frac{x^2}{s_x^2} + \frac{z^2}{s_z^2} \leq 1.

To check if a point is within the ellipse we substitute the x and z coordinate with our plane intersection coordinates.

The HLSL implementation for the ellipse thus starts the same as the plane. Only the end is different where check if the intersection point is within an ellipse, if it is we return the intersection, if not we return the max float. As input we provide the line origin, line direction, shape origin, shape forward direction, shape up direction, and the scale parameters s_x and s_z that define the ellipse. Add the following lines to our HLSL file below the intersectParaboloid function.

float intersectEllipse(float3 lineOrigin, float3 lineDir, float3 shapeOrigin, float3 shapeForwardDir, float3 shapeUpDir, float2 scaleParams)
{
    // Transform line origin and direction from world space to the shape space
    float3x3 transitionMatrix = constructTransitionMatrix(shapeForwardDir, shapeUpDir);
    float3 lO = mul(transitionMatrix, lineOrigin - shapeOrigin);
    float3 lD = mul(transitionMatrix, lineDir);

    float denominator = lD.y;
    float numerator = lO.y;

    float pIntersect = - numerator / denominator;

    float3 samplePos = lO + lD * pIntersect;
    bool isInEllipse = samplePos.x * samplePos.x / (scaleParams.x * scaleParams.x) + samplePos.z * samplePos.z / (scaleParams.y * scaleParams.y) <= 1;
    return isInEllipse ? pIntersect : maxLineDst;
}

To visualise this intersection change the intersection information in our visualise shader to the following.

// Intersect information
float intersect = intersectEllipse(lineOrigin, lineDir, shapeOrigin, shapeForwardDir, shapeUpDir, _ShapeParams.xz);
if(intersect < maxLineDst && intersect > 0)
{
    float3 pos = lineOrigin - shapeOrigin + lineDir * intersect;
    r = saturate(abs(pos.x)/_DistanceScale);
    g = saturate(abs(pos.y)/_DistanceScale);
    b = saturate(abs(pos.z)/_DistanceScale);
}
Ellipse on a 50 scale sphere object with parameters sx = 15, sz = 20 and distance scale = 20

3.2 Capped cylinder

Now that we have the ellipse function we can use it for other things like a capped cylinder, which consists of the infinite elliptical cylinder and two ellipses as caps.

To cap the cylinder we define it’s height through a scale parameter s_y, meaning that the top cap is located at y = 0.5 s_y and the bottom cap at y = -0.5 s_y. Anything above the top cap or below the bottom cap is no longer part of the shape.

In order to implement this we will first do the elliptical cylinder intersection and check if this intersection is between the caps, if not we set the distance to the max float. Next we’ll calculate the intersections for both ellipses and compare all 4 possible intersections. Luckily a capped cylinder can only have 2 intersections, the other 2 possible intersections will return the max float. Thus we can calculate and return the minimum and maximum distance intersections and return the 2 values like all our other functions.

The HLSL implementation follows the method above. As input we provide the line origin, line direction, shape origin, shape forward direction, shape up direction, and the scale parameters s_x, s_y and s_z, where s_x and s_z define the scales for both the cylinder and the ellipses and s_y the height of the cylinder. Add the following lines to our HLSL file below the intersectEllipse function.

float2 intersectCappedCylinder(float3 lineOrigin, float3 lineDir, float3 shapeOrigin, float3 shapeForwardDir, float3 shapeUpDir, float3 scaleParams)
{
    // Cylinder intersections
    float2 cylinderDst = intersectEllipticCylinder(lineOrigin, lineDir, shapeOrigin, shapeForwardDir, shapeUpDir, scaleParams.xz);

    // Cylinder intersect positions
    float3 samplePos1 = lineOrigin + lineDir * cylinderDst.x;
    float3 samplePos2 = lineOrigin + lineDir * cylinderDst.y;

    bool isMinDstBelowCap = dot(samplePos1 - shapeOrigin, shapeUpDir) <= -0.5 * scaleParams.y;
    bool isMinDstAboveCap = dot(samplePos1 - shapeOrigin, shapeUpDir) >= 0.5 * scaleParams.y;

    bool isMaxDstBelowCap = dot(samplePos2 - shapeOrigin, shapeUpDir) <= -0.5 * scaleParams.y;
    bool isMaxDstAboveCap = dot(samplePos2 - shapeOrigin, shapeUpDir) >= 0.5 * scaleParams.y;

    float minDst = cylinderDst.x;
    float maxDst = cylinderDst.y;
    if(isMinDstBelowCap || isMinDstAboveCap)
    {
        minDst = maxLineDst;
    }
    if(isMaxDstBelowCap || isMaxDstAboveCap)
    {
        maxDst = maxLineDst;
    }

    // Cap positions
    float3 ellipse1Pos = shapeOrigin - shapeUpDir * 0.5 * scaleParams.y;
    float3 ellipse2Pos = shapeOrigin + shapeUpDir * 0.5 * scaleParams.y;

    // Ellipse bottom and top cap intersections
    float ellipse1Dst = intersectEllipse(lineOrigin, lineDir, ellipse1Pos, shapeForwardDir, shapeUpDir, scaleParams.xz);
    float ellipse2Dst = intersectEllipse(lineOrigin, lineDir, ellipse2Pos, shapeForwardDir, shapeUpDir, scaleParams.xz);

    // Compare cylinder intersects to ellipse intersects for minDst
    minDst = min(min(ellipse1Dst, ellipse2Dst), minDst);
    // If value is maxLineDst make -maxLineDst for max comparison to work
    ellipse1Dst = ellipse1Dst == maxLineDst ? -maxLineDst : ellipse1Dst;
    ellipse2Dst = ellipse2Dst == maxLineDst ? -maxLineDst : ellipse2Dst;
    maxDst = maxDst == maxLineDst ? -maxLineDst : maxDst;

    // Compare cylinder intersects to ellipse intersects for maxDst
    // and make sure we don't return -maxLineDst
    maxDst = max(max(ellipse1Dst, ellipse2Dst), maxDst);
    maxDst = maxDst == -maxLineDst ? maxLineDst : maxDst;
    
    return float2(minDst, maxDst);
}

To visualise this intersection change the intersection information in our visualise shader to the following.

// Intersect information
float2 intersect = intersectCappedCylinder(lineOrigin, lineDir, shapeOrigin, shapeForwardDir, shapeUpDir, _ShapeParams.xyz);
float frontIntersection = intersect.x >= 0 ? intersect.x : intersect.y;
if(frontIntersection < maxLineDst && frontIntersection >= 0)
{
    float3 pos = lineOrigin - shapeOrigin + lineDir * frontIntersection;
    r = saturate(abs(pos.x)/_DistanceScale);
    g = saturate(abs(pos.y)/_DistanceScale);
    b = saturate(abs(pos.z)/_DistanceScale);
}
Capped cylinder on a 50 scale sphere object with parameters sx = 15, sy = 10, sz = 20 and distance scale = 20

A slightly different version of the capped cylinder was used in my supermassive black hole shader. There it was used to describe the accretion disc circling around the centre of the black hole. The only difference being an additional condition of an inner radius, where there is no cylinder or ellipse in the middle.

3.3 Capped ellipsoid

The capped ellipsoid is slightly different, as the scale of of the caps change with their heights. A cap at the centre of the ellipsoid will be greater than that of a cap at the ends. Luckily we can calculate the height dependent ellipse scales using the equation for an ellipsoid from section 2.2. This is done by setting either x or z to 0 and solve for the other as well as putting in your cap’s height as y position. From the ellipsoid equation this results in the following [4]:

\displaystyle \frac{x^2}{s_x^2} + \frac{y^2}{s_y^2} + \frac{z^2}{s_z^2} = 1,

\displaystyle s_{e,x} \rightarrow x (z = 0) = s_x \sqrt{1 - \frac{h^2}{s_y^2}},

\displaystyle s_{e,z} \rightarrow z (x=0) = s_z \sqrt{1 - \frac{h^2}{s_y^2}}.

Here s_{e,x} denotes the x-scale parameter of the ellipse, s_x the x-scale parameter of the ellipsoid and h the height position of the cap. Note that, because an ellipsoid has a finite height determined by s_y, we can remap it to a range of -1 to 1 range (-1 being the bottom of the ellipsoid and 1 the top) which is nicer to handle. In that case the equations would be of form s_{e,x} = s_x \sqrt{1 - h^2}, with the new h between -1 and 1.

The HLSL implementation is similar to the capped cylinder, meaning we first find the ellipsoid intersections and check if they are between the caps followed by a comparison with the two ellipses. As input we provide the line origin, line direction, shape origin, shape forward direction, shape up direction, the scale parameters s_x, s_y and s_z, and two cap heights which get clamped between -1 and 1. Add the following lines to our HLSL file below the intersectCappedCylinder function.

float2 intersectCappedEllipsoid(float3 lineOrigin, float3 lineDir, float3 shapeOrigin, float3 shapeForwardDir, float3 shapeUpDir, float3 scaleParams, float2 capHeight)
{
    // Ellipsoid intersections
    float2 ellipsoidDst = intersectEllipsoid(lineOrigin, lineDir, shapeOrigin, shapeForwardDir, shapeUpDir, scaleParams);

    // Cap heights in -1 to 1 range
    float cap1Height = clamp(min(capHeight.x, capHeight.y), -1, 1);
    float cap2Height = clamp(max(capHeight.x, capHeight.y), -1, 1);

    // Ellipsoid intersect positions
    float3 samplePos1 = lineOrigin + lineDir * ellipsoidDst.x;
    float3 samplePos2 = lineOrigin + lineDir * ellipsoidDst.y;

    bool isMinDstBelowCap = dot(samplePos1 - shapeOrigin, shapeUpDir) <= cap1Height * scaleParams.y;
    bool isMinDstAboveCap = dot(samplePos1 - shapeOrigin, shapeUpDir) >= cap2Height * scaleParams.y;

    bool isMaxDstBelowCap = dot(samplePos2 - shapeOrigin, shapeUpDir) <= cap1Height * scaleParams.y;
    bool isMaxDstAboveCap = dot(samplePos2 - shapeOrigin, shapeUpDir) >= cap2Height * scaleParams.y;

    float minDst = ellipsoidDst.x;
    float maxDst = ellipsoidDst.y;
    if(isMinDstBelowCap || isMinDstAboveCap)
        minDst = maxLineDst;
    if(isMaxDstBelowCap || isMaxDstAboveCap)
        maxDst = maxLineDst;
    
    // Cap positions
    float3 ellipse1Pos = shapeOrigin + cap1Height * scaleParams.y * shapeUpDir;
    float3 ellipse2Pos = shapeOrigin + cap2Height * scaleParams.y * shapeUpDir;

    // Ellipsoidal cap radius from https://keisan.casio.com/keisan/image/volume%20of%20an%20ellipsoidal%20cap.pdf
    float e1HeightSqr = cap1Height * cap1Height;
    float e1ScaleX = scaleParams.x * sqrt(1 - e1HeightSqr);
    float e1ScaleZ = scaleParams.z * sqrt(1 - e1HeightSqr);

    float e2HeightSqr = cap2Height * cap2Height;
    float e2ScaleX = scaleParams.x * sqrt(1 - e2HeightSqr);
    float e2ScaleZ = scaleParams.z * sqrt(1 - e2HeightSqr);

    // Ellipse bottom and top cap intersections
    float ellipse1Dst = intersectEllipse(lineOrigin, lineDir, ellipse1Pos, shapeForwardDir, shapeUpDir, float2(e1ScaleX, e1ScaleZ));
    float ellipse2Dst = intersectEllipse(lineOrigin, lineDir, ellipse2Pos, shapeForwardDir, shapeUpDir, float2(e2ScaleX, e2ScaleZ));

    // Compare ellipsoid intersects to ellipse intersects for minDst
    minDst = min(min(ellipse1Dst, ellipse2Dst), minDst);
    // If value is maxLineDst make -maxLineDst for max comparison to work
    ellipse1Dst = ellipse1Dst == maxLineDst ? -maxLineDst : ellipse1Dst;
    ellipse2Dst = ellipse2Dst == maxLineDst ? -maxLineDst : ellipse2Dst;
    maxDst = maxDst == maxLineDst ? -maxLineDst : maxDst;

    // Compare ellipsoid intersects to ellipse intersects for maxDst
    // and make sure we don't return -maxLineDst
    maxDst = max(max(ellipse1Dst, ellipse2Dst), maxDst);
    maxDst = maxDst == -maxLineDst ? maxLineDst : maxDst;
    
    return float2(minDst,maxDst);
}

To visualise this intersection change the intersection information in our visualise shader to the following. We use the _CapHeights float2 property to set the heights of the caps.

// Intersect information
float2 intersect = intersectCappedEllipsoid(lineOrigin, lineDir, shapeOrigin, shapeForwardDir, shapeUpDir, _ShapeParams.xyz, _CapHeights);
float frontIntersection = intersect.x >= 0 ? intersect.x : intersect.y;
if(frontIntersection < maxLineDst && frontIntersection >= 0)
{
    float3 pos = lineOrigin - shapeOrigin + lineDir * frontIntersection;
    r = saturate(abs(pos.x)/_DistanceScale);
    g = saturate(abs(pos.y)/_DistanceScale);
    b = saturate(abs(pos.z)/_DistanceScale);
}
Capped ellipsoid on a 50 scale sphere object with parameters sx = 15, sy = 10, sz = 20, distance scale = 20 and cap heights 0 and 0.75

The capped ellipsoid function was used in my potions showcase, where the volume of the potion is a sphere with a single cap. The shape (and thus cap) direction was handled by a separate c# script on the gameObject that modifies shapeUpDir.

3.4 Capped half hyperboloid

Note that the name states half a hyperboloid. This is because a capped hyperboloid can have up to 4 intersections and handling that is outside the scope of this tutorial. As such we will be looking at only the lower half (y <= 0) of the hyperboloid, which has only 2 intersections.

The caps of a hyperboloid are again made of ellipses. We can calculate the height dependent ellipse scales similarly to the ellipsoid’s. This is done by taking the equation for hyperboloids from section 2.3 and setting either x or z to 0 and solve for the other. From the hyperboloid equation this results in the following:

\displaystyle \frac{x^2}{s_x^2} - \frac{y^2}{s_y^2} + \frac{z^2}{s_z^2} = s_t,

\displaystyle s_{e,x} \rightarrow x (z=0) = s_x \sqrt{s_t + \frac{h^2}{s_y^2}},

\displaystyle s_{e,z} \rightarrow z (x=0) = s_z \sqrt{s_t + \frac{h^2}{s_y^2}}.

Here s_{e,x} denotes the x-scale parameter of the ellipse, s_x the x-scale parameter of the hyperboloid, s_t the type parameter of the hyperboloid and h the height position of the cap. Because we’re only looking at the lower half of the hyperboloid we have the extra condition that h \leq 0.

The HLSL implementation is again similar to the capped cylinder and capped ellipsoid, meaning we first find the hyperboloid intersections and check if they are between the caps followed by a comparison with the two ellipses. As input we provide the line origin, line direction, shape origin, shape forward direction, shape up direction, the scale parameters s_x, s_y and s_z, type parameter s_t, and two cap heights which get clamped between -maxLineDst and 0. Add the following lines to our HLSL file below the intersectCappedEllipsoid function.

float2 intersectCappedHalfHyperboloid(float3 lineOrigin, float3 lineDir, float3 shapeOrigin, float3 shapeForwardDir, float3 shapeUpDir, float3 scaleParams, float typeParam, float2 capHeight)
{
    // Hyperboloid intersections
    float2 hyperboloidDst = intersectHyperboloid(lineOrigin, lineDir, shapeOrigin, shapeForwardDir, shapeUpDir, scaleParams, typeParam);

    // Cap heights in -maxLineDst to 0 range
    float cap1Height = clamp(min(capHeight.x, capHeight.y), -maxLineDst, 0);
    float cap2Height = clamp(max(capHeight.x, capHeight.y), -maxLineDst, 0);

    // Ellipsoid intersect positions
    float3 samplePos1 = lineOrigin + lineDir * hyperboloidDst.x;
    float3 samplePos2 = lineOrigin + lineDir * hyperboloidDst.y;

    bool isMinDstBelowCap = dot(samplePos1 - shapeOrigin, shapeUpDir) <= cap1Height;
    bool isMinDstAboveCap = dot(samplePos1 - shapeOrigin, shapeUpDir) >= cap2Height;

    bool isMaxDstBelowCap = dot(samplePos2 - shapeOrigin, shapeUpDir) <= cap1Height;
    bool isMaxDstAboveCap = dot(samplePos2 - shapeOrigin, shapeUpDir) >= cap2Height;

    float minDst = hyperboloidDst.x;
    float maxDst = hyperboloidDst.y;
    if(isMinDstBelowCap || isMinDstAboveCap)
        minDst = maxLineDst;
    if(isMaxDstBelowCap || isMaxDstAboveCap)
        maxDst = maxLineDst;
    
    // Cap positions
    float3 ellipse1Pos = shapeOrigin + cap1Height * shapeUpDir;
    float3 ellipse2Pos = shapeOrigin + cap2Height * shapeUpDir;

    // Hyperboloid cap radius
    float sySqrR = 1 / (scaleParams.y * scaleParams.y);
    float e1HeightSqr = cap1Height * cap1Height;
    float e1ScaleX = scaleParams.x * sqrt(typeParam + e1HeightSqr * sySqrR);
    float e1ScaleZ = scaleParams.z * sqrt(typeParam + e1HeightSqr * sySqrR);

    float e2HeightSqr = cap2Height * cap2Height;
    float e2ScaleX = scaleParams.x * sqrt(typeParam + e2HeightSqr * sySqrR);
    float e2ScaleZ = scaleParams.z * sqrt(typeParam + e2HeightSqr * sySqrR);

    // Ellipse bottom and top cap intersections
    float ellipse1Dst = intersectEllipse(lineOrigin, lineDir, ellipse1Pos, shapeForwardDir, shapeUpDir, float2(e1ScaleX, e1ScaleZ));
    float ellipse2Dst = intersectEllipse(lineOrigin, lineDir, ellipse2Pos, shapeForwardDir, shapeUpDir, float2(e2ScaleX, e2ScaleZ));

    // Compare ellipsoid intersects to ellipse intersects for minDst
    minDst = min(min(ellipse1Dst, ellipse2Dst), minDst);
    // If value is maxLineDst make -maxLineDst for max comparison to work
    ellipse1Dst = ellipse1Dst == maxLineDst ? -maxLineDst : ellipse1Dst;
    ellipse2Dst = ellipse2Dst == maxLineDst ? -maxLineDst : ellipse2Dst;
    maxDst = maxDst == maxLineDst ? -maxLineDst : maxDst;

    // Compare ellipsoid intersects to ellipse intersects for maxDst
    // and make sure we don't return -maxLineDst
    maxDst = max(max(ellipse1Dst, ellipse2Dst), maxDst);
    maxDst = maxDst == -maxLineDst ? maxLineDst : maxDst;
    
    return float2(minDst,maxDst);
}

To visualise this intersection change the intersection information in our visualise shader to the following.

// Intersect information
float2 intersect = intersectCappedHalfHyperboloid(lineOrigin, lineDir, shapeOrigin, shapeForwardDir, shapeUpDir, _ShapeParams.xyz, _ShapeParams.w, _CapHeights);
float frontIntersection = intersect.x >= 0 ? intersect.x : intersect.y;
if(frontIntersection < maxLineDst && frontIntersection >= 0)
{
    float3 pos = lineOrigin - shapeOrigin + lineDir * frontIntersection;
    r = saturate(abs(pos.x)/_DistanceScale);
    g = saturate(abs(pos.y)/_DistanceScale);
    b = saturate(abs(pos.z)/_DistanceScale);
}
Capped half hyperboloid on a 50 scale sphere object with parameters sx = 5, sy = 5, sz = 5, st = 0, distance scale = 15 and cap heights -3 and -15

3.5 Rectangle

The rectangle is very similar to the ellipse, in that we do a plane intersection and check if it is within the bounds of a rectangle. The rectangular surface is bounded by the following 2 conditions:

\displaystyle |x| \leq s_x,

\displaystyle |z| \leq s_z.

Here |x| stands for the absolute value of x. To check if a point is within the rectangle we substitute the x and z coordinate with our plane intersection coordinates.

The HLSL implementation for the rectangle thus starts the same as the plane. Only the end is different where check if the intersection point is within an rectangle, if it is we return the intersection, if not we return the max float. As input we provide the line origin, line direction, shape origin, shape forward direction, shape up direction, and the scale parameters s_x and s_z that define the bounds of the rectangle. Note that the rectangle’s length is 2 s_x. Add the following lines to our HLSL file below the intersectCappedHalfHyperboloid function.

float intersectRectangle(float3 lineOrigin, float3 lineDir, float3 shapeOrigin, float3 shapeForwardDir, float3 shapeUpDir, float2 scaleParams)
{
    // Transform line origin and direction from world space to the shape space
    float3x3 transitionMatrix = constructTransitionMatrix(shapeForwardDir, shapeUpDir);
    float3 lO = mul(transitionMatrix, lineOrigin - shapeOrigin);
    float3 lD = mul(transitionMatrix, lineDir);

    float denominator = lD.y;
    float numerator = lO.y;

    float pIntersect = - numerator / denominator;

    float3 samplePos = lO + lD * pIntersect;
    bool isInRectangle = abs(samplePos.x) <= scaleParams.x && abs(samplePos.z) <= scaleParams.y;
    return isInRectangle ? pIntersect : maxLineDst;
}

To visualise this intersection change the intersection information in our visualise shader to the following.

// Intersect information
float intersect = intersectRectangle(lineOrigin, lineDir, shapeOrigin, shapeForwardDir, shapeUpDir, _ShapeParams.xz);
if(intersect < maxLineDst && intersect > 0)
{
    float3 pos = lineOrigin - shapeOrigin + lineDir * intersect;
    r = saturate(abs(pos.x)/_DistanceScale);
    g = saturate(abs(pos.y)/_DistanceScale);
    b = saturate(abs(pos.z)/_DistanceScale);
}
Rectangle on a 50 scale sphere object with parameters sx = 10, sz = 20 and distance scale = 20

3.6 Box

A box simply consists of 6 rectangles that each form a face. In this section we will find the intersection of a box through comparing 6 rectangle intersections, but I want to let you know that this is not the most optimal way of finding box intersections. If you are looking for a more optimal way I’d suggest looking at the boxIntersection function by Inigo Quilez.

For our approach there isn’t really more to it than doing a rectangle intersection for the up, down, forward, backward, right and left face of the box. You only need to remind yourself what the rectangle’s up direction and origin would be, and how the scale parameters should be ordered with respect to the forward direction.

The HLSL implementation thus consists of 6 rectangle intersections, which are then compared using min and max functions. As input we provide the line origin, line direction, shape origin, shape forward direction, shape up direction, and the scale parameters s_x, s_y and s_z that define the bounds of the box. Note that the box’s length is 2 s_x. Add the following lines to our HLSL file below the intersectCappedHalfHyperboloid function.

float2 intersectBox(float3 lineOrigin, float3 lineDir, float3 shapeOrigin, float3 shapeForwardDir, float3 shapeUpDir, float3 scaleParams)
{
    // Calculate right direction
    float3 shapeRightDir = cross(shapeForwardDir, shapeUpDir);

    // Rectangle face intersections
    float faceUp = intersectRectangle(lineOrigin, lineDir, shapeOrigin + shapeUpDir * scaleParams.y, shapeForwardDir, shapeUpDir, scaleParams.xz);
    float faceDown = intersectRectangle(lineOrigin, lineDir, shapeOrigin - shapeUpDir * scaleParams.y, shapeForwardDir, shapeUpDir, scaleParams.xz);

    float faceForward = intersectRectangle(lineOrigin, lineDir, shapeOrigin + shapeForwardDir * scaleParams.z, -shapeUpDir, shapeForwardDir, scaleParams.xy);
    float faceBackward = intersectRectangle(lineOrigin, lineDir, shapeOrigin - shapeForwardDir * scaleParams.z, -shapeUpDir, shapeForwardDir, scaleParams.xy);

    float faceRight = intersectRectangle(lineOrigin, lineDir, shapeOrigin + shapeRightDir * scaleParams.x, shapeForwardDir, shapeRightDir, scaleParams.yz);
    float faceLeft = intersectRectangle(lineOrigin, lineDir, shapeOrigin - shapeRightDir * scaleParams.x, shapeForwardDir, shapeRightDir, scaleParams.yz);

    // Min face comparison
    float minDst = min(faceUp,min(faceDown,min(faceForward,min(faceBackward, min(faceRight, faceLeft)))));
    
    // If value is maxLineDst make -maxLineDst for max comparison to work
    faceUp *= faceUp == maxLineDst ? -1 : 1;
    faceDown *= faceDown == maxLineDst ? -1 : 1;
    faceForward *= faceForward == maxLineDst ? -1 : 1;
    faceBackward *= faceBackward == maxLineDst ? -1 : 1;
    faceRight *= faceRight == maxLineDst ? -1 : 1;
    faceLeft *= faceLeft == maxLineDst ? -1 : 1;

    // Maxn face comparison and make sure we don't return -maxLineDst
    float maxDst = max(faceUp,max(faceDown,max(faceForward,max(faceBackward, max(faceRight, faceLeft)))));
    maxDst *= maxDst == -maxLineDst ? -1 : 1;

    return float2(minDst, maxDst);
}

To visualise this intersection change the intersection information in our visualise shader to the following.

// Intersect information
float2 intersect = intersectBox(lineOrigin, lineDir, shapeOrigin, shapeForwardDir, shapeUpDir, _ShapeParams.xyz);
float frontIntersection = intersect.x >= 0 ? intersect.x : intersect.y;
if(frontIntersection < maxLineDst && frontIntersection >= 0)
{
    float3 pos = lineOrigin - shapeOrigin + lineDir * frontIntersection;
    r = saturate(abs(pos.x)/_DistanceScale);
    g = saturate(abs(pos.y)/_DistanceScale);
    b = saturate(abs(pos.z)/_DistanceScale);
}
Box on a 50 scale sphere object with parameters sx = 10, sy = 5, sz = 20 and distance scale = 20

4 Conclusion

While there are still more combinations to add to this list I think 11 shape functions is a good start. This tutorial in of itself is obviously not that exiting from a visual perspective, but these functions can be used to create more interesting shaders like my black hole shader or potion shader. Be sure to show me your shaders that implement these functions by tagging me on twitter (@kelvinvanhoorn) or on reddit (u/Radagasd). Also, while not necessary I’ll very much appreciate it if you credit me when using these functions in your projects.

I do want to emphasise that while these functions give the correct answer, they aren’t necessarily the most optimal way (computation-wise) of getting there. If you’re interested you could take a look at Inigo Quilez’s functions here that covers most functions in a more optimal way.

Thank you very much for reading and I hope you learned something new. If you want to support me financially you can do so using my ko-fi page.

References

  1. Plane equation from Wikipedia’s page on planes.
  2. Images made by Sam Derbyshire for Wikipedia’s page on quadric surfaces.
  3. Equations for real quadric surfaces in Euclidean space from Wikipedia’s page on quadrics.
  4. Overview of finding the ellipses on an ellipsoid by Keisan

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: