itkSimplexMeshAdaptTopologyFilter.h

Go to the documentation of this file.

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

Program:   Insight Segmentation & Registration Toolkit
Module:    $RCSfile: itkSimplexMeshAdaptTopologyFilter.h,v $
Language:  C++
Date:      $Date: 2008-10-17 20:49:56 $
Version:   $Revision: 1.12 $

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  m_NewSimplexCellPointer;

};

} //end of namespace

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

#endif // __itkSimplexMeshAdaptTopologyFilter_h