[vtkusers] DICOM to Model Mesh

[Harold]

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.vtk.org/pipermail/vtkusers/attachments/20131126/21ded969/attachment.htm>