In a previous article I've written about how to implement a physically based rendering engine. For that engine I was using Signed Distance Fields, Constructive Solid Geometry and Ray Marching. These three components form together a convenient and powerful way to represent and manipulate many 2D/3D geometric shapes, not only for rendering applications. So I thought it would be worth the effort to make out of my previous work a standalone module, easy to reuse in various projects, and of course share it here. Even if there are tons of material about these topics on the web, I couldn't find a ready to use module in C. Now there will be mine, at least.

My module is available for download here, or you can get it from the command line as follow: wget https://baillehachepascal.dev/2023/Data/SdfCsg/sdfcsg.0.2.0.tar.gz, and decompress it using tar xvf sdfcsg.0.2.0.tar.gz. This will create a subdirectory SdfCsg, inside which executing the command make demo should compile and run the demo below (click to play). Note that it automatically download and install locally my X11Display module.

The module consists of two files: sdfcsg.h and sdfcsg.c. To use it, include the header in your project #include "sdfcsg.h", and compile for example as follow gcc -o main main.c sdfcsg.c -O3 -lm.

All the functionalities of the module are available through an opaque structure SdfCsg, which you can create with functions, one per primitive or CSG operator.

SdfCsg* SdfCsgBox( sdfCsgReal_t const width, sdfCsgReal_t const height, sdfCsgReal_t const depth); SdfCsg* SdfCsgPlane( sdfCsgReal_t const normX, sdfCsgReal_t const normY, sdfCsgReal_t const normZ, sdfCsgReal_t const dist); SdfCsg* SdfCsgSphere(sdfCsgReal_t const radius); SdfCsg* SdfCsgCone( sdfCsgReal_t const baseRadius, sdfCsgReal_t const topRadius, sdfCsgReal_t const height); SdfCsg* SdfCsgEllipsoid( sdfCsgReal_t const radiusX, sdfCsgReal_t const radiusY, sdfCsgReal_t const radiusZ); SdfCsg* SdfCsgCapsule( sdfCsgReal_t const radius, sdfCsgReal_t const length); SdfCsg* SdfCsgTorus( sdfCsgReal_t const majorRadius, sdfCsgReal_t const minorRadius); SdfCsg* SdfCsgUnion( SdfCsg* const shapeA, SdfCsg* const shapeB); SdfCsg* SdfCsgIntersection( SdfCsg* const shapeA, SdfCsg* const shapeB); SdfCsg* SdfCsgDifference( SdfCsg* const shapeA, SdfCsg* const shapeB);

A user-defined SDF function can also be used to extend the module with other shapes. The SdfCsg structure also has a void* hook in case the user needs to attach data to shape instances (extra parameters for the user-defined SDF function for example).

SdfCsg* SdfCsgUserDefined(SdfCsgFun const sdf);

To release the memory used by a shape, one can use the dedicated function:

void SdfCsgFree(SdfCsg** const that);

Transformations can be applied to shapes:

void SdfCsgTransformNone( SdfCsg* const that, unsigned int const rank); void SdfCsgTransformScale( SdfCsg* const that, unsigned int const rank, sdfCsgReal_t const scale); void SdfCsgTransformTranslate( SdfCsg* const that, unsigned int const rank, sdfCsgReal_t const x, sdfCsgReal_t const y, sdfCsgReal_t const z); void SdfCsgTransformRotate( SdfCsg* const that, unsigned int const rank, sdfCsgReal_t const quat[static 4]); void SdfCsgTransformExtrusion( SdfCsg* const that, unsigned int const rank, sdfCsgReal_t const length); void SdfCsgTransformSymmetry( SdfCsg* const that, unsigned int const rank); void SdfCsgTransformRevolution( SdfCsg* const that, unsigned int const rank); void SdfCsgTransformTwist( SdfCsg* const that, unsigned int const rank, sdfCsgReal_t const freq);

Setters and getters allow the user to access and modify shapes' properties after creation.

void SdfCsgGetSphereParams( SdfCsg const* const that, sdfCsgReal_t* const radius); void SdfCsgGetCapsuleParams( SdfCsg const* const that, sdfCsgReal_t* const radius, sdfCsgReal_t* const length); void SdfCsgGetTorusParams( SdfCsg const* const that, sdfCsgReal_t* const majorRadius, sdfCsgReal_t* const minorRadius); void SdfCsgGetEllipsoidParams( SdfCsg const* const that, sdfCsgReal_t* const radiusX, sdfCsgReal_t* const radiusY, sdfCsgReal_t* const radiusZ); void SdfCsgGetBoxParams( SdfCsg const* const that, sdfCsgReal_t* const width, sdfCsgReal_t* const height, sdfCsgReal_t* const depth); void SdfCsgGetConeParams( SdfCsg const* const that, sdfCsgReal_t* const baseRadius, sdfCsgReal_t* const topRadius, sdfCsgReal_t* const height); void SdfCsgGetPlaneParams( SdfCsg const* const that, sdfCsgReal_t* const normX, sdfCsgReal_t* const normY, sdfCsgReal_t* const normZ, sdfCsgReal_t* const dist); void SdfCsgSetSphereParams( SdfCsg* const that, sdfCsgReal_t const radius); void SdfCsgSetCapsuleParams( SdfCsg* const that, sdfCsgReal_t const radius, sdfCsgReal_t const length); void SdfCsgSetTorusParams( SdfCsg* const that, sdfCsgReal_t const majorRadius, sdfCsgReal_t const minorRadius); void SdfCsgSetEllipsoidParams( SdfCsg* const that, sdfCsgReal_t const radiusX, sdfCsgReal_t const radiusY, sdfCsgReal_t const radiusZ); void SdfCsgSetBoxParams( SdfCsg* const that, sdfCsgReal_t const width, sdfCsgReal_t const height, sdfCsgReal_t const depth); void SdfCsgSetConeParams( SdfCsg* const that, sdfCsgReal_t const baseRadius, sdfCsgReal_t const topRadius, sdfCsgReal_t const height); void SdfCsgSetPlaneParams( SdfCsg* const that, sdfCsgReal_t const normX, sdfCsgReal_t const normY, sdfCsgReal_t const normZ, sdfCsgReal_t const dist); sdfCsgReal_t SdfCsgGetThickness(SdfCsg const* const that); sdfCsgReal_t SdfCsgGetRoundness(SdfCsg const* const that); sdfCsgReal_t SdfCsgGetSmoothness(SdfCsg const* const that); void* SdfCsgGetData(SdfCsg const* const that); SdfCsg* SdfCsgGetParent(SdfCsg const* const that); SdfCsg* SdfCsgGetChildA(SdfCsg const* const that); SdfCsg* SdfCsgGetChildB(SdfCsg const* const that); void SdfCsgSetThickness( SdfCsg* const that, sdfCsgReal_t const thickness); void SdfCsgSetRoundness( SdfCsg* const that, sdfCsgReal_t const roundness); void SdfCsgSetSmoothness( SdfCsg* const that, sdfCsgReal_t const smoothness); void SdfCsgSetData( SdfCsg* const that, void* const data);

The SDF of a shape, and its gradient, can then be evaluated with:

sdfCsgReal_t SdfCsgDist( SdfCsg* const that, sdfCsgReal_t const pos[static 3]); void SdfCsgGradient( SdfCsg* const that, sdfCsgReal_t const pos[static 3], sdfCsgReal_t grad[static 3]);

And ray marching can be performed with:

void SdfCsgRayMarch( SdfCsg* const that, sdfCsgReal_t const from[static 3], sdfCsgReal_t const dir[static 3], SdfCsgRayMarchRes* const res);

The result structure of ray marching looks as follow

typedef struct SdfCsgRayMarchRes { // Flag to memorise if the ray has hit a shape bool hasHit; // Flag to memorise if the ray marching reached the maximum number of steps bool hasGaveUp; // Reference to the shape of the nearest intersection SdfCsg* shape; // Signed distance from the ray origin to the intersection sdfCsgReal_t rayDist; // Coordinates in world coordinate system of the intersection sdfCsgReal_t pos[3]; // Coordinates in the shape local coordinate system of the intersection sdfCsgReal_t posLocal[3]; // Signed distance at the intersection sdfCsgReal_t posDist; // Shortest observed signed distance between the ray and a shape during the // steps of ray marching sdfCsgReal_t minDist; } SdfCsgRayMarchRes; SdfCsgRayMarchRes SdfCsgRayMarchResCreate(sdfCsgReal_t const distMax);

The implementation of ray marching is here slightly different from the one introduced in my previous article (marching in world coordinate system instead of local coordinate systems to make things simpler and hopefully more efficient and better handling numerical imprecision problems).

Two functions to find the nearest position guaranteed to be inside/outside a shape from a given position can be conveniently used to handle the undetermination of the insideness or outsideness of the result of ray marching. This is particularly useful when computing reflection, refraction and diffraction.

bool SdfCsgFindNearestInside( SdfCsg* const that, sdfCsgReal_t const from[static 3], sdfCsgReal_t res[static 3]); bool SdfCsgFindNearestOutside( SdfCsg* const that, sdfCsgReal_t const from[static 3], sdfCsgReal_t res[static 3]);

As an example of use, the basis for a 3d renderer using this module would look like:

Render(): int const nbShape = ... shapes: array of nbShape SdfCsg shapes[0] = SdfCsgBox(...) ... loop on pixels: sdfCsgReal_t from[3] = {...} sdfCsgReal_t dir[3] = {...} ResolveRay(from, dir, shapes, nbShape, pixel.rgb) ResolveRay(from, dir, shapes, nbShape, rgb): maxDist = ... SdfCsgRayMarchRes res = SdfCsgRayMarchResCreate(maxDist) loop on shapes: SdfCsgRayMarch(shape, from, dir, res) if res.shape == NULL: rgb[0] = rgb[1] = rgb[2] = 0 else: sdfCsgReal_t dirLight[3] = posLight - res.pos sdfCsgReal_t gradient[3] SdfCsgGradient(res.shape, res.pos, gradient) sdfCsgReal_t dot = gradient dot dirLight sdfCsgReal_t intensity = dot if computing shadow: sdfCsgReal_t nextPos[3] SdfCsgFindNearestOutside(res.shape, res.pos, nextPos) SdfCsgRayMarchRes subRes = SdfCsgRayMarchResCreate(maxDist) loop on shapes: SdfCsgRayMarch(shape, nextPos, dirLight, subRes) if subRes.shape != light: intensity = 0.0 rgb = intensity * (res.shape's color at res.posLocal) if computing reflectivity: sdfCsgReal_t nextPos[3] SdfCsgFindNearestOutside(res.shape, res.pos, nextPos) sdfCsgReal_t dirReflected[3] = reflected direction given dir and gradient unsigned char nextRgb[3] ResolveRay(nextPos, dirReflected, shapes, nbShape, nextRgb) rgb = combination of rgb and nextRgb given the reflectivity if computing transparency: sdfCsgReal_t nextPos[3] SdfCsgFindNearestInside(res.shape, res.pos, nextPos) sdfCsgReal_t dirInside[3] = refracted direction given dir and gradient SdfCsgRayMarchRes resInside = SdfCsgRayMarchResCreate(distMax) SdfCsgRayMarch(res.shape, nextPos, -dirInside, resInside) SdfCsgFindNearestOutside(res.shape, resInside.pos, nextPos) sdfCsgReal_t gradientOut[3] SdfCsgGradient(res.shape, res.pos, gradientOut) sdfCsgReal_t dirOutside[3] = ray direction at exit of the shape given dirInside and gradientOut unsigned char nextRgb[3] ResolveRay(nextPos, dirOutside, shapes, nbShape, nextRgb) rgb = combination of rgb and nextRgb given the transparency

The image below gives a basic example of 3D rendering several shapes, with reflection, texturing and fake soft shadow.

Finally the complete header sdfcsg.h:

// ---------------------------- sdfcsg.h --------------------------- /* SdfCsg - a C library implementing signed distance fields, ray marching and constructive solid geometry. Copyright (C) 2023 Pascal Baillehache baillehache.pascal@gmail.com https://baillehachepascal.dev This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see <http://www.gnu.org/licenses/>. */ #ifndef SDFCSG_H #define SDFCSG_H #include <stdbool.h> // Type used for real values (1 for 'float', 2 for 'double') #define SDFCSG_REAL_TYPE 2 #if SDFCSG_REAL_TYPE==1 typedef float sdfCsgReal_t; #elif SDFCSG_REAL_TYPE==2 typedef double sdfCsgReal_t; #endif // Getter and setter for the epsilon value used during calculation. // Default is 1e-6. sdfCsgReal_t SdfCsgGetEpsilon(void); void SdfCsgSetEpsilon(sdfCsgReal_t const epsilon); // Getter and setter for the maximum number of step (to avoid infinite loop) // during ray marching. Default is 1000. unsigned int SdfCsgGetRayMarchNbMaxStep(void); void SdfCsgSetRayMarchNbMaxStep(unsigned int const nb); // Maximum number of transformation applicable to one SdfCsgShape #define SDFCSG_MAX_NB_TRANSFORM 10 // Declaration of the SdfCsg opaque structure typedef struct SdfCsg SdfCsg; // Signed distance field function type declaration // These functions return the approximated distance from a given position to // the nearest point on the shape surface. // Inputs: // that: the SdfCsg // pos: the given position, in the shape's local coordinate system // Output: // Return the distance, guaranted to be less than or equal to the exact // distance. A positive distance indicates a position outside the shape. A // negative distance indicates a position inside the shape. typedef sdfCsgReal_t (*SdfCsgFun)( SdfCsg* const that, sdfCsgReal_t const pos[static 3]); // Functions to create primitives shapes // Box centered at the origin. // Plane defined by its normal and the its distance to the origin. // Sphere centered at the origin. // Cone centered at the origin and aligned with the y axis. // Ellipsoid centered at the origin. // Capsule centered at the origin and aligned with the y axis. // Torus centered at the origin and lying in the (y=0) plane. SdfCsg* SdfCsgBox( sdfCsgReal_t const width, sdfCsgReal_t const height, sdfCsgReal_t const depth); SdfCsg* SdfCsgPlane( sdfCsgReal_t const normX, sdfCsgReal_t const normY, sdfCsgReal_t const normZ, sdfCsgReal_t const dist); SdfCsg* SdfCsgSphere(sdfCsgReal_t const radius); SdfCsg* SdfCsgCone( sdfCsgReal_t const baseRadius, sdfCsgReal_t const topRadius, sdfCsgReal_t const height); SdfCsg* SdfCsgEllipsoid( sdfCsgReal_t const radiusX, sdfCsgReal_t const radiusY, sdfCsgReal_t const radiusZ); SdfCsg* SdfCsgCapsule( sdfCsgReal_t const radius, sdfCsgReal_t const length); SdfCsg* SdfCsgTorus( sdfCsgReal_t const majorRadius, sdfCsgReal_t const minorRadius); // Functions to create CSG SdfCsg* SdfCsgUnion( SdfCsg* const shapeA, SdfCsg* const shapeB); SdfCsg* SdfCsgIntersection( SdfCsg* const shapeA, SdfCsg* const shapeB); SdfCsg* SdfCsgDifference( SdfCsg* const shapeA, SdfCsg* const shapeB); // Function to create a shape with a user defined SDF function // The argument 'sdf' is the SDF function for the canonical shape (ie after // unapplying the eventual transformations of the shape) SdfCsg* SdfCsgUserDefined(SdfCsgFun const sdf); // Free the memory used by a SdfCsg. 'that' is set to NULL after releasing the // memory. If 'that' is equal to NULL nothing happen. If the shape is a CSG, // its components are freed (recursively). void SdfCsgFree(SdfCsg** const that); // Setters and getters for the SdfCsg structure's properties // Thickness: thickness coefficient (if greater than epsilon, create a shell // out of the shape surface) // Roundness: rounding coefficient (enlarge the shape, rounding the corner) // Data: reference to optional user data // ChildA/childB: CSG components // Smoothness: smoothing coefficient for CSG (makes the CSG components blend // smoothly within each other) sdfCsgReal_t SdfCsgGetThickness(SdfCsg const* const that); sdfCsgReal_t SdfCsgGetRoundness(SdfCsg const* const that); sdfCsgReal_t SdfCsgGetSmoothness(SdfCsg const* const that); void* SdfCsgGetData(SdfCsg const* const that); SdfCsg* SdfCsgGetParent(SdfCsg const* const that); SdfCsg* SdfCsgGetChildA(SdfCsg const* const that); SdfCsg* SdfCsgGetChildB(SdfCsg const* const that); void SdfCsgSetThickness( SdfCsg* const that, sdfCsgReal_t const thickness); void SdfCsgSetRoundness( SdfCsg* const that, sdfCsgReal_t const roundness); void SdfCsgSetSmoothness( SdfCsg* const that, sdfCsgReal_t const smoothness); void SdfCsgSetData( SdfCsg* const that, void* const data); // Setters and getters for the SdfCsg primitive shapes' parameters void SdfCsgGetSphereParams( SdfCsg const* const that, sdfCsgReal_t* const radius); void SdfCsgGetCapsuleParams( SdfCsg const* const that, sdfCsgReal_t* const radius, sdfCsgReal_t* const length); void SdfCsgGetTorusParams( SdfCsg const* const that, sdfCsgReal_t* const majorRadius, sdfCsgReal_t* const minorRadius); void SdfCsgGetEllipsoidParams( SdfCsg const* const that, sdfCsgReal_t* const radiusX, sdfCsgReal_t* const radiusY, sdfCsgReal_t* const radiusZ); void SdfCsgGetBoxParams( SdfCsg const* const that, sdfCsgReal_t* const width, sdfCsgReal_t* const height, sdfCsgReal_t* const depth); void SdfCsgGetConeParams( SdfCsg const* const that, sdfCsgReal_t* const baseRadius, sdfCsgReal_t* const topRadius, sdfCsgReal_t* const height); void SdfCsgGetPlaneParams( SdfCsg const* const that, sdfCsgReal_t* const normX, sdfCsgReal_t* const normY, sdfCsgReal_t* const normZ, sdfCsgReal_t* const dist); void SdfCsgSetSphereParams( SdfCsg* const that, sdfCsgReal_t const radius); void SdfCsgSetCapsuleParams( SdfCsg* const that, sdfCsgReal_t const radius, sdfCsgReal_t const length); void SdfCsgSetTorusParams( SdfCsg* const that, sdfCsgReal_t const majorRadius, sdfCsgReal_t const minorRadius); void SdfCsgSetEllipsoidParams( SdfCsg* const that, sdfCsgReal_t const radiusX, sdfCsgReal_t const radiusY, sdfCsgReal_t const radiusZ); void SdfCsgSetBoxParams( SdfCsg* const that, sdfCsgReal_t const width, sdfCsgReal_t const height, sdfCsgReal_t const depth); void SdfCsgSetConeParams( SdfCsg* const that, sdfCsgReal_t const baseRadius, sdfCsgReal_t const topRadius, sdfCsgReal_t const height); void SdfCsgSetPlaneParams( SdfCsg* const that, sdfCsgReal_t const normX, sdfCsgReal_t const normY, sdfCsgReal_t const normZ, sdfCsgReal_t const dist); // Helper function to create a quaternion from a rotation axis and an angle (in // radians, rotation clockwise) void SdfCsgQuaternionCreate( sdfCsgReal_t const axis[static 3], sdfCsgReal_t const theta, sdfCsgReal_t quat[static 4]); // Set a transformation to a SdfCsg. The transformation are applied in // increasing order defined by 'rank'. SdfCsgTransformNone() can be used to // cancel the transformation at a given rank. It's fine to have no // transformation between two others, like: translation at rank 0, nothing at // rank 1, rotation at rank 2. // Rotations are specified using quaternion. // Extrusions occur along the z axis relative to the origin. // Symmetries are relative to the (x=0) plane. The half-space (x>0) is // duplicated in the (x<0) half-space. // Revolutions occur as follow: the intersection of the shape and (z=0) // plane is expanded radially 360 degrees around the y axis. // Twists twist the shape clockwise along the y axis at a frequency of one // rotation per 'freq' unit along the y axis. void SdfCsgTransformNone( SdfCsg* const that, unsigned int const rank); void SdfCsgTransformScale( SdfCsg* const that, unsigned int const rank, sdfCsgReal_t const scale); void SdfCsgTransformTranslate( SdfCsg* const that, unsigned int const rank, sdfCsgReal_t const x, sdfCsgReal_t const y, sdfCsgReal_t const z); void SdfCsgTransformRotate( SdfCsg* const that, unsigned int const rank, sdfCsgReal_t const quat[static 4]); void SdfCsgTransformExtrusion( SdfCsg* const that, unsigned int const rank, sdfCsgReal_t const length); void SdfCsgTransformSymmetry( SdfCsg* const that, unsigned int const rank); void SdfCsgTransformRevolution( SdfCsg* const that, unsigned int const rank); void SdfCsgTransformTwist( SdfCsg* const that, unsigned int const rank, sdfCsgReal_t const freq); // Evaluate the SDF of a shape at a given position. // Inputs: // that: the SdfCsg // pos: the given position in world coordinate system (ie the shape's // transformations are unapplied to the given position and the shape's // canonical SDF function is evaluated at the resulting position) // Output: // Return the SDF value at the requested position. sdfCsgReal_t SdfCsgDist( SdfCsg* const that, sdfCsgReal_t const pos[static 3]); // Evaluate the gradient in the SDF at a given position. // Inputs: // that: the SdfCsg // pos: the given position in world coordinate system (ie the shape's // transformations are unapplied to the given position and the shape's // canonical SDF function is evaluated at the resulting position) // grad: the resulting gradient // Output: // 'grad' is updated with the gradient. If 'pos' is on the surface of the // shape, it's equaivalent to the normal to the surface at that position, // oriented outward the surface of the shape. 'grad' is not normalised. The // symmeric difference quotient is used to estimate the gradient. void SdfCsgGradient( SdfCsg* const that, sdfCsgReal_t const pos[static 3], sdfCsgReal_t grad[static 3]); // Structure to memorise the result of SdfCsgRayMarch(). Cf that function for // details about the fields. typedef struct SdfCsgRayMarchRes { // Flag to memorise if the ray has hit a shape bool hasHit; // Flag to memorise if the ray marching reached the maximum number of steps bool hasGaveUp; // Reference to the shape of the nearest intersection SdfCsg* shape; // Signed distance from the ray origin to the intersection sdfCsgReal_t rayDist; // Coordinates in world coordinate system of the intersection sdfCsgReal_t pos[3]; // Coordinates in the shape local coordinate system of the intersection sdfCsgReal_t posLocal[3]; // Signed distance at the intersection sdfCsgReal_t posDist; // Shortest observed signed distance between the ray and a shape during the // steps of ray marching sdfCsgReal_t minDist; } SdfCsgRayMarchRes; // Create a SdfCsgRayMarchRes to be used as argument of SdfCsgRayMarch() // Input: // distMax: the maximum distance (in absolute value) from the origin of a ray // the ray marching will search for an intersection // Output: // Return a SdfCsgRayMarchRes ready to use by SdfCsgRayMarch(). // SdfCsgRayMarchRes's fields do not require memory freeing. SdfCsgRayMarchRes SdfCsgRayMarchResCreate(sdfCsgReal_t const distMax); // Ray marching for a given ray and SdfCsg // Inputs: // that: the SdfCsg // from: origin of the ray (in world coordinate system) // dir: direction of the ray (must be normalised, in world coordinate system) // res: result of the ray marching // Output: // Step iteratively along the 'dir' direction from the 'from' position until // the intersection of the ray and the surface of the shape is reached or // missed. Note that due to numerical imprecision 'res.pos' is only an // approximation, within 'sdfCsg_epsilon' in world coordinate sytem, of the // intersection. Whether it is inside or outside the shape is undefined. The // user should check the sign of 'res.distPos' and use the result // appropriately. See also the SdfCsgFindNearestInside() and // SdfCsgFindNearestOutside() functions. // The value of 'res.rayDist' is used as an upper bound for the total marched // distance during iteration. Beyond that bound, the ray is considered to // have missed the shape. When ray marching several SdfCsg, 'res' should // be created using SdfCsgRayMarchResCreate(<world size>), and // SdfCsgRayMarch() should be called with the same 'res' on all relevant // SdfCsg to ensure better performances (the algorithm doesn't search farther // than the current nearest intersection). If several ray marchings are // multithreaded, the user must use one SdfCsgRayMarchRes per thread and // combine their result when/as needed. // If the ray marching reaches the shape surface, the fields of 'res' are // updated accordingly. Else, 'res' is left unchanged except for 'res.minDist' // which is always updated as needed. 'res.minDist' can be used to // approximate soft shadows for example. If the ray marching failed to // converge (due to numerical imprecision or ill-defined SDF), // 'res.hasGaveUp' is set to true, else it is set to false. // 'from' is assumed to be outside the shape. If ray marching from inside a // shape, 'dir' must be the opposite of the actual stepping direction to ray // march backward, from inside to outside the shape. void SdfCsgRayMarch( SdfCsg* const that, sdfCsgReal_t const from[static 3], sdfCsgReal_t const dir[static 3], SdfCsgRayMarchRes* const res); // Find the nearest "inside" and "outside" position in the SDF from the given // position 'pos'. // Inputs: // that: the SdfCsg // from: the given position, in world coordinate system // res: the result position, in world coordinate system // Output: // Return true and update 'res' if the nearest position could be found, else // return false and 'res' is undefined. The search is performed using // gradient descent/ascent (for inside/outside resp.) from the given position // until a position where the SDF absolute value is greater than // 'sdfCsg_epsilon'. The search fails if such a position couldn't be found // within 'sdfCsg_rayMarchNbMaxStep' steps. bool SdfCsgFindNearestInside( SdfCsg* const that, sdfCsgReal_t const from[static 3], sdfCsgReal_t res[static 3]); bool SdfCsgFindNearestOutside( SdfCsg* const that, sdfCsgReal_t const from[static 3], sdfCsgReal_t res[static 3]); #endif

and body sdfcsg.c:

// ---------------------------- sdfcsg.c --------------------------- /* SdfCsg - a C library implementing signed distance fields, ray marching and constructive solid geometry. Copyright (C) 2023 Pascal Baillehache baillehache.pascal@gmail.com https://baillehachepascal.dev This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see <http://www.gnu.org/licenses/>. */ #include <stdlib.h> #include <math.h> #include <assert.h> #include "sdfcsg.h" // Helper macros for operations on real values #if SDFCSG_REAL_TYPE==1 #define rmin(x, y) fminf(x, y) #define rmax(x, y) fmaxf(x, y) #define rabs(x) fabsf(x) #define rsqrt(x) sqrtf(x) #define M_2PI 6.283185307179586f #elif SDFCSG_REAL_TYPE==2 #define rmin(x, y) fmin(x, y) #define rmax(x, y) fmax(x, y) #define rabs(x) fabs(x) #define rsqrt(x) sqrt(x) #define M_2PI 6.283185307179586 #endif #define rabs3(x, r) r[0] = rabs(x[0]); r[1] = rabs(x[1]); r[2] = rabs(x[2]); #define rhypot3(x) rsqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]) #define rmin3(x, y, r) \ r[0] = rmin(x[0], y[0]); r[1] = rmin(x[1], y[1]); r[2] = rmin(x[2], y[2]); #define rmax3(x, y, r) \ r[0] = rmax(x[0], y[0]); r[1] = rmax(x[1], y[1]); r[2] = rmax(x[2], y[2]); // Epsilon value used during calculation. _Thread_local static sdfCsgReal_t sdfCsg_epsilon = 1e-6; _Thread_local static sdfCsgReal_t sdfCsg_inv_2epsilon = 2e6; sdfCsgReal_t SdfCsgGetEpsilon(void) { return sdfCsg_epsilon; } void SdfCsgSetEpsilon(sdfCsgReal_t const epsilon) { sdfCsg_epsilon = epsilon; sdfCsg_inv_2epsilon = 1.0 / (2.0 * epsilon); } // Maximum number of step (to avoid infinite loop) during ray marching. _Thread_local static unsigned int sdfCsg_rayMarchNbMaxStep = 1000; unsigned int SdfCsgGetRayMarchNbMaxStep(void) { return sdfCsg_rayMarchNbMaxStep; } void SdfCsgSetRayMarchNbMaxStep(unsigned int const nb) { sdfCsg_rayMarchNbMaxStep = nb; } // Macro equivalent to a for loop iterating from 0 to ('nbIter'-1) // and storing the current iteration index in 'varName' #define loop(varName, nbIter) \ for(__typeof__((nbIter) + 0) (varName) = 0; \ (varName) < (nbIter); ++(varName)) // Types of transformations applied to a SdfCsg typedef enum SdfCsgTransformType { sdfCsg_none, sdfCsg_translate, sdfCsg_scale, sdfCsg_rotate, sdfCsg_extrusion, sdfCsg_symmetry, sdfCsg_revolution, sdfCsg_twist, } SdfCsgTransformType; // Structure to memorise a transformation applied to a SdfCsg typedef struct SdfCsgTransform { SdfCsgTransformType type; // For translation, vals[0..2] is the translation vector. // For scaling, vals[0] is the scaling coefficient. // For rotation, vals[0..2] is the rotation axis and vals[3] is the angle // in radians. // For extrusion, vals[0] is the extrusion length. // For twist, vals[0] is the frequency of the twist. sdfCsgReal_t vals[4]; } SdfCsgTransform; // Types of shape (primitives, CSGs, user defined) typedef enum SdfCsgType { sdfCsg_box, sdfCsg_sphere, sdfCsg_torus, sdfCsg_capsule, sdfCsg_ellipsoid, sdfCsg_cone, sdfCsg_plane, sdfCsg_union, sdfCsg_intersection, sdfCsg_difference, sdfCsg_user, } SdfCsgType; // Indices of primitive shape's parameters typedef enum SdfCsgParam { sdfCsg_radius = 0, sdfCsg_halfLength, sdfCsg_majorRadius = 0, sdfCsg_minorRadius, sdfCsg_radiusX = 0, sdfCsg_radiusY, sdfCsg_radiusZ, sdfCsg_halfWidth = 0, sdfCsg_halfHeight, sdfCsg_halfDepth, sdfCsg_baseRadius = 0, sdfCsg_topRadius, sdfCsg_height, sdfCsg_normX = 0, sdfCsg_normY, sdfCsg_normZ, sdfCsg_dist, sdfCsg_nbParam = 4 } SdfCsgParam; // Definition of the SdfCsg opaque structure struct SdfCsg { // Type of shape SdfCsgType type; // Transformations applied to the shape SdfCsgTransform transforms[SDFCSG_MAX_NB_TRANSFORM]; // Thickness coefficient (if greater than epsilon, create a shell out of // the shape surface) sdfCsgReal_t thickness; // Rounding coefficient (enlarge the shape, rounding the corner) sdfCsgReal_t roundness; // Paramters value for primitive sdfCsgReal_t params[sdfCsg_nbParam]; // Reference to optional user data void* data; // CSG components SdfCsg* comps[2]; // Parent shape in case of CSG SdfCsg* parent; // Smoothing coefficient for CSG (makes the CSG components blend smoothly // within each other) sdfCsgReal_t smoothness; // Canonical SDF (not accounting for the eventual transformations applied to // the shape). SdfCsgFun distCanonical; }; void SdfCsgFree(SdfCsg** const that) { if(that == NULL || *that == NULL) return; loop(i, 2) SdfCsgFree((*that)->comps + i); free(*that); *that = NULL; } sdfCsgReal_t SdfCsgGetThickness(SdfCsg const* const that) { return that->thickness; } sdfCsgReal_t SdfCsgGetRoundness(SdfCsg const* const that) { return that->roundness; } sdfCsgReal_t SdfCsgGetSmoothness(SdfCsg const* const that) { return that->smoothness; } void* SdfCsgGetData(SdfCsg const* const that) { return that->data; } SdfCsg* SdfCsgGetParent(SdfCsg const* const that) { return that->parent; } SdfCsg* SdfCsgGetChildA(SdfCsg const* const that) { return that->comps[0]; } SdfCsg* SdfCsgGetChildB(SdfCsg const* const that) { return that->comps[1]; } void SdfCsgSetThickness( SdfCsg* const that, sdfCsgReal_t const thickness) { if(thickness >= 0.0) that->thickness = thickness; } void SdfCsgSetRoundness( SdfCsg* const that, sdfCsgReal_t const roundness) { if(roundness >= 0.0) that->roundness = roundness; } void SdfCsgSetSmoothness( SdfCsg* const that, sdfCsgReal_t const smoothness) { if(smoothness >= 0.0) that->smoothness = smoothness; } void SdfCsgSetData( SdfCsg* const that, void* const data) { that->data = data; } void SdfCsgGetSphereParams( SdfCsg const* const that, sdfCsgReal_t* const radius) { assert(that->type == sdfCsg_sphere); *radius = that->params[sdfCsg_radius]; } void SdfCsgGetCapsuleParams( SdfCsg const* const that, sdfCsgReal_t* const radius, sdfCsgReal_t* const length) { assert(that->type == sdfCsg_capsule); *radius = that->params[sdfCsg_radius]; *length = 2.0 * that->params[sdfCsg_halfLength]; } void SdfCsgGetTorusParams( SdfCsg const* const that, sdfCsgReal_t* const majorRadius, sdfCsgReal_t* const minorRadius) { assert(that->type == sdfCsg_torus); *minorRadius = that->params[sdfCsg_minorRadius]; *majorRadius = that->params[sdfCsg_majorRadius]; } void SdfCsgGetEllipsoidParams( SdfCsg const* const that, sdfCsgReal_t* const radiusX, sdfCsgReal_t* const radiusY, sdfCsgReal_t* const radiusZ) { assert(that->type == sdfCsg_ellipsoid); *radiusX = that->params[sdfCsg_radiusX]; *radiusY = that->params[sdfCsg_radiusY]; *radiusZ = that->params[sdfCsg_radiusZ]; } void SdfCsgGetBoxParams( SdfCsg const* const that, sdfCsgReal_t* const width, sdfCsgReal_t* const height, sdfCsgReal_t* const depth) { assert(that->type == sdfCsg_box); *width = 2.0 * that->params[sdfCsg_halfWidth]; *height = 2.0 * that->params[sdfCsg_halfHeight]; *depth = 2.0 * that->params[sdfCsg_halfDepth]; } void SdfCsgGetConeParams( SdfCsg const* const that, sdfCsgReal_t* const baseRadius, sdfCsgReal_t* const topRadius, sdfCsgReal_t* const height) { assert(that->type == sdfCsg_cone); *baseRadius = that->params[sdfCsg_baseRadius]; *topRadius = that->params[sdfCsg_topRadius]; *height = that->params[sdfCsg_height] * 2.0; } void SdfCsgGetPlaneParams( SdfCsg const* const that, sdfCsgReal_t* const normX, sdfCsgReal_t* const normY, sdfCsgReal_t* const normZ, sdfCsgReal_t* const dist) { assert(that->type == sdfCsg_plane); *normX = that->params[sdfCsg_normX]; *normY = that->params[sdfCsg_normY]; *normZ = that->params[sdfCsg_normZ]; *dist = that->params[sdfCsg_dist]; } void SdfCsgSetSphereParams( SdfCsg* const that, sdfCsgReal_t const radius) { assert(that->type == sdfCsg_sphere); that->params[sdfCsg_radius] = radius; } void SdfCsgSetCapsuleParams( SdfCsg* const that, sdfCsgReal_t const radius, sdfCsgReal_t const length) { assert(that->type == sdfCsg_capsule); that->params[sdfCsg_radius] = radius; that->params[sdfCsg_halfLength] = length * 0.5; } void SdfCsgSetTorusParams( SdfCsg* const that, sdfCsgReal_t const majorRadius, sdfCsgReal_t const minorRadius) { assert(that->type == sdfCsg_torus); that->params[sdfCsg_majorRadius] = majorRadius; that->params[sdfCsg_minorRadius] = minorRadius; } void SdfCsgSetEllipsoidParams( SdfCsg* const that, sdfCsgReal_t const radiusX, sdfCsgReal_t const radiusY, sdfCsgReal_t const radiusZ) { assert(that->type == sdfCsg_ellipsoid); that->params[sdfCsg_radiusX] = radiusX; that->params[sdfCsg_radiusY] = radiusY; that->params[sdfCsg_radiusZ] = radiusZ; } void SdfCsgSetBoxParams( SdfCsg* const that, sdfCsgReal_t const width, sdfCsgReal_t const height, sdfCsgReal_t const depth) { assert(that->type == sdfCsg_box); that->params[sdfCsg_halfWidth] = width * 0.5; that->params[sdfCsg_halfHeight] = height * 0.5; that->params[sdfCsg_halfDepth] = depth * 0.5; } void SdfCsgSetConeParams( SdfCsg* const that, sdfCsgReal_t const baseRadius, sdfCsgReal_t const topRadius, sdfCsgReal_t const height) { assert(that->type == sdfCsg_cone); that->params[sdfCsg_baseRadius] = baseRadius; that->params[sdfCsg_topRadius] = topRadius; that->params[sdfCsg_height] = height; } void SdfCsgSetPlaneParams( SdfCsg* const that, sdfCsgReal_t const normX, sdfCsgReal_t const normY, sdfCsgReal_t const normZ, sdfCsgReal_t const dist) { assert(that->type == sdfCsg_plane); that->params[sdfCsg_normX] = normX; that->params[sdfCsg_normY] = normY; that->params[sdfCsg_normZ] = normZ; that->params[sdfCsg_dist] = dist; } void SdfCsgQuaternionCreate( sdfCsgReal_t const axis[static 3], sdfCsgReal_t const theta, sdfCsgReal_t quat[static 4]) { sdfCsgReal_t n = rhypot3(axis); sdfCsgReal_t s = sin(0.5 * theta); loop(i, 3) quat[i] = axis[i] / n * s; quat[3] = cos(0.5 * theta); } void SdfCsgTransformNone( SdfCsg* const that, unsigned int const rank) { assert(rank < SDFCSG_MAX_NB_TRANSFORM); that->transforms[rank].type = sdfCsg_none; } void SdfCsgTransformScale( SdfCsg* const that, unsigned int const rank, sdfCsgReal_t const scale) { assert(rank < SDFCSG_MAX_NB_TRANSFORM); that->transforms[rank].type = sdfCsg_scale; that->transforms[rank].vals[0] = scale; } void SdfCsgTransformTranslate( SdfCsg* const that, unsigned int const rank, sdfCsgReal_t const x, sdfCsgReal_t const y, sdfCsgReal_t const z) { assert(rank < SDFCSG_MAX_NB_TRANSFORM); that->transforms[rank].type = sdfCsg_translate; that->transforms[rank].vals[0] = x; that->transforms[rank].vals[1] = y; that->transforms[rank].vals[2] = z; } void SdfCsgTransformRotate( SdfCsg* const that, unsigned int const rank, sdfCsgReal_t const quat[static 4]) { assert(rank < SDFCSG_MAX_NB_TRANSFORM); that->transforms[rank].type = sdfCsg_rotate; loop(i, 4) that->transforms[rank].vals[i] = quat[i]; } void SdfCsgTransformExtrusion( SdfCsg* const that, unsigned int const rank, sdfCsgReal_t const length) { assert(rank < SDFCSG_MAX_NB_TRANSFORM); that->transforms[rank].type = sdfCsg_extrusion; that->transforms[rank].vals[0] = length * 0.5; } void SdfCsgTransformSymmetry( SdfCsg* const that, unsigned int const rank) { assert(rank < SDFCSG_MAX_NB_TRANSFORM); that->transforms[rank].type = sdfCsg_symmetry; } void SdfCsgTransformRevolution( SdfCsg* const that, unsigned int const rank) { assert(rank < SDFCSG_MAX_NB_TRANSFORM); that->transforms[rank].type = sdfCsg_revolution; } void SdfCsgTransformTwist( SdfCsg* const that, unsigned int const rank, sdfCsgReal_t const freq) { assert(rank < SDFCSG_MAX_NB_TRANSFORM); that->transforms[rank].type = sdfCsg_twist; that->transforms[rank].vals[0] = freq; } static sdfCsgReal_t SdfCsgDistBox( SdfCsg* const that, sdfCsgReal_t const pos[static 3]) { sdfCsgReal_t q[3]; rabs3(pos, q); loop(i, 3) q[i] -= that->params[sdfCsg_halfWidth + i]; sdfCsgReal_t maxQ[3]; sdfCsgReal_t const sdfOrig[3] = {0.0, 0.0, 0.0}; rmax3(q, sdfOrig, maxQ); sdfCsgReal_t dist = rhypot3(maxQ) + rmin(rmax(q[0], rmax(q[1], q[2])), 0.0); return dist; } static sdfCsgReal_t SdfCsgDistPlane( SdfCsg* const that, sdfCsgReal_t const pos[static 3]) { sdfCsgReal_t dist = pos[0] * that->params[sdfCsg_normX] + pos[1] * that->params[sdfCsg_normY] + pos[2] * that->params[sdfCsg_normZ] - that->params[sdfCsg_dist]; return dist; } static sdfCsgReal_t SdfCsgDistSphere( SdfCsg* const that, sdfCsgReal_t const pos[static 3]) { sdfCsgReal_t dist = rhypot3(pos) - that->params[sdfCsg_radius]; return dist; } static sdfCsgReal_t SdfCsgDistCone( SdfCsg* const that, sdfCsgReal_t const pos[static 3]) { sdfCsgReal_t q[2] = {rsqrt(pos[0] * pos[0] + pos[2] * pos[2]), pos[1]}; sdfCsgReal_t k1[2] = {that->params[sdfCsg_topRadius], that->params[sdfCsg_height]}; sdfCsgReal_t k2[2] = { that->params[sdfCsg_topRadius] - that->params[sdfCsg_baseRadius], 2.0 * that->params[sdfCsg_height] }; sdfCsgReal_t ca[2] = { q[0] - rmin( q[0], ( q[1] < 0.0 ? that->params[sdfCsg_baseRadius] : that->params[sdfCsg_topRadius] )), rabs(q[1]) - that->params[sdfCsg_height] }; sdfCsgReal_t cb[2] = {0}; loop(i, 2) { cb[i] = q[i] - k1[i]; sdfCsgReal_t k1q[2] = {k1[0] - q[0], k1[1] - q[1]}; sdfCsgReal_t c = (k1q[0] * k2[0] + k1q[1] * k2[1]) / (k2[0] * k2[0] + k2[1] * k2[1]); if(c < 0.0) c = 0.0; else if(c > 1.0) c = 1.0; cb[i] += k2[i] * c; } sdfCsgReal_t s = ((cb[0] < 0.0 && ca[1] < 0.0) ? -1.0 : 1.0); sdfCsgReal_t dist = s * rsqrt(rmin(ca[0] * ca[0] + ca[1] * ca[1], cb[0] * cb[0] + cb[1] * cb[1])); return dist; } static sdfCsgReal_t SdfCsgDistEllipsoid( SdfCsg* const that, sdfCsgReal_t const pos[static 3]) { sdfCsgReal_t q[3] = {0.0}; loop(i, 3) q[i] = pos[i] / that->params[sdfCsg_radiusX + i]; sdfCsgReal_t l = rhypot3(q); sdfCsgReal_t r[3] = {0.0}; if(l > 1e-9) { sdfCsgReal_t const a = 1.0 - 1.0 / l; loop(i, 3) r[i] = a * pos[i]; sdfCsgReal_t dist = (l > 1.0 ? rhypot3(r) : -rhypot3(r)); return dist; } else { return rmin( rmin(that->params[sdfCsg_radiusX], that->params[sdfCsg_radiusY]), that->params[sdfCsg_radiusZ]) - 1e-9; } } static sdfCsgReal_t SdfCsgDistTorus( SdfCsg* const that, sdfCsgReal_t const pos[static 3]) { sdfCsgReal_t q[3]; q[0] = rsqrt(pos[0] * pos[0] + pos[2] * pos[2]) - that->params[sdfCsg_majorRadius]; q[1] = pos[1]; q[2] = 0.0; sdfCsgReal_t dist = rhypot3(q) - that->params[sdfCsg_minorRadius]; return dist; } static sdfCsgReal_t SdfCsgDistCapsule( SdfCsg* const that, sdfCsgReal_t const pos[static 3]) { sdfCsgReal_t q[3]; q[0] = pos[0]; if(pos[1] < 0) q[1] = rmin(0, pos[1] + that->params[sdfCsg_halfLength]); else q[1] = rmax(0, pos[1] - that->params[sdfCsg_halfLength]); q[2] = pos[2]; sdfCsgReal_t dist = rhypot3(q) - that->params[sdfCsg_radius]; return dist; } static sdfCsgReal_t SdfCsgDistMerge( SdfCsg* const that, sdfCsgReal_t const pos[static 3]) { sdfCsgReal_t d0 = SdfCsgDist(that->comps[0], pos); sdfCsgReal_t d1 = SdfCsgDist(that->comps[1], pos); sdfCsgReal_t s = 0.0; if(that->smoothness > 1e-9) { sdfCsgReal_t h = rmax(that->smoothness - rabs(d0 - d1), 0.0) / that->smoothness; s = h * h * h * 0.166666666 * that->smoothness; } sdfCsgReal_t dist = rmin(d0, d1) - s; return dist; } static sdfCsgReal_t SdfCsgDistDifference( SdfCsg* const that, sdfCsgReal_t const pos[static 3]) { sdfCsgReal_t d0 = SdfCsgDist(that->comps[0], pos); sdfCsgReal_t d1 = SdfCsgDist(that->comps[1], pos); sdfCsgReal_t s = 0.0; if(that->smoothness > 1e-9) { sdfCsgReal_t h = rmax(that->smoothness - rabs(-(d0 + d1)), 0.0) / that->smoothness; s = h * h * h * 0.166666666 * that->smoothness; } sdfCsgReal_t dist = rmax(d0, -d1) + s; return dist; } static sdfCsgReal_t SdfCsgDistIntersection( SdfCsg* const that, sdfCsgReal_t const pos[static 3]) { sdfCsgReal_t d0 = SdfCsgDist(that->comps[0], pos); sdfCsgReal_t d1 = SdfCsgDist(that->comps[1], pos); sdfCsgReal_t s = 0.0; if(that->smoothness > 1e-9) { sdfCsgReal_t h = rmax(that->smoothness - rabs(d0 - d1), 0.0) / that->smoothness; s = h * h * h * 0.166666666 * that->smoothness; } sdfCsgReal_t dist = rmax(d0, d1) + s; return dist; } static SdfCsg* SdfCsgAlloc(void) { SdfCsg* that = malloc(sizeof(SdfCsg)); assert(that != NULL); loop(i, SDFCSG_MAX_NB_TRANSFORM) that->transforms[i].type = sdfCsg_none; that->thickness = 0.0; that->roundness = 0.0; loop(i, sdfCsg_nbParam) that->params[i] = 0.0; that->data = NULL; loop(i, 2) that->comps[i] = NULL; that->parent = NULL; that->smoothness = 0.0; return that; } SdfCsg* SdfCsgBox( sdfCsgReal_t const width, sdfCsgReal_t const height, sdfCsgReal_t const depth) { SdfCsg* that = SdfCsgAlloc(); that->type = sdfCsg_box; that->params[sdfCsg_halfWidth] = width * 0.5; that->params[sdfCsg_halfHeight] = height * 0.5; that->params[sdfCsg_halfDepth] = depth * 0.5; that->distCanonical = SdfCsgDistBox; return that; } SdfCsg* SdfCsgPlane( sdfCsgReal_t const normX, sdfCsgReal_t const normY, sdfCsgReal_t const normZ, sdfCsgReal_t const dist) { SdfCsg* that = SdfCsgAlloc(); that->type = sdfCsg_plane; sdfCsgReal_t l = 1.0 / rsqrt(normX * normX + normY * normY + normZ * normZ); that->params[sdfCsg_normX] = normX * l; that->params[sdfCsg_normY] = normY * l; that->params[sdfCsg_normZ] = normZ * l; that->params[sdfCsg_dist] = dist; that->distCanonical = SdfCsgDistPlane; return that; } SdfCsg* SdfCsgSphere(sdfCsgReal_t const radius) { SdfCsg* that = SdfCsgAlloc(); that->type = sdfCsg_sphere; that->params[sdfCsg_radius] = radius; that->distCanonical = SdfCsgDistSphere; return that; } SdfCsg* SdfCsgCone( sdfCsgReal_t const baseRadius, sdfCsgReal_t const topRadius, sdfCsgReal_t const height) { SdfCsg* that = SdfCsgAlloc(); that->type = sdfCsg_cone; that->params[sdfCsg_baseRadius] = baseRadius; that->params[sdfCsg_topRadius] = topRadius; that->params[sdfCsg_height] = height * 0.5; that->distCanonical = SdfCsgDistCone; return that; } SdfCsg* SdfCsgEllipsoid( sdfCsgReal_t const radiusX, sdfCsgReal_t const radiusY, sdfCsgReal_t const radiusZ) { SdfCsg* that = SdfCsgAlloc(); that->type = sdfCsg_ellipsoid; that->params[sdfCsg_radiusX] = radiusX; that->params[sdfCsg_radiusY] = radiusY; that->params[sdfCsg_radiusZ] = radiusZ; that->distCanonical = SdfCsgDistEllipsoid; return that; } SdfCsg* SdfCsgCapsule( sdfCsgReal_t const radius, sdfCsgReal_t const length) { SdfCsg* that = SdfCsgAlloc(); that->type = sdfCsg_capsule; that->params[sdfCsg_radius] = radius; that->params[sdfCsg_halfLength] = length * 0.5; that->distCanonical = SdfCsgDistCapsule; return that; } SdfCsg* SdfCsgTorus( sdfCsgReal_t const majorRadius, sdfCsgReal_t const minorRadius) { SdfCsg* that = SdfCsgAlloc(); that->type = sdfCsg_torus; that->params[sdfCsg_majorRadius] = majorRadius; that->params[sdfCsg_minorRadius] = minorRadius; that->distCanonical = SdfCsgDistTorus; return that; } SdfCsg* SdfCsgUnion( SdfCsg* const shapeA, SdfCsg* const shapeB) { SdfCsg* that = SdfCsgAlloc(); that->type = sdfCsg_union; that->comps[0] = shapeA; that->comps[1] = shapeB; shapeA->parent = that; shapeB->parent = that; that->distCanonical = SdfCsgDistMerge; return that; } SdfCsg* SdfCsgIntersection( SdfCsg* const shapeA, SdfCsg* const shapeB) { SdfCsg* that = SdfCsgAlloc(); that->type = sdfCsg_intersection; that->comps[0] = shapeA; that->comps[1] = shapeB; shapeA->parent = that; shapeB->parent = that; that->distCanonical = SdfCsgDistIntersection; return that; } SdfCsg* SdfCsgDifference( SdfCsg* const shapeA, SdfCsg* const shapeB) { SdfCsg* that = SdfCsgAlloc(); that->type = sdfCsg_difference; that->comps[0] = shapeA; that->comps[1] = shapeB; shapeA->parent = that; shapeB->parent = that; that->distCanonical = SdfCsgDistDifference; return that; } SdfCsg* SdfCsgUserDefined(SdfCsgFun const sdf) { SdfCsg* that = SdfCsgAlloc(); that->type = sdfCsg_user; that->distCanonical = sdf; return that; } static void SdfCsgApplyRotationInverse( sdfCsgReal_t const quat[4], sdfCsgReal_t const u[3], sdfCsgReal_t v[3]) { sdfCsgReal_t p[4]; loop(i, 3) p[i] = u[i]; p[3] = 0.0; sdfCsgReal_t w[4]; w[0] = p[3] * -quat[0] + p[0] * quat[3] - p[1] * -quat[2] + p[2] * -quat[1]; w[1] = p[3] * -quat[1] + p[0] * -quat[2] + p[1] * quat[3] - p[2] * -quat[0]; w[2] = p[3] * -quat[2] - p[0] * -quat[1] + p[1] * -quat[0] + p[2] * quat[3]; w[3] = p[3] * quat[3] - p[0] * -quat[0] - p[1] * -quat[1] - p[2] * -quat[2]; loop(i, 3) p[i] = quat[i]; p[3] = quat[3]; v[0] = p[3] * w[0] + p[0] * w[3] - p[1] * w[2] + p[2] * w[1]; v[1] = p[3] * w[1] + p[0] * w[2] + p[1] * w[3] - p[2] * w[0]; v[2] = p[3] * w[2] - p[0] * w[1] + p[1] * w[0] + p[2] * w[3]; } static void SdfCsgUntransform( SdfCsg const* const that, sdfCsgReal_t const pos[3], sdfCsgReal_t res[3]) { loop(i, 3) res[i] = pos[i]; for(int iTrans = SDFCSG_MAX_NB_TRANSFORM - 1; iTrans >= 0; --iTrans) { SdfCsgTransform const* trans = that->transforms + iTrans; switch(trans->type) { case sdfCsg_translate: loop(i, 3) res[i] -= trans->vals[i]; break; case sdfCsg_scale: loop(i, 3) res[i] /= trans->vals[0]; break; case sdfCsg_rotate: SdfCsgApplyRotationInverse(trans->vals, res, res); break; case sdfCsg_extrusion: if(res[2] < 0) res[2] = rmin(0, res[2] + trans->vals[0]); else res[2] = rmax(0, res[2] - trans->vals[0]); break; case sdfCsg_symmetry: res[0] = rabs(res[0]); break; case sdfCsg_revolution: res[0] = rsqrt(res[0] * res[0] + res[2] * res[2]); res[2] = 0.0; break; case sdfCsg_twist: sdfCsgReal_t theta = M_2PI * res[1] / trans->vals[0]; sdfCsgReal_t quat[4] = {0.0, sin(0.5 * theta), 0.0, cos(0.5 * theta)}; SdfCsgApplyRotationInverse(quat, res, res); break; default: break; } } } sdfCsgReal_t SdfCsgDist( SdfCsg* const that, sdfCsgReal_t const pos[static 3]) { sdfCsgReal_t posUntransformed[3]; SdfCsgUntransform(that, pos, posUntransformed); sdfCsgReal_t dist = that->distCanonical(that, posUntransformed); loop(i, SDFCSG_MAX_NB_TRANSFORM) { if(that->transforms[i].type == sdfCsg_scale) { dist *= that->transforms[i].vals[0]; } } dist -= that->roundness; if(that->thickness > 1e-9) dist = rabs(dist) - that->thickness; return dist; } void SdfCsgGradient( SdfCsg* const that, sdfCsgReal_t const pos[static 3], sdfCsgReal_t grad[static 3]) { loop(i, 3) { sdfCsgReal_t q[3]; loop(j, 3) q[j] = pos[j]; q[i] = pos[i] + sdfCsg_epsilon; grad[i] = SdfCsgDist(that, q); q[i] = pos[i] - sdfCsg_epsilon; grad[i] -= SdfCsgDist(that, q); grad[i] /= 2.0 * sdfCsg_epsilon; } } SdfCsgRayMarchRes SdfCsgRayMarchResCreate(sdfCsgReal_t const distMax) { SdfCsgRayMarchRes that = { .hasHit = false, .hasGaveUp = false, .shape = NULL, .rayDist = distMax, .pos = {0.0, 0.0, 0.0}, .posLocal = {0.0, 0.0, 0.0}, .posDist = 0.0, .minDist = distMax, }; return that; } void SdfCsgRayMarch( SdfCsg* const that, sdfCsgReal_t const from[static 3], sdfCsgReal_t const dir[static 3], SdfCsgRayMarchRes* const res) { sdfCsgReal_t pos[3]; loop(i, 3) pos[i] = from[i]; sdfCsgReal_t dist = SdfCsgDist(that, pos); sdfCsgReal_t sumDist = 0.0; res->hasGaveUp = false; unsigned int nbStep = 0; while( rabs(dist) > sdfCsg_epsilon && nbStep < sdfCsg_rayMarchNbMaxStep && rabs(sumDist) < rabs(res->rayDist) ) { ++nbStep; sumDist += dist; if(rabs(res->minDist) > rabs(dist)) res->minDist = dist; loop(i, 3) pos[i] += dist * dir[i]; dist = SdfCsgDist(that, pos); } if(rabs(dist) <= sdfCsg_epsilon) { res->hasHit = true; res->rayDist = sumDist; res->minDist = 0.0; res->shape = that; res->posDist = dist; loop(i, 3) res->pos[i] = pos[i]; SdfCsgUntransform(that, res->pos, res->posLocal); } res->hasGaveUp = (nbStep >= sdfCsg_rayMarchNbMaxStep); } static bool SdfCsgFindNearest( SdfCsg* const that, sdfCsgReal_t const from[static 3], sdfCsgReal_t res[static 3], double const dir) { unsigned int nbStep = 0; loop(i, 3) res[i] = from[i]; sdfCsgReal_t dist = SdfCsgDist(that, res); sdfCsgReal_t step = sdfCsg_epsilon; while(nbStep < sdfCsg_rayMarchNbMaxStep && rabs(dist) <= sdfCsg_epsilon) { sdfCsgReal_t gradient[3]; SdfCsgGradient(that, res, gradient); sdfCsgReal_t c = 1.0 / rsqrt(gradient[0] * gradient[0] + gradient[1] * gradient[1] + gradient[2] * gradient[2]); c *= step * dir; loop(i, 3) res[i] = from[i] + c * gradient[i]; dist = SdfCsgDist(that, res); ++nbStep; step += sdfCsg_epsilon; } return (nbStep < sdfCsg_rayMarchNbMaxStep); } bool SdfCsgFindNearestInside( SdfCsg* const that, sdfCsgReal_t const from[static 3], sdfCsgReal_t res[static 3]) { return SdfCsgFindNearest(that, from, res, -1.0); } bool SdfCsgFindNearestOutside( SdfCsg* const that, sdfCsgReal_t const from[static 3], sdfCsgReal_t res[static 3]) { return SdfCsgFindNearest(that, from, res, 1.0); }

Edit on 2023/08/27: A very interesting video about how the same idea as ray marching can be used to solve partial derivative equations.

Edit on 2024/06/11: Remove a useless calculation in SdfCsgDistIntersection; Refactor SdfCsgDistEllipsoid for faster calculation and avoid divide by zero; Update version number to 0.2.0.