How to get a reasonably good triangulated surface from a level-set volumetric field

Hello, I’m exploring some algorithms to compute curvatures on surfaces and I’m using VTK under the hoods to manage the triangulation of a level-set surface (and data structures).

The problem I’m encountering is a relatively poor quality of the triangulated surfaces once the distance function fields get triangulated (I’m using vtkMarchingCubes but I also tried the vtkExtractSurface and vtkFlyingEdges3D filters too).

Here a MWE to create a spherical surface from its level-set:

/** This example is a MWE to be posted on the VTK discourse to understand how
 * to use the vtkCleanPolyDataFilter to avoid having duplicate points when
 * triangulating a sphere with a Marching Cube on a given level-set field

#include <vtkCleanPolyData.h>
#include <vtkImageData.h>
#include <vtkMarchingCubes.h>
#include <vtkXMLImageDataWriter.h>
#include <vtkXMLPolyDataWriter.h>

#include "Surface.h"

using namespace hgve;

/** @brief The distance function for a sphere */
static double SphereDistance(double x, double y, double z, double R) {
    return x*x + y*y + z*z - R*R;

int main(int argc, char *argv[])

    // Let's define a uniform grid on which to compute our distance function
    vtkSmartPointer<vtkImageData> volume =
    volume->SetDimensions(100, 100, 100);
    volume->AllocateScalars(VTK_DOUBLE, 1);
    volume->SetSpacing(0.05, 0.05, 0.05);

    // Let's compute the value of the distance function for the sphere in all
    // the domain (i.e. for each x, y, z)
    static const double sphereRadius = 2;
    static const double sphereCenter[3] = {2.5, 2.5, 2.5};

    auto dims = volume->GetDimensions();
    for (int i = 0; i < dims[0]; i++) {
        for (int j = 0; j < dims[1]; j++) {
            for (int k = 0; k < dims[2]; k++) {
                // Let's merge coordinate index in an array
                int ijk[3] = {i, j, k};

                // Retrieve the pointId from its ijk indices
                auto pId = volume->ComputePointId(ijk);

                // Retrieve x, y, z coordinates
                auto coords = volume->GetPoint(pId);

                // Set the value of the distance at (x, y, z)
                volume->SetScalarComponentFromDouble(i, j, k, 0,
                            coords[0] - sphereCenter[0],
                            coords[1] - sphereCenter[1],
                            coords[2] - sphereCenter[2],


    vtkSmartPointer<vtkXMLImageDataWriter> iwriter =

    // Now we have the distance in the domain, let's reconstruct the surface
    // using the marching cube
    vtkSmartPointer<vtkMarchingCubes> triangulizer =

    // The surface obtained with the Marching Cubes algorithm contains
    // duplicate points, so let's clean it
    vtkSmartPointer<vtkCleanPolyData> cleaner =

    Surface surface(cleaner->GetOutput());

} // main

Is this the best I can do or not? Because in the resulting surface I get very skewed triangles. Any suggestion would be very appreciated.

Can you share an image of the poorly extracted isosurface?

Hi Bill, thank you for taking the time to reply. I didn’t get the notification so please excuse me for the late reply.

Here’s an image of the sphere with the computation of the Gauss Curvature (should be constant == 0.25 for a sphere of radius 2) made with the algorithm I’m currently evaluating [1]. As you can see I have some skewed triangles that I think belong to the intersection of the original cartesian uniform mesh cells with the sphere itself, but I don’t know in details how the vtkMarchingCubes filter works and if it does something specifically different from the reference algorithm in the literature.

[1] M. Meyer, M. Desbrun, P. Schröder, and A. H. Barr, “Discrete Differential-Geometry Operators for Triangulated 2-Manifolds,” in Visualization and Mathematics III , H.-C. Hege and K. Polthier, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg, 2003, pp. 35–57.)