itkquadedgemesheuleroperatorstest.cxx

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

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

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkQuadEdgeMeshEulerOperatorsTest.cxx,v $
  Language:  C++
  Date:      $Date: 2008-07-02 22:03:58 $
  Version:   $Revision: 1.21 $

  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 <string>

#include "itkQuadEdgeMesh.h"
#include "itkQuadEdgeMeshLineCell.h"
#include "itkQuadEdgeMeshPolygonCell.h"

#include "itkQuadEdgeMeshFunctionBase.h"

#include "itkQuadEdgeMeshEulerOperatorJoinFacetFunction.h"
#include "itkQuadEdgeMeshEulerOperatorSplitFacetFunction.h"
#include "itkQuadEdgeMeshEulerOperatorFlipEdgeFunction.h"
#include "itkQuadEdgeMeshEulerOperatorJoinVertexFunction.h"
#include "itkQuadEdgeMeshEulerOperatorSplitVertexFunction.h"
#include "itkQuadEdgeMeshEulerOperatorSplitEdgeFunction.h"
#include "itkQuadEdgeMeshEulerOperatorCreateCenterVertexFunction.h"
#include "itkQuadEdgeMeshEulerOperatorDeleteCenterVertexFunction.h"

#include "itkQuadEdgeMeshTopologyChecker.h"

typedef unsigned long IdentifierType;

template<class TMesh >
bool AssertTopologicalInvariants( TMesh *mesh, 
                                  IdentifierType NumVertices,
                                  IdentifierType NumFaces,
                                  IdentifierType NumEdges,
                                  IdentifierType NumBorders,
                                  IdentifierType Genus)
{
  typedef itk::QuadEdgeMeshTopologyChecker< TMesh > CheckerType;
  typename CheckerType::Pointer check = CheckerType::New();
  check->SetMesh( mesh );
  check->SetExpectedNumberOfPoints( NumVertices );
  check->SetExpectedNumberOfEdges( NumFaces );
  check->SetExpectedNumberOfFaces( NumEdges );
  check->SetExpectedNumberOfBoundaries( NumBorders );
  check->SetExpectedGenus( Genus );
  return( check->ValidateEulerCharacteristic( ) );
}

template< class TMesh >
void PopulateMesh( typename TMesh::Pointer mesh )
{
  typedef TMesh                         MeshType;
  typedef typename MeshType::CellType   CellType;

  typedef itk::QuadEdgeMeshPolygonCell< CellType > QEPolygonCellType;
  
  if( mesh->GetNumberOfPoints( ) )
    {
    mesh->Clear( );
    mesh->ClearFreePointAndCellIndexesLists();
    }

  /////////////////////////////////////////////////////////////
  int expectedNumPts = 25;
  int expectedNumCells = 32;
  int simpleSquareCells[96] =
  {  0,  1,  6,
     0,  6,  5,
     1,  2,  7,
     1,  7,  6,
     2,  3,  8,
     2,  8,  7,
     3,  4,  9,
     3,  9,  8,
     5,  6, 11,
     5, 11, 10,
     6,  7, 12,
     6, 12, 11,
     7,  8, 13,
     7, 13, 12,
     8,  9, 14,
     8, 14, 13,
    10, 11, 16,
    10, 16, 15,
    11, 12, 17,
    11, 17, 16,
    12, 13, 18,
    12, 18, 17,
    13, 14, 19,
    13, 19, 18,
    15, 16, 21,
    15, 21, 20,
    16, 17, 22,
    16, 22, 21,
    17, 18, 23,
    17, 23, 22,
    18, 19, 24,
    18, 24, 23 };

  typename MeshType::PointType pts[25];

  pts[ 0][0] = 0.0;  pts[ 0][1] = 0.0;  pts[ 0][2] = 0.0;
  pts[ 1][0] = 1.0;  pts[ 1][1] = 0.0;  pts[ 1][2] = 0.0;
  pts[ 2][0] = 2.0;  pts[ 2][1] = 0.0;  pts[ 2][2] = 0.0;
  pts[ 3][0] = 3.0;  pts[ 3][1] = 0.0;  pts[ 3][2] = 0.0;
  pts[ 4][0] = 4.0;  pts[ 4][1] = 0.0;  pts[ 4][2] = 0.0;
  pts[ 5][0] = 0.0;  pts[ 5][1] = 1.0;  pts[ 5][2] = 0.0;
  pts[ 6][0] = 1.0;  pts[ 6][1] = 1.0;  pts[ 6][2] = 0.0;
  pts[ 7][0] = 2.0;  pts[ 7][1] = 1.0;  pts[ 7][2] = 0.0;
  pts[ 8][0] = 3.0;  pts[ 8][1] = 1.0;  pts[ 8][2] = 0.0;
  pts[ 9][0] = 4.0;  pts[ 9][1] = 1.0;  pts[ 9][2] = 0.0;
  pts[10][0] = 0.0;  pts[10][1] = 2.0;  pts[10][2] = 0.0;
  pts[11][0] = 1.0;  pts[11][1] = 2.0;  pts[11][2] = 0.0;
  pts[12][0] = 2.0;  pts[12][1] = 2.0;  pts[12][2] = 0.0;
  pts[13][0] = 3.0;  pts[13][1] = 2.0;  pts[13][2] = 0.0;
  pts[14][0] = 4.0;  pts[14][1] = 2.0;  pts[14][2] = 0.0;
  pts[15][0] = 0.0;  pts[15][1] = 3.0;  pts[15][2] = 0.0;
  pts[16][0] = 1.0;  pts[16][1] = 3.0;  pts[16][2] = 0.0;
  pts[17][0] = 2.0;  pts[17][1] = 3.0;  pts[17][2] = 0.0;
  pts[18][0] = 3.0;  pts[18][1] = 3.0;  pts[18][2] = 0.0;
  pts[19][0] = 4.0;  pts[19][1] = 3.0;  pts[19][2] = 0.0;
  pts[20][0] = 0.0;  pts[20][1] = 4.0;  pts[20][2] = 0.0;
  pts[21][0] = 1.0;  pts[21][1] = 4.0;  pts[21][2] = 0.0;
  pts[22][0] = 2.0;  pts[22][1] = 4.0;  pts[22][2] = 0.0;
  pts[23][0] = 3.0;  pts[23][1] = 4.0;  pts[23][2] = 0.0;
  pts[24][0] = 4.0;  pts[24][1] = 4.0;  pts[24][2] = 0.0;

  for(int i=0; i<expectedNumPts; i++)
    {
    mesh->SetPoint( i, pts[i] );
    }
 
  typename CellType::CellAutoPointer cellpointer;
  QEPolygonCellType *poly;

  for(int i=0; i<expectedNumCells; i++)
    {
    poly = new QEPolygonCellType( 3 );
    cellpointer.TakeOwnership( poly );
    cellpointer->SetPointId( 0, simpleSquareCells[3*i] );
    cellpointer->SetPointId( 1, simpleSquareCells[3*i+1] );
    cellpointer->SetPointId( 2, simpleSquareCells[3*i+2] );
    mesh->SetCell( i, cellpointer );
    }
}


