Hi everybody,

I have implemented a marching tetrahedra algorithm for the isosurface extraction. For the normals, I was only able to compute them per triangle as per vertes I need information about all neighbouring triangle. However, in neighboring cells I have no information about which triangles are adjacent to the cell so I cannot compute the normals per vertex. How is this normally done? Ma marching tetrahedra code is as follows :

	
void marching_tetras::MarchingCubesSplitter(Vec3f verts[8], float interpol[8],int parity, bool wireframe, bool raytracer,vector<triangle> &triangleList) {
    //! Split each cube into 5 tetrahedra
    //! need to alternate orientation of central tetrahedron to ensure the diagonals line up
    if (parity) {
    MarchingTetras(verts[0],verts[5],verts[3],verts[6],interpol[0],interpol[5],interpol[3],interpol[6],wireframe,raytracer,triangleList);
    MarchingTetras(verts[0],verts[1],verts[3],verts[5],interpol[0],interpol[1],interpol[3],interpol[5],wireframe,raytracer,triangleList);
    MarchingTetras(verts[0],verts[4],verts[5],verts[6],interpol[0],interpol[4],interpol[5],interpol[6],wireframe,raytracer,triangleList);
    MarchingTetras(verts[0],verts[2],verts[6],verts[3],interpol[0],interpol[2],interpol[6],interpol[3],wireframe,raytracer,triangleList);
    MarchingTetras(verts[3],verts[7],verts[6],verts[5],interpol[3],interpol[7],interpol[6],interpol[5],wireframe,raytracer,triangleList);
    } else {
    MarchingTetras(verts[2],verts[4],verts[1],verts[7],interpol[2],interpol[4],interpol[1],interpol[7],wireframe,raytracer,triangleList);
    MarchingTetras(verts[5],verts[4],verts[7],verts[1],interpol[5],interpol[4],interpol[7],interpol[1],wireframe,raytracer,triangleList);
    MarchingTetras(verts[2],verts[3],verts[7],verts[1],interpol[2],interpol[3],interpol[7],interpol[1],wireframe,raytracer,triangleList);
    MarchingTetras(verts[4],verts[2],verts[1],verts[0],interpol[4],interpol[2],interpol[1],interpol[0],wireframe,raytracer,triangleList);
    MarchingTetras(verts[2],verts[7],verts[6],verts[4],interpol[2],interpol[7],interpol[6],interpol[4],wireframe,raytracer,triangleList);
    }
}
	void marching_tetras::MarchingTetras(Vec3f a, Vec3f b, Vec3f c, Vec3f d,float val_a, float val_b, float val_c, float val_d, bool wireframe, bool raytracer,vector<triangle> &triangleList) {
    //! figure out how masimulationGrid->getYSize() of the vertices are on each side of the isosurface
    //! If value is smaller than 0 -> not part of fluid. Count is the number of vertices inside the volume
    int count =
    (val_a > m_isosurface) +
    (val_b > m_isosurface) +
    (val_c > m_isosurface) +
    (val_d > m_isosurface);
    
    //! handle each case
    //! All vertices are inside the volume or outside the volume (This time of the tetrahedra, not the cube!)
    if(count==0 || count==4) return;
    
    //! Three vertices are inside the volume
    else if (count == 3) {
        //! Make sure that fourth value lies outside
        if (val_d <= m_isosurface) MarchingTetrasHelper3(a,b,c,d,val_a,val_b,val_c,val_d,wireframe,raytracer,triangleList);
        else if (val_c <= m_isosurface) MarchingTetrasHelper3(a,d,b,c,val_a,val_d,val_b,val_c,wireframe,raytracer,triangleList);
        else if (val_b <= m_isosurface) MarchingTetrasHelper3(a,c,d,b,val_a,val_c,val_d,val_b,wireframe,raytracer,triangleList);
        else MarchingTetrasHelper3(b,d,c,a,val_b,val_d,val_c,val_a,wireframe,raytracer,triangleList);
    }     
    //! Two vertices are inside the volume
    else if (count == 2) {
        //! Make sure that the last two points lie outside
        if (val_a > m_isosurface && val_b > m_isosurface) MarchingTetrasHelper2(a,b,c,d,val_a,val_b,val_c,val_d,wireframe,raytracer,triangleList);
        else if (val_a > m_isosurface && val_c > m_isosurface) MarchingTetrasHelper2(a,c,d,b,val_a,val_c,val_d,val_b,wireframe,raytracer,triangleList);
        else if (val_a > m_isosurface && val_d > m_isosurface) MarchingTetrasHelper2(a,d,b,c,val_a,val_d,val_b,val_c,wireframe,raytracer,triangleList);
        else if (val_b > m_isosurface && val_c > m_isosurface) MarchingTetrasHelper2(b,c,a,d,val_b,val_c,val_a,val_d,wireframe,raytracer,triangleList);
        else if (val_b > m_isosurface && val_d > m_isosurface) MarchingTetrasHelper2(b,d,c,a,val_b,val_d,val_c,val_a,wireframe,raytracer,triangleList);
        else MarchingTetrasHelper2(c,d,a,b,val_c,val_d,val_a,val_b,wireframe,raytracer,triangleList);
    }
    //! One vertex is inside the volume
    else if (count == 1) {
        //! Make sure that the last three points lie outside
        if (val_a > m_isosurface) MarchingTetrasHelper1(a,b,c,d,val_a,val_b,val_c,val_d,wireframe,raytracer,triangleList);
        else if (val_b > m_isosurface) MarchingTetrasHelper1(b,c,a,d,val_b,val_c,val_a,val_d,wireframe,raytracer,triangleList);
        else if (val_c > m_isosurface) MarchingTetrasHelper1(c,a,b,d,val_c,val_a,val_b,val_d,wireframe,raytracer,triangleList);
        else MarchingTetrasHelper1(d,c,b,a,val_d,val_c,val_b,val_a,wireframe,raytracer,triangleList);
    }
}
	//! when 3 vertices are within the volume
