website articles
distance functions

Intro

After having posted about the basics of distance functions in several places (pouet, my blog, shadertoy, private emails, etc), I thought it might make sense to put these together in centralized place. Here you will find the distance functions for basic primitives, plus the formulas for combining them together for building more complex shapes, as well as some distortion functions that you can use to shape your objects. Hopefully this will be usefull for those rendering scenes with raymarching. You can see some of the results you can get by using these techniques in the raymarching distance fields article. Lastly, this article doesn't include lighting tricks, nor marching acceleartion tricks or more advanced techniques as recursive primitives or fractals. In case you are looking for 2D SDF functions, you'll find them in the 2D Distance Functions page.

Primitives

All primitives are centered at the origin. You will have to transform the point to get arbitrarily rotated, translated and scaled objects (see below).

Many of these primtiives below use dot2() or ndot(), which I list here quickly before the primitives:

float dot2( in vec2 v ) { return dot(v,v); } float dot2( in vec3 v ) { return dot(v,v); } float ndot( in vec2 a, in vec2 b ) { return a.x*b.x - a.y*b.y; }

Lastly, you have working sample code of all of these primitives here: https://www.shadertoy.com/view/Xds3zN

Sphere - exact
float sdSphere( vec3 p, float s ) { return length(p)-s; }

Box - exact

float sdBox( vec3 p, vec3 b ) { vec3 q = abs(p) - b; return length(max(q,0.0)) + min(max(q.x,max(q.y,q.z)),0.0); }

Round Box - exact
float sdRoundBox( vec3 p, vec3 b, float r ) { vec3 q = abs(p) - b; return length(max(q,0.0)) + min(max(q.x,max(q.y,q.z)),0.0) - r; }

Torus - exact
float sdTorus( vec3 p, vec2 t ) { vec2 q = vec2(length(p.xz)-t.x,p.y); return length(q)-t.y; }

Capped Torus - exact
float sdCappedTorus(in vec3 p, in vec2 sc, in float ra, in float rb) { p.x = abs(p.x); float k = (sc.y*p.x>sc.x*p.y) ? dot(p.xy,sc) : length(p.xy); return sqrt( dot(p,p) + ra*ra - 2.0*ra*k ) - rb; }

float sdLink( vec3 p, float le, float r1, float r2 ) { vec3 q = vec3( p.x, max(abs(p.y)-le,0.0), p.z ); return length(vec2(length(q.xy)-r1,q.z)) - r2; }

Infinite Cylinder - exact
float sdCylinder( vec3 p, vec3 c ) { return length(p.xz-c.xy)-c.z; }

Cone - exact
float sdCone( in vec3 p, in vec2 c, float h ) { // c is the sin/cos of the angle, h is height // Alternatively pass q instead of (c,h), // which is the point at the base in 2D vec2 q = h*vec2(c.x/c.y,-1.0); vec2 w = vec2( length(p.xz), p.y ); vec2 a = w - q*clamp( dot(w,q)/dot(q,q), 0.0, 1.0 ); vec2 b = w - q*vec2( clamp( w.x/q.x, 0.0, 1.0 ), 1.0 ); float k = sign( q.y ); float d = min(dot( a, a ),dot(b, b)); float s = max( k*(w.x*q.y-w.y*q.x),k*(w.y-q.y) ); return sqrt(d)*sign(s); }

Cone - bound (not exact!)
float sdCone( vec3 p, vec2 c, float h ) { float q = length(p.xz); return max(dot(c.xy,vec2(q,p.y)),-h-p.y); }

Infinite Cone - exact
float sdCone( vec3 p, vec2 c ) { // c is the sin/cos of the angle vec2 q = vec2( length(p.xz), -p.y ); float d = length(q-c*max(dot(q,c), 0.0)); return d * ((q.x*c.y-q.y*c.x<0.0)?-1.0:1.0); }

Plane - exact
float sdPlane( vec3 p, vec4 n ) { // n must be normalized return dot(p,n.xyz) + n.w; }

Hexagonal Prism - exact
float sdHexPrism( vec3 p, vec2 h ) { const vec3 k = vec3(-0.8660254, 0.5, 0.57735); p = abs(p); p.xy -= 2.0*min(dot(k.xy, p.xy), 0.0)*k.xy; vec2 d = vec2( length(p.xy-vec2(clamp(p.x,-k.z*h.x,k.z*h.x), h.x))*sign(p.y-h.x), p.z-h.y ); return min(max(d.x,d.y),0.0) + length(max(d,0.0)); }

Triangular Prism - bound
float sdTriPrism( vec3 p, vec2 h ) { vec3 q = abs(p); return max(q.z-h.y,max(q.x*0.866025+p.y*0.5,-p.y)-h.x*0.5); }

Capsule / Line - exact
float sdCapsule( vec3 p, vec3 a, vec3 b, float r ) { vec3 pa = p - a, ba = b - a; float h = clamp( dot(pa,ba)/dot(ba,ba), 0.0, 1.0 ); return length( pa - ba*h ) - r; }

Capsule / Line - exact
float sdVerticalCapsule( vec3 p, float h, float r ) { p.y -= clamp( p.y, 0.0, h ); return length( p ) - r; }

Capped Cylinder - exact
float sdCappedCylinder( vec3 p, float h, float r ) { vec2 d = abs(vec2(length(p.xz),p.y)) - vec2(h,r); return min(max(d.x,d.y),0.0) + length(max(d,0.0)); }

Capped Cylinder - exact
float sdCappedCylinder(vec3 p, vec3 a, vec3 b, float r) { vec3 ba = b - a; vec3 pa = p - a; float baba = dot(ba,ba); float paba = dot(pa,ba); float x = length(pa*baba-ba*paba) - r*baba; float y = abs(paba-baba*0.5)-baba*0.5; float x2 = x*x; float y2 = y*y*baba; float d = (max(x,y)<0.0)?-min(x2,y2):(((x>0.0)?x2:0.0)+((y>0.0)?y2:0.0)); return sign(d)*sqrt(abs(d))/baba; }

Rounded Cylinder - exact
float sdRoundedCylinder( vec3 p, float ra, float rb, float h ) { vec2 d = vec2( length(p.xz)-2.0*ra+rb, abs(p.y) - h ); return min(max(d.x,d.y),0.0) + length(max(d,0.0)) - rb; }

Capped Cone - exact
float sdCappedCone( vec3 p, float h, float r1, float r2 ) { vec2 q = vec2( length(p.xz), p.y ); vec2 k1 = vec2(r2,h); vec2 k2 = vec2(r2-r1,2.0*h); vec2 ca = vec2(q.x-min(q.x,(q.y<0.0)?r1:r2), abs(q.y)-h); vec2 cb = q - k1 + k2*clamp( dot(k1-q,k2)/dot2(k2), 0.0, 1.0 ); float s = (cb.x<0.0 && ca.y<0.0) ? -1.0 : 1.0; return s*sqrt( min(dot2(ca),dot2(cb)) ); }

Capped Cone - exact
float sdCappedCone(vec3 p, vec3 a, vec3 b, float ra, float rb) { float rba = rb-ra; float baba = dot(b-a,b-a); float papa = dot(p-a,p-a); float paba = dot(p-a,b-a)/baba; float x = sqrt( papa - paba*paba*baba ); float cax = max(0.0,x-((paba<0.5)?ra:rb)); float cay = abs(paba-0.5)-0.5; float k = rba*rba + baba; float f = clamp( (rba*(x-ra)+paba*baba)/k, 0.0, 1.0 ); float cbx = x-ra - f*rba; float cby = paba - f; float s = (cbx < 0.0 && cay < 0.0) ? -1.0 : 1.0; return s*sqrt( min(cax*cax + cay*cay*baba, cbx*cbx + cby*cby*baba) ); }

Solid Angle - exact
float sdSolidAngle(vec3 p, vec2 c, float ra) { // c is the sin/cos of the angle vec2 q = vec2( length(p.xz), p.y ); float l = length(q) - ra; float m = length(q - c*clamp(dot(q,c),0.0,ra) ); return max(l,m*sign(c.y*q.x-c.x*q.y)); }

Round cone - exact
float sdRoundCone( vec3 p, float r1, float r2, float h ) { vec2 q = vec2( length(p.xz), p.y ); float b = (r1-r2)/h; float a = sqrt(1.0-b*b); float k = dot(q,vec2(-b,a)); if( k < 0.0 ) return length(q) - r1; if( k > a*h ) return length(q-vec2(0.0,h)) - r2; return dot(q, vec2(a,b) ) - r1; }

Round Cone - exact
float sdRoundCone(vec3 p, vec3 a, vec3 b, float r1, float r2) { // sampling independent computations (only depend on shape) vec3 ba = b - a; float l2 = dot(ba,ba); float rr = r1 - r2; float a2 = l2 - rr*rr; float il2 = 1.0/l2; // sampling dependant computations vec3 pa = p - a; float y = dot(pa,ba); float z = y - l2; float x2 = dot2( pa*l2 - ba*y ); float y2 = y*y*l2; float z2 = z*z*l2; // single square root! float k = sign(rr)*rr*rr*x2; if( sign(z)*a2*z2 > k ) return sqrt(x2 + z2) *il2 - r2; if( sign(y)*a2*y2 < k ) return sqrt(x2 + y2) *il2 - r1; return (sqrt(x2*a2*il2)+y*rr)*il2 - r1; }

Ellipsoid - bound (not exact!)
float sdEllipsoid( vec3 p, vec3 r ) { float k0 = length(p/r); float k1 = length(p/(r*r)); return k0*(k0-1.0)/k1; }

Rhombus - exact
float sdRhombus(vec3 p, float la, float lb, float h, float ra) { p = abs(p); vec2 b = vec2(la,lb); float f = clamp( (ndot(b,b-2.0*p.xz))/dot(b,b), -1.0, 1.0 ); vec2 q = vec2(length(p.xz-0.5*b*vec2(1.0-f,1.0+f))*sign(p.x*b.y+p.z*b.x-b.x*b.y)-ra, p.y-h); return min(max(q.x,q.y),0.0) + length(max(q,0.0)); }

Octahedron - exact
float sdOctahedron( vec3 p, float s) { p = abs(p); float m = p.x+p.y+p.z-s; vec3 q; if( 3.0*p.x < m ) q = p.xyz; else if( 3.0*p.y < m ) q = p.yzx; else if( 3.0*p.z < m ) q = p.zxy; else return m*0.57735027; float k = clamp(0.5*(q.z-q.y+s),0.0,s); return length(vec3(q.x,q.y-s+k,q.z-k)); }

Octahedron - bound (not exact)
float sdOctahedron( vec3 p, float s) { p = abs(p); return (p.x+p.y+p.z-s)*0.57735027; }

Pyramid - exact
float sdPyramid( vec3 p, float h) { float m2 = h*h + 0.25; p.xz = abs(p.xz); p.xz = (p.z>p.x) ? p.zx : p.xz; p.xz -= 0.5; vec3 q = vec3( p.z, h*p.y - 0.5*p.x, h*p.x + 0.5*p.y); float s = max(-q.x,0.0); float t = clamp( (q.y-0.5*p.z)/(m2+0.25), 0.0, 1.0 ); float a = m2*(q.x+s)*(q.x+s) + q.y*q.y; float b = m2*(q.x+0.5*t)*(q.x+0.5*t) + (q.y-m2*t)*(q.y-m2*t); float d2 = min(q.y,-q.x*m2-q.y*0.5) > 0.0 ? 0.0 : min(a,b); return sqrt( (d2+q.z*q.z)/m2 ) * sign(max(q.z,-p.y)); }

Triangle - exact
float udTriangle( vec3 p, vec3 a, vec3 b, vec3 c ) { vec3 ba = b - a; vec3 pa = p - a; vec3 cb = c - b; vec3 pb = p - b; vec3 ac = a - c; vec3 pc = p - c; vec3 nor = cross( ba, ac ); return sqrt( (sign(dot(cross(ba,nor),pa)) + sign(dot(cross(cb,nor),pb)) + sign(dot(cross(ac,nor),pc))<2.0) ? min( min( dot2(ba*clamp(dot(ba,pa)/dot2(ba),0.0,1.0)-pa), dot2(cb*clamp(dot(cb,pb)/dot2(cb),0.0,1.0)-pb) ), dot2(ac*clamp(dot(ac,pc)/dot2(ac),0.0,1.0)-pc) ) : dot(nor,pa)*dot(nor,pa)/dot2(nor) ); }

float udQuad( vec3 p, vec3 a, vec3 b, vec3 c, vec3 d ) { vec3 ba = b - a; vec3 pa = p - a; vec3 cb = c - b; vec3 pb = p - b; vec3 dc = d - c; vec3 pc = p - c; vec3 ad = a - d; vec3 pd = p - d; vec3 nor = cross( ba, ad ); return sqrt( (sign(dot(cross(ba,nor),pa)) + sign(dot(cross(cb,nor),pb)) + sign(dot(cross(dc,nor),pc)) + sign(dot(cross(ad,nor),pd))<3.0) ? min( min( min( dot2(ba*clamp(dot(ba,pa)/dot2(ba),0.0,1.0)-pa), dot2(cb*clamp(dot(cb,pb)/dot2(cb),0.0,1.0)-pb) ), dot2(dc*clamp(dot(dc,pc)/dot2(dc),0.0,1.0)-pc) ), dot2(ad*clamp(dot(ad,pd)/dot2(ad),0.0,1.0)-pd) ) : dot(nor,pa)*dot(nor,pa)/dot2(nor) ); }

Primitive alterations

Once we have the basic primitives, it's possible to apply some simple opearations that change their shape while still retaining exact an euclidean metric to them, which is an important property since SDFs with undistorted euclidean metric allow for faster ray marchine.

Elongation - exact

Elongating is a useful way to construct new shapes. It basically splits a primitive in two (four or eight), moves the pieces appart and and connects them. It is a perfect distance preserving operation, it does not introduce any artifacts in the SDF. Some of the basic primitives above use this technique. For exampple, the Capsule is an elongater Sphere along an axis really. You can find code here: https://www.shadertoy.com/view/Ml3fWj

float opElongate( in sdf3d primitive, in vec3 p, in vec3 h ) { vec3 q = p - clamp( p, -h, h ); return primitive( q ); } float opElongate( in sdf3d primitive, in vec3 p, in vec3 h ) { vec3 q = abs(p)-h; return primitive( max(q,0.0) ) + min(max(q.x,max(q.y,q.z)),0.0); }

The reason I provide to implementations is the following. For 1D elongations, the first function works perfecgtly and gives exact exterior and interior distances. However, the first implementation produces a small core of zero distances inside the volume for 2D and 3D elongations. Depending on your application that might be a problem. One way to create exact interior distances all the way to the very elongated core of the volume, is the following, which is in languages like GLSL that don't have function pointers or lambdas need to be implemented a bit differertly (check the code linked about in Shadertoy to see one example).

Rounding - exact

Rounding a shape is as simple as subtracting some distance (jumping to a different isosurface). The rounded box above is an example, but you can apply it to cones, heagons or any other shape like the cone in the image bellow. If you happen to be interested in preserving the overal volume of the shape, most of the times it's pretty easy to shrink the source primitive by the same amount we are rounding it by. You can find code here: https://www.shadertoy.com/view/Mt3BDj

float opRound( in sdf3d primitive, float rad ) { return primitive(p) - rad }

Onion - exact

For carving interiors or giving thickness to primitives, without performing expensive boolean operations (see below) and without distorting the distance field into a bound, one can use "onioning". You can use it multiple times to create concentric layers in your SDF. You can find code here: https://www.shadertoy.com/view/MlcBDj

float opOnion( in float sdf, in float thickness ) { return abs(sdf)-thickness; }

Revolution and extrusion from 2D - exact

Generating 3D volumes from 2D shapes has many advantages. Assuming the 2D shape defines exact distances, the resulting 3D volume is exact and way often less intensinve to evaluate than when produced from boolean operations on other volumes. Two of the most simplest way to make volumes our of flat shapes is to use extrusion and revolution (generalizations of these are easy to build, but we we'll keep simple here): You can find code here: https://www.shadertoy.com/view/4lyfzw

float opExtrusion( in vec3 p, in sdf2d primitive, in float h ) { float d = primitive(p.xy) vec2 w = vec2( d, abs(p.z) - h ); return min(max(w.x,w.y),0.0) + length(max(w,0.0)); } float opRevolution( in vec3 p, in sdf2d primitive, float o ) { vec2 q = vec2( length(p.xz) - o, p.y ); return primitive(q) }

Change of Metric - bound

Most of these functions can be modified to use other norms than the euclidean. By replacing length(p), which computes (x2+y2+z2)1/2 by (xn+yn+zn)1/n one can get variations of the basic primitives that have rounded edges rather than sharp ones. I do not recommend this technique though, since these primitives require more raymarching steps until an intersection is found than euclidean primitives. Since the they only give a bound to the real SDF, these kind of primitive alteration also doesn't play well with shadows and occlusion algorithms that rely on true SDFs for measuring distance to occluders. You can find the code here: https://www.shadertoy.com/view/ltcfDj

float length2( vec3 p ) { p=p*p; return sqrt( p.x+p.y+p.z); } float length6( vec3 p ) { p=p*p*p; p=p*p; return pow(p.x+p.y+p.z,1.0/6.0); } float length8( vec3 p ) { p=p*p; p=p*p; p=p*p; return pow(p.x+p.y+p.z,1.0/8.0); }

Primitive combinations

Sometimes you cannot simply elongate, round or onion a primitive, and you need to combine, carve or intersect basic primitives. Given the SDFs d1 and d2 of two primitives, you can use the following operators to combine together.

Union, Subtraction, Intersection - exact, bound, bound

These are the most basic combinations of pairs of primitives you can do. They correspond to the basic boolean operations. Note that opSubtraction() is not commutative and depending on the order of the operand it will produce different results.

float opUnion( float d1, float d2 ) { min(d1,d2); } float opSubtraction( float d1, float d2 ) { return max(-d1,d2); } float opIntersection( float d1, float d2 ) { return max(d1,d2); }

Smooth Union, Subtraction and Intersection - exact, bound, bound

Blendng primitives is a really powerfull tool - it allows to construct complex and organic shapes without the geometrical semas that normal boolean operations produce. There are many flavors of such operations, but the basic ones try to repalce the min() and max() functions used in the opUnion, opSubstraction and opIntersection above with smooth versions. They all accept an extra parameter called k that defines the size of the smooth transition beteewn the two primitives. It is given in actual distance units. You can find more details in the smooth minimum article article in this same site. You can code here: https://www.shadertoy.com/view/lt3BW2

float opSmoothUnion( float d1, float d2, float k ) { float h = clamp( 0.5 + 0.5*(d2-d1)/k, 0.0, 1.0 ); return mix( d2, d1, h ) - k*h*(1.0-h); } float opSmoothSubtraction( float d1, float d2, float k ) { float h = clamp( 0.5 - 0.5*(d2+d1)/k, 0.0, 1.0 ); return mix( d2, -d1, h ) + k*h*(1.0-h); } float opSmoothIntersection( float d1, float d2, float k ) { float h = clamp( 0.5 - 0.5*(d2-d1)/k, 0.0, 1.0 ); return mix( d2, d1, h ) + k*h*(1.0-h); }

Positioning

Placing primitives in different locations and orientations in space is a fundamental operation in designing SDFs. While rotations, uniform scaling and translations are exact operations, non-uniform scalling distorts the euclidean spaces and can only be bound. Therefore I do not include it here.

Rotation/Translation - exact

Since rotations and translation don't compress nor dilate space, all we need to do is simply to transform the point being sampled with the inverse of the transformation used to place an object in the scene. This code bellow assumes that transform encodes only a rotation and a translation (as a 3x4 matrix for example, or as a quaternion and a vector), and that it does not contain any scaling factors in it.

vec3 opTx( in vec3 p, in transform t, in sdf3d primitive ) { return primitive( invert(t)*p ); }

Scale - exact

Scaling an obect is slightly more tricky since that compresses/dilates spaces, so we have to take that into account on the resulting distance estimation. Still, it's not difficult to perform:

float opScale( in vec3 p, in float s, in sdf3d primitive ) { return primitive(p/s)*s; }

Symmetry - bound and exact

Symmtery is useful, since many things around us are symmetric, from humans, animals, vehicles, insturments, furniture, ... Often times, one can take shorcuts and only model half or a quarter of the desired shape, and get it duplicated automatically by using the absolute value of the domain coordinates before evaluation. For example, in the image below, there's a single object evaluation instead of two. This is a great savings in performance. You have to be aware however that the resuluting SDF might not be an exact SDF but a bound, if the object you are mirroring crosses the mirroring plane.

float opSymX( in vec3 p, in sdf3d primitive ) { p.x = abs(p.x); return primitive(p); } float opSymXZ( in vec3 p, in sdf3d primitive ) { p.xz = abs(p.xz); return primitive(p); }

Infinite Repetition

Domain repetition is a very useful operator, since it allows you to create infinitely many primitives with a single object eveluator and without increasing the memory footprint of your application. The code bellow shows how to perform the operation in the simplest way:

float opRep( in vec3 p, in vec3 c, in sdf3d primitive ) { vec3 q = mod(p+0.5*c,c)-0.5*c; return primitve( q ); }

In this code c is the repetition period (which can be different in each coordinate direction). This will work great for primitives that have a bounding box smaller than half the repetition period. If the object is big, you will need to check the 7 neighboring repetitions (in 3D, 3 in 2D) to check for closest neighbors, just as you usually do in a voronoi/Worley/cellular construction. You have an example of this in action in the following image where all of the grass field is made from a single blade which repeated infinitely in the X and Z directions with the code above. (A link to the real time animation and code for the image is right bellow the image).

Finite Repetition

Infinite domain repetition is great, but sometimes you only need a few copies or instances of a given SDF, not infinite. A frequently seen but suboptimal solution is to generate infintie copies and then clip the unwanted areas away with a box SDF. This is not ideal because the resulting SDF is not a real SDF but just a bound, since clipping through max() only prouces a bound. A much better approach is to clamp the indices of the instances instead of the SDF, and let a correct SDF emerge from the truncated/clamped indices:

vec3 opRepLim( in vec3 p, in float c, in vec3 l, in sdf3d primitive ) { vec3 q = p-c*clamp(round(p/c),-l,l); return primitive( q ); }

This procudes a rectangle of ls.x x l.y instances. Naturaly, you can create any rectanguar shape you want (or other shape) by changing the limits of the clamp or the clamp function itself. Also note that the function can be specialized to 2 or 1 dimensions easily.

Deformations and distortions

Deformations and distortions allows now to enhance the shape of primitives or even fuse different primitives together. The operations usually distort the distance field and make it non euclidean anymore, so one must be carefull when raymarching them, you will probably need to decrease your step size, if you are using a raymarcher to sample this. In principle one can compute the factor by which the step size needs to be rediced (inversely proportional to the compression of the space, which is given by the Jacobian of the deformation function). But even with dual numbers or automatic differentiation, it's usually just easier to find the constant by hand for a given primitive.

I'd say that while it is tempting to use a distortion or displacement to achieve a given shape, and I often use them myself of course, it is sometimes better to get as close to the desired shape with actual exact euclidean primitive operations (elongation, rounding, onioning, union) or tight bounded funcitons (intersection, subtraction) and then only apply as small of a distortion or displacement as possible. That way the field stays as close as possible to an actual distance field, and the raymarcher will be faster.

Displacement

The displacement example below is using sin(20*p.x)*sin(20*p.y)*sin(20*p.z) as displacement pattern, but you can of course use anything you might imagine.

float opDisplace( in sdf3d primitive, in vec3 p ) { float d1 = primitive(p); float d2 = displacement(p); return d1+d2; }

Twist

float opTwist( in sdf3d primitive, in vec3 p ) { const float k = 10.0; // or some other amount float c = cos(k*p.y); float s = sin(k*p.y); mat2 m = mat2(c,-s,s,c); vec3 q = vec3(m*p.xz,p.y); return primitive(q); }

Bend

float opCheapBend( in sdf3d primitive, in vec3 p ) { const float k = 10.0; // or some other amount float c = cos(k*p.x); float s = sin(k*p.x); mat2 m = mat2(c,-s,s,c); vec3 q = vec3(m*p.xy,p.z); return primitive(q); }

A reference implementation of most of these primitives and operators can be found here (click in the image to rotate the camera, or in the title to jump to the source code):