int itkQuadEdgeMeshEulerOperatorsTest(int argc, char * argv[])
{
  (void)argc;
  (void)argv;

  typedef itk::QuadEdgeMesh< double, 3 >                      MeshType;
  typedef MeshType::Pointer                                   MeshPointer;
  typedef MeshType::QEType                                    QEType;
  typedef MeshType::PointIdentifier                           PointIdentifier;
  typedef MeshType::PointType                                 PointType;
  typedef MeshType::CellType                                  CellType;
  typedef itk::QuadEdgeMeshLineCell< CellType >               LineCellType;
  
  typedef itk::QuadEdgeMeshEulerOperatorJoinFacetFunction< MeshType, QEType>
    JoinFacet;
  typedef itk::QuadEdgeMeshEulerOperatorSplitFacetFunction< MeshType, QEType>
    SplitFacet;
  typedef itk::QuadEdgeMeshEulerOperatorFlipEdgeFunction< MeshType, QEType>
    FlipEdge;
  typedef itk::QuadEdgeMeshEulerOperatorJoinVertexFunction< MeshType, QEType>
    JoinVertex;
  typedef itk::QuadEdgeMeshEulerOperatorSplitEdgeFunction< MeshType, QEType>
    SplitEdge;
  typedef itk::QuadEdgeMeshEulerOperatorSplitVertexFunction< MeshType, QEType>
    SplitVertex;
  typedef itk::QuadEdgeMeshEulerOperatorCreateCenterVertexFunction< MeshType, QEType>
    CreateCenterVertex;
  typedef itk::QuadEdgeMeshEulerOperatorDeleteCenterVertexFunction< MeshType, QEType>
    DeleteCenterVertex;

  // TEST TOPOLOGYCHECKER
    {
    typedef itk::QuadEdgeMeshTopologyChecker< MeshType > CheckerType;
    CheckerType::Pointer check = CheckerType::New();
    
    // test no input
    if( check->ValidateEulerCharacteristic( ) )
      {
      std::cout << "FAILED." << std::endl;
      return 1;
      }
    std::cout << "OK" << std::endl;
    
    // test with an isolated edge
    MeshPointer testmesh = MeshType::New( );  // empty mesh

    PointType dummyPoint1;
    PointType dummyPoint2;
    for( unsigned char i = 0; i < 3; i++ )
      {
      dummyPoint1[i] = 0.;
      dummyPoint2[i] = 1.;
      } 
    testmesh->SetPoint( 0, dummyPoint1 );     // add points to mesh
    testmesh->SetPoint( 1, dummyPoint2 );

    LineCellType * Line = new LineCellType;   // create an Isolated Edge
    CellType::CellAutoPointer cellpointer;    // create the corresponding AutoPointer
    cellpointer.TakeOwnership( Line );
    cellpointer->SetPointId( 0, 0 );
    cellpointer->SetPointId( 1, 1 );
    testmesh->SetCell( 0, cellpointer );      // add the cell to the mesh

    check->SetMesh( testmesh );
    if( check->ValidateEulerCharacteristic( ) )
      {
      std::cout << "FAILED." << std::endl;
      return 1;
      }
    std::cout << "OK" << std::endl;
    }

  // EULER OPERATOR TESTS
  QEType * dummy;
  MeshPointer  mesh = MeshType::New();
  PopulateMesh<MeshType>( mesh );

  // The initial configuration and numbering of simpleSquare.vtk:
  //    Vertices: 25 , Edges: 56, Faces: 32, Boundary = 1, Chi = 1
  //
  //   20 --------- 21 --------- 22 --------- 23 --------- 24
  //    |        __/ |        __/ |        __/ |        __/ |
  //    |     __/    |     __/    |     __/    |     __/    |
  //    |  __/       |  __/       |  __/       |  __/       |
  //    | /          | /          | /          | /          |
  //   15 --------- 16 --------- 17 --------- 18 --------- 19
  //    |        __/ |        __/ |        __/ |        __/ |
  //    |     __/    |     __/    |     __/    |     __/    |
  //    |  __/       |  __/       |  __/       |  __/       |
  //    | /          | /          | /          | /          |
  //   10 --------- 11 --------- 12 --------- 13 --------- 14
  //    |        __/ |        __/ |        __/ |        __/ |
  //    |     __/    |     __/    |     __/    |     __/    |
  //    |  __/       |  __/       |  __/       |  __/       |
  //    | /          | /          | /          | /          |
  //    5 ---------- 6 ---------- 7 ---------- 8 ---------  9
  //    |        __/ |        __/ |        __/ |        __/ |
  //    |     __/    |     __/    |     __/    |     __/    |
  //    |  __/       |  __/       |  __/       |  __/       |
  //    | /          | /          | /          | /          |
  //    0 ---------- 1 ---------- 2  --------- 3 ---------  4
  //
  /////////////////////////////////////////
  //
  //          Join Facet
  //
  /////////////////////////////////////////

  std::cout << "Checking JointFacet." << std::endl;  

  JoinFacet::Pointer joinFacet = JoinFacet::New( );
  std::cout << "     " << "Test No Mesh Input";
  if( joinFacet->Evaluate( (QEType*)1 ) )
    {
    std::cout << "FAILED." << std::endl;
    return 1;
    }
  std::cout << "OK" << std::endl;
  
  (void)joinFacet->GetNameOfClass(); 

  joinFacet->SetInput( mesh );
  
  std::cout << "     " << "Test QE Input not internal";
  dummy = new QEType;
  if( joinFacet->Evaluate( dummy ) )
    {
    std::cout << "FAILED." << std::endl;
    return 1;
    }
  delete dummy;
  std::cout << "OK" << std::endl;
  
  std::cout << "     " << "Test No QE Input";
  if( joinFacet->Evaluate( (QEType*)0 ) )
    {
    std::cout << "FAILED." << std::endl;
    return 1;
    }
  std::cout << "OK" << std::endl;

  //
  //   20 --------- 21 --------- 22 --------- 23 --------- 24
  //    |        __/ |        __/ |        __/ |        __/ |
  //    |     __/    |     __/    |     __/    |     __/    |
  //    |  __/       |  __/       |  __/       |  __/       |
  //    | /          | /          | /          | /          |
  //   15 --------- 16 --------- 17 --------- 18 --------- 19
  //    |        __/ |        __/ |        __/ |        __/ |
  //    |     __/    |     __/    |     __/    |     __/    |
  //    |  __/       |  __/       |  __/       |  __/       |
  //    | /          | /          | /          | /          |
  //   10 --------- 11 --------- 12 ---<-H--- 13 --------- 14
  //    |        __/ |        __/          __/ |        __/ |
  //    |     __/    |     __/ new      __/    |     __/    |
  //    |  __/       |  __/     face __/       |  __/       |
  //    | /          | /            /          | /          |
  //    5 ---------- 6 ----G->--- 7 ---------- 8 ---------  9
  //    |        __/ |        __/ |        __/ |        __/ |
  //    |     __/    |     __/    |     __/    |     __/    |
  //    |  __/       |  __/       |  __/       |  __/       |
  //    | /          | /          | /          | /          |
  //    0 ---------- 1 ---------- 2  --------- 3 ---------  4
  //
  
  std::cout << "     " << "Test Edge deletion (possible)";
  // Find an internal edge and collapse it
  QEType* DeletedEdge = mesh->FindEdge( 12, 7 );
  // Store G and H for testing the inverse with SplitFace:
  QEType* G = DeletedEdge->GetSym( )->GetLprev( );
  QEType* H = joinFacet->Evaluate( DeletedEdge );
  if( !H )
    {
    std::cout << "FAILED." << std::endl;
    return 1;
    }
   
  // Since the edge was internal we lost an edge and an face:
  if( ! AssertTopologicalInvariants< MeshType >
          ( mesh, 25, 55, 31, 1, 0 ) )