Signed Distance Function (SDF): Implicit curves or surfaces

A mathematical equation usually involves the relationship between several unknowns. For exemple:

In this case, the relation between

is implicit because none of then are expressed in terms of the other. In case we clear one of the unknowns in terms of the other, this become an explicit equation for this unknown

2D SDF: Distance to a given point

When you consider an implicit equation and you equals it to zero

the set of points

that fulfill this equation defines a curve in

(a surface in

). In our equation it corresponds to the set of points at distance 1 of the point

, that is, a circle. This is the implicit curve defined by our equation, but we can also equals the implicit equation to a value d.

Then, this is another implicit curve, now at a distance

of the point

. Or what is the same, the set of points at a distance d of the circle.

For simplicity, let's consider

the origin of coordinates and the signed distance function defined for all points in

by:

If we evaluate

in a grid of points defined in a chosen domain, we have (see the figure):

The point is outside of the circle
The point is inside of the circle
The point is on of the circle

The points on the unity cercle fulfill

x = -1.5:0.1:1.5;

y = -1.5:0.1:1.5;

[X,Y]=meshgrid(x,y);

d=sqrt(X.^2+Y.^2)-1;

dValues = [0,0]; %range of values of d to be plotted (only d=0 in this case)

figure

contour(X,Y,d,dValues)

grid on

axis equal

dValues = [-0.2,0,0.3]; %range of values of d to be plotted (three curves)

figure

contour(X,Y,d,dValues,'ShowText','on') %show level curve

grid on

axis equal

Signed Distance Function 2D: Distance to a segment

In the case of a segment we have to consider the distance according of where the point

we are measuring with respect to the segment is. In fact we have to compute which is the closest point

in the segment and compute the distance between the two points.

Let's assume the segment is defined between two points A and B, then we can parametrize the segment as the set of points

the point belongs to the segment and for

you are in the same line but outside the segment (see the figure)

In order to compute the value of

, we only need to project the vector

to the line defined by point

and direction

and from this value, the closest point is

(here <, > means dot product).

According to the value of

, you will be in region 1, 2 or 3 (see figure) of the plane. Thus, we can have

if (region 2)
if (region 3)
if (region 1)

You can code this in Matlab using some if conditions or you can obtain a compactify version of the formula according to that:

Because correspond to point A, the distance when is the same as the distance for , so we can clamp the values of using
The same argument can be applied to point B now with and the values . So we can clamp now
Combining the previous comments, we can have only one final formula for our computation:

. Then

and

A = [0,0];

B = [2,2];

x = -1:0.1:3;

y = -1:0.1:3;

[X,Y] = meshgrid(x,y);

BmA = B-A;

den = BmA*BmA';

d = zeros(size(X));

for i = 1:numel(X)

P = [X(i),Y(i)];

lamb = dot(P-A, BmA)/den;

lambt = min(1, max(0, lamb));

Q = A+lambt*BmA;

d(i) = norm(P-Q);

end

dValues = [0, 0.1, 0.2, 0.3, 0.4]; %range of values of d to be plotted (d=0 in this case)

figure

contour(X,Y,d,dValues)

grid on

axis equal

Signed Distance Function 3D: Distance to a segment

The same formulation of the case 2D can be implemented in 3D. In fact, all the formulas are vectorial formulas and are independent of the number of coordinates. You only need to consider a third coordinate for all the points, the procedure is exactly the same.

In Matlab the function to get the contours is named isosurface. Below you can see how to use it.

A = [0,0,0];

B = [2,2,2];

x = -1:0.05:3;

y = -1:0.05:3;

z = -1:0.05:3;

[X,Y,Z] = meshgrid(x,y,z);

BmA = B-A;

den = BmA*BmA';

d = zeros(size(X));

for i = 1:numel(X)

P = [X(i),Y(i),Z(i)];

lamb = dot(P-A, BmA)/den;

lambt = min(1, max(0, lamb));

Q = A+lambt*BmA;

d(i) = norm(P-Q);

end

dVal=0.5;

figure

isosurface(X,Y,Z,d,dVal);

axis equal

view([28.36 -33.00])

Signed Distance Function 3D: Distance to a Box

This problem is similar to the distance of a segment but here we have to take into account symmetries and two segment. We plot in the figure the 2D case and implement directly the 3D one.

Here we have