void marching_tetras::MarchingTetrasHelper3(Vec3f a, Vec3f b, Vec3f c, Vec3f d, float val_a, float val_b, float val_c, float val_d,bool wireframe, bool raytracer,vector<triangle> &triangleList) {
    //! Compute intersection points of isosurface with tetrahedron
    float ad_frac = (m_isosurface - val_d) / (val_a - val_d);
    float bd_frac = (m_isosurface - val_d) / (val_b - val_d);
    float cd_frac = (m_isosurface - val_d) / (val_c - val_d);
    
    //! Compute the vertices of the intersections of the isosurface with the tetrahedron
    Vec3f ad = ad_frac*a + (1-ad_frac)*d;
    Vec3f bd = bd_frac*b + (1-bd_frac)*d;
    Vec3f cd = cd_frac*c + (1-cd_frac)*d;
    
    //! Render the triangle
    storeTriangle(ad,cd,bd);    
}
	//! when 2 vertices are within the volume
void marching_tetras::MarchingTetrasHelper2(Vec3f a, Vec3f b, Vec3f c, Vec3f d,float val_a,float val_b,float val_c,float val_d,bool wireframe,bool raytracer,vector<triangle> &triangleList) {
    //! Compute the intersection points of isosurface with tetrahedron
    float ac_frac = (m_isosurface - val_c) / (val_a - val_c);
    float ad_frac = (m_isosurface - val_d) / (val_a - val_d);
    float bc_frac = (m_isosurface - val_c) / (val_b - val_c);
    float bd_frac = (m_isosurface - val_d) / (val_b - val_d);
    
    //! Compute the vertices of the TWO triangles which mark the intersection of the isosurface with the tetrahedron
    Vec3f ac = ac_frac*a + (1-ac_frac)*c;
    Vec3f ad = ad_frac*a + (1-ad_frac)*d;
    Vec3f bc = bc_frac*b + (1-bc_frac)*c;
    Vec3f bd = bd_frac*b + (1-bd_frac)*d;
    
    //! Render the two triangles
    storeTriangle(ac,bc,ad);
    storeTriangle(bc,bd,ad);    
}
	//! when 1 vertex is within the volume
void marching_tetras::MarchingTetrasHelper1(Vec3f a, Vec3f b, Vec3f c, Vec3f d,float val_a,float val_b,float val_c,float val_d, bool wireframe, bool raytracer,vector<triangle> &triangleList) {
    //! Compute the intersection points of isosurface with tetrahedron
    float ab_frac = (m_isosurface - val_b) / (val_a - val_b);
    float ac_frac = (m_isosurface - val_c) / (val_a - val_c);
    float ad_frac = (m_isosurface - val_d) / (val_a - val_d);
    
    //! Compute the vertices of the intersections of the isosurface with the tetrahedron
    Vec3f ab = ab_frac*a + (1-ab_frac)*b;
    Vec3f ac = ac_frac*a + (1-ac_frac)*c;
    Vec3f ad = ad_frac*a + (1-ad_frac)*d;
    
    //! Render the triangle
    storeTriangle(ab,ad,ac);        
}
	

First off I've never done this stuff on the GPU so probably someone else can chime in if they know better. What I do on the CPU is I have a normal accumulator in each vertex which starts out at (0,0,0). Then I go though all the triangles and add in their normals for each of their 3 vertices.  If you want, you can weight the normal you are adding in by the the angle in that corner of the triangle.  Then at the end I go through all the vertexes and normalize the normals. I'm assuming this is kind of a standard way to do it but I don't really know.

In your case on the GPU if you can feed in all your triangles on subsequent shading pass, followed by another pass for the vertexes, you should be able to do the same thing.

