Jump to content

  • Log In with Google      Sign In   
  • Create Account


A C++ code to smooth (and fix cracks in) meshes generated by the standard, original Marching Cubes algorithm


Old topic!
Guest, the last post of this topic is over 60 days old and at this point you may not reply in this topic. If you wish to continue this conversation start a new topic.

  • You cannot reply to this topic
6 replies to this topic

#1 taby   Members   -  Reputation: 336

Like
3Likes
Like

Posted 22 April 2013 - 03:00 PM

In the past, some people have asked for code or help with applying a smoothing algorithm to the meshes they get from the Marching Cubes algorithm (see: 'Polygonising a scalar field' by P. Bourke).

 

I've attached some sample C++ code to perform Taubin smoothing (see: 'Geometric Signal Processing on Polygonal Meshes' by G. Taubin and 'Implicit Fairing of Irregular Meshes using Diffusion and Curvature Flow' by M. Desbrun, et al.) on a Stereo Lithography (STL) file. The code assumes that the mesh is closed -- Attached File  mc_smooth2.zip   8.27KB   181 downloads

 

The code also eliminates the cracks in the mesh that are unfailingly generated by the standard, original Marching Cubes algorithm (Google for 'marching cubes cracks' to see that it is a common problem).

 

The cracks must be eliminated before smoothing occurs, otherwise the cracks will become holes, and the mesh will be ruined.

 

Attached is a sample mesh (including cracks) -- Attached File  fractal.zip   826.28KB   97 downloads

 

Here is an image of the sample mesh as generated by Marching Cubes, rendered using flat shading -- very blocky:

1_raw_flat_shaded.png

 

Here is an image of the cracks highlighted in red, rendered using flat shading:

3_raw_cracks.png

 

Here is an image of the sample mesh after 10 iterations of Taubin smoothing, rendered using flat shading -- the cracks have now become holes, thanks to the smoothing. This is bad:

4_taubin_smoothed_holes.png

 

Here is an image of the vertex pairs (highlighted in blue, and connected by green lines) that should have been merged -- in order to fix the cracks -- before the smoothing was applied, rendered using flat shading:

5_taubin_smoothed_merge_vertices.png

 

Here is an image of the crackless mesh after 10 iterations of Taubin smoothing, rendered using flat shading -- no holes, looks less blocky:

6_taubin_smoothed_flat_shaded.png

 

Here is an image of the crackless mesh after 10 iterations of Taubin smoothing, rendered using Gouraud shading -- no holes, looks much less blocky:

7_taubin_smoothed_gouraud_shaded.png

 

Good luck.

 

P.S. The sample mesh is an isosurface of the escape value (where isovalue = 4.0) for the quaternion Julia set generated by Z' = Z^2 + C where C = {0.3, 0.5, 0.4, 0.2} and the maximum number of iterations is 8. (see: 'Quaternion Julia Fractals' by P. Bourke). The code to generate this mesh can be found at: http://code.google.com/p/qjs-isosurface/downloads/detail?name=qjs-isosurface-v2.0.2.1.zip


Edited by taby, 23 April 2013 - 08:40 PM.


Sponsor:

#2 kloffy   Members   -  Reputation: 895

Like
1Likes
Like

Posted 22 April 2013 - 08:12 PM

Nice work! This could surely come in handy at some point.



#3 Bacterius   Crossbones+   -  Reputation: 8506

Like
4Likes
Like

Posted 22 April 2013 - 09:07 PM

Mmh, I never had cracks in Marching Cubes with constant level of detail, at least I don't remember having any. That said, cracks do occur when you try and implement LOD, in that case a smoothing or stitching algorithm is required (another alternative is the transvoxel algorithm by E. Lengyel, though I am not sure if it's patented, it is iirc)


Edited by Bacterius, 22 April 2013 - 11:54 PM.

The slowsort algorithm is a perfect illustration of the multiply and surrender paradigm, which is perhaps the single most important paradigm in the development of reluctant algorithms. The basic multiply and surrender strategy consists in replacing the problem at hand by two or more subproblems, each slightly simpler than the original, and continue multiplying subproblems and subsubproblems recursively in this fashion as long as possible. At some point the subproblems will all become so simple that their solution can no longer be postponed, and we will have to surrender. Experience shows that, in most cases, by the time this point is reached the total work will be substantially higher than what could have been wasted by a more direct approach.

 

- Pessimal Algorithms and Simplexity Analysis


#4 Eric Lengyel   Crossbones+   -  Reputation: 2260

Like
3Likes
Like

Posted 22 April 2013 - 11:51 PM

another alternative is the transvoxel algorithm by E. Lengyel, though I am not sure if it's patented, it is iirc

 

The Transvoxel algorithm is not patented. More information here:

 

http://www.terathon.com/voxels/

 

