[Insight-users] Writing itk::Mesh to disk

Thomas Boettger t.boettger at dkfz-heidelberg.de
Mon, 22 Mar 2004 09:40:13 +0100


This is a multi-part message in MIME format.
--------------020205060401050005010504
Content-Type: text/plain; charset=us-ascii; format=flowed
Content-Transfer-Encoding: 7bit

Hi Yasser,

I would convert itk to vtkPolyData or vtkUnstructuredGrid and save it to 
disk.

You find and example on how this can be done in

	InsightApplications-1.6.0/Auxiliary/vtk/vtk2itk.cxx


Furthermore I send you some code, which is based on the example, but 
where you can additionally create a vtkPolyData by calling the static method

   static vtkPolyData* meshToPolyData(typename 
MeshType::Pointer mesh, 							bool onlyTriangles = false);



Regards,

Thomas



Yasser Bashir wrote:
> Hi,
>  
> Is there a convenient way to perform file I/O for itk::Mesh? 
>  
> All of the IO filters seem to be catering to images.
>  
> Regards,
>  
> Yasser
> 
> 
> Do you Yahoo!?
> Yahoo! Finance Tax Center - File online. File on time.


-- 
Dipl.-Inform. Thomas Boettger
Deutsches Krebsforschungszentrum         (German Cancer Research Center)
Div. Medical and Biological Informatics B010    Tel: (+49) 6221-42 2328
Im Neuenheimer Feld 280                          Fax: (+49) 6221-42 2345
D-69120 Heidelberg                            e-mail: t.boettger at dkfz.de
Germany                      http://www.dkfz.de/mbi/people/thomasb.shtml


--------------020205060401050005010504
Content-Type: text/plain;
 name="MeshUtilNew.h"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="MeshUtilNew.h"

#include <itkMesh.h>
#include <itkLineCell.h>
#include <itkTriangleCell.h>
#include <itkPolygonCell.h>
#include <itkQuadrilateralCell.h>
#include <itkCellInterface.h>
#include <itkDefaultDynamicMeshTraits.h> 
#include <itkSphereMeshSource.h> 
#include <itkTransformMeshFilter.h>
#include <itkTranslationTransform.h>
#include <itkMinimumMaximumImageCalculator.h>

#include <vtkActor.h>
#include <vtkCellArray.h>
#include <vtkPolyData.h>
#include <vtkUnstructuredGrid.h>
#include <vtkPoints.h>
#include <vtkPointData.h>
#include <vtkProperty.h>
#include <vtkFloatArray.h>

//#include "itkRegularSphereMeshSource.h"

/*!
\brief The class provides mehtods for ITK - VTK mesh conversion
*
*  \todo document the inner class 
*  \todo maybe inner class should be moved out
*/

template <typename MeshType>
class MeshUtil
{

  /*!
  \brief implementation for VTK cell visitors
  \todo documentation
  */
  class VistVTKCellsClass
  {
    vtkCellArray* m_Cells;
    int* m_LastCell;
    int* m_TypeArray;
  public:
    // typedef the itk cells we are interested in
    typedef typename itk::CellInterface< typename MeshType::PixelType, 
                                         typename MeshType::CellTraits >  CellInterfaceType;
    typedef itk::LineCell<CellInterfaceType>      floatLineCell;
    typedef itk::PolygonCell<CellInterfaceType>      floatPolygonCell;
    typedef itk::TriangleCell<CellInterfaceType>      floatTriangleCell;
    typedef itk::QuadrilateralCell<CellInterfaceType> floatQuadrilateralCell;

    /*! Set the vtkCellArray that will be constructed
    */
    void SetCellArray(vtkCellArray* a) 
    {
      m_Cells = a;
    }

    /*! 
    Set the cell counter pointer
    */
    void SetCellCounter(int* i)
    {
      m_LastCell = i;
    }

    /*!
    Set the type array for storing the vtk cell types
    */
    void SetTypeArray(int* i)
    {
      m_TypeArray = i;
    }