You could calculate the normal from volume data directly by calculating the gradient, avoiding the need for a second pass. https://en.wikipedia.org/wiki/Gradient

2 minutes ago, JoeJ said:

You could calculate the normal from volume data directly by calculating the gradient, avoiding the need for a second pass. https://en.wikipedia.org/wiki/Gradient

That's true but it's sometimes a pain in the ass depending on where your data comes from. Once you start stacking functions and intersecting them it can get rather difficult.  Calculating it on the final mesh works in all cases.

Hi everyone,

thank you for your comments. Maybe I described the problem in a wrong way. I am basically going through the entire grid and extract the triangles for each cell separately. Then I can compute the normal for each triangle separately. But in order to compute the normals per vertex, I need to know for adjacent cells which triangles share a vertex on the edge. Currently I don't have this information as I am processing cell by cell and do interpolations for each cell separately. So the question is how can I find out which triangles belong to a certain vertex. Currently I have only the resulting triangles stored with 3 vertices each, but this is not the correct way I think...

Hi Thomas, Sorry, for some stupid reason I thought you were dong this on the GPU. In any case are your vertexes shared? I mean for each box that shares a vertex generated by marching cubes are they shared between adjacent cubes or do you generate them separately?

Without analyzing your code, the way I handle this is to generate all the vertexes on edge structures that are shared by the boxes.  For instance each box can have indexes to each of the 12 edge structures which can hold a triangle vertex if it exists. I need indexes because I'm using an octree of boxes so my cells aren't all of the same size. However if your boxes are all uniform you can just have a matrix of horizontal edges and one of vertical edges and index them implicitly from each box by their position.  Now you generate all try vertexes in the matrices of edge structures and they can be accessed by the boxes so a vertex is only generated once.  The box can ask the edge for it's vertex and the edge will either generate it if it needs to, or simply return the one it already has stored.  So now you build your triangles as a connected mesh.

Now as you build your triangles you have a normal accumulator included in each vertex that you created, and you add in the normal component of the triangle to each of it's three vertex accumulators. So now you have a normal at each vertex. All you need to do is go though all your vertexes in a final pass and normalize them. And that should be it.

Hi Gnollrunner,

thank you so much for your help. I tried to implement you suggestion and have now structures representing the edges. As I am working on a unified grid this is really simple. Each edge can hold an infinite amount of vertices, though in practice the number will be rather small. However, when a new vertex is computed by the marching cubes alrotihm, I need to check whether it is already contained in the edge. Currently I am comparing the x, y and z coordinates of the vertices with the == operator, but a floting point comparison will not always work. How did you solve this? Do you take a range and if the vertex lies within that range you assume it is the same vertex? For me sometimes this comparison fails although the vertex should be already in the list. Could you show me the correct code which you use to check whether the vertex is already contained in the list?

Akkk! I forgot you were doing tetrahedra and not cubes. Now I feel really stupid ?.  In any case I think it can work roughly the same way but it will be slightly more complex . I think you will have to add a few more matrices of edges.  Maybe this can work for you ...

If X, Y, and Z are the dimensions of your 3d cube, each sub cube having 6 tetrahedron. We could have:

A Matrix for tetrahedra of dimensions [X][Y][Z][6] ......... the 6 is for each of the six tetrahedra in each box
A Matrix for vertical edges of dimensions [X+1][Y+1][Z]  ....... the +1s because there is one more row and column of edges than cells
A Matrix for horizontal edges Running in the X direction of dimensions [X+1[Y][Z+1]  ..... again X+1 and Z+1 for the end plates
A Matrix for horizontal edges running in the Y direction of dimensions [X][Y+1][Z+1] ....... again for the end plates

......

Then you need to do something similar for the diagonals on each face so that's another 3 matrices for faces oriented the 3 different directions.
And you need one of dimensions [X][Y][Z] for the diagonals running though the middle of the boxes

Finally you can have a matrix of [X+1][Y+1][Z+1] for your function values if you aren't generating them on the fly.

Now given a coordinate for for a specific cell [X][Y][Z][Sub-cell] you can use a switch statement with six cases that can get you all 4 edges out of the correct tables and indexes for the values for their endpoints.  I'm not sure what you meant by "Each edge can hold an infinite amount of vertices" because all you need is one unless you are doing multiple functions over the matrix..... If you are dong things right you shouldn't need to compare any floating point values with each other . It should generate a complete mesh already put together.

Sorry if this is confusing, it's sometimes hard to describe stuff in text. I can whip out some stub data structures and functions if you want to see it in code.

Have you guys considered / tried Surface Nets as an alternative? https://0fps.net/2012/07/12/smooth-voxel-terrain-part-2/

I hope this would be easier and maybe faster but have not tried any of those algorithms. My interest is to get a robust manifold approximation from game models which may consist of multiple pieces or bad topology.