How to get normals from volume density grid?

Question

I am meshing my volume density data with marching cubes. The density data is a flat array and it's elements contain a) distance from surface, b) a normal (xyz). However this normal is for the density grid point, no the calculated vertexes. So having this data can i calculate normal for each vertex?

Answer 1

To find the surface normal for any given marched cube, just average the face normals for all of the 1-5 triangles that are generated. A face normal is calculated here:

for (vector<triangle>::const_iterator i = triangles.begin(); i != triangles.end(); i++)
{
    // Get face normal.
    vertex_3 v0 = i->vertex[1] - i->vertex[0];
    vertex_3 v1 = i->vertex[2] - i->vertex[0];
    normal = v0.cross(v1);
    normal.normalize();
    
    ...

Don't forget to sort the vertices before you perform the linear interpolation. If you use floats, which appears so, then you will experience cracks in the mesh. Taubin smoothing of the mesh turns those cracks into holes. If you sort the two vertices beforehand, then no cracks appear. If you use doubles, then no cracks appear. Just my 2 cents.

vertex_3 marching_cubes::vertex_interp(const float isovalue, vertex_3 p1, vertex_3 p2, float valp1, float valp2)
{
    // Sort the vertices so that cracks don't mess up the water-tightness of the mesh.
    // Note: the cracks don't appear if you use doubles instead of floats.
    if (p2 < p1)
    {
        vertex_3 tempv = p2;
        p2 = p1;
        p1 = tempv;

        float tempf = valp2;
        valp2 = valp1;
        valp1 = tempf;
    }

    const float epsilon = 1e-10f;

    if(fabs(isovalue - valp1) < epsilon)
        return(p1);

    if(fabs(isovalue - valp2) < epsilon)
        return(p2);

    if(fabs(valp1 - valp2) < epsilon)
        return(p1);

    float mu = (isovalue - valp1) / (valp2 - valp1);

    return p1 + (p2 - p1)*mu;
}

where:

inline bool operator<(const vertex_3 &right) const
{
    if(right.x > x)
        return true;
    else if(right.x < x)
        return false;

    if(right.y > y)
        return true;
    else if(right.y < y)
        return false;

    if(right.z > z)
        return true;
    else if(right.z < z)
        return false;

    return false;
}

Answer 2

I ended up taking the volume outside point normal and inside point normal and doing a lerp between them, where the lerp value is density float. This is how I get normals from SDF while creating volume gird:

float H = 0.001;
float dx = sdfFunc(Vector3(xPos + H, yPos, zPos)) - sdfFunc(Vector3(xPos - H, yPos, zPos));
float dy = sdfFunc(Vector3(xPos, yPos + H, zPos)) - sdfFunc(Vector3(xPos, yPos - H, zPos));
float dz = sdfFunc(Vector3(xPos, yPos, zPos + H)) - sdfFunc(Vector3(xPos, yPos, zPos - H));

float dlen = sqrt(dx * dx + dy * dy + dz * dz);

if (dlen > 0.0) {
  float mul = 1.0 / dlen;
  dx *= mul;
  dy *= mul;
  dz *= mul;
}

And when calculating vertices with marching cubes I do something like that:

Vector3 nn = Vector3();

nn.X = (1 - t) * an.X + t * bn.X;
nn.Y = (1 - t) * an.Y + t * bn.Y;
nn.Z = (1 - t) * an.Z + t * bn.Z;

where an is inside point and bn is outside point normal and t is calculated density float (vertex place along he edge).