itkmeshtest.cxx

DTMK软件开发包,此为开源软件,是一款很好的医学图像开发资源.

   * boundaries.  (The result should be 1.  Note that the quads and
   * triangles that bound the 3D cells haven't been added, so they
   * don't show up as additional neighbors.)
   */
  int numberOfNeighbors = mesh->GetCellBoundaryFeatureNeighbors(
    1,              // Topological dimension of feature.
    1,              // CellIdentifier
    1,              // CellFeatureIdentifier
    NULL);          // We don't want the neighbors themselves (yet)
  std::cout << "Number of neighbors (hex edge 1): "
            << numberOfNeighbors << ". " << std::endl;

  /**
   * Try getting the hexahedron's neighbor through its second edge.
   * This should be the test tetrahedron. (Because the boundary is
   * not defined explicitly, we use implicit relationships to determine
   * neighbors. In this case, the tetrahedron and hexahedron share
   * the two points defining the edge and are therefore considered 
   * neighbors.)
   */
  mesh->GetCellBoundaryFeatureNeighbors(
    1,              // Topological dimension of feature.
    1,              // CellIdentifier
    1,              // CellFeatureIdentifier
    &neighborSet);  // Where to put result.

  std::cout << "Neighbors (hex edge 1):" << std::endl;
  for(cell = neighborSet.begin(); cell != neighborSet.end(); ++cell)
    {
    std::cout << "Id " << *cell << ": ";
    CellAutoPointer cellPointer;

    if( mesh->GetCell(*cell, cellPointer) )
      {
      std::cout << cellPointer->GetNameOfClass();
      }
    std::cout << std::endl;
    }
 

  /**
   * Try getting the tetrahedrons's neighbor through its fourth edge.
   * This should be the test hexahedron. The boundaries are implicit
   * as in the previous example.
   */
  mesh->GetCellBoundaryFeatureNeighbors(
    1,              // Topological dimension of feature.
    0,              // CellIdentifier
    3,              // CellFeatureIdentifier
    &neighborSet); // Where to put result.

  std::cout << "Neighbors (tet edge 3):" << std::endl;
  for(cell = neighborSet.begin(); cell != neighborSet.end(); ++cell)
    {
    std::cout << "Id " << *cell << ": ";
    CellAutoPointer cellPointer;

    if( mesh->GetCell(*cell, cellPointer) )
      {
      std::cout << cellPointer->GetNameOfClass();
      }
    std::cout << std::endl;
    }

  /**
   * Create a higher order  test cell.
   */
  { // Create a local scope
    testCell.TakeOwnership( new QuadraticEdgeCellType ); // polymorphism
    unsigned long quadraticEdgePoints[3] = {0,1,2};
    testCell->SetPointIds(quadraticEdgePoints);
    mesh->SetCell(2, testCell ); // Internally transfers ownership to the mesh

    CellAutoPointer quadraticdEdgeCell;
    if(  mesh->GetCell(2, quadraticdEdgeCell) )
      {
      std::cout << "Quadratic Edge cell recovered" << std::endl;
      std::cout << "GetNameOfClass()    = " << quadraticdEdgeCell->GetNameOfClass() << std::endl;
      std::cout << "GetNumberOfPoints() = " << quadraticdEdgeCell->GetNumberOfPoints() << std::endl;
      }
    else
      {
      std::cout << "Failure: QuadraticEdge cell was not recovered" << std::endl;
      return EXIT_FAILURE;
      }

    CellAutoPointer vertexPointer;
    /**
     * Try getting one of the QuadraticEdge's vertices.
     */
    const bool vertexExists = mesh->GetCellBoundaryFeature(
                  0,    // Topological dimension of boundary.
                  2,    // CellIdentifier.
                  0,    // CellFeatureIdentifier
                  vertexPointer ); // CellPointer to return the result

    std::cout << typeid( vertexPointer ).name() << std::endl;
    std::cout << typeid( vertexPointer.GetPointer() ).name() << std::endl;
    std::cout << "GetCellBoundaryFeature() return AutoPointer owner = " << vertexPointer.IsOwner() << std::endl;
    
    QuadraticEdgeCellType::VertexType * vertex;
    try
      {
      vertex = dynamic_cast<QuadraticEdgeCellType::VertexType *>(  vertexPointer.GetPointer() );
      std::cout << "Vertex from the QuadraticEdge recovered " << std::endl;
      std::cout << vertex->GetNameOfClass() << std::endl;
      }
    catch(...)
      {
      std::cout << "CellPointer cannot be down-cast to a VertexCellType" << std::endl;
      vertex = 0;
      }
    if( vertex )
      {
      std::cout << "CellPointer was safely down-casted to a VertexCellType" << std::endl;
      }
    if( vertexExists ) 
      {
      std::cout << "Vertex number of points = " << vertexPointer->GetNumberOfPoints() << std::endl;
      std::cout << "Vertex name of class    = " << vertexPointer->GetNameOfClass() << std::endl;
      }
    else 
      {
      std::cout << "Vertex of the QuadraticEdge couldn't be extracted " << std::endl;
      }
  } // end of local scope for this part of the test. 


  /**
   * Create a higher order triangular test cell.
   */
  { // Create a local scope

    // In this block, we overwrite a cell at a particular id.
    // To avoid a memory leak, we must first grab ownership of the
    // current cell at that id so that the memory for the original
    // cell will be deleted when we leave this scope
    CellAutoPointer cellToDelete;
    mesh->GetCell(2, cellToDelete);
    cellToDelete.TakeOwnership();
    
    // Now we can construct a new cell and overwrite the id
    testCell.TakeOwnership(new QuadraticTriangleCellType); // polymorphism;
    unsigned long quadraticTrianglePoints[3] = {0,1,2};
    testCell->SetPointIds(quadraticTrianglePoints);
    mesh->SetCell(2, testCell ); // Internally transfers ownership to the mesh
    std::cout << "QuadraticTriangleCell pointer = " << (void*)testCell.GetPointer() << std::endl;
    std::cout << "QuadraticTriangleCell Owner   = " << testCell.IsOwner() << std::endl;

    CellAutoPointer quadraticdTriangleCell;

    if(  mesh->GetCell(2, quadraticdTriangleCell) )
      {
      std::cout << "Quadratic Triangle cell recovered" << std::endl;
      std::cout << "GetNameOfClass()    = " << quadraticdTriangleCell->GetNameOfClass() << std::endl;
      std::cout << "GetNumberOfPoints() = " << quadraticdTriangleCell->GetNumberOfPoints() << std::endl;
      }
    else
      {
      std::cout << "Failure: QuadraticTriangle cell was not recovered" << std::endl;
      return EXIT_FAILURE;
      }


    CellAutoPointer vertexPointer;

    /**
     * Try getting one of the QuadraticTriangle's vertices.
     */
    const bool vertexExists = mesh->GetCellBoundaryFeature(
                  0,    // Topological dimension of boundary.
                  2,    // CellIdentifier.
                  0,    // CellFeatureIdentifier
                  vertexPointer ); // CellPointer to return the result

    std::cout << typeid( vertexPointer ).name() << std::endl;
    std::cout << typeid( vertexPointer.GetPointer() ).name() << std::endl;
    std::cout << "GetCellBoundaryFeature() return AutoPointer owner = " << vertexPointer.IsOwner() << std::endl;

    QuadraticTriangleCellType::VertexType * vertex;
    try
      {
      vertex = dynamic_cast<QuadraticTriangleCellType::VertexType *>(  vertexPointer.GetPointer() );
      std::cout << "Vertex from the QuadraticTriangle recovered " << std::endl;
      std::cout << vertex->GetNameOfClass() << std::endl;
      }
    catch(...)
      {
      std::cout << "CellPointer cannot be down-cast to a VertexCellType" << std::endl;
      vertex = 0;
      }
    if( vertex )
      {
      std::cout << "CellPointer was safely down-casted to a VertexCellType" << std::endl;
      }
    if( vertexExists ) 
      {
      std::cout << "Vertex number of points = " << vertexPointer->GetNumberOfPoints() << std::endl;
      std::cout << "Vertex name of class    = " << vertexPointer->GetNameOfClass() << std::endl;
      }
    else 
      {
      std::cout << "Vertex of the QuadraticTriangle couldn't be extracted " << std::endl;
      }

    /**
     * Try getting one of the QuadraticTriangle's edges.
     */
    CellAutoPointer edgePointer;
    const bool edgeExists = mesh->GetCellBoundaryFeature(
                  1,    // Topological dimension of boundary.
                  2,    // CellIdentifier.
                  0,    // CellFeatureIdentifier
                  edgePointer ); // CellPointer to return the result

    std::cout << typeid( edgePointer ).name() << std::endl;
    std::cout << typeid( edgePointer.GetPointer() ).name() << std::endl;
    std::cout << "GetCellBoundaryFeature() return AutoPointer owner = " << edgePointer.IsOwner() << std::endl;
    
    QuadraticTriangleCellType::EdgeType * edge;
    try
      {
      edge = dynamic_cast<QuadraticTriangleCellType::EdgeType *>(  edgePointer.GetPointer() );
      std::cout << "Vertex from the QuadraticTriangle recovered " << std::endl;
      std::cout << edge->GetNameOfClass() << std::endl;
      }
    catch(...)
      {
      std::cout << "CellPointer cannot be down-cast to a VertexCellType" << std::endl;
      edge = 0;
      }
    if( edge )
      {
      std::cout << "CellPointer was safely down-casted to a VertexCellType" << std::endl;
      }
    if( edgeExists ) 
      {
      std::cout << "Edge number of points = " << edgePointer->GetNumberOfPoints() << std::endl;
      std::cout << "Edge name of class    = " << edgePointer->GetNameOfClass() << std::endl;

      // Evaluate The Shape functions for a particular parametric point
      CellType::ParametricCoordArrayType parametricCoordinates(1);
      CellType::ShapeFunctionsArrayType  weights( edgePointer->GetNumberOfPoints() );
      parametricCoordinates[0] = 0.25;
      edgePointer->EvaluateShapeFunctions( parametricCoordinates, weights );
      std::cout << "Shape Function weights = " << weights << std::endl;
      }
    else 
      {
      std::cerr << "Edge of the QuadraticTriangle couldn't be extracted " << std::endl;
      }
  } // end of local scope for this part of the test.



  /**
   * Perform some geometric operations (coordinate transformations)
   * to see if they are working.
   */
  MeshType::CoordRepType coords[MeshType::PointDimension];
  MeshType::PointIdentifier pointId;
  mesh->FindClosestPoint(coords,&pointId);

  /**
   * Compute the bounding box of the mesh
   */
  typedef itk::BoundingBox<MeshType::PointIdentifier,MeshType::PointDimension,
    MeshType::CoordRepType,MeshType::PointsContainer> BoundingBox;

  BoundingBox::Pointer bbox(BoundingBox::New());
  bbox->SetPoints(mesh->GetPoints());
  bbox->ComputeBoundingBox(); 

  /**
   * Set up some visitors 
   */
  MeshType::CellType::MultiVisitor::Pointer mv =
    MeshType::CellType::MultiVisitor::New();
  /**
   * Create a class to hold the counts of each cell type
   */
  CountClass counts;
  /**
   * Create a visitor for each cell type and set the counts class for the visitor
   */
  TetraCellVisitor::Pointer cv = TetraCellVisitor::New();
  cv->SetCountClass(&counts);
  QuadraticEdgeCellVisitor::Pointer ev = QuadraticEdgeCellVisitor::New();
  ev->SetCountClass(&counts);
  QuadraticTriangleCellVisitor::Pointer tv = QuadraticTriangleCellVisitor::New();
  tv->SetCountClass(&counts);
  mv->AddVisitor(cv);
  mv->AddVisitor(ev);
  mv->AddVisitor(tv);

  // 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);
  // print the counts found
  std::cout << "Number of TetraCellType " << counts.m_Tetra << "\n";
  std::cout << "Number of QuadraticEdgeCellType " << counts.m_QuadraticEdgeCell << "\n";
  std::cout << "Number of QuadraticTriangleCellType " << counts.m_QuadraticTriangleCellType << "\n";

  std::cout << bbox << std::endl;

  std::cout << mesh << std::endl;


  // Exercising the Graft method
  MeshType::Pointer newMesh = MeshType::New();
  newMesh->Graft( mesh );

  if( newMesh->GetNumberOfPoints() != mesh->GetNumberOfPoints() )
    {
    std::cerr << "Graft failed !, different number of points" << std::endl;
    return EXIT_FAILURE;  
    }

  if( newMesh->GetNumberOfCells() != mesh->GetNumberOfCells() )
    {
    std::cerr << "Graft failed !, different number of cells" << std::endl;
    return EXIT_FAILURE;  
    }


  std::cout << newMesh << std::endl;

  return EXIT_SUCCESS;  
}