- vtkSmartPointer<vtkPolyData> CreateMesh(vtkSmartPointer<vtkImageData> imageData, MeshingParameters parameters)
- {
- // The padding is now adjusted automatically
- double spacingImage[3];
- imageData->GetSpacing(spacingImage);
- // Compute the image padding
- double paddingImage[3] = { parameters.gaussianSmoothingActualRadiusX / spacingImage[0],
- parameters.gaussianSmoothingActualRadiusY / spacingImage[1],
- parameters.gaussianSmoothingActualRadiusZ / spacingImage[2] };
- // pad the input volume because if the segmentation touches the edge of the dataset,
- // the generated mesh will contain holes
- vtkSmartPointer<vtkImageConstantPad> padImageFilter = vtkImageConstantPad::New();
- padImageFilter->SetInputData(imageData);
- padImageFilter->SetConstant(0);
- int extent[6];
- imageData->GetExtent(extent);
- // Add one voxel for the volume boundary in order to ensure an empty space
- // pad the x dimension
- extent[0] -= (paddingImage[0] + 1);
- extent[1] += (paddingImage[0] + 1);
- // pad the y dimension
- extent[2] -= (paddingImage[1] + 1);
- extent[3] += (paddingImage[1] + 1);
- // pad the z dimension
- extent[4] -= (paddingImage[2] + 1);
- extent[5] += (paddingImage[2] + 1);
- padImageFilter->SetOutputWholeExtent(extent);
- padImageFilter->Update();
- // Convert to a floating point volume.
- // Note from proulxc: This is unfortunate, memory wise, but the VTK gaussian filter
- // on integer types performs weirdly. It appears to apply the mask in three successive
- // passes, rounding each time and causing non-symmetries and numerical errors in the
- // process.
- const double MAX_VALUE = 1.0;
- vtkSmartPointer<vtkImageCast> castImage =
- vtkSmartPointer<vtkImageCast>::New();
- castImage->SetInputConnection(padImageFilter->GetOutputPort());
- castImage->SetOutputScalarTypeToFloat();
- castImage->Update();
- vtkSmartPointer<vtkImageGaussianSmooth> gaussianSmoothingFilter =
- vtkSmartPointer<vtkImageGaussianSmooth>::New();
- gaussianSmoothingFilter->SetInputConnection(castImage->GetOutputPort());
- gaussianSmoothingFilter->SetDimensionality(3);
- double sigmaX = parameters.gaussianSmoothingFilterStdX;
- double sigmaY = parameters.gaussianSmoothingFilterStdY;
- double sigmaZ = parameters.gaussianSmoothingFilterStdZ;
- // Note that radius factor is multiplied by standard deviation to obtain
- // the actual radius. Eg: sigma = 2 and radius factor = 2 -> actual radius = 4
- // i.e. 9x9x9 kernel. We compute the radius factor so that the actual radius
- // will always be the one set in the Algo Parameters config file, regardless of sigma.
- double radiusFactorX = 1.0;
- double radiusFactorY = 1.0;
- double radiusFactorZ = 1.0;
- double actualRadiusX = parameters.gaussianSmoothingActualRadiusX;
- double actualRadiusY = parameters.gaussianSmoothingActualRadiusY;
- double actualRadiusZ = parameters.gaussianSmoothingActualRadiusZ;
- if ( sigmaX >= std::numeric_limits<double>::epsilon() )
- {
- radiusFactorX = actualRadiusX/sigmaX;
- }
- if ( sigmaY >= std::numeric_limits<double>::epsilon() )
- {
- radiusFactorY = actualRadiusY/sigmaY;
- }
- if ( sigmaZ >= std::numeric_limits<double>::epsilon() )
- {
- radiusFactorZ = actualRadiusZ/sigmaZ;
- }
- gaussianSmoothingFilter->SetRadiusFactors(radiusFactorX,radiusFactorY,radiusFactorZ);
- gaussianSmoothingFilter->SetStandardDeviations(sigmaX,sigmaY,sigmaZ);
- gaussianSmoothingFilter->Update();
- vtkSmartPointer<vtkMarchingCubes> surface = vtkSmartPointer<vtkMarchingCubes>::New();
- surface->SetInputConnection(gaussianSmoothingFilter->GetOutputPort());
- surface->ComputeNormalsOn();
- surface->SetValue(0,parameters.isoSurfaceThreshold);
- surface->Update();
- vtkAlgorithmOutput* marchingCubesOutput = surface->GetOutputPort();
- // Verify if the output is topologically correct using the Euler-Poincare formula
- // for a manifold mesh topologically equivalent to a sphere.
- // V - E + F = 2 where V is the number of vertices, E is the
- // number of edges (i.e. nbTriangles*3/2 since every edge is shared by two triangles)
- // and F is the number of faces.
- // This is necessary because of a poorly documented bug in VTK whereas some volumes
- // are split into distinct isosurfaces (due to memory constraints?) and incorrectly
- // stitched back together.
- // See http://public.kitware.com/pipermail/vtkusers/2009-July/053115.html
- int nbVertices = surface->GetOutput()->GetNumberOfPoints();
- int nbTriangles = surface->GetOutput()->GetNumberOfPolys();
- vtkSmartPointer<vtkCleanPolyData> cleanSurface = vtkSmartPointer<vtkCleanPolyData>::New();
- if (nbVertices - nbTriangles/2 != 2)
- {
- ORTHO_LOG_INFO("Meshing: Output of Marching Cubes is "
- "topologically incorrect. Attempting to clean the mesh. ");
- cleanSurface->SetInputConnection(surface->GetOutputPort());
- cleanSurface->Update();
- nbVertices = cleanSurface->GetOutput()->GetNumberOfPoints();
- nbTriangles = cleanSurface->GetOutput()->GetNumberOfPolys();
- if (nbVertices - nbTriangles/2 != 2)
- {
- ORTHO_LOG_ERROR("Meshing: Output of Marching Cubes could not "
- "be cleaned into a topologically correct mesh.")
- }
- marchingCubesOutput = cleanSurface->GetOutputPort();
- }
- double ratioPoly = static_cast<double>(parameters.targetTriangleCount)/nbTriangles;
- double targetReduction = 1.0-ratioPoly;
- vtkSmartPointer<vtkQuadricDecimation> decimatedSurface =
- vtkSmartPointer<vtkQuadricDecimation>::New();
- decimatedSurface->SetInputConnection(marchingCubesOutput);
- decimatedSurface->SetTargetReduction(targetReduction);
- decimatedSurface->Update();
- vtkSmartPointer<vtkPolyData> resultingMesh = decimatedSurface->GetOutput();
- return resultingMesh;
- }
[text] Meshing code
Viewer
*** This page was generated with the meta tag "noindex, nofollow". This happened because you selected this option before saving or the system detected it as spam. This means that this page will never get into the search engines and the search bot will not crawl it. There is nothing to worry about, you can still share it with anyone.
Editor
You can edit this paste and save as new: