itkmeshtest.cxx

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

/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkMeshTest.cxx,v $
  Language:  C++
  Date:      $Date: 2008-01-18 18:53:14 $
  Version:   $Revision: 1.49 $

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif

#include "itkMesh.h"
#include "itkTetrahedronCell.h"
#include "itkHexahedronCell.h"
#include "itkQuadraticEdgeCell.h"
#include "itkQuadraticTriangleCell.h"
#include "itkBoundingBox.h"
#include "itkFileOutputWindow.h"

#include <iostream>
#include <string>

namespace itkMeshTestTypes { 
  // this namespace helps to isolate the types defined blow
  // when all the code is included in the test driver.
 

  
/**
 * Some typedefs to make things easier.
 */

/**
 * Define a mesh type that stores a PixelType of "int".  Use the defaults
 * for the other template parameters.
 */
typedef itk::Mesh<int>  MeshType;
typedef MeshType::CellTraits  CellTraits;


/**
 * Define a few cell types which uses a PixelType of "int".  Again,
 * use the defaults for the other parameters.  Note that a cell's template
 * parameters must match those of the mesh into which it is inserted.
 */
typedef itk::CellInterface< int, CellTraits >           CellInterfaceType;
typedef itk::LineCell<CellInterfaceType>                LineCellType;
typedef itk::TetrahedronCell<CellInterfaceType>         TetraCellType;
typedef itk::HexahedronCell<CellInterfaceType>          HexaCellType;
typedef itk::QuadraticEdgeCell<CellInterfaceType>       QuadraticEdgeCellType;
typedef itk::QuadraticTriangleCell<CellInterfaceType>   QuadraticTriangleCellType;



/**
 * Typedef the generic cell type for the mesh.  It is an abstract class,
 * so we can only use information from it, like get its pointer type.
 */
typedef MeshType::CellType              CellType;
typedef CellType::CellAutoPointer       CellAutoPointer;

/**
 * The type of point stored in the mesh. Because mesh was instantiated
 * with defaults (itkDefaultStaticMeshTraits), the point dimension is 3 and
 * the coordinate representation is float.
 */
typedef MeshType::PointType  PointType;

/**
 * The mesh that is created consists of a single hexahedron and a single
 * tetrahedron. (The tetra is inside of the hex.)
 */


// create a class to store counts of cells found in the visit pass
class CountClass
{
public:
  CountClass()
    {
      m_Tetra =0;
      m_QuadraticEdgeCell =0;
      m_QuadraticTriangleCellType =0;
    }
  int m_Tetra;
  int m_QuadraticEdgeCell;
  int m_QuadraticTriangleCellType;
};


// Create a class that can be used to visit cells of
// different types via overloading the Visit method
class VisitCells
{
public:  
  void SetCountClass(CountClass* c) 
    {
      m_CountClass = c;
    }

  void Visit(unsigned long , TetraCellType*)
    {
      m_CountClass->m_Tetra++;
    }
  void Visit(unsigned long , QuadraticEdgeCellType*)
    {
      m_CountClass->m_QuadraticEdgeCell++;
    }

  void Visit(unsigned long , QuadraticTriangleCellType*)
    {
      m_CountClass->m_QuadraticTriangleCellType++;
    }
  CountClass* m_CountClass;
  VisitCells()
    {
      m_CountClass = 0;
    };
};

typedef itk::CellInterfaceVisitorImplementation<
  int, MeshType::CellTraits,
  TetraCellType,
  VisitCells> TetraCellVisitor;

typedef itk::CellInterfaceVisitorImplementation<
  int, MeshType::CellTraits,
  QuadraticEdgeCellType,
  VisitCells> QuadraticEdgeCellVisitor;

typedef itk::CellInterfaceVisitorImplementation<
  int, MeshType::CellTraits,
  QuadraticTriangleCellType,
  VisitCells> QuadraticTriangleCellVisitor;

}