In general, a correctly implemented Marching Cubes algorithm generating a mesh with a single LOD will only produce cracks if it doesn't have a consistent way of choosing polarity for the so-called "ambiguous cases". This can be solved by using a fixed polarity based on corner states or some face-level choice function as used in the MC33 algorithm. See Section 3.1.2 of my dissertation at the above link for some discussion of these. Using fixed polarity is easy, and it never generates any holes in the resulting mesh.

 

A good MC implementation will generate a smooth mesh to begin with if the data at each voxel location has a range of values instead of just a binary in or out state. The ugly stair-stepping only shows up if you're forced to put each vertex right in middle of each isosurface-crossing edge because you don't have enough information to do anything more intelligent.



#5 taby   Members   -  Reputation: 336

Like
0Likes
Like

Posted 23 April 2013 - 10:36 AM

Yes, I am using the plain vanilla, ancient, original MC that produces topological inconsistencies (writing that makes it seem worse than it is, because, really, the cracks generally seem to be just that; cracks, not gaping holes; each boundary consists of an even number of edges where the majority of the angle is distributed amongst all but two of the vertices) -- I've once tried the version with topological guarantees / case ambiguity resolution by Lewiner, et al. back in the mid-200x's, but I went with the version from Bourke's site in the end because it does a relatively decent job with relatively minimal headache.

 

I updated the first post to indicate that I'm using the standard, original MC algorithm (see: P. Bourke's 'Polygonising a scalar field').

 

If anyone has a full, cost-free, public domain C++ implementation of an adaptive MC algorithm to share, I'd be more than happy to absorb it into my toolkit -- perhaps GD could consider the possibility of adding in a 'recipes' section of the site to gather all of these kinds of things into one nicely organized spot. Apparently I implemented one when I was experimenting back in the mid-200x's and I forgot.

 

