# How to get normals from volume density grid?

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?

• What kind of results do you get when you simply interpolate the grid normals according to the vertex position? Apr 8, 2021 at 11:42
• I'm not sure how to do that. Apr 8, 2021 at 11:58
• Just a standard tricubic interpolation? Apr 8, 2021 at 12:02
• In MC the vertex is placed on the edge, so I can relate only to the points that make up the edge. Apr 8, 2021 at 12:15
• You can try that. What results do you get? Apr 8, 2021 at 12:16

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;
}

• This does not appear to answer the question "How to get normals" Apr 10, 2021 at 18:06
• Apologies. I have not enough reputation to leave comments in general. I'll add to it. Apr 11, 2021 at 2:10

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