for points in region 1
for points in region 3
for points in region 1

The formula to unify regions 1, 2 3 can be expressed as:

, where means the absolute value of its coordinates.

The general formula can be derived to include also internal points in region 4 and is expressed as: (this is not as intuitive as the segment example)

where

is defined as the absolut value of the maximum component of a vector.

We can apply this formulas directly to 3D vectors

% Centered 3D Box

Rx=0.5;

Ry=0.5;

Rz=0.5;

R=[Rx,Ry,Rz];

x = -2:0.05:2;

y = -2:0.05:2;

z = -2:0.05:2;

[X,Y,Z] = meshgrid(x,y,z);

d = zeros(size(X));

for i = 1:numel(X)

P = [X(i),Y(i),Z(i)];

Q = abs(P)-R;

d(i) = norm(max(Q-R,0))-min(norm(abs(P)-R,"inf"),0);

end

dVal=0;

figure

isosurface(X,Y,Z,d,dVal);

axis equal

view([61.72 21.60])

dVal = 0.2; % now plot a little bit out of the zero surface

figure

isosurface(X,Y,Z,d,dVal);

axis equal

view([61.72 21.60])

Mixture of several objects

We now consider plotting the surface corresponding to many objects. When using signed distance functions

to two different bodies

and

the following properties apply (see the figure):

gives the union of the two bodies
gives the intersection of the two bodies
gives the difference of the two bodies

Below we show an example using spheres and segments.

First sphere:

Let's now consider the distance to a point in 3D (in fact we consider square distance), that is a sphere. Like in the 2D case we can compute. We also introduce some lighting properties in order to obtain better visualizations.

% here we define the implicit function directly

x = -3:0.1:3;

y=x;

z=x;

[X,Y,Z]=meshgrid(x,y,z);

c = [0,0,0.5]; %center of the sphere

r = 1; %radius

d1 = zeros(size(X));

for i=1:numel(X)

p=[X(i),Y(i),Z(i)];

d1(i)=(p(1)-c(1)).^2 +(p(2)-c(2)).^2 +(p(3)-c(3)).^2 - r^2;

end

dVal=0.0;

figure

patch(isosurface(X,Y,Z,d1,dVal,X),...

'FaceColor','interp','EdgeColor','none')

view(30,-15)

axis equal

grid

colormap gray

Second sphere:

Now we add a second small sphere and we consider the union of the two spheres using the new distance

% Add a new small sphere

c = [0,0,-0.5]; %center of the sphere

r = 0.5; %radius

d2 = zeros(size(X));

for i=1:numel(X)

p=[X(i),Y(i),Z(i)];

d2(i)=(p(1)-c(1)).^2 +(p(2)-c(2)).^2 +(p(3)-c(3)).^2 - r^2;

end

dd=min(d1,d2);

dVal=0.2;

figure

patch(isosurface(X,Y,Z,dd,dVal,X),...

'FaceColor','interp','EdgeColor','none')

view(30,-15)

axis equal

grid

colormap gray

Two Spheres and a Segment :

Finally we add a third object from the distance to a given segment using the new distance

A = [-1,0,-1.1];

B = [1,0,-1.1];

BmA = B-A;

den = BmA*BmA';

d3 = zeros(size(X));

for i = 1:numel(X)

P = [X(i),Y(i),Z(i)];

lamb = dot(P-A, BmA)/den;

lambt = min(1, max(0, lamb));

Q = A+lambt*BmA;

d3(i) = norm(P-Q)-0.2;

end

dd=min(d1,min(d2,d3));

dVal=0.2;

figure

patch(isosurface(X,Y,Z,dd,dVal,X),...

'FaceColor','interp','EdgeColor','none')

view(30,-15)

axis equal

grid

colormap gray

view([-51.53 19.80])

Substractting two bodies:

You can also subtract two bodies if you consider the distance function

% Add a new small sphere

dd=max(d1,-d2);

dVal=0.02;

figure

patch(isosurface(X,Y,Z,dd,dVal,X),...

'FaceColor','interp','EdgeColor','none')

view(30,-15)

axis equal

grid

colormap gray

%note: the boundary artifacts are due to the size of the discretization in

%the volume points (try a high resolution for improve).

If you want to have fun, I suggest you watch the video Painting a Selfie Girl, with Maths from Iñigo Quilez. It is an amazing application of these techniques.