Simple SDF Collisions
31/10/2026
Here is a super simple algorithm for computing when two SDFs are colliding/intersecting.
While there is a variety of research online on computing collisions between SDFs (and some hint at the method I describe here), most papers contain pretty complex algorithms with the goal of maximising performance using different heuristics. But one of the beautiful things about SDFs are their simplicity - so I think there is value in presenting an SDF collision algorithm with the same complexity as something like SDF ray tracing, which is what I've tried to provide here.
We start with two shapes defined by SDFs. And we assume we already have bounding boxes for these shapes:
The first step is to compute the intersection of these two bounding boxes. This is going to be our region of the space where we are going to search for a collision.
Once we have an initial bounding box we are going to sample each of our shapes' SDFs at the center of this bounding box. We then compute the maximum of these two distances. If this maximum is less than zero then we've found a collision - we've found a point inside both shapes:
bool sdIntersects(Bound b)
{
// Compute the bound center
Vector2 p = BoundCenter(b);
// Sample each of the SDFs at the bound center
float d0 = sdShape0(p);
float d1 = sdShape1(p);
// If the sampled point contains an intersection return true
if (maxf(d0, d1) <= 0.0f)
{
return true;
}
...
If the value is greater than zero then we need to keep looking. But where should we sample next?
Well, the trick is to compute the smallest possible combined SDF value contained within this bounding box. This smallest possible combined SDF value will be the combined value we sampled at the center, minus the radius of the bounding box. This is because the greatest amount by which our combined SDF value can change over a given area is exactly the euclidean distance of the furthest point in that area from the point we sampled. See this fantastic blog post or this paper for a little more detail on why that is the case.
So if the smallest possible value is still greater than zero, then no intersection is possible anywhere within this region of the space and we can return that there is no collision.
...
// Compute the bound radius
float r = BoundRadius(b);
// If there is no possible intersection return false
if (maxf(d0, d1) - r > 0.0f)
{
return false;
}
...
Otherwise we simply split the bounding box into two and perform exactly the same process on each of the two smaller sub-bounding boxes:
...
// Otherwise split the Bounding Box in two and recurse
Bound s0, s1;
BoundSplit(&s0, &s1, b);
return sdIntersects(s0) || sdIntersects(s1);
}
As this is a recursive process we need to ensure it doesn't go on forever. A simple way to do this is to stop when the bounding box radius is below some threshold.
At this point we need to decide what to do. If we want to be conservative we can return true if the bounding box may contain a collision using our previous method. Otherwise we can return true only if we find a collision right in the center of the bound itself. For something like frustum culling the conservative answer is probably the better option. For something like physics simulation, some minor under-estimation of the collision might be better.
Here is everything all in one function:
bool sdIntersects(Bound b)
{
// Compute the bound center and radius
Vector2 p = BoundCenter(b);
float r = BoundRadius(b);
// Sample each of the SDFs at the bound center
float d0 = sdShape0(p);
float d1 = sdShape1(p);
// If the sampled point contains an intersection return true
if (maxf(d0, d1) <= 0.0f)
{
return true;
}
// If there is no possible intersection or the radius is too small return false
float minrad = 0.02f;
if (maxf(d0, d1) - r > margin || r < minrad)
{
return false;
}
// Otherwise split the Bounding Box in two and recurse
Bound s0, s1;
BoundSplit(&s0, &s1, b);
return sdIntersects(s0) || sdIntersects(s1);
}
And here is a visualization of it in action:
This algorithm is basically a depth-first binary search, so the complexity is \( \log(n) \) with respect to the size of the initial bounding box.
A breadth-first search tends to perform better since it searches the space in a coarse-to-fine way which means it can early-out more quickly for simple shapes. We can implement this with a queue instead of recursion:
bool sdIntersects2(Bound b)
{
// Insert the initial bound into the queue
Bound queue[INTERSECTS_QUEUE_MAX];
int queue_num = 1;
int queue_start = 0;
queue[queue_start] = b;
// Iterate over the queue
while (queue_num)
{
// Process the next bound in the queue
Bound c = queue[queue_start];
queue_start = (queue_start + 1) % INTERSECTS_QUEUE_MAX;
queue_num--;
// Compute the bound center and radius
Vector2 p = BoundCenter(c);
float r = BoundRadius(c);
// Sample the SDFs at the bound center
float d0 = sdShape0(p);
float d1 = sdShape1(p);
// If the sampled point contains an intersection return true
if (maxf(d0, d1) <= 0.0f)
{
return true;
}
// If the bound contains a possible intersection,
// the radius is large enough, and the queue is not full
float minrad = 0.02f;
if (maxf(d0, d1) - r > 0.0f
|| r < minrad
|| queue_num + 2 > INTERSECTS_QUEUE_MAX)
{
continue;
}
// Split the bound into two add to the queue
Bound s0, s1;
BoundSplit(&s0, &s1, c);
queue[(queue_start + queue_num + 0) % INTERSECTS_QUEUE_MAX] = s0;
queue[(queue_start + queue_num + 1) % INTERSECTS_QUEUE_MAX] = s1;
queue_num += 2;
}
// No bounds left to search, return false
return false;
}
If we want to use this method for collision response in something like a physics engine then we also need to have the contact points on the surface of each SDF shape, as well as the collision normal(s).
One idea for this is to try to find not just if there is a collision, but also the point within the intersection region which has the greatest combined depth to both surfaces. This can be done with a slight adaption to our algorithm. Rather than early-outing when we find a collision we're instead going to update our best candidate point (and it's combined depth), and only ignore regions of the space which are either entirely outside the intersection region, or cannot contain a better candidate.
void sdPenetration(
Vector2 * point,
float * depth,
Bound b)
{
// Compute the bound center and radius
Vector2 p = BoundCenter(b);
float r = BoundRadius(b);
// Sample each of the SDFs at the bound center
float d0 = sdShape0(p);
float d1 = sdShape1(p);
// Update best depth found
if (maxf(d0, d1) <= 0.0f && d0 + d1 <= *depth)
{
*depth = d0 + d1;
*point = p;
}
// There is no possible intersection or improvement
float minrad = 0.02f;
if (maxf(d0, d1) - r > 0.0f || d0 + d1 - 2 * r > *depth || r < minrad)
{
return;
}
Bound s0, s1;
BoundSplit(&s0, &s1, b);
sdPenetration(point, depth, s0);
sdPenetration(point, depth, s1);
}
Like before, we can also do a breadth-first non-recursive version of this:
void sdPenetration2(
Vector2 * point,
float * depth,
Bound b)
{
// Insert the initial bound into the queue
Bound queue[INTERSECTS_QUEUE_MAX];
int queue_num = 1;
int queue_start = 0;
queue[queue_start] = b;
// Iterate over the queue
while (queue_num)
{
// Process the next bound in the queue
Bound c = queue[queue_start];
queue_start = (queue_start + 1) % INTERSECTS_QUEUE_MAX;
queue_num--;
// Compute the bound center and radius
Vector2 p = BoundCenter(c);
float r = BoundRadius(c);
// Sample each of the SDFs at the bound center
float d0 = sdShape0(p);
float d1 = sdShape1(p);
// Update best depth found
if (maxf(d0, d1) <= 0.0f && d0 + d1 <= *depth)
{
*depth = d0 + d1;
*point = p;
}
// There is no possible intersection or improvement
float minrad = 0.02f;
if (maxf(d0, d1) - r > 0.0f
|| d0 + d1 - 2 * r > *depth
|| r < minrad
|| queue_num + 2 > INTERSECTS_QUEUE_MAX)
{
continue;
}
Bound s0, s1;
BoundSplit(&s0, &s1, c);
queue[(queue_start + queue_num + 0) % INTERSECTS_QUEUE_MAX] = s0;
queue[(queue_start + queue_num + 1) % INTERSECTS_QUEUE_MAX] = s1;
queue_num += 2;
}
}
This naturally requires more evaluations, but gives us a point in the intersection region which we can then project back onto the surface of each SDF using the distance and the normal:
float depth = 0.0f;
Vector2 point = (Vector2){ 0.0f, 0.0f };
sdPenetration(&point, &depth, intersection);
if (depth < 0.0f)
{
float d0 = sdShape0(point);
float d1 = sdShape1(point);
Vector2 n0 = normalShape0(point); // Normal shape0
Vector2 n1 = normalShape1(point); // Normal shape1
Vector2 c0 = Vector2Add(point, Vector2Scale(n0, -d0)); // Surface point on shape0
Vector2 c1 = Vector2Add(point, Vector2Scale(n1, -d1)); // Surface point on shape1
}
Here is how it looks in action:
This is probably not the 100% correct way (if there even is a correct way) to generate contact points, but it seems to work reasonably well in practice.
There is something else to keep in mind here, which is that because the SDFs are non-convex there are actually multiple collision regions from which we might want to generate multiple contact points. If we want this we're going to need a more complex search algorithm which uses a heuristic to determine when there are multiple contact regions. That's a bit out-of-the-scope of this post.
Like with other collision detection algorithms, it can also help to introduce a "margin" and expand all of our shapes by some small fixed distance. This makes contact points and normals smoother and more stable, and also deals with the fact that certain SDF constructions don't have accurate interior distances.
Luckily because we're dealing with SDFs, this is easy - we simply subtract our margins from our SDF distance functions.
Here is a little visualisation:
If we increase our margin arbitrarily we can effectively compute the distance between two SDFs like how we might use the GJK algorithm - although naturally computation cost grows large as the overlapping region grows large.
Here is an adaption to the previous algorithm to finding the minimum distance between two SDFs:
bool sdDistance(
Vector2 * point,
float * dist,
Bound b)
{
// Compute the bound center and radius
Vector2 p = BoundCenter(b);
float r = BoundRadius(b);
// Sample each of the SDFs at the bound center
float d0 = sdShape0(p);
float d1 = sdShape1(p);
// If the sampled point contains an intersection return true
if (maxf(d0, d1) <= 0.0f)
{
*dist = 0.0f;
*point = p;
return true;
}
// Update smallest combined dist found
if (minf(d0, d1) > 0.0f && d0 + d1 <= *dist)
{
*dist = d0 + d1;
*point = p;
}
// There is no possible improvement
float minrad = 0.02f;
if (minf(d0, d1) + r <= 0.0f || d0 + d1 - 2 * r > *dist || r < minrad)
{
return false;
}
Bound s0, s1;
BoundSplit(&s0, &s1, b);
return sdDistance(point, dist, s0) || sdDistance(point, dist, s1);
}
And here it is in action:
I hope that gives you a sense of how this basic algorithm can be adapted for different things. For example, it would not be a difficult extension to the above to make it returns bounds on the distance instead of a single estimate.
It's also worth noting that this same algorithm can be adapted to use interval arithmetic to compute the bounds of the SDF functions within the bounding box. This might work better when you're bounding boxes are long and thin, I'm not sure.
If you want to have a play around with this you can find a web demo here.
And all the source code for this article is here.
I hope you've found that useful!