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

**URL:** https://discourse.vtk.org/t/how-to-get-a-reasonably-good-triangulated-surface-from-a-level-set-volumetric-field/765
**Category:** Support
**Created:** [April 16, 2019, 4:17pm UTC](https://discourse.vtk.org/t/how-to-get-a-reasonably-good-triangulated-surface-from-a-level-set-volumetric-field/765 "2019-04-16T16:17:51Z")
**Posts on this page:** 3
**Page:** 1

<div class="post-metadata">

### Author: ![rdbisme](https://discourse.vtk.org/user_avatar/discourse.vtk.org/rdbisme/32/309_2.png) [@rdbisme](https://discourse.vtk.org/u/rdbisme)
#### Post date: [April 16, 2019, 4:17pm UTC](https://discourse.vtk.org/t/how-to-get-a-reasonably-good-triangulated-surface-from-a-level-set-volumetric-field/765/1 "2019-04-16T16:17:51Z")

</div>

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:

```auto
/** 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 =
        vtkSmartPointer<vtkImageData>::New();
    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,
                        SphereDistance(
                            coords[0] - sphereCenter[0],
                            coords[1] - sphereCenter[1],
                            coords[2] - sphereCenter[2],
                            sphereRadius)
                        );

            }
        }
    }

    vtkSmartPointer<vtkXMLImageDataWriter> iwriter =
        vtkSmartPointer<vtkXMLImageDataWriter>::New();
    iwriter->SetInputData(volume);
    iwriter->SetFileName("sphere.vti");
    iwriter->Write();

    // Now we have the distance in the domain, let's reconstruct the surface
    // using the marching cube
    vtkSmartPointer<vtkMarchingCubes> triangulizer =
        vtkSmartPointer<vtkMarchingCubes>::New();
    triangulizer->SetInputData(volume);
    triangulizer->Update();

    // The surface obtained with the Marching Cubes algorithm contains
    // duplicate points, so let's clean it
    vtkSmartPointer<vtkCleanPolyData> cleaner =
        vtkSmartPointer<vtkCleanPolyData>::New();
    cleaner->SetInputData(triangulizer->GetOutput());
    cleaner->SetAbsoluteTolerance(1E-3);
    cleaner->Update();

    Surface surface(cleaner->GetOutput());
    surface.ComputeCurvatures();
    surface.Write("sphere.vtp");

} // 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.

---

<div class="post-metadata">

### Author: ![lorensen](https://discourse.vtk.org/user_avatar/discourse.vtk.org/lorensen/32/20_2.png) [@lorensen](https://discourse.vtk.org/u/lorensen)
#### Post date: [April 16, 2019, 9:25pm UTC](https://discourse.vtk.org/t/how-to-get-a-reasonably-good-triangulated-surface-from-a-level-set-volumetric-field/765/2 "2019-04-16T21:25:08Z")

</div>

Can you share an image of the poorly extracted isosurface?

---

<div class="post-metadata">

### Author: ![rdbisme](https://discourse.vtk.org/user_avatar/discourse.vtk.org/rdbisme/32/309_2.png) [@rdbisme](https://discourse.vtk.org/u/rdbisme)
#### Post date: [April 19, 2019, 11:59am UTC](https://discourse.vtk.org/t/how-to-get-a-reasonably-good-triangulated-surface-from-a-level-set-volumetric-field/765/3 "2019-04-19T11:59:28Z")

</div>

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.

 ![sphere-levelset](https://discourse.vtk.org/uploads/default/original/1X/f35b491afbdeedf71ccfdcb69a3704f76c733956.jpeg)

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