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 -- [attachment=15206:mc_smooth2.zip]

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) -- [attachment=15174:fractal.zip]

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

[attachment=15175:1_raw_flat_shaded.png]

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

[attachment=15177: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:

[attachment=15178: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:

[attachment=15179: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:

[attachment=15180: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:

[attachment=15181: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

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

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

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.size())
            continue;

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

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

            float edge_length = vertices.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.size(); j++)
        {
            size_t neighbour_j = vertex_to_vertex_indices[j];
            displacements += (vertices[neighbour_j] - vertices)*weights[j];
        }
    }

    // Apply per-vertex displacement.
    for(size_t i = 0; i < vertices.size(); i++)
        vertices += displacements*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.size())
            continue;

        vector<float> weights(vertex_to_vertex_indices.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.size(); j++)
        {
            size_t angle_count = 0;

            size_t neighbour_j = vertex_to_vertex_indices[j];

            // Find out which two triangles are shared by the edge.
            for(size_t k = 0; k < vertex_to_triangle_indices.size(); k++)
            {
                for(size_t l = 0; l < vertex_to_triangle_indices[neighbour_j].size(); l++)
                {
                    size_t tri0_index = vertex_to_triangle_indices[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 - 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.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.size(); j++)
        {
            size_t neighbour_j = vertex_to_vertex_indices[j];

            displacements += (vertices[neighbour_j] - vertices)*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 += displacements*scale;
}

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 += displacement*scale*((2pi - total_angle)/(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.