int itkMeshTest(int, char* [] )
{

  using namespace itkMeshTestTypes; // open the namespace here.
                                    // this is safe because only happens localy.
  
  itk::FileOutputWindow::Pointer fow = itk::FileOutputWindow::New();
  fow->SetInstance(fow);
  
  /**
   * Define the 3d geometric positions for 8 points in a cube.
   */
  MeshType::CoordRepType testPointCoords[8][3]
    = { {0,0,0}, {9,0,0}, {9,0,9}, {0,0,9},
        {0,9,0}, {9,9,0}, {9,9,9}, {0,9,9} };
  
  /**
   * List the points that the tetrahedron will use from the mesh.
   */
  unsigned long tetraPoints[4] = {0,1,2,4};
  
  /**
   * List the points that the hexahedron will use from the mesh.
   */
  unsigned long hexaPoints[8] = {0,1,2,3,4,5,6,7};
  
  
  /**
   * Create the mesh through its object factory.
   */
  MeshType::Pointer mesh = MeshType::New();  
  mesh->DebugOn();

  /**
   * Add our test points to the mesh.
   * mesh->SetPoint(pointId, point)
   * Note that the constructor for Point is public, and takes an array
   * of coordinates for the point.
   */
  for(int i=0; i < 8 ; ++i)
    {
    mesh->SetPoint(i, PointType(testPointCoords[i]));
    }

  /** 
   * Specify the method used for allocating cells
   */
   mesh->SetCellsAllocationMethod( MeshType::CellsAllocatedDynamicallyCellByCell );

  /**
   * Create the test cell. Note that testCell is a generic auto
   * pointer to a cell; in this example it ends up pointing to
   * different types of cells.
   */
  CellAutoPointer testCell; 
  testCell.TakeOwnership(  new TetraCellType ); // polymorphism

  /**
   * Assign the points to the tetrahedron through their identifiers.
   */
  testCell->SetPointIds(tetraPoints);

  /**
   * Add the test cell to the mesh.
   * mesh->SetCell(cellId, cell)
   */
  mesh->SetCell(0, testCell ); // Transfer ownership to the mesh
  std::cout << "TetrahedronCell pointer = " << (void const *)testCell.GetPointer() << std::endl;
  std::cout << "TetrahedronCell Owner   = " << testCell.IsOwner() << std::endl;
  
  /**
   * Create another test cell.
   */
  testCell.TakeOwnership( new HexaCellType ); // polymorphism
  testCell->SetPointIds(hexaPoints);
  mesh->SetCell(1, testCell ); // Internally transfers ownership to the mesh


  CellAutoPointer hexaCell;
  if(  mesh->GetCell(1, hexaCell) )
    {
    std::cout << "Hexahedron cell recovered" << std::endl;
    std::cout << "GetNameOfClass()    = " << hexaCell->GetNameOfClass() << std::endl;
    std::cout << "GetNumberOfPoints() = " << hexaCell->GetNumberOfPoints() << std::endl;
    }
  else
    {
    std::cout << "Failure: hexahedron cell was not recovered" << std::endl;
    return EXIT_FAILURE;
    }

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

  std::cout << typeid( cellPointer0 ).name() << std::endl;
  std::cout << typeid( cellPointer0.GetPointer() ).name() << std::endl;
  std::cout << "GetCellBoundaryFeature() return AutoPointer owner = " << cellPointer0.IsOwner() << std::endl;
  
  HexaCellType::FaceType * quad;
  try
    {
    quad = dynamic_cast<HexaCellType::FaceType *>(  cellPointer0.GetPointer() );
    std::cout << "Quad face recovered " << std::endl;
    std::cout << quad->GetNameOfClass() << std::endl;
    }
  catch(...)
    {
    std::cout << "CellPointer cannot be down-cast to a QuadCellType" << std::endl;
    quad = 0;
    }
  if( quad )
    {
    std::cout << "CellPointer was safely down-casted to a QuadCellType" << std::endl;
    }

  if( faceExists ) 
    {
    std::cout << cellPointer0->GetNumberOfPoints() << std::endl;
    std::cout << cellPointer0->GetNameOfClass() << std::endl;
    }
  else 
    {
    std::cout << "Hexahedron face couldn't be extracted " << std::endl;
    }
  
  /**
   * Allocate an explicit boundary line.
   */
  CellAutoPointer boundLine; 
  boundLine.TakeOwnership(  new LineCellType );
  
  /**
   * We don't want the hexahedron to consider the tetrahedron a neighbor
   * across its first edge, so don't add the tetrahedron as a using cell.
   */
  boundLine->AddUsingCell(1);
  boundLine->SetPointId(0,0);
  boundLine->SetPointId(1,1);

  mesh->SetCell( 2,             // New ID for new cell
                 boundLine);    // New cell being added
  
  mesh->SetBoundaryAssignment(1,  // Topological dimension.
            1,                    // CellIdentifier
            0,                    // CellFeatureIdentifier
            2);                   // Cell ID of boundary
  std::cout << "boundLine.IsOwner() = " << boundLine.IsOwner() << std::endl;
  std::cout << "boundLine.GetPointer() = " << boundLine.GetPointer() << std::endl;
  
  /**
   * Try getting the hexahedron's neighbor through its first edge.
   * This should be the test tetrahedron, except that we have done an
   * explicit assignment which removes this.
   */
  std::set<MeshType::CellIdentifier>  neighborSet;
  std::set<MeshType::CellIdentifier>::iterator cell;
  mesh->GetCellBoundaryFeatureNeighbors(
    1,              // Topological dimension of feature.
    1,              // CellIdentifier
    0,              // CellFeatureIdentifier
    &neighborSet); // Where to put result.

  std::cout << "Neighbors (hex edge 0):" << 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 querying the for the number of neighbors the hexahedron has
   * through its second edge, without getting the actual set of