1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkExtractCells.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 Copyright (c) Sandia Corporation
17 See Copyright.txt or http://www.paraview.org/HTML/Copyright.html for details.
18 ----------------------------------------------------------------------------*/
19
20 #include "vtkExtractCells.h"
21
22 #include "vtkCellArray.h"
23 #include "vtkIdTypeArray.h"
24 #include "vtkIntArray.h"
25 #include "vtkUnsignedCharArray.h"
26 #include "vtkUnstructuredGrid.h"
27 #include "vtkCell.h"
28 #include "vtkNew.h"
29 #include "vtkPoints.h"
30 #include "vtkPointData.h"
31 #include "vtkPointSet.h"
32 #include "vtkCellData.h"
33 #include "vtkIntArray.h"
34 #include "vtkInformation.h"
35 #include "vtkInformationVector.h"
36 #include "vtkObjectFactory.h"
37 #include "vtkTimeStamp.h"
38 #include "vtkSMPTools.h"
39
40 vtkStandardNewMacro(vtkExtractCells);
41
42 #include <algorithm>
43 #include <numeric>
44 #include <vector>
45
46 namespace {
47 struct FastPointMap
48 {
49 using ConstIteratorType = const vtkIdType*;
50
51 vtkNew<vtkIdList> Map;
52 vtkIdType LastInput;
53 vtkIdType LastOutput;
54
CBegin__anona6f042220111::FastPointMap55 ConstIteratorType CBegin() const
56 {
57 return this->Map->GetPointer(0);
58 }
59
CEnd__anona6f042220111::FastPointMap60 ConstIteratorType CEnd() const
61 {
62 return this->Map->GetPointer(this->Map->GetNumberOfIds());
63 }
64
Reset__anona6f042220111::FastPointMap65 vtkIdType* Reset(vtkIdType numValues)
66 {
67 this->LastInput = -1;
68 this->LastOutput = -1;
69 this->Map->SetNumberOfIds(numValues);
70 return this->Map->GetPointer(0);
71 }
72
73 // Map inputId to the new PointId. If inputId is invalid, return -1.
LookUp__anona6f042220111::FastPointMap74 vtkIdType LookUp(vtkIdType inputId)
75 {
76 vtkIdType outputId = -1;
77 ConstIteratorType first;
78 ConstIteratorType last;
79
80 if (this->LastOutput >= 0)
81 {
82 // Here's the optimization: since the point ids are usually requested
83 // with some locality, we can reduce the search range by caching the
84 // results of the last lookup. This reduces the number of lookups and
85 // improves CPU cache behavior.
86
87 // Offset is the distance (in input space) between the last lookup and
88 // the current id. Since the point map is sorted and unique, this is the
89 // maximum distance that the current ID can be from the previous one.
90 vtkIdType offset = inputId - this->LastInput;
91
92 // Our search range is from the last output location
93 first = this->CBegin() + this->LastOutput;
94 last = first + offset;
95
96 // Ensure these are correctly ordered (offset may be < 0):
97 if (last < first)
98 {
99 std::swap(first, last);
100 }
101
102 // Adjust last to be past-the-end:
103 ++last;
104
105 // Clamp to map bounds:
106 first = std::max(first, this->CBegin());
107 last = std::min(last, this->CEnd());
108 }
109 else
110 { // First run, use full range:
111 first = this->CBegin();
112 last = this->CEnd();
113 }
114
115 outputId = this->BinaryFind(first, last, inputId);
116 if (outputId >= 0)
117 {
118 this->LastInput = inputId;
119 this->LastOutput = outputId;
120 }
121
122 return outputId;
123 }
124
125 private:
126 // Modified version of std::lower_bound that returns as soon as a value is
127 // found (rather than finding the beginning of a sequence). Returns the
128 // position in the list, or -1 if not found.
BinaryFind__anona6f042220111::FastPointMap129 vtkIdType BinaryFind(ConstIteratorType first, ConstIteratorType last,
130 vtkIdType val) const
131 {
132 vtkIdType len = last - first;
133
134 while (len > 0)
135 {
136 // Select median
137 vtkIdType half = len / 2;
138 ConstIteratorType middle = first + half;
139
140 const vtkIdType &mVal = *middle;
141 if (mVal < val)
142 { // This soup is too cold.
143 first = middle;
144 ++first;
145 len = len - half - 1;
146 }
147 else if (val < mVal)
148 { // This soup is too hot!
149 len = half;
150 }
151 else
152 { // This soup is juuuust right.
153 return middle - this->Map->GetPointer(0);
154 }
155 }
156
157 return -1;
158 }
159 };
160 } // end anon namespace
161
162 class vtkExtractCellsSTLCloak
163 {
164 public:
165 std::vector<vtkIdType> CellIds;
166 vtkTimeStamp ModifiedTime;
167 vtkTimeStamp SortTime;
168 FastPointMap PointMap;
169
Modified()170 void Modified()
171 {
172 this->ModifiedTime.Modified();
173 }
174
IsPrepared() const175 inline bool IsPrepared() const
176 {
177 return this->ModifiedTime.GetMTime() < this->SortTime.GetMTime();
178 }
179
Prepare(vtkIdType numInputCells)180 void Prepare(vtkIdType numInputCells)
181 {
182 if (!this->IsPrepared())
183 {
184 vtkSMPTools::Sort(this->CellIds.begin(), this->CellIds.end());
185 auto last = std::unique(this->CellIds.begin(), this->CellIds.end());
186 auto first=this->CellIds.begin();
187
188 // These lines clamp the ids to the number of cells in the
189 // dataset. Otherwise segfaults occur when cellIds are greater than the
190 // number of input cells, in particular when cellId==numInputCells.
191 auto clampLast = std::find(first, last, numInputCells);
192 last = ( clampLast != last ? clampLast : last );
193
194 this->CellIds.resize(std::distance(first, last));
195 this->SortTime.Modified();
196 }
197 }
198 };
199
200 //----------------------------------------------------------------------------
vtkExtractCells()201 vtkExtractCells::vtkExtractCells()
202 {
203 this->SubSetUGridCellArraySize = 0;
204 this->InputIsUgrid = 0;
205 this->CellList = new vtkExtractCellsSTLCloak;
206 }
207
208 //----------------------------------------------------------------------------
~vtkExtractCells()209 vtkExtractCells::~vtkExtractCells()
210 {
211 delete this->CellList;
212 }
213
214 //----------------------------------------------------------------------------
SetCellList(vtkIdList * l)215 void vtkExtractCells::SetCellList(vtkIdList *l)
216 {
217 delete this->CellList;
218 this->CellList = new vtkExtractCellsSTLCloak;
219
220 if (l != nullptr)
221 {
222 this->AddCellList(l);
223 }
224 }
225
226 //----------------------------------------------------------------------------
AddCellList(vtkIdList * l)227 void vtkExtractCells::AddCellList(vtkIdList *l)
228 {
229 const vtkIdType inputSize = l ? l->GetNumberOfIds() : 0;
230 if (inputSize == 0)
231 {
232 return;
233 }
234
235 const vtkIdType *inputBegin = l->GetPointer(0);
236 const vtkIdType *inputEnd = inputBegin + inputSize;
237
238 const std::size_t oldSize = this->CellList->CellIds.size();
239 const std::size_t newSize = oldSize + static_cast<std::size_t>(inputSize);
240 this->CellList->CellIds.resize(newSize);
241
242 auto outputBegin = this->CellList->CellIds.begin();
243 std::advance(outputBegin, oldSize);
244
245 std::copy(inputBegin, inputEnd, outputBegin);
246
247 this->CellList->Modified();
248 }
249
250 //----------------------------------------------------------------------------
AddCellRange(vtkIdType from,vtkIdType to)251 void vtkExtractCells::AddCellRange(vtkIdType from, vtkIdType to)
252 {
253 if (to < from || to < 0 )
254 {
255 vtkWarningMacro("Bad cell range: (" << to << "," << from << ")");
256 return;
257 }
258
259 // This range specification is inconsistent with C++. Left for backward
260 // compatibility reasons.
261 const vtkIdType inputSize = to - from + 1; // +1 to include 'to'
262 const std::size_t oldSize = this->CellList->CellIds.size();
263 const std::size_t newSize = oldSize + static_cast<std::size_t>(inputSize);
264
265 this->CellList->CellIds.resize(newSize);
266
267 auto outputBegin = this->CellList->CellIds.begin() + oldSize;
268 auto outputEnd = this->CellList->CellIds.begin() + newSize;
269
270 std::iota(outputBegin, outputEnd, from);
271
272 this->CellList->Modified();
273 }
274
275 //----------------------------------------------------------------------------
GetMTime()276 vtkMTimeType vtkExtractCells::GetMTime()
277 {
278 vtkMTimeType mTime = this->Superclass::GetMTime();
279 mTime = std::max(mTime, this->CellList->ModifiedTime.GetMTime());
280 mTime = std::max(mTime, this->CellList->SortTime.GetMTime());
281 return mTime;
282 }
283
284 //----------------------------------------------------------------------------
RequestData(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)285 int vtkExtractCells::RequestData(
286 vtkInformation *vtkNotUsed(request),
287 vtkInformationVector **inputVector,
288 vtkInformationVector *outputVector)
289 {
290 // get the info objects
291 vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
292 vtkInformation *outInfo = outputVector->GetInformationObject(0);
293
294 // get the input and output
295 vtkDataSet *input = vtkDataSet::SafeDownCast(
296 inInfo->Get(vtkDataObject::DATA_OBJECT()));
297 vtkUnstructuredGrid *output = vtkUnstructuredGrid::SafeDownCast(
298 outInfo->Get(vtkDataObject::DATA_OBJECT()));
299
300 // Sort/uniquify the cell ids if needed.
301 vtkIdType numCellsInput = input->GetNumberOfCells();
302 this->CellList->Prepare(numCellsInput);
303
304 this->InputIsUgrid =
305 ((vtkUnstructuredGrid::SafeDownCast(input)) != nullptr);
306
307 vtkIdType numCells = static_cast<vtkIdType>(this->CellList->CellIds.size());
308
309 if (numCells == numCellsInput)
310 {
311 this->Copy(input, output);
312 return 1;
313 }
314
315 vtkPointData *PD = input->GetPointData();
316 vtkCellData *CD = input->GetCellData();
317
318 if (numCells == 0)
319 {
320 // set up a ugrid with same data arrays as input, but
321 // no points, cells or data.
322
323 output->Allocate(1);
324
325 output->GetPointData()->CopyGlobalIdsOn();
326 output->GetPointData()->CopyAllocate(PD, VTK_CELL_SIZE);
327 output->GetCellData()->CopyGlobalIdsOn();
328 output->GetCellData()->CopyAllocate(CD, 1);
329
330 vtkPoints *pts = vtkPoints::New();
331 pts->SetNumberOfPoints(0);
332
333 output->SetPoints(pts);
334
335 pts->Delete();
336
337 return 1;
338 }
339
340 vtkPointData *newPD = output->GetPointData();
341 vtkCellData *newCD = output->GetCellData();
342
343 vtkIdType numPoints = reMapPointIds(input);
344
345 newPD->CopyGlobalIdsOn();
346 newPD->CopyAllocate(PD, numPoints);
347
348 newCD->CopyGlobalIdsOn();
349 newCD->CopyAllocate(CD, numCells);
350
351 vtkPoints *pts = vtkPoints::New();
352 if(vtkPointSet* inputPS = vtkPointSet::SafeDownCast(input))
353 {
354 // preserve input datatype
355 pts->SetDataType(inputPS->GetPoints()->GetDataType());
356 }
357 pts->SetNumberOfPoints(numPoints);
358
359 // Copy points and point data:
360 vtkPointSet *pointSet;
361 if ((pointSet = vtkPointSet::SafeDownCast(input)))
362 { // Optimize when a vtkPoints object exists in the input:
363 vtkNew<vtkIdList> dstIds; // contiguous range [0, numPoints)
364 dstIds->SetNumberOfIds(numPoints);
365 std::iota(dstIds->GetPointer(0), dstIds->GetPointer(numPoints), 0);
366
367 pts->InsertPoints(dstIds, this->CellList->PointMap.Map, pointSet->GetPoints());
368 newPD->CopyData(PD, this->CellList->PointMap.Map, dstIds);
369 }
370 else
371 { // Slow path if we have to query the dataset:
372 for (vtkIdType newId = 0; newId < numPoints; ++newId)
373 {
374 vtkIdType oldId = this->CellList->PointMap.Map->GetId(newId);
375 pts->SetPoint(newId, input->GetPoint(oldId));
376 newPD->CopyData(PD, oldId, newId);
377 }
378 }
379
380 output->SetPoints(pts);
381 pts->Delete();
382
383 if (this->InputIsUgrid)
384 {
385 this->CopyCellsUnstructuredGrid(input, output);
386 }
387 else
388 {
389 this->CopyCellsDataSet(input, output);
390 }
391
392 this->CellList->PointMap.Reset(0);
393 output->Squeeze();
394
395 return 1;
396 }
397
398 //----------------------------------------------------------------------------
Copy(vtkDataSet * input,vtkUnstructuredGrid * output)399 void vtkExtractCells::Copy(vtkDataSet *input, vtkUnstructuredGrid *output)
400 {
401 // If input is unstructured grid just deep copy through
402 if (this->InputIsUgrid)
403 {
404 output->DeepCopy(vtkUnstructuredGrid::SafeDownCast(input));
405 return;
406 }
407
408 vtkIdType numPoints = input->GetNumberOfPoints();
409 vtkIdType numCells = input->GetNumberOfCells();
410
411 vtkPointData *PD = input->GetPointData();
412 vtkCellData *CD = input->GetCellData();
413 vtkPointData *newPD = output->GetPointData();
414 vtkCellData *newCD = output->GetCellData();
415 newPD->CopyAllocate(PD, numPoints);
416 newCD->CopyAllocate(CD, numCells);
417
418 output->Allocate(numCells);
419
420 vtkPoints *pts = vtkPoints::New();
421 pts->SetNumberOfPoints(numPoints);
422 output->SetPoints(pts);
423
424 for (vtkIdType i=0; i<numPoints; i++)
425 {
426 pts->SetPoint(i, input->GetPoint(i));
427 }
428 newPD->DeepCopy(PD);
429 pts->Delete();
430
431 vtkIdList *cellPoints = vtkIdList::New();
432 for (vtkIdType cellId=0; cellId < numCells; cellId++)
433 {
434 input->GetCellPoints(cellId, cellPoints);
435 output->InsertNextCell(input->GetCellType(cellId), cellPoints);
436 }
437 newCD->DeepCopy(CD);
438 cellPoints->Delete();
439
440 output->Squeeze();
441 }
442
443 //----------------------------------------------------------------------------
reMapPointIds(vtkDataSet * grid)444 vtkIdType vtkExtractCells::reMapPointIds(vtkDataSet *grid)
445 {
446 vtkIdType totalPoints = grid->GetNumberOfPoints();
447
448 char *temp = new char [totalPoints];
449
450 if (!temp)
451 {
452 vtkErrorMacro(<< "vtkExtractCells::reMapPointIds memory allocation");
453 return 0;
454 }
455 memset(temp, 0, totalPoints);
456
457 vtkIdType numberOfIds = 0;
458 vtkIdType i;
459 vtkIdType id;
460 vtkIdList *ptIds = vtkIdList::New();
461 std::vector<vtkIdType>::const_iterator cellPtr;
462
463 if (!this->InputIsUgrid)
464 {
465 for (cellPtr = this->CellList->CellIds.cbegin();
466 cellPtr != this->CellList->CellIds.cend();
467 ++cellPtr)
468 {
469 grid->GetCellPoints(*cellPtr, ptIds);
470
471 vtkIdType nIds = ptIds->GetNumberOfIds();
472
473 vtkIdType *ptId = ptIds->GetPointer(0);
474
475 for (i=0; i<nIds; ++i)
476 {
477 id = *ptId++;
478
479 if (temp[id] == 0)
480 {
481 ++numberOfIds;
482 temp[id] = 1;
483 }
484 }
485 }
486 }
487 else
488 {
489 vtkUnstructuredGrid *ugrid = vtkUnstructuredGrid::SafeDownCast(grid);
490
491 this->SubSetUGridCellArraySize = 0;
492
493 vtkIdType *cellArray = ugrid->GetCells()->GetPointer();
494 vtkIdType *locs = ugrid->GetCellLocationsArray()->GetPointer(0);
495
496 this->SubSetUGridCellArraySize = 0;
497 vtkIdType maxid = ugrid->GetCellLocationsArray()->GetMaxId();
498
499 for (cellPtr = this->CellList->CellIds.cbegin();
500 cellPtr != this->CellList->CellIds.cend();
501 ++cellPtr)
502 {
503 if (*cellPtr > maxid) continue;
504
505 vtkIdType loc = locs[*cellPtr];
506
507 vtkIdType nIds = cellArray[loc++];
508
509 this->SubSetUGridCellArraySize += (1 + nIds);
510
511 for (i=0; i<nIds; ++i)
512 {
513 id = cellArray[loc++];
514
515 if (temp[id] == 0)
516 {
517 ++numberOfIds;
518 temp[id] = 1;
519 }
520 }
521 }
522 }
523 ptIds->Delete();
524 ptIds = nullptr;
525
526 vtkIdType *pointMap = this->CellList->PointMap.Reset(numberOfIds);
527
528 for (id=0; id<totalPoints; id++)
529 {
530 if (temp[id])
531 {
532 (*pointMap++) = id;
533 }
534 }
535
536 delete [] temp;
537
538 return numberOfIds;
539 }
540
541 //----------------------------------------------------------------------------
CopyCellsDataSet(vtkDataSet * input,vtkUnstructuredGrid * output)542 void vtkExtractCells::CopyCellsDataSet(vtkDataSet *input,
543 vtkUnstructuredGrid *output)
544 {
545 output->Allocate(static_cast<vtkIdType>(this->CellList->CellIds.size()));
546
547 vtkCellData *oldCD = input->GetCellData();
548 vtkCellData *newCD = output->GetCellData();
549
550 // We only create vtkOriginalCellIds for the output data set if it does not
551 // exist in the input data set. If it is in the input data set then we
552 // let CopyData() take care of copying it over.
553 vtkIdTypeArray *origMap = nullptr;
554 if(oldCD->GetArray("vtkOriginalCellIds") == nullptr)
555 {
556 origMap = vtkIdTypeArray::New();
557 origMap->SetNumberOfComponents(1);
558 origMap->SetName("vtkOriginalCellIds");
559 newCD->AddArray(origMap);
560 origMap->Delete();
561 }
562
563 vtkIdList *cellPoints = vtkIdList::New();
564
565 std::vector<vtkIdType>::const_iterator cellPtr;
566
567 for (cellPtr = this->CellList->CellIds.cbegin();
568 cellPtr != this->CellList->CellIds.cend();
569 ++cellPtr)
570 {
571 vtkIdType cellId = *cellPtr;
572
573 input->GetCellPoints(cellId, cellPoints);
574
575 for (int i=0; i < cellPoints->GetNumberOfIds(); i++)
576 {
577 vtkIdType oldId = cellPoints->GetId(i);
578
579 vtkIdType newId = this->CellList->PointMap.LookUp(oldId);
580 assert("Old id exists in map." && newId >= 0);
581
582 cellPoints->SetId(i, newId);
583 }
584 vtkIdType newId = output->InsertNextCell(input->GetCellType(cellId), cellPoints);
585
586 newCD->CopyData(oldCD, cellId, newId);
587 if(origMap)
588 {
589 origMap->InsertNextValue(cellId);
590 }
591 }
592
593 cellPoints->Delete();
594 }
595
596 //----------------------------------------------------------------------------
CopyCellsUnstructuredGrid(vtkDataSet * input,vtkUnstructuredGrid * output)597 void vtkExtractCells::CopyCellsUnstructuredGrid(vtkDataSet *input,
598 vtkUnstructuredGrid *output)
599 {
600 vtkUnstructuredGrid *ugrid = vtkUnstructuredGrid::SafeDownCast(input);
601 if (ugrid == nullptr)
602 {
603 this->CopyCellsDataSet(input, output);
604 return;
605 }
606
607 vtkCellData *oldCD = input->GetCellData();
608 vtkCellData *newCD = output->GetCellData();
609
610 // We only create vtkOriginalCellIds for the output data set if it does not
611 // exist in the input data set. If it is in the input data set then we
612 // let CopyData() take care of copying it over.
613 vtkIdTypeArray *origMap = nullptr;
614 if(oldCD->GetArray("vtkOriginalCellIds") == nullptr)
615 {
616 origMap = vtkIdTypeArray::New();
617 origMap->SetNumberOfComponents(1);
618 origMap->SetName("vtkOriginalCellIds");
619 newCD->AddArray(origMap);
620 origMap->Delete();
621 }
622
623 vtkIdType numCells = static_cast<vtkIdType>(this->CellList->CellIds.size());
624
625 vtkCellArray *cellArray = vtkCellArray::New(); // output
626 vtkIdTypeArray *newcells = vtkIdTypeArray::New();
627 newcells->SetNumberOfValues(this->SubSetUGridCellArraySize);
628 cellArray->SetCells(numCells, newcells);
629 vtkIdType cellArrayIdx = 0;
630
631 vtkIdTypeArray *locationArray = vtkIdTypeArray::New();
632 locationArray->SetNumberOfValues(numCells);
633
634 vtkUnsignedCharArray *typeArray = vtkUnsignedCharArray::New();
635 typeArray->SetNumberOfValues(numCells);
636
637 vtkIdType nextCellId = 0;
638
639 std::vector<vtkIdType>::const_iterator cellPtr; // input
640 vtkIdType *cells = ugrid->GetCells()->GetPointer();
641 vtkIdType maxid = ugrid->GetCellLocationsArray()->GetMaxId();
642 vtkIdType *locs = ugrid->GetCellLocationsArray()->GetPointer(0);
643 vtkUnsignedCharArray *types = ugrid->GetCellTypesArray();
644
645 for (cellPtr = this->CellList->CellIds.cbegin();
646 cellPtr != this->CellList->CellIds.cend();
647 ++cellPtr)
648 {
649 if (*cellPtr > maxid) continue;
650
651 vtkIdType oldCellId = *cellPtr;
652
653 vtkIdType loc = locs[oldCellId];
654 int size = static_cast<int>(cells[loc]);
655 vtkIdType *pts = cells + loc + 1;
656 unsigned char type = types->GetValue(oldCellId);
657
658 locationArray->SetValue(nextCellId, cellArrayIdx);
659 typeArray->SetValue(nextCellId, type);
660
661 newcells->SetValue(cellArrayIdx++, size);
662
663 for (int i=0; i<size; i++)
664 {
665 vtkIdType oldId = *pts++;
666 vtkIdType newId = this->CellList->PointMap.LookUp(oldId);
667 assert("Old id exists in map." && newId >= 0);
668
669 newcells->SetValue(cellArrayIdx++, newId);
670 }
671
672 newCD->CopyData(oldCD, oldCellId, nextCellId);
673 if(origMap)
674 {
675 origMap->InsertNextValue(oldCellId);
676 }
677 nextCellId++;
678 }
679
680 output->SetCells(typeArray, locationArray, cellArray);
681
682 typeArray->Delete();
683 locationArray->Delete();
684 newcells->Delete();
685 cellArray->Delete();
686 }
687
688 //----------------------------------------------------------------------------
FillInputPortInformation(int,vtkInformation * info)689 int vtkExtractCells::FillInputPortInformation(int, vtkInformation *info)
690 {
691 info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataSet");
692 return 1;
693 }
694
695 //----------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)696 void vtkExtractCells::PrintSelf(ostream& os, vtkIndent indent)
697 {
698 this->Superclass::PrintSelf(os,indent);
699 }
700