[Insight-users] Contrib: TriangleCell with centroid and normal

Mon Jul 25 05:46:52 EDT 2005

Hello all,

I needed to have the centroid and normal of the TriangleCells in a Mesh
available for a metric.  To support this, and to avoid adding the extra
baggage to the normal TriangleCell (where people might not want to pay for
the extra storage) I've created a new class, imaginatively called
TriangleCellNC.  (Suggestions for a better name welcome!)

All the TriangleCellNC class does is keep track of the extra two attributes,
and allow them to be set and queried.  Then I've written a visitor (that
gets attached to the mesh) and traverses all the faces, calculating the
centroid and normal.  The normal calculation assumes the counter-clockwise
winding rule (as in OpenGL) so the point indices need to be in CCW order for
"outside" to be inferred correctly.

There's also a test program that creates a mesh, does the above
calculations, then prints out the result.

I'd like to contribute this. All feedback welcome!

  :: Gavin

-- 
Gavin Baker                                      Complex Systems Group
http://www.cs.mu.oz.au/~gavinb             The University of Melbourne
-------------- next part --------------
/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: $
  Language:  C++
  Date:      $Date: $
  Version:   $Revision: $

  Copyright (c) 2005 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.

=========================================================================*/
#ifndef __itkTriangleCellNC_h
#define __itkTriangleCellNC_h

#include "itkTriangleCell.h"
#include "itkCovariantVector.h"

namespace itk
{

/** \class TriangleCellNC
 * TriangleCell represents a triangle for a Mesh, with a centroid and normal.
 *
 * Template parameters for TriangleCell:
 *
 * TCellInterface =
 *     The cell interface from the mesh
 *
 * This behaves identically to the TriangleCell, except extra storage is
 * added to track the centroid and the normal.  The calculation of the
 * centroid and normal are carried out by
 * TriangleCellCentroidAndNormalVisitor, which needs to be applied to the
 * mesh to update these attributes.
 *
 * \ingroup MeshObjects
 *
 * \author Gavin Baker, Complex Systems Group, The University of Melbourne
 */

template < typename TCellInterface >
class ITK_EXPORT TriangleCellNC: public TriangleCell<TCellInterface>
{
public:
  itkTypeMacro(TriangleCellNC, TriangleCell);

  /** Standard cell class typedefs */
  itkCellCommonTypedefs(TriangleCellNC);
  itkCellInheritedTypedefs(TriangleCell<TCellInterface>);

  /** The triangle surface normal */
  typedef CovariantVector<PixelType> CovariantType;

  TriangleCellNC() {}

  /* */

  void SetCentroid( const PointType& p )
  {
    m_Centroid = p;
  }

  void GetCentroid( PointType* p ) const
  {
    *p = m_Centroid;
  }

  void SetNormal( const CovariantType& n )
  {
    m_Normal = n;
  }

  void GetNormal( CovariantType* n ) const
  {
    *n = m_Normal;
  }

protected:

  /** The centroid of the cell, average position in each dimension */
  PointType m_Centroid;

  /** The surface normal */
  CovariantType m_Normal;

private:
  TriangleCellNC(const Self&); //purposely not implemented
  void operator=(const Self&); //purposely not implemented

};

} // end namespace itk

#endif
-------------- next part --------------
/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: $
  Language:  C++
  Date:      $Date: $
  Version:   $Revision: $

  Copyright (c) 2005 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.

=========================================================================*/
#ifndef __itkTriangleCellCentroidAndNormalVisitor_h
#define __itkTriangleCellCentroidAndNormalVisitor_h

#include "itkTriangleCellNC.h"

#include <vnl/vnl_vector.h>
#include <vnl/vnl_cross.h>

namespace itk
{

/** \class TriangleCellCentroidAndNormalVisitor
 * Visits a mesh, calculating the centroid and surface normal for each face.
 *
 * Specifically designed to work with the TriangleCellNC, since it has the
 * added support for keeping track of its own centroid and normal.
 *
 * NB: Assumes the counter-clockwise (CCW) winding rule for determining the
 * inside and outside of a face (the default in OpenGL).
 *
 * \ingroup MeshObjects
 *
 * \author Gavin Baker, Complex Systems Group, The University of Melbourne
 */

template<typename TMeshType>
class ITK_EXPORT TriangleCellCentroidAndNormalVisitor
{
public:

  typedef TMeshType                         MeshType;
  typedef typename MeshType::CellType       CellType;
  typedef typename MeshType::PointType      PointType;
  typedef itk::TriangleCellNC<CellType>     TriType;
  typedef typename TriType::CovariantType   CovariantType;
  typedef vnl_vector<float>                 vectorf;

  /**
   * The mesh must be supplied so we can retrieve the points.
   */
  void SetMesh( MeshType* mesh )
  {
    m_Mesh = mesh;
  }

  /**
   * Visit each TriangleCellNC, retrieve its points from the mesh, and
   * calculate its centroid (by averaging in each dimension) and its normal
   * (by taking the cross-product of vectors constructed according to the
   * CCW winding rule).
   */
  void Visit( unsigned long cellid, TriType* tri )
  {
    if ( ! m_Mesh )
      {
      throw "a hissy fit";
      }

    typename TriType::PointIdIterator pointId = tri->PointIdsBegin();

    // Centroid

    PointType p0;
    PointType p1;
    PointType p2;

    m_Mesh->GetPoint( *pointId++, &p0 );
    m_Mesh->GetPoint( *pointId++, &p1 );
    m_Mesh->GetPoint( *pointId++, &p2 );

    PointType centroid;

    centroid[0] = (p0[0]+p1[0]+p2[0])/3.0;
    centroid[1] = (p0[1]+p1[1]+p2[1])/3.0;
    centroid[2] = (p0[2]+p1[2]+p2[2])/3.0;

    tri->SetCentroid( centroid );

    // Face Normal

    vectorf v0 = p1.GetVnlVector() - p0.GetVnlVector();
    vectorf v1 = p2.GetVnlVector() - p0.GetVnlVector();

    vectorf normal = vnl_cross_3d( v0, v1 );

    CovariantType cov_normal;
    cov_normal.Set_vnl_vector( normal );
    cov_normal.Normalize();

    tri->SetNormal( cov_normal );
  }

protected:

  typename MeshType::Pointer m_Mesh;
};

} /* namespace itk */

#endif
-------------- next part --------------
/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: $
  Language:  C++
  Date:      $Date: $
  Version:   $Revision: $

     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.

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

#include "itkMesh.h"
#include "itkTriangleCell.h"
#include "itkTriangleCellCentroidAndNormalVisitor.h"

typedef float                             PixelType;
typedef itk::Mesh< PixelType, 3 >         MeshType;

typedef MeshType::CellType                CellType;

typedef itk::TriangleCellNC< CellType >   TriangleType;

// Simple printing visitor

template <typename TMesh>
class PrintingVisitor
  {
      typedef TMesh MeshType;
      typedef typename TriangleType::PointType PointType;
      typedef typename TriangleType::CovariantType CovariantType;

  public:

      void SetMesh( MeshType* mesh )
          {
              m_Mesh = mesh;
          }

      void Visit(unsigned long cellId, TriangleType * tri )
          {
              std::cout << "Cell " << cellId << std::endl;

              PointType p;
              tri->GetCentroid( &p );
              std::cout << "    Centroid: " << p << std::endl;

              CovariantType n;
              tri->GetNormal( &n );
              std::cout << "    Normal: " << n << std::endl;
          }

  protected:

      typename MeshType::Pointer m_Mesh;
};

int main(int, char *[])
{
  MeshType::Pointer  mesh = MeshType::New();

  // Creating the points and inserting them in the mesh
  //
  MeshType::PointType   point0;
  MeshType::PointType   point1;
  MeshType::PointType   point2;
  MeshType::PointType   point3;

  point0[0] = -1; point0[1] = -1; point0[2] = -1; 
  point1[0] =  1; point1[1] =  1; point1[2] = -1; 
  point2[0] =  1; point2[1] = -1; point2[2] =  1; 
  point3[0] = -1; point3[1] =  1; point3[2] =  1; 

  mesh->SetPoint( 0, point0 );
  mesh->SetPoint( 1, point1 );
  mesh->SetPoint( 2, point2 );
  mesh->SetPoint( 3, point3 );

  // Creating and associating the Tetrahedron
  //
  CellType::CellAutoPointer cellpointer;

  // Creating and associating the Triangles
  //
  cellpointer.TakeOwnership( new TriangleType );
  cellpointer->SetPointId( 0, 0 );
  cellpointer->SetPointId( 1, 1 );
  cellpointer->SetPointId( 2, 2 );
  mesh->SetCell( 1, cellpointer );

  cellpointer.TakeOwnership( new TriangleType );
  cellpointer->SetPointId( 0, 0 );
  cellpointer->SetPointId( 1, 2 );
  cellpointer->SetPointId( 2, 3 );
  mesh->SetCell( 2, cellpointer );

  cellpointer.TakeOwnership( new TriangleType );
  cellpointer->SetPointId( 0, 0 );
  cellpointer->SetPointId( 1, 3 );
  cellpointer->SetPointId( 2, 1 );
  mesh->SetCell( 3, cellpointer );

  cellpointer.TakeOwnership( new TriangleType );
  cellpointer->SetPointId( 0, 3 );
  cellpointer->SetPointId( 1, 2 );
  cellpointer->SetPointId( 2, 1 );
  mesh->SetCell( 4, cellpointer );

  // Simple verification of the number of points and cells inserted
  //
  std::cout << "# Points= " << mesh->GetNumberOfPoints() << std::endl;
  std::cout << "# Cell  = " << mesh->GetNumberOfCells() << std::endl;

  // Visitor - have one to calculate and the other to print

  typedef itk::TriangleCellCentroidAndNormalVisitor<MeshType> NCVisitorType;

  typedef PrintingVisitor<MeshType> PrintingVisitorType;

  typedef itk::CellInterfaceVisitorImplementation<
      PixelType, MeshType::CellTraits, TriangleType,
      NCVisitorType> NCTriangleVisitorType;

  typedef itk::CellInterfaceVisitorImplementation<
      PixelType, MeshType::CellTraits, TriangleType,
      PrintingVisitorType> PrintingTriangleVisitorType;

  NCTriangleVisitorType::Pointer triangleVisitor =
      NCTriangleVisitorType::New();

  PrintingTriangleVisitorType::Pointer printingVisitor =
      PrintingTriangleVisitorType::New();

  triangleVisitor->SetMesh( mesh ); 
  printingVisitor->SetMesh( mesh ); 

  // Multi-Visitor

  typedef CellType::MultiVisitor CellMultiVisitorType;
  CellMultiVisitorType::Pointer multiVisitor = CellMultiVisitorType::New();  

  // Process the mesh...
  multiVisitor->AddVisitor( triangleVisitor );
  mesh->Accept( multiVisitor );

  // Now print it (We can reuse the same multiVisitor because only one
  // visitor for each topology is kept at a time.  It should probably use a
  // hash map or something, so multiple visitors can actually be used.)
  multiVisitor->AddVisitor( printingVisitor );
  mesh->Accept( multiVisitor );

  return 0;
}