Problems about 3D construction with tetrahedral mesh from DICOM series

classic Classic list List threaded Threaded
1 message Options
Reply | Threaded
Open this post in threaded view
|

Problems about 3D construction with tetrahedral mesh from DICOM series

smallping
I want to read the DICOM series, and then create tetrahedral mesh for volume modeling.

I do like this:
1) Read the 3D image data from a series of DICOM files
2) Use CurvatureFlowImageFilter to smooth the image
3) Run an image segmentation method
4) Use the vtkContour filter to generate the surface mesh
5) Save the surface mesh data to a .vtk file
6) Using TetGen, to generate the tetrahedral mesh from the surface mesh, and save into .vtk file

But I have problems at step 6, when I try to use TetGen to generate the tetrahedral mesh from the surface mesh, it failed.
Then I use ParaView to visualize my surface mesh, I find maybe there are some problems with the surface mesh.
Because in the surface mesh, there are many holes, and a surface circle around the object, I do not know what is this.
The surface mesh is shown as the following picture.



Can somebody know how to fix this ?
Or give me an example to do this ?
I need help from you.
Thanks

The following is my codes for step 2 to step 6 to generate the surface mesh
/*=========================================================================
#include "itkCommand.h"
#include "itkImage.h"
#include "itkVTKImageExport.h"
#include "itkVTKImageImport.h"
#include "itkConfidenceConnectedImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkRGBPixel.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"

#include "vtkImageImport.h"
#include "vtkImageExport.h"
#include "vtkRenderer.h"
#include "vtkRenderWindow.h"
#include "vtkRenderWindowInteractor.h"
#include "vtkActor.h"
#include "vtkPolyData.h"
#include "vtkPolyDataMapper.h"
#include "vtkContourFilter.h"
#include "vtkImageData.h"
#include "vtkDataSet.h"
#include "vtkProperty.h"
#include "vtkImagePlaneWidget.h"
#include "vtkCellPicker.h"
#include "vtkPolyDataWriter.h"

#include "itkCurvatureFlowImageFilter.h"
#include "vtkSmartPointer.h"
#include "vtkSTLWriter.h"


template <typename ITK_Exporter, typename VTK_Importer>
void ConnectPipelines(ITK_Exporter exporter, VTK_Importer* importer)
{
  importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback());
  importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback());
  importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback());
  importer->SetSpacingCallback(exporter->GetSpacingCallback());
  importer->SetOriginCallback(exporter->GetOriginCallback());
  importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback());
  importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback());
  importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback());
  importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback());
  importer->SetDataExtentCallback(exporter->GetDataExtentCallback());
  importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback());
  importer->SetCallbackUserData(exporter->GetCallbackUserData());
}

template <typename VTK_Exporter, typename ITK_Importer>
void ConnectPipelines(VTK_Exporter* exporter, ITK_Importer importer)
{
  importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback());
  importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback());
  importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback());
  importer->SetSpacingCallback(exporter->GetSpacingCallback());
  importer->SetOriginCallback(exporter->GetOriginCallback());
  importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback());
  importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback());
  importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback());
  importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback());
  importer->SetDataExtentCallback(exporter->GetDataExtentCallback());
  importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback());
  importer->SetCallbackUserData(exporter->GetCallbackUserData());
}


int main(int argc, char * argv [] )
{
  //DICOM-IMG.vtk is read from DICOM series
  argv[1]="/home/smallping/Codes/ITK/Test/MeshSurface/DICOM-IMG.vtk";
 
  try
  {
    typedef float        InputPixelType;
    typedef unsigned char MaskPixelType;

    const unsigned int Dimension = 3;
    typedef itk::Image< InputPixelType, Dimension > InputImageType;
    typedef itk::Image< MaskPixelType,  Dimension > MaskImageType ;

    typedef itk::ImageFileReader< InputImageType >  ReaderType    ;

    ReaderType::Pointer reader  = ReaderType::New();
    reader->SetFileName( argv[1] );
    reader->Update();

    ////////////////////////////////////////////////////////////// Smooth
    typedef itk::CurvatureFlowImageFilter< InputImageType, InputImageType > CurvatureFlowImageFilterType;
    CurvatureFlowImageFilterType::Pointer smoothing = CurvatureFlowImageFilterType::New();

    smoothing->SetInput( reader->GetOutput() );
    smoothing->SetNumberOfIterations( 10 );
    smoothing->SetTimeStep( 0.0625 );

    //////////////////////////////////////////////////////////////Segmentation
    typedef itk::ConfidenceConnectedImageFilter<InputImageType, MaskImageType> ConnectedFilterType;
    ConnectedFilterType::Pointer filter = ConnectedFilterType::New();

    filter->SetInput( smoothing->GetOutput() );

    filter->SetNumberOfIterations(2);
    filter->SetReplaceValue(255);
    filter->SetMultiplier(2.5);
     
    InputImageType::Pointer inputImage = reader->GetOutput();
    InputImageType::SizeType  size  = inputImage->GetBufferedRegion().GetSize();
    InputImageType::IndexType start = inputImage->GetBufferedRegion().GetIndex();

    InputImageType::IndexType seed;
    seed[0] = start[0] + size[0] / 2;
    seed[1] = start[1] + size[1] / 2;
    seed[2] = start[2] + size[2] / 2;

    filter->SetSeed( seed );

    //////////////////////////////////////////////////////////////
    typedef itk::VTKImageExport< InputImageType > ExportFilter1Type;
    typedef itk::VTKImageExport< MaskImageType  > ExportFilter2Type;

    ExportFilter1Type::Pointer itkExporter1 = ExportFilter1Type::New();
    ExportFilter2Type::Pointer itkExporter2 = ExportFilter2Type::New();

    itkExporter1->SetInput( reader->GetOutput() );
    itkExporter2->SetInput( filter->GetOutput() );

    vtkImageImport* vtkImporter1 = vtkImageImport::New();  
    ConnectPipelines(itkExporter1, vtkImporter1);
    vtkImageImport* vtkImporter2 = vtkImageImport::New();  
    ConnectPipelines(itkExporter2, vtkImporter2);
   
    vtkImporter1->Update();
     
    //////////////////////////////////////////////////////////////Surface
    vtkContourFilter * contour = vtkContourFilter::New();
    contour->SetInput( vtkImporter2->GetOutput() );
    contour->SetValue(0, 128);

    //////////////////////////////////////////////////////////////
    vtkPolyDataWriter * vtkWriter = vtkPolyDataWriter::New();
    vtkWriter->SetFileName("/home/smallping/Codes/ITK/Test/MeshSurface/DICOM-SURF.vtk"); //save as vtk
    vtkWriter->SetInput( contour->GetOutput() );
    vtkWriter->Write();

    vtkSmartPointer< vtkSTLWriter > stlWriter = vtkSmartPointer< vtkSTLWriter >::New();  //save as stl
    stlWriter->SetInput( contour->GetOutput() );
    stlWriter->SetFileName( "/home/smallping/Codes/ITK/Test/MeshSurface/DICOM-SURF.stl" );
    stlWriter->Update();
 
    //////////////////////////////////////////////////////////////
    vtkImporter1->Delete();
    vtkImporter2->Delete();
    contour->Delete();
  }
  catch( itk::ExceptionObject & e )
  {
    std::cerr << "Exception catched !! " << e << std::endl;
  }

  return 0;
}