1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtk3DLinearGridCrinkleExtractor.cxx
5
6 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7 All rights reserved.
8 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9
10 This software is distributed WITHOUT ANY WARRANTY; without even
11 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12 PURPOSE. See the above copyright notice for more information.
13
14 =========================================================================*/
15
16 // Hide VTK_DEPRECATED_IN_9_0_0() warnings for this class.
17 #define VTK_DEPRECATION_LEVEL 0
18
19 #include "vtk3DLinearGridCrinkleExtractor.h"
20
21 #include "vtk3DLinearGridInternal.h"
22 #include "vtkArrayListTemplate.h" // For processing attribute data
23 #include "vtkCellArray.h"
24 #include "vtkCellData.h"
25 #include "vtkCellTypes.h"
26 #include "vtkCompositeDataIterator.h"
27 #include "vtkCompositeDataSet.h"
28 #include "vtkImplicitFunction.h"
29 #include "vtkInformation.h"
30 #include "vtkInformationVector.h"
31 #include "vtkLogger.h"
32 #include "vtkMultiBlockDataSet.h"
33 #include "vtkObjectFactory.h"
34 #include "vtkPlane.h"
35 #include "vtkPointData.h"
36 #include "vtkSMPThreadLocalObject.h"
37 #include "vtkSMPTools.h"
38 #include "vtkStaticCellLinksTemplate.h"
39 #include "vtkStreamingDemandDrivenPipeline.h"
40 #include "vtkUnsignedCharArray.h"
41 #include "vtkUnstructuredGrid.h"
42
43 #include <atomic>
44 #include <vector>
45
46 vtkStandardNewMacro(vtk3DLinearGridCrinkleExtractor);
47 vtkCxxSetObjectMacro(vtk3DLinearGridCrinkleExtractor, ImplicitFunction, vtkImplicitFunction);
48
49 //------------------------------------------------------------------------------
50 // Macros immediately below are just used to make code easier to
51 // read. Invokes functor _op _num times depending on serial (_seq==1) or
52 // parallel processing mode. The _REDUCE_ version is used to called functors
53 // with a Reduce() method).
54 #define EXECUTE_SMPFOR(_seq, _num, _op) \
55 if (!_seq) \
56 { \
57 vtkSMPTools::For(0, _num, _op); \
58 } \
59 else \
60 { \
61 _op(0, _num); \
62 }
63
64 #define EXECUTE_REDUCED_SMPFOR(_seq, _num, _op, _nt) \
65 if (!_seq) \
66 { \
67 vtkSMPTools::For(0, _num, _op); \
68 } \
69 else \
70 { \
71 _op.Initialize(); \
72 _op(0, _num); \
73 _op.Reduce(); \
74 } \
75 _nt = _op.NumThreadsUsed;
76
77 namespace
78 { // anonymous
79
80 //========================= Quick implicit function cell selection =============
81
82 // Compute an array that classifies each point with respect to the current
83 // implicit function (i.e. above the function(=2), below the function(=1), on
84 // the function(=0)). InOutArray is allocated here and should be deleted by
85 // the invoking code. InOutArray is an unsigned char array to simplify bit
86 // fiddling later on (i.e., Intersects() method). A fast path is available
87 // for vtkPlane implicit functions.
88 //
89 // The reason we compute this unsigned char array as compared to an array of
90 // function values is to reduce the amount of memory used, and writing to
91 // memory, since these are significant costs for large data.
92
93 // Templated for explicit point representations of real type
94 template <typename TP>
95 struct PlaneClassifyPoints;
96 template <typename TP>
97 struct FunctionClassifyPoints;
98
99 // General classification capability
100 struct Classify
101 {
102 unsigned char* InOutArray;
103
Classify__anon2aa6e4590111::Classify104 Classify(vtkPoints* pts) { this->InOutArray = new unsigned char[pts->GetNumberOfPoints()]; }
105
106 // Check if a list of points intersects the plane
Intersects__anon2aa6e4590111::Classify107 static bool Intersects(const unsigned char* inout, vtkIdType npts, const vtkIdType* pts)
108 {
109 unsigned char onOneSideOfPlane = inout[pts[0]];
110 for (vtkIdType i = 1; onOneSideOfPlane && i < npts; ++i)
111 {
112 onOneSideOfPlane &= inout[pts[i]];
113 }
114 return (!onOneSideOfPlane);
115 }
116 };
117
118 // Faster path for vtkPlane
119 template <typename TP>
120 struct PlaneClassifyPoints : public Classify
121 {
122 TP* Points;
123 double Origin[3];
124 double Normal[3];
125
PlaneClassifyPoints__anon2aa6e4590111::PlaneClassifyPoints126 PlaneClassifyPoints(vtkPoints* pts, vtkPlane* plane)
127 : Classify(pts)
128 {
129 this->Points = static_cast<TP*>(pts->GetVoidPointer(0));
130 plane->GetOrigin(this->Origin);
131 plane->GetNormal(this->Normal);
132 }
133
operator ()__anon2aa6e4590111::PlaneClassifyPoints134 void operator()(vtkIdType ptId, vtkIdType endPtId)
135 {
136 double p[3], zero = double(0), eval;
137 double *n = this->Normal, *o = this->Origin;
138 TP* pts = this->Points + 3 * ptId;
139 unsigned char* ioa = this->InOutArray + ptId;
140 for (; ptId < endPtId; ++ptId)
141 {
142 // Access each point
143 p[0] = static_cast<double>(*pts);
144 ++pts;
145 p[1] = static_cast<double>(*pts);
146 ++pts;
147 p[2] = static_cast<double>(*pts);
148 ++pts;
149
150 // Evaluate position of the point with the plane. Invoke inline,
151 // non-virtual version of evaluate method.
152 eval = vtkPlane::Evaluate(n, o, p);
153
154 // Point is either above(=2), below(=1), or on(=0) the plane.
155 *ioa++ = (eval > zero ? 2 : (eval < zero ? 1 : 0));
156 }
157 }
158 };
159
160 // General path for vtkImplicitFunction
161 template <typename TP>
162 struct FunctionClassifyPoints : public Classify
163 {
164 TP* Points;
165 vtkImplicitFunction* Function;
166
FunctionClassifyPoints__anon2aa6e4590111::FunctionClassifyPoints167 FunctionClassifyPoints(vtkPoints* pts, vtkImplicitFunction* f)
168 : Classify(pts)
169 , Function(f)
170 {
171 this->Points = static_cast<TP*>(pts->GetVoidPointer(0));
172 }
173
operator ()__anon2aa6e4590111::FunctionClassifyPoints174 void operator()(vtkIdType ptId, vtkIdType endPtId)
175 {
176 double p[3], zero = double(0), eval;
177 TP* pts = this->Points + 3 * ptId;
178 unsigned char* ioa = this->InOutArray + ptId;
179 for (; ptId < endPtId; ++ptId)
180 {
181 // Access each point
182 p[0] = static_cast<double>(*pts);
183 ++pts;
184 p[1] = static_cast<double>(*pts);
185 ++pts;
186 p[2] = static_cast<double>(*pts);
187 ++pts;
188
189 // Evaluate position of the point wrt the implicit function. This
190 // call better be thread safe.
191 eval = this->Function->FunctionValue(p);
192
193 // Point is either above(=2), below(=1), or on(=0) the plane.
194 *ioa++ = (eval > zero ? 2 : (eval < zero ? 1 : 0));
195 }
196 }
197 };
198
199 // Base class for extracting cells and points from the input
200 // vtkUnstructuredGrid.
201 struct ExtractCellsBase
202 {
203 typedef std::vector<vtkIdType> CellArrayType;
204 typedef std::vector<vtkIdType> OriginCellType;
205 typedef std::vector<unsigned char> CellTypesType;
206
207 // Track local data on a per-thread basis. In the Reduce() method this
208 // information will be used to composite the data from each thread.
209 struct LocalDataType
210 {
211 CellArrayType LocalCells;
212 OriginCellType LocalOrigins;
213 CellTypesType LocalTypes;
214 vtkIdType LocalNumCells;
215 CellIter LocalCellIter;
216
LocalDataType__anon2aa6e4590111::ExtractCellsBase::LocalDataType217 LocalDataType()
218 : LocalNumCells(0)
219 {
220 }
221 };
222
223 const unsigned char* InOut;
224 CellIter* Iter;
225 vtkIdType InputNumPts;
226 vtkIdType OutputNumPts;
227 vtkIdType OutputNumCells;
228 vtkIdType TotalSize;
229 vtkUnstructuredGrid* Grid;
230 vtkCellArray* Cells;
231 bool CopyPointData;
232 bool CopyCellData;
233 vtkIdType* PointMap;
234 vtkIdType* CellMap;
235 int NumThreadsUsed;
236 vtkSMPThreadLocal<LocalDataType> LocalData;
237
ExtractCellsBase__anon2aa6e4590111::ExtractCellsBase238 ExtractCellsBase(vtkIdType inNumPts, CellIter* c, unsigned char* inout, vtkUnstructuredGrid* grid,
239 vtkCellArray* cells, bool copyPtData, bool copyCellData)
240 : InOut(inout)
241 , Iter(c)
242 , InputNumPts(inNumPts)
243 , OutputNumPts(0)
244 , OutputNumCells(0)
245 , TotalSize(0)
246 , Grid(grid)
247 , Cells(cells)
248 , CopyPointData(copyPtData)
249 , CopyCellData(copyCellData)
250 , PointMap(nullptr)
251 , CellMap(nullptr)
252 , NumThreadsUsed(0)
253 {
254 }
255
256 // Set up the iteration process
Initialize__anon2aa6e4590111::ExtractCellsBase257 void Initialize()
258 {
259 auto& localData = this->LocalData.Local();
260 localData.LocalCellIter = *(this->Iter);
261 }
262
263 }; // ExtractCellsBase
264
265 // Traverse all cells and extract intersected cells
266 struct ExtractCells : public ExtractCellsBase
267 {
ExtractCells__anon2aa6e4590111::ExtractCells268 ExtractCells(vtkIdType inNumPts, CellIter* c, unsigned char* inout, vtkUnstructuredGrid* grid,
269 vtkCellArray* cells, bool copyPtData, bool copyCellData)
270 : ExtractCellsBase(inNumPts, c, inout, grid, cells, copyPtData, copyCellData)
271 {
272 }
273
Initialize__anon2aa6e4590111::ExtractCells274 void Initialize() { this->ExtractCellsBase::Initialize(); }
275
276 // operator() method extracts cells
operator ()__anon2aa6e4590111::ExtractCells277 void operator()(vtkIdType cellId, vtkIdType endCellId)
278 {
279 auto& localData = this->LocalData.Local();
280 auto& lCells = localData.LocalCells;
281 auto& lOrigins = localData.LocalOrigins;
282 auto& lTypes = localData.LocalTypes;
283 CellIter* cellIter = &localData.LocalCellIter;
284 const vtkIdType* c = cellIter->Initialize(cellId); // connectivity array
285 const unsigned char* inout = this->InOut;
286 vtkIdType& lNumCells = localData.LocalNumCells;
287 vtkIdType npts;
288
289 for (; cellId < endCellId; ++cellId)
290 {
291 // Does the implicit function cut this cell?
292 npts = cellIter->NumVerts;
293 if (Classify::Intersects(inout, npts, c))
294 {
295 ++lNumCells;
296 lTypes.emplace_back(cellIter->GetCellType(cellId));
297 lCells.emplace_back(npts);
298 const vtkIdType* pts = cellIter->GetCellIds(cellId);
299 for (auto i = 0; i < npts; ++i)
300 {
301 lCells.emplace_back(pts[i]);
302 }
303 if (this->CopyCellData)
304 {
305 lOrigins.emplace_back(cellId); // to support cell data copying
306 }
307 } // if implicit function intersects
308 c = cellIter->Next(); // move to the next cell
309 } // for all cells in this batch
310 }
311
312 // Composite local thread data. Basically build the output unstructured grid.
Reduce__anon2aa6e4590111::ExtractCells313 void Reduce()
314 {
315 // Count the number of cells, and the number of threads used. Figure out the total
316 // length of the cell array.
317 vtkIdType numCells = 0;
318 vtkIdType size = 0;
319 for (const auto& threadData : this->LocalData)
320 {
321 numCells += threadData.LocalNumCells;
322 size += static_cast<vtkIdType>(threadData.LocalCells.size());
323 this->NumThreadsUsed++;
324 }
325 this->OutputNumCells = numCells;
326 this->TotalSize = size;
327
328 // Now allocate the cell array, offset array, and cell types array.
329 this->Cells->AllocateExact(numCells, size - numCells);
330 vtkNew<vtkUnsignedCharArray> cellTypes;
331 unsigned char* ctptr = static_cast<unsigned char*>(cellTypes->WriteVoidPointer(0, numCells));
332
333 // If cell data is requested, roll up generating cell ids
334 vtkIdType* cellMap = nullptr;
335 if (this->CopyCellData)
336 {
337 this->CellMap = cellMap = new vtkIdType[numCells];
338 }
339
340 // Now composite the cell-related information
341 for (const auto& threadData : this->LocalData)
342 {
343 const auto& lCells = threadData.LocalCells;
344 const auto& lTypes = threadData.LocalTypes;
345 const auto& lOrigins = threadData.LocalOrigins;
346 numCells = threadData.LocalNumCells;
347
348 this->Cells->AppendLegacyFormat(lCells.data(), static_cast<vtkIdType>(lCells.size()));
349 ctptr = std::copy_n(lTypes.cbegin(), numCells, ctptr);
350 if (this->CopyCellData)
351 {
352 cellMap = std::copy_n(lOrigins.cbegin(), numCells, cellMap);
353 }
354 } // for this thread
355
356 // Define the grid
357 this->Grid->SetCells(cellTypes, this->Cells);
358
359 } // Reduce
360
361 }; // ExtractCells
362
363 // Traverse all cells to extract intersected cells and remapped points
364 struct ExtractPointsAndCells : public ExtractCellsBase
365 {
ExtractPointsAndCells__anon2aa6e4590111::ExtractPointsAndCells366 ExtractPointsAndCells(vtkIdType inNumPts, CellIter* c, unsigned char* inout,
367 vtkUnstructuredGrid* grid, vtkCellArray* cells, bool copyPtData, bool copyCellData)
368 : ExtractCellsBase(inNumPts, c, inout, grid, cells, copyPtData, copyCellData)
369 {
370 this->PointMap = new vtkIdType[inNumPts];
371 std::fill_n(this->PointMap, inNumPts, (-1));
372 }
373
Initialize__anon2aa6e4590111::ExtractPointsAndCells374 void Initialize() { this->ExtractCellsBase::Initialize(); }
375
376 // operator() method identifies cells and points to extract
operator ()__anon2aa6e4590111::ExtractPointsAndCells377 void operator()(vtkIdType cellId, vtkIdType endCellId)
378 {
379 auto& localData = this->LocalData.Local();
380 auto& lCells = localData.LocalCells;
381 auto& lOrigins = localData.LocalOrigins;
382 auto& lTypes = localData.LocalTypes;
383 CellIter* cellIter = &localData.LocalCellIter;
384 const vtkIdType* c = cellIter->Initialize(cellId); // connectivity array
385 const unsigned char* inout = this->InOut;
386 vtkIdType& lNumCells = localData.LocalNumCells;
387 vtkIdType npts;
388 vtkIdType* pointMap = this->PointMap;
389
390 for (; cellId < endCellId; ++cellId)
391 {
392 // Does the implicit function cut this cell?
393 npts = cellIter->NumVerts;
394 if (Classify::Intersects(inout, npts, c))
395 {
396 ++lNumCells;
397 lTypes.emplace_back(cellIter->GetCellType(cellId));
398 lCells.emplace_back(npts);
399 const vtkIdType* pts = cellIter->GetCellIds(cellId);
400 for (auto i = 0; i < npts; ++i)
401 {
402 pointMap[pts[i]] = 1; // this point is used
403 lCells.emplace_back(pts[i]);
404 }
405 if (this->CopyCellData)
406 {
407 lOrigins.emplace_back(cellId); // to support cell data copying
408 }
409 } // if implicit function intersects
410 c = cellIter->Next(); // move to the next cell
411 } // for all cells in this batch
412 }
413
414 // Composite local thread data. Basically build the output unstructured grid.
Reduce__anon2aa6e4590111::ExtractPointsAndCells415 void Reduce()
416 {
417 // Generate point map
418 vtkIdType globalPtId = 0;
419 vtkIdType* ptMap = this->PointMap;
420 for (auto ptId = 0; ptId < this->InputNumPts; ++ptId)
421 {
422 if (this->PointMap[ptId] > 0)
423 {
424 ptMap[ptId] = globalPtId++;
425 }
426 }
427 this->OutputNumPts = globalPtId;
428
429 // Count the number of cells, and the number of threads used. Figure out the total
430 // length of the cell array.
431 vtkIdType numCells = 0;
432 vtkIdType size = 0;
433 for (const auto& threadData : this->LocalData)
434 {
435 numCells += threadData.LocalNumCells;
436 size += static_cast<vtkIdType>(threadData.LocalCells.size());
437 this->NumThreadsUsed++;
438 }
439 this->OutputNumCells = numCells;
440 this->TotalSize = size;
441
442 // Now allocate the cell array, offset array, and cell types array.
443 this->Cells->AllocateExact(numCells, size - numCells);
444 vtkNew<vtkUnsignedCharArray> cellTypes;
445 unsigned char* ctptr = static_cast<unsigned char*>(cellTypes->WriteVoidPointer(0, numCells));
446
447 // If cell data is requested, roll up generating cell ids
448 vtkIdType* cellMap = nullptr;
449 this->CellMap = cellMap = new vtkIdType[numCells];
450
451 // Now composite the cell-related information
452 for (const auto& threadData : this->LocalData)
453 {
454 const auto& lCells = threadData.LocalCells;
455 const auto& lTypes = threadData.LocalTypes;
456 const auto& lOrigins = threadData.LocalOrigins;
457 numCells = threadData.LocalNumCells;
458
459 ctptr = std::copy_n(lTypes.cbegin(), numCells, ctptr);
460 if (this->CopyCellData)
461 {
462 cellMap = std::copy_n(lOrigins.cbegin(), numCells, cellMap);
463 }
464
465 // Need to do this in a loop since the pointIds are mapped through ptMap:
466 auto threadCells = lCells.cbegin();
467 for (auto i = 0; i < numCells; ++i)
468 {
469 const vtkIdType npts = *threadCells++;
470 this->Cells->InsertNextCell(static_cast<int>(npts));
471 for (auto j = 0; j < npts; ++j)
472 {
473 this->Cells->InsertCellPoint(ptMap[*threadCells++]);
474 }
475 } // over all the cells in this thread
476 } // for this thread
477
478 // Define the grid
479 this->Grid->SetCells(cellTypes, this->Cells);
480
481 } // Reduce
482
483 }; // ExtractPointsAndCells
484
485 // Copy cell data from input to output
486 struct CopyCellAttributes
487 {
488 ArrayList* Arrays;
489 const vtkIdType* CellMap;
490
CopyCellAttributes__anon2aa6e4590111::CopyCellAttributes491 CopyCellAttributes(ArrayList* arrays, const vtkIdType* cellMap)
492 : Arrays(arrays)
493 , CellMap(cellMap)
494 {
495 }
496
operator ()__anon2aa6e4590111::CopyCellAttributes497 void operator()(vtkIdType cellId, vtkIdType endCellId)
498 {
499 for (; cellId < endCellId; ++cellId)
500 {
501 this->Arrays->Copy(this->CellMap[cellId], cellId);
502 }
503 }
504 };
505
506 // Generate point coordinates
507 template <typename TPIn, typename TPOut>
508 struct GeneratePoints
509 {
510 const TPIn* InPts;
511 const vtkIdType* PointMap;
512 TPOut* OutPts;
513
GeneratePoints__anon2aa6e4590111::GeneratePoints514 GeneratePoints(TPIn* inPts, vtkIdType* ptMap, TPOut* outPts)
515 : InPts(inPts)
516 , PointMap(ptMap)
517 , OutPts(outPts)
518 {
519 }
520
operator ()__anon2aa6e4590111::GeneratePoints521 void operator()(vtkIdType ptId, vtkIdType endPtId)
522 {
523 const TPIn* p = this->InPts + 3 * ptId;
524 const vtkIdType* ptMap = this->PointMap;
525 TPOut *outPts = this->OutPts, *x;
526
527 for (; ptId < endPtId; ++ptId, p += 3)
528 {
529 if (ptMap[ptId] >= 0)
530 {
531 x = outPts + 3 * ptMap[ptId];
532 *x++ = static_cast<TPOut>(p[0]);
533 *x++ = static_cast<TPOut>(p[1]);
534 *x = static_cast<TPOut>(p[2]);
535 }
536 }
537 }
538 };
539
540 // Copy point data from input to output.
541 struct CopyPointAttributes
542 {
543 ArrayList* Arrays;
544 const vtkIdType* PointMap;
545
CopyPointAttributes__anon2aa6e4590111::CopyPointAttributes546 CopyPointAttributes(ArrayList* arrays, const vtkIdType* ptMap)
547 : Arrays(arrays)
548 , PointMap(ptMap)
549 {
550 }
551
operator ()__anon2aa6e4590111::CopyPointAttributes552 void operator()(vtkIdType ptId, vtkIdType endPtId)
553 {
554 const vtkIdType* ptMap = this->PointMap;
555 for (; ptId < endPtId; ++ptId)
556 {
557 if (ptMap[ptId] >= 0)
558 {
559 this->Arrays->Copy(ptId, ptMap[ptId]);
560 }
561 }
562 }
563 };
564
565 } // anonymous namespace
566
567 //------------------------------------------------------------------------------
568 // Construct an instance of the class.
vtk3DLinearGridCrinkleExtractor()569 vtk3DLinearGridCrinkleExtractor::vtk3DLinearGridCrinkleExtractor()
570 {
571 this->ImplicitFunction = nullptr;
572 this->CopyPointData = true;
573 this->CopyCellData = false;
574 this->RemoveUnusedPoints = false;
575 this->OutputPointsPrecision = DEFAULT_PRECISION;
576 this->SequentialProcessing = false;
577 this->NumberOfThreadsUsed = 0;
578 }
579
580 //------------------------------------------------------------------------------
~vtk3DLinearGridCrinkleExtractor()581 vtk3DLinearGridCrinkleExtractor::~vtk3DLinearGridCrinkleExtractor()
582 {
583 this->SetImplicitFunction(nullptr);
584 }
585
586 //------------------------------------------------------------------------------
587 // Overload standard modified time function. If the implicit function
588 // definition is modified, then this object is modified as well.
GetMTime()589 vtkMTimeType vtk3DLinearGridCrinkleExtractor::GetMTime()
590 {
591 vtkMTimeType mTime = this->Superclass::GetMTime();
592 if (this->ImplicitFunction != nullptr)
593 {
594 vtkMTimeType mTime2 = this->ImplicitFunction->GetMTime();
595 return (mTime2 > mTime ? mTime2 : mTime);
596 }
597 else
598 {
599 return mTime;
600 }
601 }
602
603 //------------------------------------------------------------------------------
604 // Specialized implicit function extraction filter to handle unstructured
605 // grids with 3D linear cells (tetrahedras, hexes, wedges, pyradmids, voxels)
606 //
ProcessPiece(vtkUnstructuredGrid * input,vtkImplicitFunction * f,vtkUnstructuredGrid * grid)607 int vtk3DLinearGridCrinkleExtractor::ProcessPiece(
608 vtkUnstructuredGrid* input, vtkImplicitFunction* f, vtkUnstructuredGrid* grid)
609 {
610 if (input == nullptr || f == nullptr || grid == nullptr)
611 {
612 // Not really an error
613 return 1;
614 }
615
616 // Make sure there is input data to process
617 vtkPoints* inPts = input->GetPoints();
618 vtkIdType numPts = 0;
619 if (inPts)
620 {
621 numPts = inPts->GetNumberOfPoints();
622 }
623 vtkCellArray* cells = input->GetCells();
624 vtkIdType numCells = 0;
625 if (cells)
626 {
627 numCells = cells->GetNumberOfCells();
628 }
629 if (numPts <= 0 || numCells <= 0)
630 {
631 vtkLog(INFO, "Empty input");
632 return 0;
633 }
634
635 // Check the input point type. Only real types are supported.
636 int inPtsType = inPts->GetDataType();
637 if ((inPtsType != VTK_FLOAT && inPtsType != VTK_DOUBLE))
638 {
639 vtkLog(ERROR, "Input point type not supported");
640 return 0;
641 }
642
643 // Output cells go here.
644 vtkCellArray* newCells = vtkCellArray::New();
645
646 // Set up the cells for processing. A specialized iterator is used to traverse the cells.
647 unsigned char* cellTypes =
648 static_cast<unsigned char*>(input->GetCellTypesArray()->GetVoidPointer(0));
649 CellIter* cellIter = new CellIter(numCells, cellTypes, cells);
650
651 // Classify the cell points based on the specified implicit function. A
652 // fast path is available for planes.
653 unsigned char* inout = nullptr;
654 int ptsType = inPts->GetDataType();
655 if (vtkPlane::SafeDownCast(f) != nullptr)
656 { // plane fast path
657 if (ptsType == VTK_FLOAT)
658 {
659 PlaneClassifyPoints<float> classify(inPts, static_cast<vtkPlane*>(f));
660 EXECUTE_SMPFOR(this->SequentialProcessing, numPts, classify);
661 inout = classify.InOutArray;
662 }
663 else if (ptsType == VTK_DOUBLE)
664 {
665 PlaneClassifyPoints<double> classify(inPts, static_cast<vtkPlane*>(f));
666 EXECUTE_SMPFOR(this->SequentialProcessing, numPts, classify);
667 inout = classify.InOutArray;
668 }
669 }
670 else
671 { // general implicit function fast path
672 if (ptsType == VTK_FLOAT)
673 {
674 FunctionClassifyPoints<float> classify(inPts, f);
675 EXECUTE_SMPFOR(this->SequentialProcessing, numPts, classify);
676 inout = classify.InOutArray;
677 }
678 else if (ptsType == VTK_DOUBLE)
679 {
680 FunctionClassifyPoints<double> classify(inPts, f);
681 EXECUTE_SMPFOR(this->SequentialProcessing, numPts, classify);
682 inout = classify.InOutArray;
683 }
684 }
685
686 // Depending on whether we are going to eliminate unused points, use
687 // different extraction techniques. There is a large performance
688 // difference if points are compacted.
689 vtkIdType outNumCells = 0;
690 vtkIdType* cellMap = nullptr;
691 vtkIdType outNumPts = 0;
692 vtkIdType* ptMap = nullptr;
693 vtkPointData* inPD = input->GetPointData();
694 vtkCellData* inCD = input->GetCellData();
695 if (!this->RemoveUnusedPoints)
696 {
697 ExtractCells extract(
698 numPts, cellIter, inout, grid, newCells, this->CopyPointData, this->CopyCellData);
699 EXECUTE_REDUCED_SMPFOR(
700 this->SequentialProcessing, numCells, extract, this->NumberOfThreadsUsed);
701
702 outNumCells = extract.OutputNumCells;
703 cellMap = extract.CellMap;
704
705 grid->SetPoints(inPts);
706 if (this->CopyPointData)
707 {
708 vtkPointData* outPD = grid->GetPointData();
709 outPD->PassData(inPD);
710 }
711 }
712
713 else
714 {
715 ExtractPointsAndCells extract(
716 numPts, cellIter, inout, grid, newCells, this->CopyPointData, this->CopyCellData);
717 EXECUTE_REDUCED_SMPFOR(
718 this->SequentialProcessing, numCells, extract, this->NumberOfThreadsUsed);
719
720 outNumPts = extract.OutputNumPts;
721 ptMap = extract.PointMap;
722
723 outNumCells = extract.OutputNumCells;
724 cellMap = extract.CellMap;
725 }
726
727 // Copy cell data if requested
728 if (this->CopyCellData)
729 {
730 vtkCellData* outCD = grid->GetCellData();
731 ArrayList arrays;
732 outCD->CopyAllocate(inCD, outNumCells);
733 arrays.AddArrays(outNumCells, inCD, outCD);
734 CopyCellAttributes copyCellData(&arrays, cellMap);
735 EXECUTE_SMPFOR(this->SequentialProcessing, outNumCells, copyCellData);
736 delete[] cellMap;
737 }
738
739 if (this->RemoveUnusedPoints)
740 {
741 // Create the output points if not passing through. Only real types are
742 // supported. Use the point map to create them.
743 int inType = inPts->GetDataType(), outType;
744 void *inPtr, *outPtr;
745 vtkPoints* outPts = vtkPoints::New();
746 if (this->OutputPointsPrecision == vtkAlgorithm::DEFAULT_PRECISION)
747 {
748 outType = inType;
749 }
750 else if (this->OutputPointsPrecision == vtkAlgorithm::SINGLE_PRECISION)
751 {
752 outType = VTK_FLOAT;
753 }
754 else // if(this->OutputPointsPrecision == vtkAlgorithm::DOUBLE_PRECISION)
755 {
756 outType = VTK_DOUBLE;
757 }
758 outPts->SetDataType(outType);
759 outPts->SetNumberOfPoints(outNumPts);
760
761 // Generate points using the point map
762 inPtr = inPts->GetData()->GetVoidPointer(0);
763 outPtr = outPts->GetData()->GetVoidPointer(0);
764 if (inType == VTK_DOUBLE && outType == VTK_DOUBLE)
765 {
766 GeneratePoints<double, double> generatePts((double*)inPtr, ptMap, (double*)outPtr);
767 EXECUTE_SMPFOR(this->SequentialProcessing, numPts, generatePts);
768 }
769 else if (inType == VTK_FLOAT && outType == VTK_FLOAT)
770 {
771 GeneratePoints<float, float> generatePts((float*)inPtr, ptMap, (float*)outPtr);
772 EXECUTE_SMPFOR(this->SequentialProcessing, numPts, generatePts);
773 }
774 else if (inType == VTK_DOUBLE && outType == VTK_FLOAT)
775 {
776 GeneratePoints<double, float> generatePts((double*)inPtr, ptMap, (float*)outPtr);
777 EXECUTE_SMPFOR(this->SequentialProcessing, numPts, generatePts);
778 }
779 else // if ( inType == VTK_FLOAT && outType == VTK_DOUBLE )
780 {
781 GeneratePoints<float, double> generatePts((float*)inPtr, ptMap, (double*)outPtr);
782 EXECUTE_SMPFOR(this->SequentialProcessing, numPts, generatePts);
783 }
784 grid->SetPoints(outPts);
785 outPts->Delete();
786
787 // Use the point map to copy point data if desired
788 if (this->CopyPointData)
789 {
790 vtkPointData* outPD = grid->GetPointData();
791 ArrayList arrays;
792 outPD->CopyAllocate(inPD, outNumPts);
793 arrays.AddArrays(outNumPts, inPD, outPD);
794 CopyPointAttributes copyPointData(&arrays, ptMap);
795 EXECUTE_SMPFOR(this->SequentialProcessing, numPts, copyPointData);
796 delete[] ptMap;
797 }
798 }
799
800 // Report the results of execution
801 vtkLog(INFO,
802 "Extracted: " << grid->GetNumberOfPoints() << " points, " << grid->GetNumberOfCells()
803 << " cells");
804
805 // Clean up
806 delete[] inout;
807 delete cellIter;
808 newCells->Delete();
809
810 return 1;
811 }
812
813 //------------------------------------------------------------------------------
814 // The output dataset type varies depending on the input type.
RequestDataObject(vtkInformation *,vtkInformationVector ** inputVector,vtkInformationVector * outputVector)815 int vtk3DLinearGridCrinkleExtractor::RequestDataObject(
816 vtkInformation*, vtkInformationVector** inputVector, vtkInformationVector* outputVector)
817 {
818 vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
819 if (!inInfo)
820 {
821 return 0;
822 }
823
824 vtkDataObject* inputDO = vtkDataObject::GetData(inputVector[0], 0);
825 vtkDataObject* outputDO = vtkDataObject::GetData(outputVector, 0);
826 assert(inputDO != nullptr);
827
828 vtkInformation* outInfo = outputVector->GetInformationObject(0);
829
830 if (vtkUnstructuredGrid::SafeDownCast(inputDO))
831 {
832 if (vtkUnstructuredGrid::SafeDownCast(outputDO) == nullptr)
833 {
834 outputDO = vtkUnstructuredGrid::New();
835 outInfo->Set(vtkDataObject::DATA_OBJECT(), outputDO);
836 outputDO->Delete();
837 }
838 return 1;
839 }
840
841 if (vtkCompositeDataSet::SafeDownCast(inputDO))
842 {
843 // For any composite dataset, we're create a vtkMultiBlockDataSet as output;
844 if (vtkMultiBlockDataSet::SafeDownCast(outputDO) == nullptr)
845 {
846 outputDO = vtkMultiBlockDataSet::New();
847 outInfo->Set(vtkDataObject::DATA_OBJECT(), outputDO);
848 outputDO->Delete();
849 }
850 return 1;
851 }
852
853 vtkLog(ERROR, "Not sure what type of output to create!");
854 return 0;
855 }
856
857 //------------------------------------------------------------------------------
858 // Specialized extraction filter to handle unstructured grids with 3D
859 // linear cells (tetrahedras, hexes, wedges, pyradmids, voxels)
860 //
RequestData(vtkInformation *,vtkInformationVector ** inputVector,vtkInformationVector * outputVector)861 int vtk3DLinearGridCrinkleExtractor::RequestData(
862 vtkInformation*, vtkInformationVector** inputVector, vtkInformationVector* outputVector)
863 {
864 // Get the input and output
865 vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
866 vtkInformation* outInfo = outputVector->GetInformationObject(0);
867
868 vtkUnstructuredGrid* inputGrid =
869 vtkUnstructuredGrid::SafeDownCast(inInfo->Get(vtkDataObject::DATA_OBJECT()));
870 vtkUnstructuredGrid* outputGrid =
871 vtkUnstructuredGrid::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));
872
873 vtkCompositeDataSet* inputCDS =
874 vtkCompositeDataSet::SafeDownCast(inInfo->Get(vtkDataObject::DATA_OBJECT()));
875 vtkMultiBlockDataSet* outputMBDS =
876 vtkMultiBlockDataSet::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));
877
878 // Make sure we have valid input and output of some form
879 if ((inputGrid == nullptr || outputGrid == nullptr) &&
880 (inputCDS == nullptr || outputMBDS == nullptr))
881 {
882 return 0;
883 }
884
885 // Need an implicit function to do the cutting
886 vtkImplicitFunction* f = this->ImplicitFunction;
887 if (!f)
888 {
889 vtkLog(ERROR, "Implicit function not defined");
890 return 0;
891 }
892
893 // If the input is an unstructured grid, then simply process this single
894 // grid producing a single output vtkUnstructuredGrid.
895 if (inputGrid)
896 {
897 this->ProcessPiece(inputGrid, f, outputGrid);
898 }
899
900 // Otherwise it is an input composite data set and each unstructured grid
901 // contained in it is processed, producing a vtkGrid that is added to
902 // the output multiblock dataset.
903 else
904 {
905 vtkUnstructuredGrid* grid;
906 vtkUnstructuredGrid* output;
907 outputMBDS->CopyStructure(inputCDS);
908 vtkSmartPointer<vtkCompositeDataIterator> inIter;
909 inIter.TakeReference(inputCDS->NewIterator());
910 for (inIter->InitTraversal(); !inIter->IsDoneWithTraversal(); inIter->GoToNextItem())
911 {
912 auto ds = inIter->GetCurrentDataObject();
913 if ((grid = vtkUnstructuredGrid::SafeDownCast(ds)))
914 {
915 output = vtkUnstructuredGrid::New();
916 this->ProcessPiece(grid, f, output);
917 outputMBDS->SetDataSet(inIter, output);
918 output->Delete();
919 }
920 else
921 {
922 vtkLog(INFO, << "This filter only processes unstructured grids");
923 }
924 }
925 }
926
927 return 1;
928 }
929
930 //------------------------------------------------------------------------------
SetOutputPointsPrecision(int precision)931 void vtk3DLinearGridCrinkleExtractor::SetOutputPointsPrecision(int precision)
932 {
933 this->OutputPointsPrecision = precision;
934 this->Modified();
935 }
936
937 //------------------------------------------------------------------------------
GetOutputPointsPrecision() const938 int vtk3DLinearGridCrinkleExtractor::GetOutputPointsPrecision() const
939 {
940 return this->OutputPointsPrecision;
941 }
942
943 //------------------------------------------------------------------------------
CanFullyProcessDataObject(vtkDataObject * object)944 bool vtk3DLinearGridCrinkleExtractor::CanFullyProcessDataObject(vtkDataObject* object)
945 {
946 auto ug = vtkUnstructuredGrid::SafeDownCast(object);
947 auto cd = vtkCompositeDataSet::SafeDownCast(object);
948
949 if (ug)
950 {
951 // Get list of cell types in the unstructured grid
952 vtkNew<vtkCellTypes> cellTypes;
953 ug->GetCellTypes(cellTypes);
954 for (vtkIdType i = 0; i < cellTypes->GetNumberOfTypes(); ++i)
955 {
956 unsigned char cellType = cellTypes->GetCellType(i);
957 if (cellType != VTK_VOXEL && cellType != VTK_TETRA && cellType != VTK_HEXAHEDRON &&
958 cellType != VTK_WEDGE && cellType != VTK_PYRAMID)
959 {
960 // Unsupported cell type, can't process data
961 return false;
962 }
963 }
964
965 // All cell types are supported, can process data.
966 return true;
967 }
968 else if (cd)
969 {
970 bool supported = true;
971 vtkSmartPointer<vtkCompositeDataIterator> iter;
972 iter.TakeReference(cd->NewIterator());
973 iter->SkipEmptyNodesOn();
974 for (iter->InitTraversal(); !iter->IsDoneWithTraversal(); iter->GoToNextItem())
975 {
976 auto leafDS = iter->GetCurrentDataObject();
977 if (!CanFullyProcessDataObject(leafDS))
978 {
979 supported = false;
980 break;
981 }
982 }
983 return supported;
984 }
985
986 return false; // not a vtkUnstructuredGrid nor a composite dataset
987 }
988
989 //------------------------------------------------------------------------------
FillInputPortInformation(int,vtkInformation * info)990 int vtk3DLinearGridCrinkleExtractor::FillInputPortInformation(int, vtkInformation* info)
991 {
992 info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkUnstructuredGrid");
993 info->Append(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkCompositeDataSet");
994 return 1;
995 }
996
997 //------------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)998 void vtk3DLinearGridCrinkleExtractor::PrintSelf(ostream& os, vtkIndent indent)
999 {
1000 this->Superclass::PrintSelf(os, indent);
1001
1002 os << indent << "Implicit Function: " << this->ImplicitFunction << "\n";
1003
1004 os << indent << "Copy Point Data: " << (this->CopyPointData ? "true\n" : "false\n");
1005 os << indent << "Copy Cell Data: " << (this->CopyCellData ? "true\n" : "false\n");
1006
1007 os << indent << "RemoveUnusedPoints: " << (this->RemoveUnusedPoints ? "true\n" : "false\n");
1008
1009 os << indent << "Precision of the output points: " << this->OutputPointsPrecision << "\n";
1010
1011 os << indent << "Sequential Processing: " << (this->SequentialProcessing ? "true\n" : "false\n");
1012 }
1013
1014 #undef EXECUTE_SMPFOR
1015 #undef EXECUTE_REDUCED_SMPFOR
1016 #undef MAX_CELL_VERTS
1017