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 itkDelaunayConformingQuadEdgeMeshFilter_h
19 #define itkDelaunayConformingQuadEdgeMeshFilter_h
20 
21 #include "itkIntTypes.h"
22 #include "itkPriorityQueueContainer.h"
23 #include "itkQuadEdgeMeshToQuadEdgeMeshFilter.h"
24 #include "itkQuadEdgeMeshEulerOperatorFlipEdgeFunction.h"
25 #include "itkMath.h"
26 
27 namespace itk
28 {
29 /**
30  *  \class DelaunayConformingQuadEdgeMeshFilter
31  *
32  *  \brief FIXME Add documentation
33  *
34  * \ingroup ITKQuadEdgeMeshFiltering
35  */
36 template< typename TInputMesh, typename TOutputMesh=TInputMesh >
37 class ITK_TEMPLATE_EXPORT DelaunayConformingQuadEdgeMeshFilter:
38   public QuadEdgeMeshToQuadEdgeMeshFilter< TInputMesh, TOutputMesh >
39 {
40 public:
41   /** Basic types. */
42   using Self = DelaunayConformingQuadEdgeMeshFilter;
43   using Superclass = QuadEdgeMeshToQuadEdgeMeshFilter< TInputMesh,
44                                             TOutputMesh >;
45   using Pointer = SmartPointer< Self >;
46   using ConstPointer = SmartPointer< const Self >;
47 
48   /** Input types. */
49   using InputMeshType = TInputMesh;
50   using InputMeshPointer = typename InputMeshType::Pointer;
51   using InputCoordRepType = typename InputMeshType::CoordRepType;
52   using InputPointType = typename InputMeshType::PointType;
53   using InputPointVectorType = typename InputPointType::VectorType;
54   using InputPointIdentifier = typename InputMeshType::PointIdentifier;
55   using InputQEType = typename InputMeshType::QEType;
56   using InputVectorType = typename InputMeshType::VectorType;
57   using InputEdgeListType = typename InputMeshType::EdgeListType;
58   using InputPixelType = typename InputMeshType::PixelType;
59   using InputTraits = typename InputMeshType::Traits;
60 
61   static constexpr unsigned int InputVDimension = InputMeshType::PointDimension;
62 
63   using InputPointsContainer = typename InputMeshType::PointsContainer;
64   using InputPointsContainerConstIterator = typename InputMeshType::PointsContainerConstIterator;
65   using InputCellsContainerConstIterator = typename InputMeshType::CellsContainerConstIterator;
66   using InputEdgeCellType = typename InputMeshType::EdgeCellType;
67   using InputPolygonCellType = typename InputMeshType::PolygonCellType;
68   using InputPointIdList = typename InputMeshType::PointIdList;
69 
70   using InputQEIterator = typename InputQEType::IteratorGeom;
71 
72   /** Output types. */
73   using OutputMeshType = TOutputMesh;
74   using OutputMeshPointer = typename OutputMeshType::Pointer;
75   using OutputCoordRepType = typename OutputMeshType::CoordRepType;
76   using OutputPointType = typename OutputMeshType::PointType;
77   using OutputPointIdentifier = typename OutputMeshType::PointIdentifier;
78   using OutputCellType = typename OutputMeshType::CellType;
79   using OutputCellIdentifier = typename OutputMeshType::CellIdentifier;
80   using OutputEdgeCellType = typename OutputMeshType::EdgeCellType;
81   using OutputQEType = typename OutputMeshType::QEType;
82   using OutputLineCellIdentifier = typename OutputQEType::LineCellIdentifier;
83   using OutputVectorType = typename OutputMeshType::VectorType;
84   using OutputQEIterator = typename OutputQEType::IteratorGeom;
85   using OutputPointsContainerPointer = typename OutputMeshType::PointsContainerPointer;
86   using OutputPointsContainerIterator = typename OutputMeshType::PointsContainerIterator;
87   using OutputCellsContainer = typename OutputMeshType::CellsContainer;
88   using OutputCellsContainerIterator = typename OutputMeshType::CellsContainerIterator;
89 
90   static constexpr unsigned int OutputVDimension = OutputMeshType::PointDimension;
91 
92   itkNewMacro(Self);
93   itkTypeMacro(DelaunayConformingQuadEdgeMeshFilter, QuadEdgeMeshToQuadEdgeMeshFilter);
94 
95   itkGetConstMacro(NumberOfEdgeFlips, SizeValueType);
96 
97 public:
98   using OutputEdgeCellListType = std::list< OutputEdgeCellType * >;
99   using OutputEdgeCellListIterator = typename OutputEdgeCellListType::iterator;
100 
101   using CriterionValueType = double;
102   using PriorityType = std::pair< bool, CriterionValueType >;
103 
104   using PriorityQueueItemType = MaxPriorityQueueElementWrapper<
105     OutputEdgeCellType *, PriorityType, long >;
106 
107   using PriorityQueueType = PriorityQueueContainer< PriorityQueueItemType *,
108                                   ElementWrapperPointerInterface< PriorityQueueItemType * >,
109                                   PriorityType,
110                                   long >;
111 
112   using PriorityQueuePointer = typename PriorityQueueType::Pointer;
113   using QueueMapType = std::map< OutputEdgeCellType *, PriorityQueueItemType * >;
114   using QueueMapIterator = typename QueueMapType::iterator;
115 
116   using FlipEdgeFunctionType = QuadEdgeMeshEulerOperatorFlipEdgeFunction<
117     OutputMeshType, OutputQEType >;
118   using FlipEdgeFunctionPointer = typename FlipEdgeFunctionType::Pointer;
119 
SetListOfConstrainedEdges(const OutputEdgeCellListType & iList)120   void SetListOfConstrainedEdges(const OutputEdgeCellListType & iList)
121   {
122     m_ListOfConstrainedEdges = iList;
123   }
124 
125 protected:
126   DelaunayConformingQuadEdgeMeshFilter();
127   ~DelaunayConformingQuadEdgeMeshFilter() override;
128   void PrintSelf(std::ostream & os, Indent indent) const override;
129 
130   OutputEdgeCellListType m_ListOfConstrainedEdges;
131   PriorityQueuePointer   m_PriorityQueue;
132   QueueMapType           m_QueueMapper;
133 
134   SizeValueType           m_NumberOfEdgeFlips;
135   FlipEdgeFunctionPointer m_FlipEdge;
136 
137   void GenerateData() override;
138 
139   void InitializePriorityQueue();
140 
141   void Process();
142 
143   inline CriterionValueType
Dyer07Criterion(OutputMeshType * iMesh,OutputQEType * iEdge)144   Dyer07Criterion(OutputMeshType *iMesh, OutputQEType *iEdge) const
145   {
146     OutputPointIdentifier id1 = iEdge->GetOrigin();
147     OutputPointIdentifier id2 = iEdge->GetDestination();
148 
149     OutputPointIdentifier idA = iEdge->GetLnext()->GetDestination();
150     OutputPointIdentifier idB = iEdge->GetRnext()->GetOrigin();
151 
152     OutputPointType pt1 = iMesh->GetPoint(id1);
153     OutputPointType pt2 = iMesh->GetPoint(id2);
154     OutputPointType ptA = iMesh->GetPoint(idA);
155     OutputPointType ptB = iMesh->GetPoint(idB);
156 
157     OutputVectorType v1A = ptA - pt1;
158     OutputVectorType v1B = ptB - pt1;
159     OutputVectorType v2A = ptA - pt2;
160     OutputVectorType v2B = ptB - pt2;
161 
162     OutputCoordRepType sq_norm1A = v1A * v1A;
163     OutputCoordRepType sq_norm1B = v1B * v1B;
164     OutputCoordRepType sq_norm2A = v2A * v2A;
165     OutputCoordRepType sq_norm2B = v2B * v2B;
166 
167     auto dotA = static_cast< CriterionValueType >( v1A * v2A );
168     auto dotB = static_cast< CriterionValueType >( v1B * v2B );
169     auto den  = static_cast< CriterionValueType >( sq_norm1A * sq_norm2A );
170 
171     if ( den != 0. )
172       {
173       dotA /= std::sqrt(den);
174       }
175 
176     if ( dotA > 1. )
177       {
178       dotA = 1.;
179       }
180 
181     if ( dotA < -1. )
182       {
183       dotA = -1.;
184       }
185 
186     den = static_cast< CriterionValueType >( sq_norm1B * sq_norm2B );
187 
188     if ( den != 0. )
189       {
190       dotB /= std::sqrt(den);
191       }
192 
193     if ( dotB > 1. )
194       {
195       dotB = 1.;
196       }
197 
198     if ( dotB < -1. )
199       {
200       dotB = -1.;
201       }
202 
203     return ( std::acos(dotA) + std::acos(dotB) - itk::Math::pi );
204   }
205 
206 private:
207 
208   DelaunayConformingQuadEdgeMeshFilter(const Self &) = delete;
209   void operator=(const Self &) = delete;
210 };
211 } // end namespace itk
212 
213 #include "itkDelaunayConformingQuadEdgeMeshFilter.hxx"
214 
215 #endif
216