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

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

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

If the point belongs to the segment and for or 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

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

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

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.

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

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

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

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.

© Numerical Factory 2021 (Agustí Roig in memoriam, who suggested me to implement this script for our students).