[vtkusers] DICOM to Model Mesh

[Marco Nawijn]

Would it be acceptable to step out of VTK for meshing?
You could try to do something like the following:
1. Import DICOM data
2. Select region of interest
3. Export region of interest to STL
4. Use dedicated mesh generator like netgen

http://sourceforge.net/apps/mediawiki/netgen-mesher/index.php?title=Main_Page
    to mesh the geometry
5. Import mesh into VTK

Export/Import can ofcourse also be replaced by
direct calls from C++ into the netgen library.
I am not sure if direct VTK support is available in
netgen.


Regards,

Marco



On Wed, Nov 27, 2013 at 5:20 AM, Harold <boogiedoll at yahoo.com> wrote:

> Hi Miro,
>
> Thank you for the guide. I'm new to vtk so while messing around, I managed
> to apply the followings to a set of dicom images (in a very crude way) and
> save it as a .vtk file.
> 1. Marching Cubes - get mesh
> 2. Decimation - reduce the number of triangles
> 3. Smoothing - smooth the mesh
> 4. Cropping - select the largest region using connectivity filter
>
> I guess my next step is to find out how to do pre-processing and refining
> of the defected mesh, and also making statistical models (if there are
> enough training images...).
>
> By the way, here's my crude code based on the the  MarchingCubes Example
> at http://www.vtk.org/Wiki/VTK/Examples/Cxx/Modelling/MarchingCubes".
> Hopefully it can be a small help to other novice like me out there.
>
> /////////////////////////////////
> int main(int argc, char *argv[])
> {
>   vtkSmartPointer<vtkImageData> volume =
>     vtkSmartPointer<vtkImageData>::New();
>   double isoValue;
>
>   if (argc < 3)
>     {
> // Create a sphere if there's no dicom input
>     vtkSmartPointer<vtkSphereSource> sphereSource =
>       vtkSmartPointer<vtkSphereSource>::New();
>     sphereSource->SetPhiResolution(20);
>     sphereSource->SetThetaResolution(20);
>     sphereSource->Update();
>
>     double bounds[6];
>     sphereSource->GetOutput()->GetBounds(bounds);
>     for (unsigned int i = 0; i < 6; i += 2)
>       {
>       double range = bounds[i+1] - bounds[i];
>       bounds[i]   = bounds[i] - .1 * range;
>       bounds[i+1] = bounds[i+1] + .1 * range;
>       }
>     vtkSmartPointer<vtkVoxelModeller> voxelModeller =
>       vtkSmartPointer<vtkVoxelModeller>::New();
>     voxelModeller->SetSampleDimensions(50,50,50);
>     voxelModeller->SetModelBounds(bounds);
>     voxelModeller->SetScalarTypeToFloat();
>     voxelModeller->SetMaximumDistance(.1);
>
>     voxelModeller->SetInputConnection(sphereSource->GetOutputPort());
>     voxelModeller->Update();
>     isoValue = 0.5;
>     volume->DeepCopy(voxelModeller->GetOutput());
>     }
>   else
>     {
>     // Read dicom image when dicom image folder path is specified
>     vtkSmartPointer<vtkDICOMImageReader> reader =
>       vtkSmartPointer<vtkDICOMImageReader>::New();
>     reader->SetDirectoryName(argv[1]);
>     reader->Update();
>     volume->DeepCopy(reader->GetOutput());
>     isoValue = atof(argv[2]);
>     }
>
>   // Use MarchingCubes to generate iso-surface
>   vtkSmartPointer<vtkMarchingCubes> surface =
>     vtkSmartPointer<vtkMarchingCubes>::New();
>
>   #if VTK_MAJOR_VERSION <= 5
>     surface->SetInput(volume);
>   #else
>     surface->SetInputData(volume);
>   #endif
>   surface->ComputeNormalsOn();
>   surface->ComputeScalarsOn();
>   surface->SetValue(0, isoValue);
>   std::cout << "Marching Cube finished..." << std::endl;
>
>   // Create polydata from iso-surface
>   vtkSmartPointer<vtkPolyData> marched =
>     vtkSmartPointer<vtkPolyData>::New();
>   surface->Update();
>   marched->DeepCopy(surface->GetOutput());
>   std::cout<<"Number of points: " << marched->GetNumberOfPoints() <<
> std::endl;
>
>   // Decimation to reduce the number of triangles
>   vtkSmartPointer<vtkDecimatePro> decimator =
>     vtkDecimatePro::New();
>   decimator->SetInputData(marched);
>   decimator->SetTargetReduction(0.5);
>   decimator->SetPreserveTopology(1);
>   decimator->Update();
>   std::cout << "Decimation finished...." << std::endl;
>   // Smoothing
>   vtkSmartPointer<vtkSmoothPolyDataFilter> smoother =
> vtkSmartPointer<vtkSmoothPolyDataFilter>::New();
>   smoother->SetInputData(decimator->GetOutput());
>   smoother->SetNumberOfIterations(5);
>   smoother->SetFeatureAngle(60);
>   smoother->SetRelaxationFactor(0.05);
>   smoother->FeatureEdgeSmoothingOff();
>   std::cout << "Smoothing mesh finished...." << std::endl;
>
>  // Select the largest region
>   vtkSmartPointer<vtkPolyDataConnectivityFilter> connectivityFilter =
> vtkSmartPointer<vtkPolyDataConnectivityFilter>::New();
>   connectivityFilter->SetInputConnection(smoother->GetOutputPort());
>   connectivityFilter->ScalarConnectivityOff();
>   connectivityFilter->SetExtractionModeToLargestRegion();
>   connectivityFilter->Update();
>
>  // Create final polygon mesh
>   vtkSmartPointer<vtkPolyData> mesh =
>     vtkSmartPointer<vtkPolyData>::New();
>   mesh->ShallowCopy(connectivityFilter->GetOutput());
>   std::cout<<"Number of points in the final polygon: " <<
> mesh->GetNumberOfPoints() << std::endl;
>
>   // Visualization
>   vtkSmartPointer<vtkRenderer> renderer =
>     vtkSmartPointer<vtkRenderer>::New();
>   renderer->SetBackground(.1, .2, .3);
>
>   vtkSmartPointer<vtkRenderWindow> renderWindow =
>     vtkSmartPointer<vtkRenderWindow>::New();
>   renderWindow->AddRenderer(renderer);
>   vtkSmartPointer<vtkRenderWindowInteractor> interactor =
>     vtkSmartPointer<vtkRenderWindowInteractor>::New();
>   interactor->SetRenderWindow(renderWindow);
>
>   vtkSmartPointer<vtkPolyDataMapper> mapper =
>     vtkSmartPointer<vtkPolyDataMapper>::New();
>   mapper->SetInputData(mesh);
>   mapper->ScalarVisibilityOff();
>
>   vtkSmartPointer<vtkActor> actor =
>     vtkSmartPointer<vtkActor>::New();
>   actor->SetMapper(mapper);
>
>   renderer->AddActor(actor);
>
>   renderWindow->Render();
>   interactor->Start();
>
>   // Save the mesh as a .vtk file
>   vtkSmartPointer<vtkPolyDataWriter> writer = vtkPolyDataWriter::New();
>   writer->SetFileName("mesh.vtk");
>   writer->SetInputData(mesh);
>   writer->SetFileTypeToASCII();
>   writer->Write();
>
>   return EXIT_SUCCESS;
> }
>
>
> Best Regards,
> Harold
>
>
>   On Wednesday, November 27, 2013 12:06 PM, Miro Drahos <
> mdrahos at robodoc.com> wrote:
>  This is basically a segmentation problem, which in general is a
> difficult problem.
> One such pipeline could be as follows:
> 1.) image preprocessing (cropping, smoothing, etc).
> 2.) thresholding to get a filled binary mask of your object (vertebra)
> 3.) running Marching Cubes to get a mesh.
>
> This is a very crude and basic approach, but something to get you
> started. Note that marching cubes might give you a mesh with topological
> defects (holes, handles), from the noise in your data.
>
> There are a variety of other sophisticated algorithms, e.g. using
> statistical methods if you have some already segmented data (like
> creating a mean shape from the previously segmented data, and using PCA
> to do shape analysis, or some machine learning approach).
>
> Another approach would be to use deformable models, i.e. define forces
> that would (iteratively) drive an initial mesh to fit the image. The
> forces would be based on image intensity, gradient and/or laplacian, and
> internal forces to keep the mesh regular and relatively smooth.
>
> There is a lot that has been done in this domain, so look at some
> articles on segmentation.
>
> HTH,
> Miro
>
>
> On 11/25/2013 06:20 AM, Boogiedoll wrote:
> > Hi,
> > I would like to ask the experts about common way in vtk to create a
> spine bone (e.g L1) mesh from a CT DICOM stack, including the
> pre-processing.
> >
> > Thanks & Regards,
> > Harold
> > _______________________________________________
> > Powered by www.kitware.com
> >
> > Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
> >
> > Please keep messages on-topic and check the VTK FAQ at:
> http://www.vtk.org/Wiki/VTK_FAQ
> >
> > Follow this link to subscribe/unsubscribe:
> > http://www.vtk.org/mailman/listinfo/vtkusers
>
>
>
>
> _______________________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Please keep messages on-topic and check the VTK FAQ at:
> http://www.vtk.org/Wiki/VTK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://www.vtk.org/mailman/listinfo/vtkusers
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.vtk.org/pipermail/vtkusers/attachments/20131127/751f1ddc/attachment.htm>