Exact signed distance function of a rhombic dodecahedron

Question

I'm playing around with a ray marcher engine with voxels shaped as rhombic dodecahedra since they can tessallate the 3D space perfectly like hexagons do for the 2D space. The thing is that I only managed to come up with an approximation of the signed distance function for that solid. It's okay to render the voxels but it gives really sharp edges. I would like to make them look better with some rounded edges using an exact SDF.

Does anyone know the exact SDF of a rhombic dodecaheron ?

Here's my (not exact) SDF :

float Sdf(Voxel voxel, vec3 position) {
    position = position - voxel.Position;
    // Exploit the solid's symmetries
    position = vec3(abs(position.x), sign(position.z) * position.y, abs(position.z));

    // distance to each face
    float a = dot(vec3(1.  ,   0.           , 0.          ), position) - voxelSize;
    float b = dot(vec3(0.5 ,   0.           , sqrt3 / 2.  ), position) - voxelSize;
    float c = dot(vec3(0.  , - sqrt2 / sqrt3, 1. / sqrt3  ), position) - voxelSize;
    float d = dot(vec3(0.5 ,   sqrt2 / sqrt3, sqrt3 / 6.  ), position) - voxelSize;
    float e = dot(vec3(0.5 , - sqrt2 / sqrt3, - sqrt3 / 6.), position) - voxelSize;

    return max(a, max(b, max(c, max(d, e))));
}

And this is how it looks like:

EDIT :

So I had some new ideas on how to approach this problem. In this new technique, I try to find the different areas that would have different distance formulas (depending on whether the closest point to the solid is on a face, an edge or a vertex) and compute their distance separately. So far, I only have a correct separation and distance for the areas where the closest point is on a face or on an edge. The code looks like this :

float Sdf(Voxel voxel, vec3 position){
    position = position - voxel.Position;

    vec3 normals[6];
    float dists[6];
    float signs[6];

    // The rhombic dodecahedron has 6 pairs of opposite faces.
    // Assign a normal to each pair of faces.
    normals[0] = vec3(1.  ,   0.           , 0.          );
    normals[1] = vec3(0.5 ,   0.           , sqrt3 / 2.  );
    normals[2] = vec3(0.5 ,   0.           , - sqrt3 / 2.);
    normals[3] = vec3(0.  , - sqrt2 / sqrt3, 1. / sqrt3  );
    normals[4] = vec3(0.5 ,   sqrt2 / sqrt3, sqrt3 / 6.  );
    normals[5] = vec3(0.5 , - sqrt2 / sqrt3, - sqrt3 / 6.);

    // Compute the distance to each face (the sign tells which face of the pair of faces is closest)
    for (int i = 0; i < 6; i++){
        dists[i] = dot(normals[i], position);
        signs[i] = sign(dists[i]);
        dists[i] = max(0., abs(dists[i]) - voxelSize);
    }

    bool sorted = false;
    while(!sorted){
        sorted = true;
        for (int i = 0; i < 5; i++){
            if (dists[i] < dists[i+1]){
                vec3 n = normals[i];
                float d = dists[i];
                float s = signs[i];
                normals[i] = normals[i+1];
                dists[i] = dists[i+1];
                signs[i] = signs[i+1];
                normals[i+1] = n;
                dists[i+1] = d;
                signs[i+1] = s;
                sorted = false;
            }
        }
    }
    
    if (dists[1] < dists[0] * 0.5){ //0.5 comes from sin(pi/6)
        // case where the closest point is on a face
        return dists[0];
    }
    else if (dists[2] == 0.){
        // case where the closest point is on an edge
        float a = 0.5 * dists[0];
        float b = 0.5 * (dists[1] - a);
        return length((dists[0] - b) * signs[0] * normals[0]
                    + (dists[1] - a) * signs[1] * normals[1]);
    }
    else {
        //dunno
    }
}

Here's a quick sketch to show how my approach

The remaining cases get crazier. I'll update the post if any breakthroughs are make.

Answer 1

Preface

Assume the rhombic dodecahedron is oriented as seen here.

This orientation makes it be symmetric with respect to mirroring along all three coordinate planes.

What you then need, is an sdf that works for this shape in the (+, +, +) octant.

Alternative approximation

The shape can be seen as the union of three planes. Utilising the manhattan distance, it it possible to define an alternative to the approximated sdf:

float Sdf(float3 position, float size)
{
    position = abs(position);

    float halfsqrt2 = sqrt(2) / 2;

    float side1 = (position.x + position.y - size) * halfsqrt2;
    float side2 = (position.x + position.z - size) * halfsqrt2;
    float side3 = (position.y + position.z - size) * halfsqrt2;

    return max(max(side1, side2), side3);
}

This can be abbreviated to:

float Sdf(float3 position, float size)
{
    position = abs(position);

    return (max(max(position.x + position.y, position.x + position.z), position.y + position.z) - size) * sqrt(2) / 2;
}

Exact sdf

The shape can also be seen as a transformed cube. The matrices to transform from one to the other are

0 1 1
1 0 1
1 1 0

and it's inverse

-0.5  0.5  0.5
 0.5 -0.5  0.5
 0.5  0.5 -0.5

(There is multiple matrices that work, but this is the one I used)

After transforming, it is possible to use the sdf of a rounded box. Here is a code example:

float SdfRoundRhombicDodecahedron(float3 position, float size, float edgeSize)
{
    float3x3 w2r = float3x3(
        0, 1, 1,
        1, 0, 1,
        1, 1, 0
    );

    position = abs(position);
    position = mul(w2r, position);
    position -= clamp(position, 0, size - edgeSize);
    return length(position) - edgeSize;
}

The transformed and mirrored rounded cube then has the shape of a rhombic dodecahedron. Here is a visualisation of it made in Unity