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