1 /*========================================================================= 2 * 3 * Copyright Insight Software Consortium 4 * 5 * Licensed under the Apache License, Version 2.0 (the "License"); 6 * you may not use this file except in compliance with the License. 7 * You may obtain a copy of the License at 8 * 9 * http://www.apache.org/licenses/LICENSE-2.0.txt 10 * 11 * Unless required by applicable law or agreed to in writing, software 12 * distributed under the License is distributed on an "AS IS" BASIS, 13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 * See the License for the specific language governing permissions and 15 * limitations under the License. 16 * 17 *=========================================================================*/ 18 #ifndef itkDiscreteGaussianCurvatureQuadEdgeMeshFilter_h 19 #define itkDiscreteGaussianCurvatureQuadEdgeMeshFilter_h 20 21 #include "itkDiscreteCurvatureQuadEdgeMeshFilter.h" 22 #include "itkMath.h" 23 24 25 namespace itk 26 { 27 /** 28 * \class DiscreteGaussianCurvatureQuadEdgeMeshFilter 29 * \brief see the following paper 30 * title: Discrete Differential-Geometry Operators for Triangulated 2-Manifolds 31 * authors: Mark Meyer, Mathieu Desbrun, Peter Schroder, Alan H. Barr 32 * conference: VisMath '02 33 * location: Berlin (Germany) 34 * \author: Arnaud Gelas, Alexandre Gouaillard 35 * \ingroup ITKQuadEdgeMeshFiltering 36 */ 37 template< typename TInputMesh, typename TOutputMesh=TInputMesh > 38 class DiscreteGaussianCurvatureQuadEdgeMeshFilter: 39 public DiscreteCurvatureQuadEdgeMeshFilter< TInputMesh, TOutputMesh > 40 { 41 public: 42 ITK_DISALLOW_COPY_AND_ASSIGN(DiscreteGaussianCurvatureQuadEdgeMeshFilter); 43 44 using Self = DiscreteGaussianCurvatureQuadEdgeMeshFilter; 45 using Pointer = SmartPointer< Self >; 46 using ConstPointer = SmartPointer< const Self >; 47 using Superclass = DiscreteCurvatureQuadEdgeMeshFilter< 48 TInputMesh, TOutputMesh >; 49 50 using InputMeshType = typename Superclass::InputMeshType; 51 using InputMeshPointer = typename Superclass::InputMeshPointer; 52 53 using OutputMeshType = typename Superclass::OutputMeshType; 54 using OutputMeshPointer = typename Superclass::OutputMeshPointer; 55 using OutputPointsContainerPointer = typename Superclass::OutputPointsContainerPointer; 56 using OutputPointsContainerIterator = typename Superclass::OutputPointsContainerIterator; 57 using OutputPointType = typename Superclass::OutputPointType; 58 using OutputVectorType = typename Superclass::OutputVectorType; 59 using OutputCoordType = typename Superclass::OutputCoordType; 60 using OutputPointIdentifier = typename Superclass::OutputPointIdentifier; 61 using OutputCellIdentifier = typename Superclass::OutputCellIdentifier; 62 using OutputQEType = typename Superclass::OutputQEType; 63 using OutputMeshTraits = typename Superclass::OutputMeshTraits; 64 using OutputCurvatureType = typename Superclass::OutputCurvatureType; 65 using TriangleType = typename Superclass::TriangleType; 66 67 /** Run-time type information (and related methods). */ 68 itkTypeMacro(DiscreteGaussianCurvatureQuadEdgeMeshFilter, DiscreteCurvatureQuadEdgeMeshFilter); 69 70 /** New macro for creation of through a Smart Pointer */ 71 itkNewMacro(Self); 72 73 #ifdef ITK_USE_CONCEPT_CHECKING 74 // Begin concept checking 75 itkConceptMacro( OutputIsFloatingPointCheck, 76 ( Concept::IsFloatingPoint< OutputCurvatureType > ) ); 77 // End concept checking 78 #endif 79 80 protected: 81 DiscreteGaussianCurvatureQuadEdgeMeshFilter() = default; 82 ~DiscreteGaussianCurvatureQuadEdgeMeshFilter() override = default; 83 EstimateCurvature(const OutputPointType & iP)84 OutputCurvatureType EstimateCurvature(const OutputPointType & iP) override 85 { 86 OutputMeshPointer output = this->GetOutput(); 87 88 OutputQEType *qe = iP.GetEdge(); 89 90 if ( qe != nullptr ) 91 { 92 OutputQEType *qe_it = qe; 93 OutputQEType *qe_it2; 94 95 OutputPointType q0, q1; 96 97 OutputCurvatureType sum_theta = 0.; 98 OutputCurvatureType area = 0.; 99 100 do 101 { 102 // cell_id = qe_it->GetLeft(); 103 qe_it2 = qe_it->GetOnext(); 104 q0 = output->GetPoint( qe_it->GetDestination() ); 105 q1 = output->GetPoint( qe_it2->GetDestination() ); 106 107 // Compute Angle; 108 sum_theta += static_cast< OutputCurvatureType >( 109 TriangleType::ComputeAngle(q0, iP, q1) ); 110 area += this->ComputeMixedArea(qe_it, qe_it2); 111 qe_it = qe_it2; 112 } 113 while ( qe_it != qe ); 114 115 return ( 2.0 * itk::Math::pi - sum_theta ) / area; 116 } 117 118 return 0.; 119 } 120 }; 121 } 122 123 #endif 124