ITK  4.1.0
Insight Segmentation and Registration Toolkit
itkSimplexMeshAdaptTopologyFilter.h
Go to the documentation of this file.
00001 /*=========================================================================
00002  *
00003  *  Copyright Insight Software Consortium
00004  *
00005  *  Licensed under the Apache License, Version 2.0 (the "License");
00006  *  you may not use this file except in compliance with the License.
00007  *  You may obtain a copy of the License at
00008  *
00009  *         http://www.apache.org/licenses/LICENSE-2.0.txt
00010  *
00011  *  Unless required by applicable law or agreed to in writing, software
00012  *  distributed under the License is distributed on an "AS IS" BASIS,
00013  *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
00014  *  See the License for the specific language governing permissions and
00015  *  limitations under the License.
00016  *
00017  *=========================================================================*/
00018 #ifndef __itkSimplexMeshAdaptTopologyFilter_h
00019 #define __itkSimplexMeshAdaptTopologyFilter_h
00020 
00021 #include "itkPolygonCell.h"
00022 #include "itkCellInterfaceVisitor.h"
00023 
00024 #include "itkSimplexMesh.h"
00025 #include "itkMeshToMeshFilter.h"
00026 #include "itkVectorContainer.h"
00027 
00028 #include "vxl_version.h"
00029 #if VXL_VERSION_DATE_FULL > 20040406
00030 #include "vnl/vnl_cross.h"
00031 #define itk_cross_3d vnl_cross_3d
00032 #else
00033 #define itk_cross_3d cross_3d
00034 #endif
00035 
00036 namespace itk
00037 {
00049 template< class TInputMesh, class TOutputMesh >
00050 class ITK_EXPORT SimplexMeshAdaptTopologyFilter:public MeshToMeshFilter< TInputMesh, TOutputMesh >
00051 {
00052 public:
00054   typedef SimplexMeshAdaptTopologyFilter Self;
00055 
00057   typedef MeshToMeshFilter< TInputMesh, TOutputMesh > Superclass;
00058 
00060   typedef SmartPointer< Self > Pointer;
00061 
00063   typedef SmartPointer< const Self > ConstPointer;
00064 
00066   itkNewMacro(Self);
00067 
00069   itkTypeMacro(SimplexMeshAdaptTopologyFilter, MeshToMeshFilter);
00070 
00071   typedef TInputMesh                                                InputMeshType;
00072   typedef typename InputMeshType::Pointer                           InputMeshPointer;
00073   typedef typename InputMeshType::PointType                         InputPointType;
00074   typedef typename InputMeshType::VectorType                        InputVectorType;
00075   typedef typename InputMeshType::PixelType                         InputPixelType;
00076   typedef typename InputMeshType::MeshTraits::CellTraits            InputCellTraitsType;
00077   typedef typename InputMeshType::CellType                          InputCellType;
00078   typedef typename InputMeshType::PointIdentifier                   PointIdentifier;
00079   typedef typename InputMeshType::CellIdentifier                    CellIdentifier;
00080   typedef typename InputCellType::PointIdIterator                   InputCellPointIdIterator;
00081   typedef typename InputCellType::CellAutoPointer                   InputCellAutoPointer;
00082   typedef typename InputMeshType::CellAutoPointer                   CellAutoPointer;
00083   typedef          itk::PolygonCell< InputCellType >                InputPolygonType;
00084   typedef typename InputPolygonType::PointIdIterator                InputPolygonPointIdIterator;
00085   typedef CovariantVector< typename InputVectorType::ValueType, 3 > CovariantVectorType;
00086   typedef  TOutputMesh                                              OutputMeshType;
00087   typedef typename OutputMeshType::Pointer                          OutputMeshPointer;
00088   typedef typename OutputMeshType::CellType                         OutputCellType;
00089   typedef          itk::PolygonCell< OutputCellType >               OutputPolygonType;
00090 
00091   typedef typename itk::MapContainer< CellIdentifier, double > DoubleValueMapType;
00092   typedef typename DoubleValueMapType::Iterator                DoubleContainerIterator;
00093 
00101   class SimplexCellVisitor
00102   {
00103 public:
00104     InputMeshPointer            mesh;
00105     double                      totalArea;
00106     double                      totalCurvature;
00107     double                      minCellSize;
00108     double                      maxCellSize;
00109     typename DoubleValueMapType::Pointer areaMap;
00110     typename DoubleValueMapType::Pointer curvatureMap;
00111 
00112     double minCurvature;
00113     double maxCurvature;
00114 
00115     SimplexCellVisitor()
00116     {
00117       areaMap = DoubleValueMapType::New();
00118       curvatureMap = DoubleValueMapType::New();
00119       totalArea = 0;
00120       totalCurvature = 0;
00121       minCellSize = NumericTraits< double >::max();
00122       maxCellSize = 0;
00123       minCurvature = NumericTraits< double >::max();
00124       maxCurvature = 0;
00125     }
00126 
00130     void Visit(CellIdentifier cellId, InputPolygonType *poly)
00131     {
00132       typename InputPolygonType::PointIdIterator it =  poly->PointIdsBegin();
00133 
00134       double        meanCurvature = 0;
00135       PointIdentifier refPoint = *it;
00136       double        val = mesh->GetMeanCurvature(*it++);
00137       meanCurvature += vcl_abs(val);
00138 
00139       PointIdentifier id1 = *it;
00140       val = mesh->GetMeanCurvature(*it++);
00141       meanCurvature += vcl_abs(val);
00142 
00143       PointIdentifier id2;
00144 
00145       double area = 0;
00146 
00147       int cnt = 0;
00148 
00149       while ( it != poly->PointIdsEnd() )
00150         {
00151         id2 = *it;
00152         area += ComputeArea(refPoint, id1, id2);
00153         id1 = id2;
00154         val = mesh->GetMeanCurvature(*it);
00155         meanCurvature += vcl_abs(val);
00156         cnt++;
00157         it++;
00158         }
00159 
00160       meanCurvature /= (double)cnt;
00161       totalArea += area;
00162       totalCurvature += meanCurvature;
00163 
00164       areaMap->InsertElement(cellId, area);
00165       curvatureMap->InsertElement(cellId, meanCurvature);
00166 
00167       if ( area > maxCellSize ) { maxCellSize = area; }
00168       if ( area < minCellSize ) { minCellSize = area; }
00169       if ( meanCurvature > maxCurvature ) { maxCurvature = meanCurvature; }
00170       if ( meanCurvature < minCurvature ) { minCurvature = meanCurvature; }
00171     }
00172 
00173     double ComputeArea(PointIdentifier p1, PointIdentifier p2, PointIdentifier p3)
00174     {
00175       InputPointType v1, v2, v3;
00176 
00177       v1.Fill(0);
00178       v2.Fill(0);
00179       v3.Fill(0);
00180 
00181       mesh->GetPoint(p1, &v1);
00182       mesh->GetPoint(p2, &v2);
00183       mesh->GetPoint(p3, &v3);
00184       return vcl_abs (itk_cross_3d( ( v2 - v1 ).GetVnlVector(), ( v3 - v1 ).GetVnlVector() ).two_norm() / 2.0);
00185     }
00186 
00187     typename DoubleValueMapType::Pointer GetAreaMap()
00188     {
00189       return areaMap;
00190     }
00191 
00192     typename DoubleValueMapType::Pointer GetCurvatureMap()
00193     {
00194       return curvatureMap;
00195     }
00196 
00197     double GetTotalMeshArea()
00198     {
00199       return totalArea;
00200     }
00201 
00202     double GetTotalMeanCurvature()
00203     {
00204       return totalCurvature / ( curvatureMap->Size() );
00205     }
00206 
00207     double GetMaximumCellSize()
00208     {
00209       return maxCellSize;
00210     }
00211 
00212     double GetMinimumCellSize()
00213     {
00214       return minCellSize;
00215     }
00216 
00217     double GetMaximumCurvature()
00218     {
00219       return maxCurvature;
00220     }
00221 
00222     double GetMinimumCurvature()
00223     {
00224       return minCurvature;
00225     }
00226   };
00227 
00228   // cell visitor stuff
00229   typedef itk::CellInterfaceVisitorImplementation< InputPixelType,
00230                                                    InputCellTraitsType,
00231                                                    InputPolygonType,
00232                                                    SimplexCellVisitor >
00233   SimplexVisitorInterfaceType;
00234 
00235   typedef typename SimplexVisitorInterfaceType::Pointer SimplexVisitorInterfacePointer;
00236   typedef typename InputCellType::MultiVisitor          CellMultiVisitorType;
00237   typedef typename CellMultiVisitorType::Pointer        CellMultiVisitorPointer;
00238 
00239   itkSetMacro(Threshold, double);
00240   itkGetConstMacro(Threshold, double);
00241 
00242   itkSetMacro(SelectionMethod, int);
00243   itkGetConstMacro(SelectionMethod, int);
00244 
00245   itkGetConstMacro(ModifiedCount, int);
00246 protected:
00247 
00248   SimplexMeshAdaptTopologyFilter();
00249   ~SimplexMeshAdaptTopologyFilter();
00250   SimplexMeshAdaptTopologyFilter(const Self &) {}
00251 
00252   void operator=(const Self &) {}
00253 
00254   void PrintSelf(std::ostream & os, Indent indent) const;
00255 
00256   virtual void GenerateData();
00257 
00261   void Initialize();
00262 
00268   void ComputeCellParameters();
00269 
00271   void CopyInputMeshToOutputMeshGeometryData();
00272 
00278   void ModifyNeighborCells(CellIdentifier id1, CellIdentifier id2, PointIdentifier insertPointId);
00279 
00283   InputPointType ComputeCellCenter(InputCellAutoPointer & simplexCell);
00284 
00288   CellIdentifier m_IdOffset;
00289 
00294   double m_Threshold;
00295 
00299   int m_SelectionMethod;
00300 
00305   int m_ModifiedCount;
00306 
00311   OutputMeshPointer m_Output;
00312 
00313   InputCellAutoPointer m_NewSimplexCellPointer;
00314 };
00315 } //end of namespace
00316 
00317 #ifndef ITK_MANUAL_INSTANTIATION
00318 #include "itkSimplexMeshAdaptTopologyFilter.hxx"
00319 #endif
00320 
00321 #endif // __itkSimplexMeshAdaptTopologyFilter_h
00322