    /*!
    Visit a line and create the VTK_LINE cell   
    */
    void Visit(unsigned long , floatLineCell* t)
    {
      m_Cells->InsertNextCell(2,  (vtkIdType*)t->PointIdsBegin());
      m_TypeArray[*m_LastCell] = VTK_LINE;
      (*m_LastCell)++;
    }

    /*!
    Visit a line and create the VTK_POLYGON cell   
    */
    void Visit(unsigned long , floatPolygonCell* t)
    {
 
      unsigned long num = t->GetNumberOfVertices();
      if (num > 3) {
        m_Cells->InsertNextCell(num, (vtkIdType*)t->PointIdsBegin());
        m_TypeArray[*m_LastCell] = VTK_POLYGON;
        (*m_LastCell)++;
      }
    }

    /*!
    Visit a triangle and create the VTK_TRIANGLE cell   
    */
    void Visit(unsigned long , floatTriangleCell* t)
    {
      m_Cells->InsertNextCell(3,  (vtkIdType*)t->PointIdsBegin());
      m_TypeArray[*m_LastCell] = VTK_TRIANGLE;
      (*m_LastCell)++;
    }

    /*! 
    Visit a triangle and create the VTK_QUAD cell 
    */
    void Visit(unsigned long , floatQuadrilateralCell* t)
    {
      m_Cells->InsertNextCell(4,  (vtkIdType*)t->PointIdsBegin());
      m_TypeArray[*m_LastCell] = VTK_QUAD;
      (*m_LastCell)++;
    }
  };


  typedef typename MeshType::CellType                CellType;
  typedef typename itk::LineCell< CellType >         LineType;
  typedef typename itk::PolygonCell< CellType >         PolygonType;
  typedef typename itk::TriangleCell< CellType >     TriangleType;


  /*! 
  default line cell visitior definition 
  */
  typedef typename itk::CellInterfaceVisitorImplementation<
                                      typename MeshType::PixelType,
                                      typename MeshType::CellTraits,
                                      LineType, 
                                      VistVTKCellsClass>			LineVisitor;

  /*! 
  default line cell visitior definition 
  */
  typedef typename itk::CellInterfaceVisitorImplementation<
                                      typename MeshType::PixelType,
                                      typename MeshType::CellTraits,
                                      PolygonType, 
                                      VistVTKCellsClass>			PolygonVisitor;

  /*! 
  default triangle cell visitior definition 
  */
  typedef typename itk::CellInterfaceVisitorImplementation<typename MeshType::PixelType,
    typename MeshType::CellTraits,
    itk::TriangleCell<itk::CellInterface<typename MeshType::PixelType, typename MeshType::CellTraits >
    >, VistVTKCellsClass>			TriangleVisitor;


  /*! 
  default quad cell visitior definition 
  */
  typedef typename itk::CellInterfaceVisitorImplementation<typename MeshType::PixelType, typename MeshType::CellTraits,
    itk::QuadrilateralCell< itk::CellInterface<typename MeshType::PixelType, typename MeshType::CellTraits > >, 
    VistVTKCellsClass> QuadrilateralVisitor;


public:


  /*!
  create an itkMesh object from a vtkPolyData
  */
  static typename MeshType::Pointer meshFromPolyData(vtkPolyData* poly)
  {
    // Create a new mesh
    typename MeshType::Pointer output = MeshType::New();
    output->SetCellsAllocationMethod( MeshType::CellsAllocatedDynamicallyCellByCell );  

    typedef typename MeshType::CellDataContainer MeshCellDataContainerType;

    output->SetCellData(MeshCellDataContainerType::New());

    // Get the points from vtk
    vtkPoints* vtkpoints = poly->GetPoints();
    int numPoints = poly->GetNumberOfPoints();

    // Create a compatible point container for the mesh
    // the mesh is created with a null points container
    // MeshType::PointsContainer::Pointer points = 
    //   MeshType::PointsContainer::New();
    // // Resize the point container to be able to fit the vtk points
    // points->Reserve(numPoints);
    // // Set the point container on the mesh
    //output->SetPoints(points);
    for(int i =0; i < numPoints; i++)
    {
      typename MeshType::PixelType* apoint = (typename MeshType::PixelType*) vtkpoints->GetPoint(i);
      typename MeshType::PointType p;
      p[0] = apoint[0];
      p[1] = apoint[1];
      p[2] = apoint[2];
      output->SetPoint( i, p );
    }

    vtkCellArray* vtkcells = poly->GetPolys();
    // extract the cell id's from the vtkUnstructuredGrid
    int numcells = vtkcells->GetNumberOfCells();
    int* vtkCellTypes = new int[numcells];
    int cellId =0;
    for(; cellId < numcells; cellId++)
    {
      vtkCellTypes[cellId] = poly->GetCellType(cellId);
    }

    // cells->Reserve(numcells);
    vtkIdType npts;
    vtkIdType* pts;
    cellId = 0;

    typedef typename MeshType::MeshTraits   OMeshTraits;
    typedef typename OMeshTraits::PixelType       OPixelType; 
    typedef typename MeshType::CellTraits					CellTraits;
    typedef typename itk::CellInterface<OPixelType, CellTraits>   CellInterfaceType;
    typedef typename itk::TriangleCell<CellInterfaceType>         TriCellType;
    typedef typename TriCellType::CellAutoPointer	           TriCellPointer;

    TriCellPointer newCell;

    for(vtkcells->InitTraversal(); vtkcells->GetNextCell(npts, pts); cellId++)
    {
      switch(vtkCellTypes[cellId])
      {
      case VTK_TRIANGLE:
        {
          newCell.TakeOwnership( new TriCellType );
          newCell->SetPointIds((unsigned long*)pts);
          output->SetCell(cellId, newCell );
          output->SetCellData(cellId, (typename MeshType::PixelType)3);
          break;    
        }  

      case VTK_QUAD:
      case VTK_EMPTY_CELL:
      case VTK_VERTEX:
      case VTK_POLY_VERTEX:
      case VTK_LINE:
      case VTK_POLY_LINE:
      case VTK_TRIANGLE_STRIP:
      case VTK_POLYGON:
      case VTK_PIXEL:
      case VTK_TETRA:
      case VTK_VOXEL:
      case VTK_HEXAHEDRON:
      case VTK_WEDGE:
      case VTK_PYRAMID:
      case VTK_PARAMETRIC_CURVE:
      case VTK_PARAMETRIC_SURFACE:
      default:
        std::cerr << "Warning unhandled cell type " 
          << vtkCellTypes[cellId] << std::endl;
      }
    }
    //output->Print(std::cout);
    return output;
  }


  /*!
    create an vtkUnstructuredGrid object from an itkMesh
  */
  static vtkUnstructuredGrid* meshToUnstructuredGrid(typename MeshType::Pointer mesh)
  {


    // Get the number of points in the mesh
    int numPoints = mesh->GetNumberOfPoints();
    if(numPoints == 0)
    {
      mesh->Print(std::cerr);
      std::cerr << "no points in Grid " << std::endl;
      exit(-1);
    }
    // Create a vtkUnstructuredGrid
    vtkUnstructuredGrid* vgrid = vtkUnstructuredGrid::New();
    // Create the vtkPoints object and set the number of points
    vtkPoints* vpoints = vtkPoints::New();
    vpoints->SetNumberOfPoints(numPoints);
    // iterate over all the points in the itk mesh filling in
    // the vtkPoints object as we go
    typename MeshType::PointsContainer::Pointer points = mesh->GetPoints();

    for(typename MeshType::PointsContainer::Iterator i = points->Begin();
      i != points->End(); ++i)
    {
      // Get the point index from the point container iterator
      int idx = i->Index();


      // Set the vtk point at the index with the the coord array from itk
      // itk returns a const pointer, but vtk is not const correct, so
      // we have to use a const cast to get rid of the const
      //      vpoints->SetPoint(idx, const_cast<DATATYPE*>(i->Value().GetDataPointer()));
      vpoints->SetPoint(idx, (typename MeshType::PixelType*)(i->Value().GetDataPointer()));
    }
    // Set the points on the vtk grid
    vgrid->SetPoints(vpoints);
    // Now create the cells using the MulitVisitor
    // 1. Create a MultiVisitor
    typename MeshType::CellType::MultiVisitor::Pointer mv =
      typename MeshType::CellType::MultiVisitor::New();
    // 2. Create a triangle and quadrilateral visitor
    typename TriangleVisitor::Pointer tv = TriangleVisitor::New();
    typename QuadrilateralVisitor::Pointer qv =  QuadrilateralVisitor::New();
    // 3. Set up the visitors
    int vtkCellCount = 0; // running counter for current cell being inserted into vtk
    int numCells = mesh->GetNumberOfCells();
    int *types = new int[numCells]; // type array for vtk 
    // create vtk cells and estimate the size
    vtkCellArray* cells = vtkCellArray::New();
    cells->EstimateSize(numCells, 4);
    // Set the TypeArray CellCount and CellArray for both visitors
    tv->SetTypeArray(types);
    tv->SetCellCounter(&vtkCellCount);
    tv->SetCellArray(cells);
    qv->SetTypeArray(types);
    qv->SetCellCounter(&vtkCellCount);
    qv->SetCellArray(cells);
    // add the visitors to the multivisitor
    mv->AddVisitor(tv);
    mv->AddVisitor(qv);
    // Now ask the mesh to accept the multivisitor which
    // will Call Visit for each cell in the mesh that matches the
    // cell types of the visitors added to the MultiVisitor
    mesh->Accept(mv);
    // Now set the cells on the vtk grid with the type array and cell array
    vgrid->SetCells(types, cells);  
    // Clean up vtk objects (no vtkSmartPointer ... )
    cells->Delete();
    vpoints->Delete();

    //std::cout << "meshToUnstructuredGrid end" << std::endl;
    return vgrid;
  }


