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