[text] Meshing code

Viewer

copydownloadembedprintName: Meshing code
  1. vtkSmartPointer<vtkPolyData> CreateMesh(vtkSmartPointer<vtkImageData> imageData, MeshingParameters parameters)
  2. {
  3.    // The padding is now adjusted automatically 
  4.    double spacingImage[3];
  5.    imageData->GetSpacing(spacingImage);
  6.  
  7.    // Compute the image padding
  8.    double paddingImage[3] = { parameters.gaussianSmoothingActualRadiusX  / spacingImage[0],
  9.       parameters.gaussianSmoothingActualRadiusY  / spacingImage[1],
  10.       parameters.gaussianSmoothingActualRadiusZ  / spacingImage[2] };
  11.  
  12.    // pad the input volume because if the segmentation touches the edge of the dataset,
  13.    // the generated mesh will contain holes
  14.    vtkSmartPointer<vtkImageConstantPad> padImageFilter = vtkImageConstantPad::New();
  15.    padImageFilter->SetInputData(imageData);
  16.    padImageFilter->SetConstant(0);
  17.    int extent[6];
  18.    imageData->GetExtent(extent);
  19.  
  20.    // Add one voxel for the volume boundary in order to ensure an empty space
  21.    // pad the x dimension
  22.    extent[0] -= (paddingImage[0] + 1); 
  23.    extent[1] += (paddingImage[0] + 1);
  24.    // pad the y dimension
  25.    extent[2] -= (paddingImage[1] + 1);
  26.    extent[3] += (paddingImage[1] + 1);
  27.    // pad the z dimension
  28.    extent[4] -= (paddingImage[2] + 1); 
  29.    extent[5] += (paddingImage[2] + 1); 
  30.    padImageFilter->SetOutputWholeExtent(extent);
  31.    padImageFilter->Update();
  32.    
  33.    // Convert to a floating point volume. 
  34.    // Note from proulxc: This is unfortunate, memory wise, but the VTK gaussian filter
  35.    // on integer types performs weirdly. It appears to apply the mask in three successive
  36.    // passes, rounding each time and causing non-symmetries and numerical errors in the
  37.    // process. 
  38.    const double MAX_VALUE = 1.0;
  39.    vtkSmartPointer<vtkImageCast> castImage = 
  40.       vtkSmartPointer<vtkImageCast>::New();
  41.    castImage->SetInputConnection(padImageFilter->GetOutputPort());
  42.    castImage->SetOutputScalarTypeToFloat();
  43.    castImage->Update();
  44.    
  45.    vtkSmartPointer<vtkImageGaussianSmooth> gaussianSmoothingFilter = 
  46.       vtkSmartPointer<vtkImageGaussianSmooth>::New();
  47.    gaussianSmoothingFilter->SetInputConnection(castImage->GetOutputPort());
  48.  
  49.    gaussianSmoothingFilter->SetDimensionality(3);
  50.  
  51.    double sigmaX = parameters.gaussianSmoothingFilterStdX;
  52.    double sigmaY = parameters.gaussianSmoothingFilterStdY;
  53.    double sigmaZ = parameters.gaussianSmoothingFilterStdZ;
  54.  
  55.    // Note that radius factor is multiplied by standard deviation to obtain
  56.    // the actual radius. Eg: sigma = 2 and radius factor = 2 -> actual radius = 4
  57.    // i.e. 9x9x9 kernel. We compute the radius factor so that the actual radius
  58.    // will always be the one set in the Algo Parameters config file, regardless of sigma.
  59.    double radiusFactorX = 1.0;
  60.    double radiusFactorY = 1.0;
  61.    double radiusFactorZ = 1.0;
  62.    double actualRadiusX = parameters.gaussianSmoothingActualRadiusX;
  63.    double actualRadiusY = parameters.gaussianSmoothingActualRadiusY;
  64.    double actualRadiusZ = parameters.gaussianSmoothingActualRadiusZ;
  65.    if ( sigmaX >= std::numeric_limits<double>::epsilon() )
  66.    {
  67.       radiusFactorX = actualRadiusX/sigmaX;
  68.    }
  69.    if ( sigmaY >= std::numeric_limits<double>::epsilon() )
  70.    {
  71.       radiusFactorY = actualRadiusY/sigmaY;
  72.    }
  73.    if ( sigmaZ >= std::numeric_limits<double>::epsilon() )
  74.    {
  75.       radiusFactorZ = actualRadiusZ/sigmaZ;
  76.    }
  77.  
  78.    gaussianSmoothingFilter->SetRadiusFactors(radiusFactorX,radiusFactorY,radiusFactorZ);
  79.    gaussianSmoothingFilter->SetStandardDeviations(sigmaX,sigmaY,sigmaZ);
  80.    gaussianSmoothingFilter->Update();
  81.  
  82.  
  83.    vtkSmartPointer<vtkMarchingCubes> surface = vtkSmartPointer<vtkMarchingCubes>::New();
  84.    surface->SetInputConnection(gaussianSmoothingFilter->GetOutputPort());
  85.    surface->ComputeNormalsOn();
  86.  
  87.    surface->SetValue(0,parameters.isoSurfaceThreshold);
  88.    surface->Update();
  89.  
  90.    vtkAlgorithmOutput* marchingCubesOutput = surface->GetOutputPort();
  91.  
  92.    // Verify if the output is topologically correct using the Euler-Poincare formula 
  93.    // for a manifold mesh topologically equivalent to a sphere. 
  94.    // V - E + F = 2 where V is the number of vertices, E is the 
  95.    // number of edges (i.e. nbTriangles*3/2 since every edge is shared by two triangles)
  96.    // and F is the number of faces.
  97.    // This is necessary because of a poorly documented bug in VTK whereas some volumes
  98.    // are split into distinct isosurfaces (due to memory constraints?) and incorrectly
  99.    // stitched back together.  
  100.    // See http://public.kitware.com/pipermail/vtkusers/2009-July/053115.html
  101.    int nbVertices = surface->GetOutput()->GetNumberOfPoints();
  102.    int nbTriangles = surface->GetOutput()->GetNumberOfPolys();
  103.    vtkSmartPointer<vtkCleanPolyData> cleanSurface = vtkSmartPointer<vtkCleanPolyData>::New();
  104.    if (nbVertices - nbTriangles/2 != 2)
  105.    {
  106.       ORTHO_LOG_INFO("Meshing: Output of Marching Cubes is "
  107.          "topologically incorrect. Attempting to clean the mesh. ");
  108.       cleanSurface->SetInputConnection(surface->GetOutputPort());
  109.       cleanSurface->Update();
  110.       nbVertices = cleanSurface->GetOutput()->GetNumberOfPoints();
  111.       nbTriangles = cleanSurface->GetOutput()->GetNumberOfPolys();
  112.       if (nbVertices - nbTriangles/2 != 2)
  113.       {
  114.          ORTHO_LOG_ERROR("Meshing: Output of Marching Cubes could not "
  115.             "be cleaned into a topologically correct mesh.")
  116.       }
  117.       marchingCubesOutput = cleanSurface->GetOutputPort();
  118.    }
  119.  
  120.    double ratioPoly = static_cast<double>(parameters.targetTriangleCount)/nbTriangles;
  121.    double targetReduction = 1.0-ratioPoly;
  122.  
  123.    vtkSmartPointer<vtkQuadricDecimation> decimatedSurface = 
  124.       vtkSmartPointer<vtkQuadricDecimation>::New();
  125.    decimatedSurface->SetInputConnection(marchingCubesOutput);
  126.    decimatedSurface->SetTargetReduction(targetReduction);
  127.    decimatedSurface->Update();
  128.    vtkSmartPointer<vtkPolyData> resultingMesh = decimatedSurface->GetOutput();
  129.  
  130.    return resultingMesh;
  131. }

Editor

You can edit this paste and save as new:


File Description
  • Meshing code
  • Paste Code
  • 30 Jan-2023
  • 5.85 Kb
You can Share it: