[Insight-developers] Change suggestions for itkVTKPolyDataWriter and Simplex Meshes

[Thomas Boettger]

Dear Alex,

I modified the code of the itkVTKPolyDataWriter class slightly for two reasons:

1. I did not have you CVS patch for the point id squeezing. So I
introduced a runnning index at take the empty points into account
during save time. This work for


2. I wanted to use the writer for storing simplex meshes and found it
was working for triangle meshes only as you assumed all cells being
triangles when writing the total number of polygon ids:

  outputFile << "POLYGONS " << numberOfCells << " " << 4 *
numberOfCells << std::endl;

I changed it. I count the number of cell points for all cells
(numOfEntries) and store this number:

  outputFile << "POLYGONS " << numberOfCells << " " << numberOfCells +
numOfEntries << std::endl;

Find the complete GenerateData method below. You might want to
consider adding the changes. I think this fix is better than calling
the SqueezePointID method as the mesh is intact without squeezed
PointIds and one should be able to store any polydata mesh with the
polydata writer.

2.Simplex Meshes:
I furthermore adapted my deformation classes to run one a TriangleMesh
without explicit conversion. The deformation filter will take care of
creating all simplex mesh relationships necessary for the deformation
algorithm internally and update the original triangle mesh in
ComputeOutput. I did this as a proof of concept and it works. I will
send you the code sometime next week.

I am still not sure what would be the best way of porting the simplex
mesh deformation algo.

In general I do believe it is better to have a stric separation, i.e.
convert to simplex mesh -> deform -> store simplex mesh or -> convert
to triangle mesh -> store triangle mesh.

On the other hand I showed with my latest implementation tests one can
completely ignore the simplex meshes and only use the deformation
algorithm hiding the dual simplex mesh representation from the user
completely.

Furthermore those simplex meshes look so nice I think one should not
try hide them from a user in clincal applications visualization parts
:-)

Well there is loads of thing which need to be discussed. Maybe we
should start putting together a document discussing the different
options and pros & cons before a decision can be made.


  Kind regards,

    Thomas


template<class TInputMesh>
void
VTKPolyDataWriter<TInputMesh>
::GenerateData()
{
  std::cout << __LINE__ << " GenerateData" << std::endl;

  this->m_Input->SetCellsAllocationMethod(
      InputMeshType::CellsAllocatedDynamicallyCellByCell );

  if( this->m_FileName == "" )
    {
    itkExceptionMacro("No FileName");
    return;
    }

  //
  // Write to output file
  //
  std::ofstream outputFile( this->m_FileName.c_str() );

  if( !outputFile.is_open() )
    {
    itkExceptionMacro("Unable to open file\n"
        "outputFilename= " << this->m_FileName );
    return;
    }

  outputFile << "# vtk DataFile Version 2.0" << std::endl;
  outputFile << "File written by itkVTKPolyDataWriter" << std::endl;
  outputFile << "ASCII" << std::endl;
  outputFile << "DATASET POLYDATA" << std::endl;

  unsigned int numberOfPoints = this->m_Input->GetNumberOfPoints();
  outputFile << "POINTS " << numberOfPoints << " float" << std::endl;


  std::map<unsigned long,unsigned long> pointIDMap;

  unsigned long runningIndex = 0;

  PointIterator pointIterator = this->m_Input->GetPoints()->Begin();
  PointIterator pointEnd = this->m_Input->GetPoints()->End();

  while( pointIterator != pointEnd )
  {
    PointType point = pointIterator.Value();
    pointIDMap.insert(std::make_pair(pointIterator.Index(), runningIndex) );

    outputFile << point[0] << " " << point[1] << " " << point[2] << std::endl;

    runningIndex++;
    pointIterator++;
  }


  // here we are taking the linecells out of the count
  unsigned int numberOfCells = 0;
  CellIterator cellIterator = this->m_Input->GetCells()->Begin();
  CellIterator cellEnd = this->m_Input->GetCells()->End();

  unsigned int numOfEntries = 0;
  while( cellIterator != cellEnd )
    {
    if( !(cellIterator.Value()->GetType() == 1 ) ) // LINE_CELL
      {
      numOfEntries += cellIterator.Value()->GetNumberOfPoints();
      numberOfCells++;
      }
    cellIterator++;
    }
  outputFile << "POLYGONS " << numberOfCells << " " << numberOfCells +
numOfEntries << std::endl;

  // here we should do a multipass algorithms
  // one pass per type in a vtkPolyData
  // here i rule out the line cells,
  // but triangle strips and others shoudl also be ruled out
  //  or handled nicely
  cellIterator = this->m_Input->GetCells()->Begin();
  while( cellIterator != cellEnd )
    {
    CellType * cellPointer = cellIterator.Value();
    if( !(cellIterator.Value()->GetType() == 1) ) // LINE_CELL
      {
      PointIdIterator pointIdIterator = cellPointer->PointIdsBegin();
      PointIdIterator pointIdEnd = cellPointer->PointIdsEnd();

      outputFile << cellPointer->GetNumberOfPoints();

      while( pointIdIterator != pointIdEnd )
      {
        unsigned long nextId = (*pointIDMap.find(*pointIdIterator)).second;
        outputFile << " " << nextId;
        pointIdIterator++;
      }

      outputFile << std::endl;
      }
    cellIterator++;
    }

  outputFile.close();
}