In any case, the major point of the post was mesh smoothing -- if anyone has any kind of implementation of Taubin smoothing that uses 'Fujiwara' (scale-dependent) or curvature normal (cotan) weighting, but doesn't destroy the mesh, that would be nice too. I've had no luck with these; the literature is abundant (to the point where it's conflicting).


Edited by taby, 24 April 2013 - 09:18 PM.


#6 taby   Members   -  Reputation: 336

Like
0Likes
Like

Posted 23 April 2013 - 10:51 AM

Nice work! This could surely come in handy at some point.

 

Thank you. I hope it comes it handy.

 

The paper 'Geometric Signal Processing on Polygonal Meshes' by G. Taubin mentions a few edge weight schemes:

- constant unit weight: w_ij = 1

- inverse edge length weight: w_ij = 1 / |e_ij|

- curvature normal-esque (cotan) weight: w_ij = 1 / tan(theta_ij0) + 1 / tan(theta_ij1)

 

The code in the zip file from the first post uses the constant unit weight scheme (see indexed_mesh::laplace_smooth()).

 

I didn't put the other schemes into the code zip file because I'm not sure if I'm implementing them correctly, but I'll include them here in case anyone wishes to play with them to see if/where I'm going wrong.

 

Here's the code to replace the laplace_smooth function with the inverse edge length weight scheme:

 

void indexed_mesh::inverse_edge_length_smooth(const float scale)
{
    vector<vertex_3> displacements(vertices.size(), vertex_3(0, 0, 0));

    // Get per-vertex displacement.
    for(size_t i = 0; i < vertices.size(); i++)
    {
        // Skip rogue vertices (which were probably made rogue during a previous
        // attempt to fix mesh cracks).
        if(0 == vertex_to_vertex_indices[i].size())
            continue;

        vector<float> weights(vertex_to_vertex_indices[i].size(), 0.0f);

        // Calculate weights based on inverse edge lengths.
        for(size_t j = 0; j < vertex_to_vertex_indices[i].size(); j++)
        {
            size_t neighbour_j = vertex_to_vertex_indices[i][j];

            float edge_length = vertices[i].distance(vertices[neighbour_j]);

            if(0 == edge_length)
                edge_length = numeric_limits<float>::epsilon();

            weights[j] = 1.0f / edge_length;
        }

        // Normalize the weights so that they sum up to 1.
        float s = 0;

        for(size_t j = 0; j < weights.size(); j++)
            s += weights[j];

        if(0 == s)
            s = numeric_limits<float>::epsilon();

        for(size_t j = 0; j < weights.size(); j++)
            weights[j] /= s;

        // Sum the displacements.
        for(size_t j = 0; j < vertex_to_vertex_indices[i].size(); j++)
        {
            size_t neighbour_j = vertex_to_vertex_indices[i][j];
            displacements[i] += (vertices[neighbour_j] - vertices[i])*weights[j];
        }
    }

    // Apply per-vertex displacement.
    for(size_t i = 0; i < vertices.size(); i++)
        vertices[i] += displacements[i]*scale;
}

 

Here's the code to replace the laplace_smooth function with the curvature normal-esque weight scheme (it works most times, but no guarantees):

 

void indexed_mesh::curvature_normal_smooth(const float scale)
{
    vector<vertex_3> displacements(vertices.size(), vertex_3(0, 0, 0));

    // Get per-vertex displacement.
    for(size_t i = 0; i < vertices.size(); i++)
    {
        if(0 == vertex_to_vertex_indices[i].size())
            continue;

        vector<float> weights(vertex_to_vertex_indices[i].size(), 0.0f);

        size_t angle_error = 0;

        // For each vertex pair (ie. each edge),
        // calculate weight based on the two opposing angles (ie. curvature normal scheme).
        for(size_t j = 0; j < vertex_to_vertex_indices[i].size(); j++)
        {
            size_t angle_count = 0;

            size_t neighbour_j = vertex_to_vertex_indices[i][j];

            // Find out which two triangles are shared by the edge.
            for(size_t k = 0; k < vertex_to_triangle_indices[i].size(); k++)
            {
                for(size_t l = 0; l < vertex_to_triangle_indices[neighbour_j].size(); l++)
                {
                    size_t tri0_index = vertex_to_triangle_indices[i][k];
                    size_t tri1_index = vertex_to_triangle_indices[neighbour_j][l];

                    // This will occur twice per edge.
                    if(tri0_index == tri1_index)
                    {
                        // Find the third vertex in this triangle (the vertex that doesn't belong to the edge).
                        for(size_t m = 0; m < 3; m++)
                        {
                            // This will occur once per triangle.
                            if(triangles[tri0_index].vertex_indices[m] != i && triangles[tri0_index].vertex_indices[m] != neighbour_j)
                            {
                                size_t opp_vert_index = triangles[tri0_index].vertex_indices[m];

                                // Get the angle opposite of the edge.
                                vertex_3 a = vertices[i] - vertices[opp_vert_index];
                                vertex_3 b = vertices[neighbour_j] - vertices[opp_vert_index];
                                a.normalize();
                                b.normalize();

                                float dot = a.dot(b);

                                if(-1 > dot)
                                    dot = -1;
                                else if(1 < dot)
                                    dot = 1;

                                float angle = acosf(dot);

                                // Curvature normal weighting.
                                float slope = tanf(angle);

                                if(0 == slope)
                                    slope = numeric_limits<float>::epsilon();

                                // Note: Some weights will be negative, due to obtuse triangles.
                                // You may wish to do weights[j] += fabsf(1.0f / slope); here.
                                weights[j] += 1.0f / slope;

                                angle_count++;

                                break;
                            }
                        }

                        // Since we found a triangle match, we can skip to the first vertex's next triangle.
                        break;
                    }
                }
            } // End of: Find out which two triangles are shared by the vertex pair.

            if(angle_count != 2)
                angle_error++;

        } // End of: For each vertex pair (ie. each edge).

        if(angle_error != 0)
        {
            cout << "Warning: Vertex " << i << " belongs to " << angle_error << " edges that do not belong to two triangles (" << vertex_to_vertex_indices[i].size() - angle_error << " edges were OK)." << endl;
            cout << "Your mesh probably has cracks or holes in it." << endl;
        }

        // Normalize the weights so that they sum up to 1.
        float s = 0;

        // Note: Some weights will be negative, due to obtuse triangles.
        // You may wish to do s += fabsf(weights[j]); here.
        for(size_t j = 0; j < weights.size(); j++)
            s += weights[j];

        if(0 == s)
            s = numeric_limits<float>::epsilon();

        for(size_t j = 0; j < weights.size(); j++)
            weights[j] /= s;

        // Sum the displacements.
        for(size_t j = 0; j < vertex_to_vertex_indices[i].size(); j++)
        {
            size_t neighbour_j = vertex_to_vertex_indices[i][j];

            displacements[i] += (vertices[neighbour_j] - vertices[i])*weights[j];
        }
    }

    // To do: Find out why there are cases where displacement is much, much, much larger than all edge lengths put together.

    // Apply per-vertex displacement.
    for(size_t i = 0; i < vertices.size(); i++)
        vertices[i] += displacements[i]*scale;
}

Edited by taby, 24 April 2013 - 09:36 AM.


#7 taby   Members   -  Reputation: 336

Like
1Likes
Like

Posted 29 April 2013 - 11:12 AM

Another simple smoothing algorithm to try is where the scale is multiplied by a factor related to each vertex's angle deficit (curvature), ie.:

 

vertex[i] += displacement[i]*scale*((2pi - total_angle[i])/(2pi)).

 

This should smooth out spikes and pits, but leave ridges, valleys and flat regions relatively untouched -- which seems ideal for the task at hand.

 

Not sure what the algorithm's called, nor how much it would differ from the curvature normal weighting, but will try later today.

 

Update:  implemented as suggested, it doesn't work entirely as expected, but will try other things related to it.


Edited by taby, 01 May 2013 - 03:53 PM.





Old topic!
Guest, the last post of this topic is over 60 days old and at this point you may not reply in this topic. If you wish to continue this conversation start a new topic.



PARTNERS