itkSimplexMeshAdaptTopologyFilter.h

Go to the documentation of this file.

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

Program:   Insight Segmentation & Registration Toolkit
Module:    $RCSfile: itkSimplexMeshAdaptTopologyFilter.h,v $
Language:  C++
Date:      $Date: 2007/02/01 18:06:01 $
Version:   $Revision: 1.11 $

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.

=========================================================================*/
#ifndef __itkSimplexMeshAdaptTopologyFilter_h
#define __itkSimplexMeshAdaptTopologyFilter_h

#include "itkMesh.h"
#include "itkPolygonCell.h"
#include "itkMapContainer.h"
#include "itkCellInterfaceVisitor.h"

#include "itkSimplexMesh.h"
#include "itkSimplexMeshGeometry.h"
#include "itkMeshToMeshFilter.h"
#include "itkVectorContainer.h"
#include <itkCovariantVector.h>

#include <vxl_version.h>
#if VXL_VERSION_DATE_FULL > 20040406
# include <vnl/vnl_cross.h>
# define itk_cross_3d vnl_cross_3d
#else
# define itk_cross_3d cross_3d
#endif

namespace itk
  {
  template <class TInputMesh, class TOutputMesh>
class SimplexMeshAdaptTopologyFilter : public MeshToMeshFilter<TInputMesh, TOutputMesh>
  {

  public:
    typedef SimplexMeshAdaptTopologyFilter Self;

    typedef MeshToMeshFilter<TInputMesh, TOutputMesh> Superclass;

    typedef SmartPointer<Self>  Pointer;

    typedef SmartPointer<const Self>  ConstPointer;

    itkNewMacro(Self);

    itkTypeMacro(SimplexMeshAdaptTopologyFilter,MeshToMeshFilter);

    typedef TInputMesh                                              InputMeshType;
    typedef typename InputMeshType::Pointer                         InputMeshPointer;
    typedef typename InputMeshType::PointType                       InputPointType;
    typedef typename InputMeshType::VectorType                      InputVectorType;
    typedef typename InputMeshType::PixelType                       InputPixelType;
    typedef typename InputMeshType::MeshTraits::CellTraits          InputCellTraitsType;
    typedef typename InputMeshType::CellType                        InputCellType;
    typedef typename InputCellType::PointIdIterator                 InputCellPointIdIterator;
    typedef typename InputCellType::CellAutoPointer                 InputCellAutoPointer;
    typedef typename InputMeshType::CellAutoPointer                 CellAutoPointer;
    typedef          itk::PolygonCell<InputCellType>                InputPolygonType;
    typedef typename InputPolygonType::PointIdIterator              InputPolygonPointIdIterator;

    typedef CovariantVector< 
                    typename InputVectorType::ValueType, 3 >        CovariantVectorType;
    typedef  TOutputMesh                                            OutputMeshType;
    typedef typename OutputMeshType::Pointer                        OutputMeshPointer;
    typedef typename OutputMeshType::CellType                       OutputCellType;
    typedef          itk::PolygonCell<OutputCellType>               OutputPolygonType;

    typedef          itk::MapContainer<unsigned long, double>       DoubleValueMapType;
    typedef typename DoubleValueMapType::Iterator                   DoubleContainerIterator;


    class SimplexCellVisitor
      {

      public:
        InputMeshPointer mesh;
        double totalArea;
        double totalCurvature;
        double minCellSize;
        double maxCellSize;
        DoubleValueMapType::Pointer areaMap;
        DoubleValueMapType::Pointer curvatureMap;

        double minCurvature;
        double maxCurvature;

        SimplexCellVisitor()
          {
          areaMap = DoubleValueMapType::New();
          curvatureMap = DoubleValueMapType::New();
          totalArea = 0;
          totalCurvature = 0;
          minCellSize = NumericTraits<double>::max();
          maxCellSize = 0;
          minCurvature = NumericTraits<double>::max();
          maxCurvature = 0;
          }


        void Visit(unsigned long cellId, InputPolygonType * poly)
          {
          typename InputPolygonType::PointIdIterator it =  poly->PointIdsBegin();

          double meanCurvature = 0;
          unsigned long refPoint = *it;
          double val = mesh->GetMeanCurvature(*it++);
          meanCurvature += vcl_abs(val);

          unsigned long id1 = *it;
          val = mesh->GetMeanCurvature(*it++);
          meanCurvature += vcl_abs(val);

          unsigned long id2;

          double area = 0;

          int cnt = 0; 

          while ( it != poly->PointIdsEnd() )
            {
            id2 = *it;
            area += ComputeArea(refPoint,id1,id2);
            id1 = id2;
            val = mesh->GetMeanCurvature(*it);
            meanCurvature += vcl_abs(val);
            cnt++;
            it++;
            }

          meanCurvature /= (double)cnt;
          totalArea += area;
          totalCurvature += meanCurvature;

          areaMap->InsertElement(cellId, area);
          curvatureMap->InsertElement(cellId, meanCurvature);

          if (area > maxCellSize ) maxCellSize = area;
          if (area < minCellSize ) minCellSize = area;
          if (meanCurvature > maxCurvature ) maxCurvature = meanCurvature;
          if (meanCurvature < minCurvature ) minCurvature = meanCurvature;
          }

        double ComputeArea(unsigned long p1,unsigned long p2, unsigned long p3)
          {
          InputPointType v1,v2,v3;
          mesh->GetPoint(p1, &v1);
          mesh->GetPoint(p2, &v2);
          mesh->GetPoint(p3, &v3);
          return vcl_abs (itk_cross_3d((v2-v1).GetVnlVector(), (v3-v1).GetVnlVector()).two_norm() /2.0);
          }

        DoubleValueMapType::Pointer GetAreaMap()
          {
          return areaMap;
          }

        DoubleValueMapType::Pointer GetCurvatureMap()
          {
          return curvatureMap;
          }

        double GetTotalMeshArea()
          {
          return totalArea;
          }

        double GetTotalMeanCurvature()
          {
          return totalCurvature/(curvatureMap->Size());
          }

        double GetMaximumCellSize()
          {
          return maxCellSize;
          }

        double GetMinimumCellSize()
          {
          return minCellSize;
          }
        double GetMaximumCurvature()
          {
          return maxCurvature;
          }

        double GetMinimumCurvature()
          {
          return minCurvature;
          }

      };

    // cell visitor stuff
    typedef itk::CellInterfaceVisitorImplementation<InputPixelType,
      InputCellTraitsType,
      InputPolygonType,
      SimplexCellVisitor>
      SimplexVisitorInterfaceType;

    typedef typename SimplexVisitorInterfaceType::Pointer   SimplexVisitorInterfacePointer;
    typedef typename InputCellType::MultiVisitor            CellMultiVisitorType;
    typedef typename CellMultiVisitorType::Pointer          CellMultiVisitorPointer;


    itkSetMacro(Threshold, double);
    itkGetMacro(Threshold, double);

    itkSetMacro(SelectionMethod, int);
    itkGetMacro(SelectionMethod, int);

    itkGetMacro(ModifiedCount, int);


  protected:

    SimplexMeshAdaptTopologyFilter();

    ~SimplexMeshAdaptTopologyFilter();

    SimplexMeshAdaptTopologyFilter(const Self&) 
      {
      }

    void operator=(const Self&)
      {
      }

    void PrintSelf(std::ostream& os, Indent indent) const;

    virtual void GenerateData();


    void Initialize();

    void ComputeCellParameters();

    void InsertNewCells();

    void ModifyNeighborCells(unsigned long id1, unsigned long id2, unsigned long insertPointId);


    InputPointType ComputeCellCenter(InputCellAutoPointer &simplexCell);


    unsigned long m_IdOffset;


    double m_Threshold;

    int  m_SelectionMethod;

    int m_ModifiedCount;

    OutputMeshPointer m_Output;

    InputCellAutoPointer  NewSimplexCellPointer;

  };

  } //end of namespace

#ifndef ITK_MANUAL_INSTANTIATION
#include "itkSimplexMeshAdaptTopologyFilter.txx"
#endif

#endif // __itkSimplexMeshAdaptTopologyFilter_h

---