  /*!
    create an vtkUnstructuredGrid object from an itkMesh
  */
  static vtkPolyData* meshToPolyData(typename MeshType::Pointer mesh, bool onlyTriangles = false)
  {


    // Get the number of points in the mesh
    int numPoints = mesh->GetNumberOfPoints();
    if(numPoints == 0)
    {
      mesh->Print(std::cerr);
      std::cerr << "no points in Grid " << std::endl;
    }
    // Create a vtkUnstructuredGrid
    vtkPolyData* vgrid = vtkPolyData::New();
    // Create the vtkPoints object and set the number of points
    vtkPoints* vpoints = vtkPoints::New();
    vpoints->SetNumberOfPoints(numPoints);


    vtkFloatArray * scalars = vtkFloatArray::New();
    scalars->SetNumberOfComponents(1);
    scalars->SetNumberOfTuples(numPoints);

    // iterate over all the points in the itk mesh filling in
    // the vtkPoints object as we go
    typename MeshType::PointsContainer::Pointer points = mesh->GetPoints();

    for(typename MeshType::PointsContainer::Iterator i = points->Begin();
      i != points->End(); ++i)
    {
      // Get the point index from the point container iterator
      int idx = i->Index();
      
      typename MeshType::PointType p;
      p = i->Value();
      double  pd[3];
      pd[0]=p[0];
      pd[1]=p[1];
      pd[2]=p[2];

      typename MeshType::PixelType dis;

      mesh->GetPointData(i->Index(), &dis);

     
      // Set the vtk point at the index with the the coord array from itk
      // itk returns a const pointer, but vtk is not const correct, so
      // we have to use a const cast to get rid of the const
      //      vpoints->SetPoint(idx, const_cast<DATATYPE*>(i->Value().GetDataPointer()));
      vpoints->SetPoint(idx, pd);
      scalars->InsertTuple1( idx, dis );
    }


    // Set the points on the vtk grid
    vgrid->SetPoints(vpoints);
    vgrid->GetPointData()->SetScalars(scalars);
    vgrid->GetPointData()->CopyAllOn();


    // Now create the cells using the MulitVisitor
    // 1. Create a MultiVisitor
    typedef typename MeshType::CellType::MultiVisitor MeshMV;
    typename MeshMV::Pointer mv = MeshMV::New();
    // 2. Create a triangle and quadrilateral visitor
    typename LineVisitor::Pointer lv = LineVisitor::New();
    typename PolygonVisitor::Pointer pv = PolygonVisitor::New();
    typename TriangleVisitor::Pointer tv = TriangleVisitor::New();
    typename QuadrilateralVisitor::Pointer qv =  QuadrilateralVisitor::New();
    // 3. Set up the visitors
    int vtkCellCount = 0; // running counter for current cell being inserted into vtk
    int numCells = mesh->GetNumberOfCells();
    int *types = new int[numCells]; // type array for vtk 
    // create vtk cells and estimate the size
    vtkCellArray* cells = vtkCellArray::New();
    cells->EstimateSize(numCells, 4);
    // Set the TypeArray CellCount and CellArray for both visitors

    lv->SetTypeArray(types);
    lv->SetCellCounter(&vtkCellCount);
    lv->SetCellArray(cells);

    pv->SetTypeArray(types);
    pv->SetCellCounter(&vtkCellCount);
    pv->SetCellArray(cells);

    tv->SetTypeArray(types);
    tv->SetCellCounter(&vtkCellCount);
    tv->SetCellArray(cells);
    qv->SetTypeArray(types);
    qv->SetCellCounter(&vtkCellCount);
    qv->SetCellArray(cells);
    // add the visitors to the multivisitor

    if (onlyTriangles) {
      mv->AddVisitor(tv);
      mesh->Accept(mv);
      vgrid->SetStrips(cells);  
    }
    else 
    {
      mv->AddVisitor(tv);
      mv->AddVisitor(lv);
      mv->AddVisitor(pv);
  //    mv->AddVisitor(qv);
      // Now ask the mesh to accept the multivisitor which
      // will Call Visit for each cell in the mesh that matches the
      // cell types of the visitors added to the MultiVisitor
      mesh->Accept(mv);
      // Now set the cells on the vtk grid with the type array and cell array
      vgrid->SetPolys(cells);  
  //    vgrid->SetTriangles(cells);  
      vgrid->SetLines(cells);
    // Clean up vtk objects (no vtkSmartPointer ... )
    }
    cells->Delete();
    vpoints->Delete();

    //std::cout << "meshToUnstructuredGrid end" << std::endl;
    return vgrid;
  }


 static typename MeshType::Pointer createSphereMesh(typename MeshType::PointType center, typename MeshType::PointType scale, int* resolution)
  {
    	typedef typename itk::SphereMeshSource<MeshType> SphereSource;

      typename SphereSource::Pointer mySphereSource = SphereSource::New();

			mySphereSource->SetCenter(center);
			mySphereSource->SetScale(scale);
      mySphereSource->SetResolutionX(resolution[0]);
			mySphereSource->SetResolutionY(resolution[1]);
			mySphereSource->SetSquareness1(1);
			mySphereSource->SetSquareness2(1);
			mySphereSource->Update();
			mySphereSource->GetOutput();

      typename MeshType::Pointer resultMesh = mySphereSource->GetOutput();
      resultMesh->Register();

		  return resultMesh;
  }



  static typename MeshType::Pointer translateMesh(typename MeshType::PointType vec, typename MeshType::Pointer inputMesh)
  {
    typedef typename itk::TranslationTransform< typename MeshType::PixelType, 3>  TransformType;
    typedef typename itk::TransformMeshFilter<MeshType, MeshType, TransformType > MeshFilterType;

    typename TransformType::Pointer translation = TransformType::New();
    translation->SetIdentity();
    
    typename TransformType::ParametersType params = translation->GetParameters();
    unsigned int i;

    for (i=0; i < 3; i++) {
      params[i]= vec[i];
    }
    translation->SetParameters(params);

    typename MeshFilterType::Pointer filter = MeshFilterType::New();
    filter->SetTransform( translation );
    filter->SetInput( inputMesh );
    filter->Update();

    typename MeshType::Pointer resultMesh = filter->GetOutput();
    resultMesh->Register();
    
    return resultMesh;
  }

};


--------------020205060401050005010504--