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 /*=========================================================================
19  *
20  *  Portions of this file are subject to the VTK Toolkit Version 3 copyright.
21  *
22  *  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
23  *
24  *  For complete copyright, license and disclaimer of warranty information
25  *  please refer to the NOTICE file at the top of the ITK source tree.
26  *
27  *=========================================================================*/
28 #ifndef itkSimplexMeshAdaptTopologyFilter_hxx
29 #define itkSimplexMeshAdaptTopologyFilter_hxx
30 
31 #include "itkSimplexMeshAdaptTopologyFilter.h"
32 
33 namespace itk
34 {
35 template< typename TInputMesh, typename TOutputMesh >
SimplexMeshAdaptTopologyFilter()36 SimplexMeshAdaptTopologyFilter< TInputMesh, TOutputMesh >::SimplexMeshAdaptTopologyFilter() :
37   m_IdOffset(0),
38   m_Output(TOutputMesh::New())
39 {
40   this->ProcessObject::SetNumberOfRequiredOutputs(1);
41   this->ProcessObject::SetNthOutput( 0, m_Output.GetPointer() );
42 }
43 
44 template< typename TInputMesh, typename TOutputMesh >
45 void SimplexMeshAdaptTopologyFilter< TInputMesh, TOutputMesh >
GenerateData()46 ::GenerateData()
47 {
48   this->Initialize();
49   this->ComputeCellParameters();
50 }
51 
52 //
53 template< typename TInputMesh, typename TOutputMesh >
54 void SimplexMeshAdaptTopologyFilter< TInputMesh, TOutputMesh >
Initialize()55 ::Initialize()
56 {
57   m_ModifiedCount = 0;
58   this->CopyInputMeshToOutputMeshPoints();
59   this->CopyInputMeshToOutputMeshPointData();
60   this->CopyInputMeshToOutputMeshCellData();
61   this->CopyInputMeshToOutputMeshCells();
62   this->CopyInputMeshToOutputMeshGeometryData();
63 }
64 
65 
66 template< typename TInputMesh, typename TOutputMesh >
67 void SimplexMeshAdaptTopologyFilter< TInputMesh, TOutputMesh >
CopyInputMeshToOutputMeshGeometryData()68 ::CopyInputMeshToOutputMeshGeometryData()
69 {
70   const InputMeshType *inputMesh   =  this->GetInput();
71   OutputMeshPointer    outputMesh  =  this->GetOutput();
72 
73   const PointIdentifier numberOfPoints = inputMesh->GetNumberOfPoints();
74 
75   using GeometryMapType = typename InputMeshType::GeometryMapType;
76   using GeometryMapPointer = typename InputMeshType::GeometryMapPointer;
77   using GeometryMapConstIterator = typename InputMeshType::GeometryMapConstIterator;
78 
79   GeometryMapPointer inputGeometryData = inputMesh->GetGeometryData();
80 
81   GeometryMapPointer outputGeometryData = GeometryMapType::New();
82 
83   outputGeometryData->Reserve( numberOfPoints );
84 
85   GeometryMapConstIterator inputGeometryItr = inputGeometryData->Begin();
86 
87   for( PointIdentifier pointId = 0; pointId < numberOfPoints; pointId++ )
88     {
89     auto * outputGeometryDataItem = new SimplexMeshGeometry;
90     SimplexMeshGeometry * inputGeometryDataItem = inputGeometryItr.Value();
91     outputGeometryDataItem->CopyFrom( *inputGeometryDataItem );
92     outputGeometryData->InsertElement( pointId, outputGeometryDataItem );
93     ++inputGeometryItr;
94     }
95 
96   outputMesh->SetGeometryData( outputGeometryData );
97   outputMesh->SetLastCellId( inputMesh->GetLastCellId() );
98 }
99 
100 
101 template< typename TInputMesh, typename TOutputMesh >
102 void SimplexMeshAdaptTopologyFilter< TInputMesh, TOutputMesh >
ComputeCellParameters()103 ::ComputeCellParameters()
104 {
105   OutputMeshPointer    outputMesh  =  this->GetOutput();
106 
107   // Ensure that cells will be deallocated by the Mesh.
108   outputMesh->SetCellsAllocationMethod( TInputMesh::CellsAllocatedDynamicallyCellByCell );
109 
110   SimplexVisitorInterfacePointer simplexVisitor = SimplexVisitorInterfaceType::New();
111   simplexVisitor->mesh = outputMesh;
112   CellMultiVisitorPointer mv = CellMultiVisitorType::New();
113   mv->AddVisitor(simplexVisitor);
114   outputMesh->Accept(mv);
115 
116   typename DoubleValueMapType::Pointer areas = simplexVisitor->GetAreaMap();
117   DoubleContainerIterator     areaIt = areas->Begin();
118   typename DoubleValueMapType::Pointer curvatures = simplexVisitor->GetCurvatureMap();
119   DoubleContainerIterator     curvatureIt = curvatures->Begin();
120 
121   double averageCurvature = simplexVisitor->GetTotalMeanCurvature();
122 
123   double rangeCellSize = simplexVisitor->GetMaximumCellSize() - simplexVisitor->GetMinimumCellSize();
124   double rangeCurvature  = simplexVisitor->GetMaximumCurvature() - simplexVisitor->GetMinimumCurvature();
125 
126   while ( curvatureIt != curvatures->End() )
127     {
128     bool doRefinement = false;
129 
130     const bool conditionA1 = ( m_SelectionMethod == 0 );
131     const bool conditionA2 = ( curvatureIt.Value() > averageCurvature );
132 
133     const double limit1 =        0.05 * rangeCellSize + simplexVisitor->GetMinimumCellSize();
134     const double limit2 = m_Threshold * rangeCellSize + simplexVisitor->GetMinimumCellSize();
135 
136     const bool conditionA3 = ( areaIt.Value() > limit1 );
137     const bool conditionA4 = ( areaIt.Value() > limit2 );
138 
139     if ( conditionA1 && ( ( conditionA2 && conditionA3 ) || conditionA4 ) )
140       {
141       doRefinement = true;
142       }
143     else
144       {
145       const bool conditionB1 = ( m_SelectionMethod == 1 );
146       const bool conditionB2 = curvatureIt.Value() > m_Threshold * rangeCurvature;
147       const bool conditionB3 = areaIt.Value() > 0.05 * rangeCellSize;
148       const bool conditionB4 = areaIt.Value() > m_Threshold * rangeCellSize;
149 
150       if ( conditionB1  && ( ( conditionB2 && conditionB3  ) || conditionB4 ) )
151         {
152         doRefinement = true;
153         }
154       }
155 
156     if ( doRefinement )
157       {
158       m_ModifiedCount++;
159 
160       InputCellAutoPointer poly;
161       outputMesh->GetCell(curvatureIt.Index(), poly);
162 
163       InputPointType cellCenter = this->ComputeCellCenter(poly);
164 
165       typename InputPolygonType::PointIdIterator pointIds = poly->PointIdsBegin();
166 
167       PointIdentifier lineOneFirstIdx = *pointIds;
168       pointIds++;
169       PointIdentifier lineOneSecondIdx = *pointIds;
170 
171       unsigned short cnt = 0;
172 
173       while ( cnt < poly->GetNumberOfPoints() / 2 - 1 )
174         {
175         pointIds++;
176         cnt++;
177         }
178       PointIdentifier lineTwoFirstIdx = *pointIds;
179       pointIds++;
180       PointIdentifier lineTwoSecondIdx = *pointIds;
181 
182       PointIdentifier newPointId = outputMesh->GetNumberOfPoints();
183       PointIdentifier firstNewIndex = newPointId;
184       PointIdentifier secondNewIndex = newPointId + 1;
185 
186       //create first new point
187       InputPointType newMidPoint, helperPoint;
188       InputPointType p1, p2;
189       p1.Fill(0);
190       p2.Fill(0);
191       outputMesh->GetPoint(lineOneFirstIdx, &p1);
192       outputMesh->GetPoint(lineOneSecondIdx, &p2);
193 
194       helperPoint.SetToMidPoint(p1, p2);
195       newMidPoint.SetToMidPoint(helperPoint, cellCenter);
196 
197       outputMesh->SetPoint(firstNewIndex, newMidPoint);
198       outputMesh->SetGeometryData( firstNewIndex, new itk::SimplexMeshGeometry() );
199 
200       outputMesh->ReplaceNeighbor(lineOneFirstIdx, lineOneSecondIdx, firstNewIndex);
201       outputMesh->ReplaceNeighbor(lineOneSecondIdx, lineOneFirstIdx, firstNewIndex);
202 
203       //create second new point
204       outputMesh->GetPoint(lineTwoFirstIdx, &p1);
205       outputMesh->GetPoint(lineTwoSecondIdx, &p2);
206 
207       helperPoint.SetToMidPoint(p1, p2);
208       newMidPoint.SetToMidPoint(helperPoint, cellCenter);
209 
210       outputMesh->SetPoint(secondNewIndex, newMidPoint);
211       outputMesh->SetGeometryData( secondNewIndex, new itk::SimplexMeshGeometry() );
212 
213       outputMesh->ReplaceNeighbor(lineTwoFirstIdx, lineTwoSecondIdx, secondNewIndex);
214       outputMesh->ReplaceNeighbor(lineTwoSecondIdx, lineTwoFirstIdx, secondNewIndex);
215 
216       outputMesh->AddNeighbor(firstNewIndex, secondNewIndex);
217       outputMesh->AddNeighbor(firstNewIndex, lineOneFirstIdx);
218       outputMesh->AddNeighbor(firstNewIndex, lineOneSecondIdx);
219 
220       outputMesh->AddNeighbor(secondNewIndex, lineTwoSecondIdx);
221       outputMesh->AddNeighbor(secondNewIndex, firstNewIndex);
222       outputMesh->AddNeighbor(secondNewIndex, lineTwoFirstIdx);
223 
224       CovariantVectorType lineOneFirstNormal = outputMesh->ComputeNormal(lineOneFirstIdx);
225       CovariantVectorType firstNewNormal = outputMesh->ComputeNormal(firstNewIndex);
226 
227       CovariantVectorType lineTwoFirstNormal = outputMesh->ComputeNormal(lineTwoFirstIdx);
228       CovariantVectorType secondNewNormal = outputMesh->ComputeNormal(secondNewIndex);
229 
230       double prod;
231 
232       prod = dot_product( firstNewNormal.GetVnlVector(), lineOneFirstNormal.GetVnlVector() );
233 
234       if ( prod < 0 )
235         {
236         outputMesh->SwapNeighbors(firstNewIndex, lineOneFirstIdx, lineOneSecondIdx);
237         firstNewNormal = outputMesh->ComputeNormal(firstNewIndex);
238         }
239 
240       prod = dot_product( secondNewNormal.GetVnlVector(), lineTwoFirstNormal.GetVnlVector() );
241       if ( prod < 0 )
242         {
243         outputMesh->SwapNeighbors(secondNewIndex, lineTwoFirstIdx, lineTwoSecondIdx);
244         secondNewNormal = outputMesh->ComputeNormal(secondNewIndex);
245         }
246 
247       outputMesh->AddEdge(firstNewIndex, secondNewIndex);
248 
249       // splitting cell
250       PointIdentifier    newPointIndex = NumericTraits< PointIdentifier >::ZeroValue();
251       auto * polygon = new OutputPolygonType;
252       InputCellAutoPointer newPolygonPointer1;
253       newPolygonPointer1.TakeOwnership(polygon);
254 
255       pointIds = poly->PointIdsBegin();
256 
257       PointIdentifier firstPointId = *pointIds++;
258 
259       while ( *pointIds != lineTwoSecondIdx )
260         {
261         newPolygonPointer1->SetPointId(newPointIndex++, *pointIds++);
262         }
263 
264       newPolygonPointer1->SetPointId(newPointIndex++, secondNewIndex);
265       newPolygonPointer1->SetPointId(newPointIndex++, firstNewIndex);
266 
267       auto * polygon2 = new OutputPolygonType;
268       newPointIndex = 0;
269 
270       while ( pointIds != poly->PointIdsEnd() )
271         {
272         polygon2->SetPointId(newPointIndex++, *pointIds++);
273         }
274 
275       InputCellAutoPointer newPolygonPointer2;
276       newPolygonPointer2.TakeOwnership( polygon2 );
277 
278       newPolygonPointer2->SetPointId(newPointIndex++, firstPointId);
279       newPolygonPointer2->SetPointId(newPointIndex++, firstNewIndex);
280       newPolygonPointer2->SetPointId(newPointIndex++, secondNewIndex);
281 
282       outputMesh->ReplaceFace(curvatureIt.Index(), newPolygonPointer1);
283       outputMesh->AddFace(newPolygonPointer2);
284 
285       outputMesh->BuildCellLinks();
286 
287       this->ModifyNeighborCells(lineOneFirstIdx, lineOneSecondIdx, firstNewIndex);
288       this->ModifyNeighborCells(lineTwoFirstIdx, lineTwoSecondIdx, secondNewIndex);
289 
290       } // end if cell must be modified
291     ++areaIt;
292     ++curvatureIt;
293     }
294 }
295 
296 
297 template< typename TInputMesh, typename TOutputMesh >
298 void
299 SimplexMeshAdaptTopologyFilter< TInputMesh, TOutputMesh >
ModifyNeighborCells(CellIdentifier id1,CellIdentifier id2,PointIdentifier insertPointId)300 ::ModifyNeighborCells(CellIdentifier id1, CellIdentifier id2, PointIdentifier insertPointId)
301 {
302   OutputMeshPointer    outputMesh  =  this->GetOutput();
303 
304   std::set< PointIdentifier >           cells1 =   outputMesh->GetCellLinks()->GetElement(id1);
305   std::set< PointIdentifier >           cells2 =   outputMesh->GetCellLinks()->GetElement(id2);
306   auto cellIt = cells1.begin();
307 
308   std::set< PointIdentifier > result;
309 
310   while ( cellIt != cells1.end() )
311     {
312     auto found = std::find(cells2.begin(), cells2.end(), *cellIt);
313     if ( found != cells2.end() )
314       {
315       result.insert(*cellIt);
316       }
317     cellIt++;
318     }
319 
320   cellIt = result.begin();
321 
322   while ( cellIt != result.end() )
323     {
324     InputCellAutoPointer nextCell;
325     outputMesh->GetCell(*cellIt, nextCell);
326 
327     if ( nextCell->GetNumberOfPoints() == 2 )
328       {
329       InputCellPointIdIterator lineIt =  nextCell->PointIdsBegin();
330       PointIdentifier          first = *lineIt++;
331       PointIdentifier          second = *lineIt;
332 
333       outputMesh->AddEdge(first, insertPointId);
334       outputMesh->AddEdge(insertPointId, second);
335 
336       // Take over the cell and release its memory
337       outputMesh->GetCells()->DeleteIndex(*cellIt);
338       nextCell.TakeOwnership();
339       nextCell.Reset();
340       }
341     else if ( nextCell->GetNumberOfPoints() > 3 )
342       {
343       m_NewSimplexCellPointer.TakeOwnership(new OutputPolygonType);
344       InputPolygonPointIdIterator             pointIt =  nextCell->PointIdsBegin();
345       PointIdentifier cnt = NumericTraits< PointIdentifier >::ZeroValue();
346       PointIdentifier first = *pointIt++;
347       PointIdentifier startId = first;
348 
349       PointIdentifier second = 0;
350 
351       while ( pointIt != nextCell->PointIdsEnd() )
352         {
353         m_NewSimplexCellPointer->SetPointId(cnt++,  first);
354         second = *pointIt;
355 
356         if ( ( id1 == first && id2 == second ) || ( id2 == first && id1 == second ) )
357           {
358           m_NewSimplexCellPointer->SetPointId(cnt++,  insertPointId);
359           }
360         first = second;
361         pointIt++;
362         }
363 
364       m_NewSimplexCellPointer->SetPointId(cnt++,  second);
365       if ( ( id1 == second && id2 == startId ) || ( id2 ==  second && id1 == startId ) )
366         {
367         m_NewSimplexCellPointer->SetPointId(cnt++,  insertPointId);
368         }
369 
370       outputMesh->ReplaceFace(*cellIt, m_NewSimplexCellPointer);
371       }
372     cellIt++;
373     }
374 
375   outputMesh->BuildCellLinks();
376 }
377 
378 /* PrintSelf. */
379 template< typename TInputMesh, typename TOutputMesh >
380 void
381 SimplexMeshAdaptTopologyFilter< TInputMesh, TOutputMesh >
PrintSelf(std::ostream & os,Indent indent) const382 ::PrintSelf(std::ostream & os, Indent indent) const
383 {
384   Superclass::PrintSelf(os, indent);
385   os << indent << "Threshold: " << m_Threshold << std::endl;
386   os << indent << "SelectionMethod: " << m_SelectionMethod << std::endl;
387   os << indent << "ModifiedCount: " << m_ModifiedCount << std::endl;
388 }
389 
390 template< typename TInputMesh, typename TOutputMesh >
391 typename SimplexMeshAdaptTopologyFilter< TInputMesh, TOutputMesh >::InputPointType
392 SimplexMeshAdaptTopologyFilter< TInputMesh, TOutputMesh >
ComputeCellCenter(InputCellAutoPointer & simplexCell)393 ::ComputeCellCenter(InputCellAutoPointer & simplexCell)
394 {
395   OutputMeshPointer    outputMesh  =  this->GetOutput();
396   InputPolygonPointIdIterator pointIt =  simplexCell->PointIdsBegin();
397 
398   InputVectorType tmp;
399   InputPointType  p1, cellCenter;
400 
401   p1.Fill(0);
402   cellCenter.Fill(0);
403 
404   // compute the cell center first
405   while ( pointIt != simplexCell->PointIdsEnd() )
406     {
407     outputMesh->GetPoint(*pointIt, &p1);
408     cellCenter += p1.GetVectorFromOrigin();
409     pointIt++;
410     }
411 
412   tmp.SetVnlVector( cellCenter.GetVnlVector() / simplexCell->GetNumberOfPoints() );
413   cellCenter.Fill(0.0);
414   cellCenter += tmp;
415 
416   return cellCenter;
417 }
418 } // end of namspace itk
419 
420 #endif // itkSimplexMeshAdaptTopologyFilter_hxx
421