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