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 itkQuadEdgeMeshEulerOperatorCreateCenterVertexFunction_hxx
19 #define itkQuadEdgeMeshEulerOperatorCreateCenterVertexFunction_hxx
20
21 #include "itkQuadEdgeMeshEulerOperatorCreateCenterVertexFunction.h"
22
23 namespace itk
24 {
25 template< typename TMesh, typename TQEType >
26 typename QuadEdgeMeshEulerOperatorCreateCenterVertexFunction< TMesh, TQEType >::OutputType
Evaluate(QEType * e)27 QuadEdgeMeshEulerOperatorCreateCenterVertexFunction< TMesh, TQEType >::Evaluate(QEType *e)
28 {
29 // Is there any input ?
30 #ifndef NDEBUG
31 if ( !e )
32 {
33 itkDebugMacro("Input is not an edge.");
34 return ( (QEType *)nullptr );
35 }
36
37 if ( !this->m_Mesh )
38 {
39 itkDebugMacro("No mesh present.");
40 return ( (OutputType)nullptr );
41 }
42
43 // Is left face set ?
44 if ( !e->IsLeftSet() )
45 {
46 itkDebugMacro("Argument edge has no left face.");
47 return ( (OutputType)nullptr );
48 }
49 #endif
50
51 // remove left face
52 this->m_Mesh->DeleteFace( e->GetLeft() );
53
54 // create new point geometry
55 unsigned int sum = 0;
56 VectorType vec;
57 vec.Fill(0);
58 PointIdentifier pid = this->m_Mesh->FindFirstUnusedPointIndex();
59 using AssociatedBarycenters = std::map< QEType *, PointIdentifier >;
60 AssociatedBarycenters m_AssocBary;
61 using QEIterator = typename QEType::IteratorGeom;
62 QEIterator lit = e->BeginGeomLnext();
63 while ( lit != e->EndGeomLnext() )
64 {
65 QEType *g = lit.Value();
66 vec += this->m_Mesh->GetVector( g->GetOrigin() );
67 sum++;
68 m_AssocBary[g] = pid;
69 lit++;
70 } // rof
71 vec /= CoordRepType(sum);
72 PointType p;
73 for ( unsigned int i = 0; i < 3; i++ )
74 {
75 p[i] = vec[i];
76 }
77
78 // add new point to mesh
79 this->m_NewPointID = this->m_Mesh->AddPoint(p);
80 PointIdentifier tempPoint;
81
82 // create edges and faces
83 tempPoint = e->GetDestination();
84 this->m_Mesh->AddFaceTriangle(this->m_NewPointID,
85 e->GetOrigin(),
86 tempPoint);
87 QEType *edgeRef = this->m_Mesh->FindEdge(this->m_NewPointID, tempPoint);
88 while ( !edgeRef->IsLeftSet() )
89 {
90 tempPoint = edgeRef->GetLnext()->GetDestination();
91 this->m_Mesh->AddFaceTriangle(this->m_NewPointID,
92 edgeRef->GetLnext()->GetOrigin(),
93 tempPoint);
94 edgeRef = this->m_Mesh->FindEdge(this->m_NewPointID, tempPoint);
95 }
96
97 return ( e->GetLnext() );
98 }
99
100 } // end namespace itk
101
102 #endif
103