Hi,
I am using VTK for medical applications, and encounter a possible bug on vtkPolyConnectivityFilter. I use the function to generate an STL from only the largest connected region. I doubt the spacing parameter of VTK data causes this phenomenon. Your help is greatly appreciated.
Following code defines a cylinder on 3D space, and generates two STL files from the cylinder. First STL with ordinary axis pitches has only a part of a cylinder, and the second STL with axis pitch of small deviations, difference between float value and double value, has full cylinder.
The program uses VTK-8.1.1 and runs on Ubuntu 18.04.02 LTS.
#include “vtkAutoInit.h”
VTK_MODULE_INIT(vtkRenderingOpenGL2); // VTK was built with vtkRenderingOpenGL2
VTK_MODULE_INIT(vtkInteractionStyle);
// vtk headers
#include <vtkImageImport.h>
#include <vtkMarchingCubes.h>
#include <vtkPolyDataConnectivityFilter.h>
#include <vtkProperty.h>
#include <vtkSTLReader.h>
#include <vtkSTLWriter.h>
// local window display
#include <vtkActor.h>
#include <vtkPolyDataMapper.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include
#include <unistd.h>
using namespace std;
//----------------------------------------------------//
// show stl file //
//----------------------------------------------------//
void showStl(string fnstl) {
// define renderer
vtkRenderer *renderer = vtkRenderer::New();
renderer->SetBackground(0.5, 0.5, 0.5); // Background color
vtkRenderWindow *renderWindow = vtkRenderWindow::New();
renderWindow->AddRenderer(renderer);
vtkRenderWindowInteractor *renderWindowInteractor = vtkRenderWindowInteractor::New();
renderWindowInteractor->SetRenderWindow(renderWindow);
// show stl
vtkSTLReader *stlreader = vtkSTLReader::New();
stlreader->SetFileName(fnstl.c_str());
stlreader->Update();
// setup the visualization pipeline
vtkPolyDataMapper *stlmapper = vtkPolyDataMapper::New();
stlmapper->SetInputConnection(stlreader->GetOutputPort());
// actor
vtkActor *stlactor = vtkActor::New();
stlactor->GetProperty()->SetColor(0.5,0.5,0.5);
stlactor->GetProperty()->SetOpacity(0.5);
stlactor->SetMapper(stlmapper);
// set stl to rendere
renderer->AddActor(stlactor);
stlreader->Delete();
stlactor->Delete();
renderWindow->Render();
renderWindowInteractor->Start();
// delete vtk objects
renderer->Delete();
renderWindowInteractor->Delete();
renderWindow->Delete();
stlreader->Delete();
stlmapper->Delete();
stlactor->Delete();
}
//----------------------------------------------------//
// make stl file from voxel data //
//----------------------------------------------------//
void makeStl(float *fpvox, string fnstl, int pix[3], double spc[3]) {
// Convert the c-style image to a vtkImageData
vtkImageImport *imp = vtkImageImport::New();
imp->SetDataSpacing(spc[0], spc[1], spc[2]);
imp->SetDataOrigin(0, 0, 0);
imp->SetWholeExtent(0, pix[0]-1, 0, pix[1]-1, 0, pix[2]-1);
imp->SetDataExtentToWholeExtent();
imp->SetDataScalarTypeToFloat();
imp->SetNumberOfScalarComponents(1);
imp->SetImportVoidPointer(fpvox);
imp->Update();
// apply marching cubes
vtkMarchingCubes *sfc = vtkMarchingCubes::New();
sfc->SetInputConnection(imp->GetOutputPort());
sfc->GenerateValues(1, 1, 1);
sfc->Update();
// to remain largest region
vtkPolyDataConnectivityFilter *flt = vtkPolyDataConnectivityFilter::New();
flt->SetInputConnection(sfc->GetOutputPort());
flt->SetExtractionModeToLargestRegion();
// generate stl file (binary)
vtkSTLWriter *stw = vtkSTLWriter::New();
stw->SetInputConnection(flt->GetOutputPort());
stw->SetFileName(fnstl.c_str());
stw->SetFileTypeToBinary();
stw->Write();
// delete vtk classes
imp->Delete();
sfc->Delete();
flt->Delete();
stw->Delete();
}
//----------------------------------------------------//
// entry point //
//----------------------------------------------------//
int main(int argc, char *argv[]) {
// define sizes
int pix[3] = {400, 400, 150};
double spc[3] = {0.626953, 0.626953, 3};
// allocate memory
int n = pix[2] * pix[1] * pix[0]; // sizes on axes
float *fpvox = (float *)malloc(n * sizeof(float)); // pitches on axes
// make test data: cylinder
int cntr[2] = {pix[0] / 2, pix[1] / 2}; // cylinder center on x-y plane
int rds = 100; // radius
for (int k=0; k<pix[2]; k++) {
for (int j=0; j<pix[1]; j++) {
for (int i=0; i<pix[0]; i++) {
int inx = (k * pix[1] + j) * pix[0] + i;
float r = sqrt((i - cntr[0]) * (i - cntr[0]) + (j - cntr[1]) * (j - cntr[1]));
if (r <= rds) {
fpvox[inx] = 1;
} else {
fpvox[inx] = 0;
}
}
}
}
// make stl: erroneous
string fnstl0 = "stl0.stl";
makeStl(fpvox, fnstl0, pix, spc);
// show stl
if (fork() == 0) {
showStl(fnstl0);
}
// make small error on axis pitch
spc[0] = (float)spc[0];
spc[1] = (float)spc[1];
spc[2] = (float)spc[2];
// make stl: correct
string fnstl1 = "stl1.stl";
makeStl(fpvox, fnstl1, pix, spc);
// release memory
free(fpvox);
// show stl
showStl(fnstl1);
return EXIT_SUCCESS;
}