Inigo Quilez   ::     ::  


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 useful 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 or in this video.

I have other article where I explain how to compute SDFs and raymarch more complex objects such as recursive primitives, fractals or acceleration struture based meshes. You can also often construct 3D SDFs by extruding or doing revolutions on 2D SDFs, so have a look to those too. And if you carea bout how I derived these formulas, you can explore some of the video tutorials I have on the topic.

In this index below, each primitives, modifier and operator function in this page you'll find a "exact" or "bound" note. This refers to the properties of the SDF that is generated or returned by the function. An "exact" SDF is one that retains all the qualities of a true SDF in Euclidean space - it really measures a distance exactly, meaning and its gradient always has length one. A "bound" SDF is no longer a true SDF (being pedantic) and only returns a lower bound to the real SDF, which can still be useful in certain scenarios. SDFs that are "exact" are generally desired over the "bound" ones because they work better with a bigger number of algorithms and techniques, and produce higher quality results. However some primitives (like the ellipsoid) or operators (like the smooth minimum here) cannot be "exact" because the very mathematics that describe them prevent it. In those cases, I propose a "bound" version and mark it so to prevent confusion. If a primitive that has an "exact" implementation can also be approximated by a "bound" version and sometimes I sometimes document that one too, but not too often because while locally faster than their "exact" counterpart they often backfire in the form of a global performance hit (ironically) due to their poorer quality in measuring distances.

All primitives here 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 primtives 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: and


Sphere - exact   (

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

Box - exact   (Youtube Tutorial with derivation:

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; }

Box Frame - exact   (

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

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; }

Link - exact   (

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, vec3 n, float h ) { // n must be normalized return dot(p,n) + h; }

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; }

Vertical Capped Cylinder - exact   (

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

Arbitrary 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)); }

Cut Sphere - exact   (

float sdCutSphere( vec3 p, float r, float h ) { // sampling independent computations (only depend on shape) float w = sqrt(r*r-h*h); // sampling dependant computations vec2 q = vec2( length(p.xz), p.y ); float s = max( (h-r)*q.x*q.x+w*w*(h+r-2.0*q.y), h*q.x-w*q.y ); return (s<0.0) ? length(q)-r : (q.x<w) ? h - q.y : length(q-vec2(w,h)); }

Cut Hollow Sphere - exact   (

float sdCutHollowSphere( vec3 p, float r, float h, float t ) { // sampling independent computations (only depend on shape) float w = sqrt(r*r-h*h); // sampling dependant computations vec2 q = vec2( length(p.xz), p.y ); return ((h*q.x<w*q.y) ? length(q-vec2(w,h)) : abs(length(q)-r) ) - t; }

Death Star - exact   (

float sdDeathStar( in vec3 p2, in float ra, float rb, in float d ) { // sampling independent computations (only depend on shape) float a = (ra*ra - rb*rb + d*d)/(2.0*d); float b = sqrt(max(ra*ra-a*a,0.0)); // sampling dependant computations vec2 p = vec2( p2.x, length(p2.yz) ); if( p.x*b-p.y*a > d*max(b-p.y,0.0) ) return length(p-vec2(a,b)); else return max( (length(p )-ra), -(length(p-vec2(d,0))-rb)); }

Round cone - exact

float sdRoundCone( vec3 p, float r1, float r2, float h ) { // sampling independent computations (only depend on shape) float b = (r1-r2)/h; float a = sqrt(1.0-b*b); // sampling dependant computations vec2 q = vec2( length(p.xz), p.y ); 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 =; 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) ); }

Quad - exact   (

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 operations 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 apart 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 example, the Capsule is an elongated Sphere along an axis really. You can find code here:

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 perfectly 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 differently (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, hexagons or any other shape like the cone in the image below. If you happen to be interested in preserving the overall 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:

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:

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 intensive to evaluate than when produced from boolean operations on other volumes. Two of the most simplest way to make volumes out 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:

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 they only give a bound to the real SDF, this 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:

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, bound

These are the most basic combinations of pairs of primitives you can do. They correspond to the basic boolean operations. Please note that only the Union of two SDFs returns a true SDF, not the Subtraction or Intersection. To make it more subtle, this is only true in the exterior of the SDF (where distances are positive) and not in the interior. You can learn more about this and how to work around it in the article "Interior Distances". Also 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 ) { return 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 - bound, bound, bound

Blending primitives is a really powerful 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 replace 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 between 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:

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); }


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 scaling 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 below 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

Symmetry is useful, since many things around us are symmetric, from humans, animals, vehicles, instruments, furniture, ... Oftentimes, one can take shortcuts 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 evaluator and without increasing the memory footprint of your application. The code below 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 primitive( 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 repeates 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 below 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 infinite 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 produces 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. You can see this in action here as a 2D shader, although it works in any number of dimensions: Code:

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 produces a rectangle of l.x times l.y instances. Naturally, you can create any rectangular 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 allow 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 careful 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 reduced (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 functions (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.


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; }


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); }


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):