GSaha567/seq_level_training_data
052
1text,length,is_long_context,metric_val,label_metric2"/*=========================================================================3 4 Program: Visualization Toolkit5 Module: vtkBinnedDecimation.cxx6 7 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen8 All rights reserved.9 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.10 11 This software is distributed WITHOUT ANY WARRANTY; without even12 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR13 PURPOSE. See the above copyright notice for more information.14 15=========================================================================*/16#include ""vtkBinnedDecimation.h""17 18#include ""vtkArrayDispatch.h""19#include ""vtkArrayListTemplate.h"" // For processing attribute data20#include ""vtkBoundingBox.h""21#include ""vtkCellArray.h""22#include ""vtkCellArrayIterator.h""23#include ""vtkCellData.h""24#include ""vtkDataArray.h""25#include ""vtkDataArrayRange.h""26#include ""vtkExecutive.h""27#include ""vtkFloatArray.h""28#include ""vtkIdTypeArray.h""29#include ""vtkInformation.h""30#include ""vtkInformationVector.h""31#include ""vtkLogger.h""32#include ""vtkMath.h""33#include ""vtkObjectFactory.h""34#include ""vtkPointData.h""35#include ""vtkPolyData.h""36#include ""vtkSMPThreadLocal.h""37#include ""vtkSMPTools.h""38#include ""vtkTriangle.h""39 40#include <atomic>41#include <vector>42 43vtkStandardNewMacro(vtkBinnedDecimation);44 45//----------------------------------------------------------------------------46// Core algorithms and functors to implement threading and type dispatch. Note47// that four different algorithms are implemented, partially for demonstrative48// purposes but also for the fun of it. If I had to pick one, I'd pick49// algorithm 4) BIN_CENTERS which scales better for large numbers of bins, and50// generally produces the best results.51 52namespace53{54 55// There are essentially four different algorithms implemented here depending56// on how the output points are generated in each bin: 1) reuse the input57// points; 2) generate new points by selecting one of the points falling58// into each bin; 3) generate new points at bin centers; and 4) average the59// points (and attribute data) contained within the bins.60 61// Functor to bin points to generate bin ids. The binning functor is common62// to all the algorithms / output options.63template <typename PointsT, typename TIds>64struct BinPoints65{66 PointsT* Points;67 TIds* BinIds;68 int Divisions[3];69 double Bounds[6];70 double H[3];71 // For performance72 double hX, hY, hZ;73 double fX, fY, fZ, bX, bY, bZ;74 vtkIdType xD, yD, zD, xyD;75 76 BinPoints(77 PointsT* pts, TIds* binIds, const int* dims, const double* bounds, const double* spacing)78 : Points(pts)79 , BinIds(binIds)80 {81 for (auto i = 0; i < 3; ++i)82 {83 this->Divisions[i] = dims[i];84 this->Bounds[2 * i] = bounds[2 * i];85 this->Bounds[2 * i + 1] = bounds[2 * i + 1];86 this->H[i] = spacing[i];87 }88 89 // Used for performance90 this->hX = this->H[0] = spacing[0];91 this->hY = this->H[1] = spacing[1];92 this->hZ = this->H[2] = spacing[2];93 this->fX = 1.0 / spacing[0];94 this->fY = 1.0 / spacing[1];95 this->fZ = 1.0 / spacing[2];96 this->bX = this->Bounds[0] = bounds[0];97 this->Bounds[1] = bounds[1];98 this->bY = this->Bounds[2] = bounds[2];99 this->Bounds[3] = bounds[3];100 this->bZ = this->Bounds[4] = bounds[4];101 this->Bounds[5] = bounds[5];102 this->xD = this->Divisions[0];103 this->yD = this->Divisions[1];104 this->zD = this->Divisions[2];105 this->xyD = this->Divisions[0] * this->Divisions[1];106 }107 108 // Functions to obtain bin index inlined for performance.109 void GetBinIndices(const double* x, int ijk[3]) const110 {111 // Compute point index. Make sure it lies within range of locator.112 TIds tmp0 = static_cast<TIds>(((x[0] - bX) * fX));113 TIds tmp1 = static_cast<TIds>(((x[1] - bY) * fY));114 TIds tmp2 = static_cast<TIds>(((x[2] - bZ) * fZ));115 116 ijk[0] = tmp0 < 0 ? 0 : (tmp0 >= xD ? xD - 1 : tmp0);117 ijk[1] = tmp1 < 0 ? 0 : (tmp1 >= yD ? yD - 1 : tmp1);118 ijk[2] = tmp2 < 0 ? 0 : (tmp2 >= zD ? zD - 1 : tmp2);119 }120 121 TIds GetBinIndex(const double* x) const122 {123 int ijk[3];124 this->GetBinIndices(x, ijk);125 return (ijk[0] + ijk[1] * xD + ijk[2] * xyD);126 }127 128 void operator()(vtkIdType ptId, vtkIdType endPtId)129 {130 const auto points = vtk::DataArrayTupleRange<3>(this->Points, ptId, endPtId);131 double x[3];132 TIds* bins = this->BinIds + ptId;133 for (const auto tuple : points)134 {135 x[0] = static_cast<double>(tuple[0]);136 x[1] = static_cast<double>(tuple[1]);137 x[2] = static_cast<double>(tuple[2]);138 *bins++ = this->GetBinIndex(x);139 }140 }141};142 143// Generate the output triangles. This functor is common to144// three of the algorithms #1-3.145template <typename TIds, typename TPtMap>146struct GenerateTriangles147{148 const TIds* BinIds;149 const TPtMap* PointMap;150 vtkCellArray* Tris;151 vtkSMPThreadLocal<vtkSmartPointer<vtkCellArrayIterator>> CellIterator;152 const TIds* TriMap;153 vtkIdType* OutTris;154 vtkIdType* OutTriOffsets;155 ArrayList* Arrays;156 157 GenerateTriangles(TIds* bins, TPtMap* ptMap, vtkCellArray* tris, TIds* triMap, vtkIdType* outTris,158 vtkIdType* outTriOffsets, ArrayList* arrays)159 : BinIds(bins)160 , PointMap(ptMap)161 , Tris(tris)162 , TriMap(triMap)163 , OutTris(outTris)164 , OutTriOffsets(outTriOffsets)165 , Arrays(arrays)166 {167 }168 169 void Initialize() { this->CellIterator.Local().TakeReference(this->Tris->NewIterator()); }170 171 void operator()(vtkIdType triId, vtkIdType endTriId)172 {173 const TIds* binIds = this->BinIds;174 const TPtMap* ptMap = this->PointMap;175 vtkIdType npts;176 const vtkIdType* tri;177 vtkCellArrayIterator* cellIter = this->CellIterator.Local();178 const TIds* triMap = this->TriMap;179 vtkIdType *outTris = this->OutTris, *outTri;180 vtkIdType *outTriOffsets = this->OutTriOffsets, *outOffsets;181 182 for (; triId < endTriId; ++triId)183 {184 if ((triMap[triId + 1] - triMap[triId]) > 0) // spit out triangle185 {186 cellIter->GetCellAtId(triId, npts, tri);187 outOffsets = outTriOffsets + triMap[triId];188 *outOffsets = triMap[triId] * 3;189 outTri = outTris + *outOffsets;190 outTri[0] = ptMap[binIds[tri[0]]];191 outTri[1] = ptMap[binIds[tri[1]]];192 outTri[2] = ptMap[binIds[tri[2]]];193 if (this->Arrays) // copy cell data if requested194 {195 this->Arrays->Copy(triId, triMap[triId]);196 }197 }198 }199 }200 201 void Reduce() {}202};203 204//============================ Implement algorithms ==================205 206//============================ 1) Reuse INPUT_POINTS =================207// Traverse cells and mark points and cells that are included in the output208template <typename TIds>209struct SelectOutput210{211 const TIds* BinIds;212 unsigned char* PointUses;213 vtkCellArray* Tris;214 TIds* TriMap;215 vtkSMPThreadLocal<vtkSmartPointer<vtkCellArrayIterator>> CellIterator;216 217 SelectOutput(TIds* bins, unsigned char* ptUses, vtkCellArray* tris, TIds* triMap)218 : BinIds(bins)219 , PointUses(ptUses)220 , Tris(tris)221 , TriMap(triMap)222 {223 }224 225 void Initialize() { this->CellIterator.Local().TakeReference(this->Tris->NewIterator()); }226 227 void operator()(vtkIdType triId, vtkIdType endTriId)228 {229 vtkIdType npts;230 const vtkIdType* tri;231 vtkCellArrayIterator* cellIter = this->CellIterator.Local();232 TIds* triMap = this->TriMap + triId;233 unsigned char* ptUses = this->PointUses;234 235 for (; triId < endTriId; ++triId, ++triMap)236 {237 cellIter->GetCellAtId(triId, npts, tri);238 // All three points have to be in different bins239 if (this->BinIds[tri[0]] != this->BinIds[tri[1]] &&240 this->BinIds[tri[0]] != this->BinIds[tri[2]] &&241 this->BinIds[tri[1]] != this->BinIds[tri[2]])242 {243 *triMap = 1;244 ptUses[tri[0]] = 1;245 ptUses[tri[1]] = 1;246 ptUses[tri[2]] = 1;247 }248 else // mark excluded from output249 {250 *triMap = 0;251 }252 }253 }254 255 void Reduce() {}256};257 258// Only initialize bins that actually contain a point. We are doing259// this to avoid using std::atomic which the other algorithms use.260template <typename TIds>261struct InitializePointMap262{263 const TIds* BinIds;264 const unsigned char* PointUses;265 TIds* PointMap;266 267 InitializePointMap(TIds* binIds, unsigned char* ptUses, TIds* ptMap)268 : BinIds(binIds)269 , PointUses(ptUses)270 , PointMap(ptMap)271 {272 }273 274 void operator()(vtkIdType ptId, vtkIdType endPtId)275 {276 const TIds* binIds = this->BinIds;277 const unsigned char* ptUses = this->PointUses + ptId;278 TIds* ptMap = this->PointMap;279 280 for (; ptId < endPtId; ++ptId)281 {282 if (*ptUses++ > 0)283 {284 ptMap[binIds[ptId]] = (-1); // mark unvisited285 }286 }287 }288};289 290//----------------------------------------------------------------------------291// Decimation algorithm #1 which reuses the input points.292template <typename PointsT, typename TIds>293void ReuseDecimate(vtkIdType numPts, PointsT* pts, vtkIdType numTris, vtkCellArray* tris,294 vtkCellData* inCD, vtkCellData* outCD, vtkIdType numBins, const int* dims, const double* bounds,295 const double* spacing, vtkPolyData* output)296{297 // Setup execution. Several arrays are used to transform the data.298 // The bin id of each point299 TIds* binIds = new TIds[numPts];300 // Is the point used in the output ?301 unsigned char* ptUses = new unsigned char[numPts];302 std::fill_n(ptUses, numPts, 0);303 // The output point id assigned to each bin (if bin contains an output point).304 TIds* ptMap = new TIds[numBins];305 // Which triangle cells are output? And later, the offsets into the306 // output cell array.307 TIds* triMap = new TIds[numTris + 1];308 309 // Bin points to generate a bin index for each point310 BinPoints<PointsT, TIds> binPoints(pts, binIds, dims, bounds, spacing);311 vtkSMPTools::For(0, numPts, binPoints);312 313 // Select which triangles and points are sent to the output.314 SelectOutput<TIds> selectOutput(binIds, ptUses, tris, triMap);315 vtkSMPTools::For(0, numTris, selectOutput);316 317 // Initialize the point map, only the bins that contain something318 InitializePointMap<TIds> initPtMap(binIds, ptUses, ptMap);319 vtkSMPTools::For(0, numPts, initPtMap);320 321 // Prefix sums to roll up the points and cells, and setup offsets for322 // subsequent threading. This could be threaded, although the gains323 // are likely modest.324 vtkIdType numOutPts = 0;325 for (vtkIdType ptId = 0; ptId < numPts; ++ptId)326 {327 if (ptUses[ptId] > 0)328 {329 if (ptMap[binIds[ptId]] < 0)330 {331 ptMap[binIds[ptId]] = ptId;332 ++numOutPts;333 }334 }335 }336 vtkIdType mark, numOutTris = 0;337 for (vtkIdType triId = 0; triId < numTris; ++triId)338 {339 mark = triMap[triId];340 triMap[triId] = numOutTris;341 numOutTris += mark;342 }343 triMap[numTris] = numOutTris;344 345 // Produce the decimated output346 vtkCellArray* outTrisArray = output->GetPolys();347 vtkNew<vtkIdTypeArray> outConn;348 vtkIdType* outTris = outConn->WritePointer(0, numOutTris * 3);349 vtkNew<vtkIdTypeArray> outOffsets;350 vtkIdType* outTriOffsets = outOffsets->WritePointer(0, numOutTris + 1);351 outTriOffsets[numOutTris] = 3 * numOutTris;352 353 ArrayList arrays;354 if (outCD) // copy cell data if requested355 {356 outCD->CopyAllocate(inCD, numOutTris);357 arrays.AddArrays(numOutTris, inCD, outCD);358 }359 360 // Produce output triangles.361 GenerateTriangles<TIds, TIds> genTris(362 binIds, ptMap, tris, triMap, outTris, outTriOffsets, (outCD != nullptr ? (&arrays) : nullptr));363 vtkSMPTools::For(0, numTris, genTris);364 outTrisArray->SetData(outOffsets, outConn);365 366 // Clean up367 delete[] triMap;368 delete[] ptMap;369 delete[] ptUses;370 delete[] binIds;371}372 373//----------------------------------------------------------------------------374// Dispatch to the decimate algorithm #1 which reuses input points.375struct PointReuseWorker376{377 template <typename DataT>378 void operator()(DataT* pts, bool largeIds, vtkCellArray* tris, vtkCellData* inCD,379 vtkCellData* outCD, int* divs, double bounds[6], double spacing[3], vtkPolyData* output)380 {381 vtkIdType numPts = pts->GetNumberOfTuples();382 vtkIdType numTris = tris->GetNumberOfCells();383 vtkIdType numBins = divs[0] * divs[1] * divs[2];384 385 // Use the appropriate id type for memory and performance reasons.386 if (!largeIds)387 {388 ReuseDecimate<DataT, int>(389 numPts, pts, numTris, tris, inCD, outCD, numBins, divs, bounds, spacing, output);390 }391 else392 {393 ReuseDecimate<DataT, vtkIdType>(394 numPts, pts, numTris, tris, inCD, outCD, numBins, divs, bounds, spacing, output);395 }396 }397};398 399//==========2) Generate new points from selected points BIN_POINTS =============400//==========3) Generate new points at BIN_CENTERS ==============================401// These algorithms #2 and #3 are essentially the same with the difference402// being how the output points are created.403 404// Traverse cells and map input points and cells to output points and cells405template <typename TIds>406struct MapOutput407{408 const TIds* BinIds;409 std::atomic<TIds>* PointMap;410 vtkCellArray* Tris;411 TIds* TriMap;412 vtkSMPThreadLocal<vtkSmartPointer<vtkCellArrayIterator>> CellIterator;413 414 MapOutput(TIds* bins, std::atomic<TIds>* ptMap, vtkCellArray* tris, TIds* triMap)415 : BinIds(bins)416 , PointMap(ptMap)417 , Tris(tris)418 , TriMap(triMap)419 {420 }421 422 // This method is used to select a point id within a bin from a potential423 // set of contributing point ids, which becomes the single point id424 // associated with the bin. Since there are possibly multiple, simultaneous425 // writes to a bin, an atomic is used to prevent data races etc.426 inline void WritePtId(std::atomic<TIds>& binId, vtkIdType ptId)427 {428 // Because of zero initialization, a negative ptId is written. The end result is429 // that we select the smallest point id from the set of points within a bin.430 TIds currentId, targetId = (-(ptId + 1));431 do432 {433 currentId = binId;434 if (currentId < targetId)435 {436 return;437 }438 } while (!atomic_compare_exchange_weak(&binId, ¤tId, targetId));439 }440 441 void Initialize() { this->CellIterator.Local().TakeReference(this->Tris->NewIterator()); }442 443 void operator()(vtkIdType triId, vtkIdType endTriId)444 {445 vtkIdType npts;446 const vtkIdType* tri;447 vtkCellArrayIterator* cellIter = this->CellIterator.Local();448 TIds* triMap = this->TriMap + triId;449 std::atomic<TIds>* ptMap = this->PointMap;450 TIds binIds[3];451 452 for (; triId < endTriId; ++triId, ++triMap)453 {454 cellIter->GetCellAtId(triId, npts, tri);455 binIds[0] = this->BinIds[tri[0]];456 binIds[1] = this->BinIds[tri[1]];457 binIds[2] = this->BinIds[tri[2]];458 459 // All three points have to be in different bins for triangle insertion460 if (binIds[0] != binIds[1] && binIds[0] != binIds[2] && binIds[1] != binIds[2])461 {462 *triMap = 1;463 this->WritePtId(ptMap[binIds[0]], tri[0]);464 this->WritePtId(ptMap[binIds[1]], tri[1]);465 this->WritePtId(ptMap[binIds[2]], tri[2]);466 }467 else // mark excluded from output468 {469 *triMap = 0;470 }471 }472 }473 474 // Reduce is required if Initialize() defined.475 void Reduce() {}476};477 478// Count the number of points in each z-slice of the binning volume. The479// resulting slice offsets are used to thread the generation of the output480// points.481template <typename TIds>482struct CountPoints483{484 const int* Dims;485 std::atomic<TIds>* PointMap;486 int* SliceOffsets;487 488 CountPoints(const int* dims, std::atomic<TIds>* ptMap, int* sliceOffsets)489 : Dims(dims)490 , PointMap(ptMap)491 , SliceOffsets(sliceOffsets)492 {493 }494 495 void Initialize() {}496 497 void operator()(vtkIdType slice, vtkIdType endSlice)498 {499 int binOffset = slice * this->Dims[0] * this->Dims[1];500 501 for (; slice < endSlice; ++slice)502 {503 vtkIdType numSlicePts = 0;504 for (auto j = 0; j < this->Dims[1]; ++j)505 {506 for (auto i = 0; i < this->Dims[0]; ++i)507 {508 if (this->PointMap[binOffset] != 0)509 {510 ++numSlicePts;511 }512 ++binOffset;513 }514 }515 this->SliceOffsets[slice] = numSlicePts;516 } // for all slices in this batch517 }518 519 void Reduce()520 {521 // Prefix sum to roll up total point count across all of the slices.522 TIds numSlicePts, numNewPts = 0;523 for (auto i = 0; i < this->Dims[2]; ++i)524 {525 numSlicePts = this->SliceOffsets[i];526 this->SliceOffsets[i] = numNewPts;527 numNewPts += numSlicePts;528 }529 this->SliceOffsets[this->Dims[2]] = numNewPts;530 }531};532 533// Generate output points; either from bin centers, or from534// selecting one of the points in the bin.535template <typename PointsT, typename TIds>536struct GenerateBinPoints537{538 int PointGenerationMode;539 const double* Bounds;540 const double* Spacing;541 const int* Dims;542 int* SliceOffsets;543 std::atomic<TIds>* PointMap;544 PointsT* InPoints;545 ArrayList* Arrays;546 float* OutPoints;547 548 GenerateBinPoints(int genMode, const double* bounds, const double* spacing, const int* dims,549 int* sliceOffsets, std::atomic<TIds>* ptMap, PointsT* inPts, ArrayList* arrays, float* outPts)550 : PointGenerationMode(genMode)551 , Bounds(bounds)552 , Spacing(spacing)553 , Dims(dims)554 , SliceOffsets(sliceOffsets)555 , PointMap(ptMap)556 , InPoints(inPts)557 , Arrays(arrays)558 , OutPoints(outPts)559 {560 }561 562 void operator()(vtkIdType slice, vtkIdType endSlice)563 {564 int binOffset = slice * this->Dims[0] * this->Dims[1];565 vtkIdType oldPtId;566 vtkIdType newPtId = this->SliceOffsets[slice];567 float* xOut;568 double xIn[3];569 const auto pts = vtk::DataArrayTupleRange<3>(this->InPoints);570 571 for (; slice < endSlice; ++slice)572 {573 for (auto j = 0; j < this->Dims[1]; ++j)574 {575 for (auto i = 0; i < this->Dims[0]; ++i)576 {577 oldPtId = this->PointMap[binOffset];578 if (oldPtId != 0)579 {580 oldPtId = -(oldPtId + 1); // transform back to non-negative point id581 xOut = this->OutPoints + 3 * newPtId;582 if (this->PointGenerationMode == vtkBinnedDecimation::BIN_CENTERS)583 {584 xIn[0] = this->Bounds[0] + ((0.5 + static_cast<double>(i)) * this->Spacing[0]);585 xIn[1] = this->Bounds[2] + ((0.5 + static_cast<double>(j)) * this->Spacing[1]);586 xIn[2] = this->Bounds[4] + ((0.5 + static_cast<double>(slice)) * this->Spacing[2]);587 }588 else // genMode == vtkBinnedDecimation::BIN_POINTS589 {590 const auto xp = pts[oldPtId];591 xIn[0] = xp[0];592 xIn[1] = xp[1];593 xIn[2] = xp[2];594 }595 xOut[0] = xIn[0];596 xOut[1] = xIn[1];597 xOut[2] = xIn[2];598 this->PointMap[binOffset] = newPtId; // update to new point id599 if (this->Arrays) // copy point data if requested600 {601 this->Arrays->Copy(oldPtId, newPtId);602 }603 newPtId++;604 }605 ++binOffset;606 }607 }608 } // for all slices in this batch609 }610};611 612//----------------------------------------------------------------------------613// Decimation algorithms #2-3 which generates new points for each bin. Either614// a bin center point is generated, or one of the points contained in the bin615// is selected and copied to the output.616template <typename PointsT, typename TIds>617void BinPointsDecimate(int genMode, vtkIdType numPts, PointsT* pts, vtkPointData* inPD,618 vtkPointData* outPD, vtkIdType numTris, vtkCellArray* tris, vtkCellData* inCD, vtkCellData* outCD,619 vtkIdType numBins, const int* dims, const double* bounds, const double* spacing,620 vtkPolyData* output)621{622 // Setup execution. Several arrays are used to transform the data.623 // The bin id of each point.624 TIds* binIds = new TIds[numPts];625 626 // Now bin points to generate a bin index for each point.627 BinPoints<PointsT, TIds> binPoints(pts, binIds, dims, bounds, spacing);628 vtkSMPTools::For(0, numPts, binPoints);629 630 // The ptMap is the output point id assigned to each bin (if the bin631 // contains an output point). Note that multiple, simultaneous writes can632 // occur to a bin hence the use of atomics. Initialize to zero via {}. Zero633 // is a problem because a ptId can == zero; as a workaround, we'll634 // initially use negative ids, and convert to positive ids in the635 // final composition.636 std::atomic<TIds>* ptMap = new std::atomic<TIds>[numBins] {};637 638 // Is the triangle output? And eventually the offset into the output cell array.639 TIds* triMap = new TIds[numTris + 1];640 641 // Begin to construct mappings of input points and cells, to output points642 // and cells.643 MapOutput<TIds> mapOutput(binIds, ptMap, tris, triMap);644 vtkSMPTools::For(0, numTris, mapOutput);645 646 // Now generate the new points. First generate new point ids, and then647 // produce the actual points.648 int* sliceOffsets = new int[dims[2] + 1];649 CountPoints<TIds> countPts(dims, ptMap, sliceOffsets);650 vtkSMPTools::For(0, dims[2], countPts);651 int numNewPts = sliceOffsets[dims[2]];652 653 vtkNew<vtkPoints> newPts;654 newPts->SetDataType(VTK_FLOAT); // could be the same type as the input point type655 newPts->SetNumberOfPoints(numNewPts);656 ArrayList ptArrays;657 if (outPD) // copy point data if requested658 {659 outPD->CopyAllocate(inPD, numNewPts);660 ptArrays.AddArrays(numNewPts, inPD, outPD);661 }662 663 GenerateBinPoints<PointsT, TIds> genPts(genMode, bounds, spacing, dims, sliceOffsets, ptMap, pts,664 (outPD != nullptr ? (&ptArrays) : nullptr),665 vtkFloatArray::FastDownCast(newPts->GetData())->GetPointer(0));666 vtkSMPTools::For(0, dims[2], genPts);667 output->SetPoints(newPts);668 669 // Create a mapping of the input triangles to the output triangles.670 vtkIdType mark, numOutTris = 0;671 for (vtkIdType triId = 0; triId < numTris; ++triId)672 {673 mark = triMap[triId];674 triMap[triId] = numOutTris;675 numOutTris += mark;676 }677 triMap[numTris] = numOutTris;678 679 // Produce the decimated output. We'll directly create the offset680 // and connectivity arrays for the output polydata.681 vtkCellArray* outTrisArray = output->GetPolys();682 vtkNew<vtkIdTypeArray> outConn;683 vtkIdType* outTris = outConn->WritePointer(0, numOutTris * 3);684 vtkNew<vtkIdTypeArray> outOffsets;685 vtkIdType* outTriOffsets = outOffsets->WritePointer(0, numOutTris + 1);686 outTriOffsets[numOutTris] = 3 * numOutTris;687 688 ArrayList arrays;689 if (outCD) // copy cell data if requested690 {691 outCD->CopyAllocate(inCD, numOutTris);692 arrays.AddArrays(numOutTris, inCD, outCD);693 }694 695 GenerateTriangles<TIds, std::atomic<TIds>> genTris(696 binIds, ptMap, tris, triMap, outTris, outTriOffsets, (outCD != nullptr ? (&arrays) : nullptr));697 vtkSMPTools::For(0, numTris, genTris);698 outTrisArray->SetData(outOffsets, outConn);699 700 // Clean up701 delete[] sliceOffsets;702 delete[] triMap;703 delete[] ptMap;704 delete[] binIds;705}706 707//----------------------------------------------------------------------------708// Invoke the decimate algorithm which generates either a selected bin point,709// or bin centered point. Depending on the id size, use large 64-bit ids or710// 32-bit ids (enhances performance and reduces memory usage).711struct BinPointsWorker712{713 template <typename DataT>714 void operator()(DataT* pts, vtkPointData* inPD, vtkPointData* outPD, bool largeIds, int genMode,715 vtkCellArray* tris, vtkCellData* inCD, vtkCellData* outCD, int* divs, double bounds[6],716 double spacing[3], vtkPolyData* output)717 {718 vtkIdType numPts = pts->GetNumberOfTuples();719 vtkIdType numTris = tris->GetNumberOfCells();720 vtkIdType numBins = divs[0] * divs[1] * divs[2];721 722 if (!largeIds)723 {724 BinPointsDecimate<DataT, int>(genMode, numPts, pts, inPD, outPD, numTris, tris, inCD, outCD,725 numBins, divs, bounds, spacing, output);726 }727 else728 {729 BinPointsDecimate<DataT, vtkIdType>(genMode, numPts, pts, inPD, outPD, numTris, tris, inCD,730 outCD, numBins, divs, bounds, spacing, output);731 }732 }733};734 735//========================== 4) Generate points from BIN_AVERAGES ============736// Traverse cells and mark output triangles. Sort points based on the bins737// they fall into. We need to keep track of the inserted points in each bin,738// which are combined later to produce an average position.739 740//------------------------------------------------------------------------------741// Sort bins with associated bin ids. This creates runs of points for each742// bin, which are later averaged to create a new point in each bin.743template <typename TIds>744struct BinTuple745{746 TIds PtId; // originating point id747 TIds Bin; // i-j-k index into bin space748 749 // Operator< used to support the sort operation. Just sort on bin750 // id; points within a bin can be in any order.751 bool operator<(const BinTuple& tuple) const { return (Bin < tuple.Bin ? true : false); }752};753 754template <typename PointsT, typename TIds>755struct BinPointTuples : public BinPoints<PointsT, TIds>756{757 BinTuple<TIds>* BinTuples;758 759 BinPointTuples(PointsT* pts, BinTuple<TIds>* binTuples, const int* dims, const double* bounds,760 const double* spacing)761 : BinPoints<PointsT, TIds>(pts, nullptr, dims, bounds, spacing)762 , BinTuples(binTuples)763 {764 }765 766 void operator()(vtkIdType ptId, vtkIdType endPtId)767 {768 const auto points = vtk::DataArrayTupleRange<3>(this->Points, ptId, endPtId);769 double x[3];770 BinTuple<TIds>* bins = this->BinTuples + ptId;771 for (const auto tuple : points)772 {773 (*bins).PtId = ptId++;774 x[0] = static_cast<double>(tuple[0]);775 x[1] = static_cast<double>(tuple[1]);776 x[2] = static_cast<double>(tuple[2]);777 (*bins).Bin = this->GetBinIndex(x);778 ++bins;779 }780 }781};782 783template <typename TIds>784struct MarkBinnedTris785{786 const BinTuple<TIds>* BinTuples;787 vtkCellArray* Tris;788 TIds* TriMap;789 vtkSMPThreadLocal<vtkSmartPointer<vtkCellArrayIterator>> CellIterator;790 791 MarkBinnedTris(BinTuple<TIds>* bt, vtkCellArray* tris, TIds* triMap)792 : BinTuples(bt)793 , Tris(tris)794 , TriMap(triMap)795 {796 }797 798 void Initialize() { this->CellIterator.Local().TakeReference(this->Tris->NewIterator()); }799 800 void operator()(vtkIdType triId, vtkIdType endTriId)801 {802 vtkIdType npts;803 const vtkIdType* tri;804 vtkCellArrayIterator* cellIter = this->CellIterator.Local();805 TIds* triMap = this->TriMap + triId;806 TIds binIds[3];807 808 for (; triId < endTriId; ++triId, ++triMap)809 {810 cellIter->GetCellAtId(triId, npts, tri);811 binIds[0] = this->BinTuples[tri[0]].Bin;812 binIds[1] = this->BinTuples[tri[1]].Bin;813 binIds[2] = this->BinTuples[tri[2]].Bin;814 815 // All three points have to be in different bins for triangle insertion816 if (binIds[0] != binIds[1] && binIds[0] != binIds[2] && binIds[1] != binIds[2])817 {818 *triMap = 1;819 }820 else // mark excluded from output821 {822 *triMap = 0;823 }824 }825 }826 827 void Reduce() {}828};829 830// Produce the output triangles from the bin tuples.831template <typename TIds>832struct BinAveTriangles833{834 const BinTuple<TIds>* BinTuples;835 vtkCellArray* Tris;836 vtkSMPThreadLocal<vtkSmartPointer<vtkCellArrayIterator>> CellIterator;837 const TIds* TriMap;838 vtkIdType* OutTris;839 vtkIdType* OutTriOffsets;840 ArrayList* Arrays;841 842 BinAveTriangles(const BinTuple<TIds>* bt, vtkCellArray* tris, TIds* triMap, vtkIdType* outTris,843 vtkIdType* outTriOffsets, ArrayList* arrays)844 : BinTuples(bt)845 , Tris(tris)846 , TriMap(triMap)847 , OutTris(outTris)848 , OutTriOffsets(outTriOffsets)849 , Arrays(arrays)850 {851 }852 853 void Initialize() { this->CellIterator.Local().TakeReference(this->Tris->NewIterator()); }854 855 void operator()(vtkIdType triId, vtkIdType endTriId)856 {857 const BinTuple<TIds>* binTuples = this->BinTuples;858 vtkIdType npts;859 const vtkIdType* tri;860 vtkCellArrayIterator* cellIter = this->CellIterator.Local();861 const TIds* triMap = this->TriMap;862 vtkIdType *outTris = this->OutTris, *outTri;863 vtkIdType *outTriOffsets = this->OutTriOffsets, *outOffsets;864 865 for (; triId < endTriId; ++triId)866 {867 if ((triMap[triId + 1] - triMap[triId]) > 0) // spit out triangle868 {869 cellIter->GetCellAtId(triId, npts, tri);870 871 outOffsets = outTriOffsets + triMap[triId];872 *outOffsets = triMap[triId] * 3;873 outTri = outTris + *outOffsets;874 875 outTri[0] = binTuples[tri[0]].Bin; // set the bin id876 outTri[1] = binTuples[tri[1]].Bin;877 outTri[2] = binTuples[tri[2]].Bin;878 if (this->Arrays) // copy cell data if requested879 {880 this->Arrays->Copy(triId, triMap[triId]);881 }882 }883 }884 }885 886 void Reduce() {}887};888 889// Generate the output triangles by rewriting bin ids into new point ids890template <typename TIds>891struct GenerateAveTriangles892{893 const BinTuple<TIds>* BinTuples;894 const TIds* Offsets;895 vtkIdType* OutTris;896 897 GenerateAveTriangles(const BinTuple<TIds>* bt, TIds* offsets, vtkIdType* outTris)898 : BinTuples(bt)899 , Offsets(offsets)900 , OutTris(outTris)901 {902 }903 904 void operator()(vtkIdType triId, vtkIdType endTriId)905 {906 const BinTuple<TIds>* binTuples = this->BinTuples;907 const TIds* offsets = this->Offsets;908 vtkIdType* outTri = this->OutTris + 3 * triId;909 910 for (; triId < endTriId; ++triId, outTri += 3)911 {912 outTri[0] = (*(binTuples + offsets[outTri[0]])).PtId;913 outTri[1] = (*(binTuples + offsets[outTri[1]])).PtId;914 outTri[2] = (*(binTuples + offsets[outTri[2]])).PtId;915 }916 }917};918 919// A clever way to build offsets in parallel. Basically each thread builds920// offsets across a range of the sorted map.921template <typename TIds>922struct MapOffsets923{924 const BinTuple<TIds>* BinTuples;925 TIds* Offsets;926 TIds NumPts;927 TIds NumBins;928 TIds BatchSize;929 930 MapOffsets(BinTuple<TIds>* bt, TIds* offsets, TIds numPts, TIds numBins, TIds numBatches)931 : BinTuples(bt)932 , Offsets(offsets)933 , NumPts(numPts)934 , NumBins(numBins)935 {936 this->BatchSize = static_cast<int>(ceil(static_cast<double>(numPts) / numBatches));937 }938 939 // Traverse sorted points (i.e., tuples) and update bin offsets.940 void operator()(vtkIdType batch, vtkIdType batchEnd)941 {942 TIds* offsets = this->Offsets;943 const BinTuple<TIds>* curPt = this->BinTuples + batch * this->BatchSize;944 const BinTuple<TIds>* endBatchPt = this->BinTuples + batchEnd * this->BatchSize;945 const BinTuple<TIds>* endPt = this->BinTuples + this->NumPts;946 const BinTuple<TIds>* prevPt;947 endBatchPt = (endBatchPt > endPt ? endPt : endBatchPt);948 949 // Special case at the very beginning of the bin tuples array. If950 // the first point is in bin# N, then all bins up and including951 // N must refer to the first point.952 if (curPt == this->BinTuples)953 {954 prevPt = this->BinTuples;955 std::fill_n(offsets, curPt->Bin + 1, 0); // offset to the first points956 } // at the very beginning of the map (sorted points array)957 958 // We are entering this functor somewhere in the interior of the959 // mapped points array. All we need to do is point to the entry960 // position because we are interested only in prevPt->Bin.961 else962 {963 prevPt = curPt;964 } // else in the middle of a batch965 966 // Okay we have a starting point for a bin run. Now we can begin967 // filling in the offsets in this batch. A previous thread should968 // have/will have completed the previous and subsequent runs outside969 // of the [batch,batchEnd) range970 for (curPt = prevPt; curPt < endBatchPt;)971 {972 for (; curPt->Bin == prevPt->Bin && curPt <= endBatchPt; ++curPt)973 {974 ; // advance975 }976 // Fill in any gaps in the offset array977 if (curPt < endPt) // still within range of points978 {979 std::fill_n(offsets + prevPt->Bin + 1, curPt->Bin - prevPt->Bin, curPt - this->BinTuples);980 prevPt = curPt;981 }982 else // at the end of the points983 {984 std::fill_n(985 offsets + prevPt->Bin + 1, this->NumBins - prevPt->Bin - 1, curPt - this->BinTuples);986 }987 } // for all batches in this range988 } // operator()989};990 991// Count the number of averaged points in each z-slice of the binning volume.992// The resulting offsets are used to thread the generation of the output993// points.994template <typename TIds>995struct CountAvePts996{997 const int* Dims;998 const TIds* Offsets;999 int* SliceOffsets;1000 1001 CountAvePts(const int* dims, const TIds* offsets, int* sliceOffsets)1002 : Dims(dims)1003 , Offsets(offsets)1004 , SliceOffsets(sliceOffsets)1005 {1006 }1007 1008 void Initialize() {}1009 1010 void operator()(vtkIdType slice, vtkIdType endSlice)1011 {1012 int binNum = slice * this->Dims[0] * this->Dims[1];1013 1014 for (; slice < endSlice; ++slice)1015 {1016 vtkIdType numSlicePts = 0;1017 for (auto j = 0; j < this->Dims[1]; ++j)1018 {1019 for (auto i = 0; i < this->Dims[0]; ++i)1020 {1021 if ((this->Offsets[binNum + 1] - this->Offsets[binNum]) > 0)1022 {1023 ++numSlicePts;1024 }1025 ++binNum;1026 }1027 }1028 this->SliceOffsets[slice] = numSlicePts;1029 } // for all slices in this batch1030 }1031 1032 void Reduce()1033 {1034 // Prefix sum to roll up total point count in each slice1035 TIds numSlicePts, numNewPts = 0;1036 for (auto i = 0; i < this->Dims[2]; ++i)1037 {1038 numSlicePts = this->SliceOffsets[i];1039 this->SliceOffsets[i] = numNewPts;1040 numNewPts += numSlicePts;1041 }1042 this->SliceOffsets[this->Dims[2]] = numNewPts;1043 }1044};1045 1046// Generate points from the binning -- in this case from the average1047// position of all points in each bin.1048template <typename PointsT, typename TIds>1049struct GenerateAveBinPoints1050{1051 const int* Dims;1052 PointsT* InPoints;1053 const int* SliceOffsets;1054 BinTuple<TIds>* BinTuples;1055 const TIds* Offsets;1056 ArrayList* Arrays;1057 float* OutPoints;1058 vtkSMPThreadLocal<std::vector<vtkIdType>> PtIds;1059 1060 GenerateAveBinPoints(const int* dims, PointsT* inPts, const int* sliceOffsets,1061 BinTuple<TIds>* binTuples, const TIds* offsets, ArrayList* arrays, float* outPts)1062 : Dims(dims)1063 , InPoints(inPts)1064 , SliceOffsets(sliceOffsets)1065 , BinTuples(binTuples)1066 , Offsets(offsets)1067 , Arrays(arrays)1068 , OutPoints(outPts)1069 {1070 }1071 1072 void operator()(vtkIdType slice, vtkIdType endSlice)1073 {1074 int binNum = slice * this->Dims[0] * this->Dims[1];1075 vtkIdType newPtId = this->SliceOffsets[slice];1076 float* xOut;1077 double xAve[3];1078 const auto pts = vtk::DataArrayTupleRange<3>(this->InPoints);1079 BinTuple<TIds>* binTuples = this->BinTuples;1080 const TIds* offsets = this->Offsets;1081 BinTuple<TIds>* pIds;1082 TIds pId;1083 std::vector<vtkIdType> v = this->PtIds.Local();1084 1085 for (; slice < endSlice; ++slice)1086 {1087 for (auto j = 0; j < this->Dims[1]; ++j)1088 {1089 for (auto i = 0; i < this->Dims[0]; ++i)1090 {1091 TIds npts = offsets[binNum + 1] - offsets[binNum];1092 if (npts > 0)1093 {1094 // Average the points in the bin1095 xAve[0] = xAve[1] = xAve[2] = 0.0;1096 pIds = binTuples + offsets[binNum];1097 v.resize(npts);1098 for (auto idx = 0; idx < npts; ++idx)1099 {1100 pId = (*(pIds + idx)).PtId;1101 v[idx] = pId;1102 const auto p = pts[pId];1103 xAve[0] += p[0];1104 xAve[1] += p[1];1105 xAve[2] += p[2];1106 }1107 xAve[0] /= static_cast<double>(npts);1108 xAve[1] /= static_cast<double>(npts);1109 xAve[2] /= static_cast<double>(npts);1110 1111 xOut = this->OutPoints + 3 * newPtId;1112 xOut[0] = xAve[0];1113 xOut[1] = xAve[1];1114 xOut[2] = xAve[2];1115 if (this->Arrays) // average point data if requested1116 {1117 this->Arrays->Average(npts, v.data(), newPtId);1118 }1119 (*pIds).PtId = newPtId; // update to new point id1120 newPtId++;1121 }1122 ++binNum;1123 } // for i1124 } // for j1125 } // for all slices in this batch1126 }1127};1128 1129//----------------------------------------------------------------------------1130// Decimation algorithm which generates new points for each bin by averaging1131// the point coordinates and point attributes in each bin. This algorithm1132// typically produces the best results. For small numbers of bins it can be1133// a tad slower, but generally scales better as the number of bins is1134// increased.1135template <typename PointsT, typename TIds>1136void AvePointsDecimate(vtkIdType numPts, PointsT* pts, vtkPointData* inPD, vtkPointData* outPD,1137 vtkIdType numTris, vtkCellArray* tris, vtkCellData* inCD, vtkCellData* outCD, vtkIdType numBins,1138 const int* dims, const double* bounds, const double* spacing, vtkPolyData* output)1139{1140 // Setup execution. Several arrays are used to transform the data.1141 // Define the bin id and associated point id of each point.1142 BinTuple<TIds>* binTuples = new BinTuple<TIds>[numPts];1143 1144 // Now bin points to generate a bin index for each point.1145 BinPointTuples<PointsT, TIds> binPoints(pts, binTuples, dims, bounds, spacing);1146 vtkSMPTools::For(0, numPts, binPoints);1147 1148 // Initially, triMap indicates which triangles are output. And then1149 // contains the offsets into the output triangles array.1150 TIds* triMap = new TIds[numTris + 1];1151 1152 // Begin to construct mappings of input points and cells, to output points1153 // and cells. First identify the triangles to be sent to the output.1154 MarkBinnedTris<TIds> markBinnedTris(binTuples, tris, triMap);1155 vtkSMPTools::For(0, numTris, markBinnedTris);1156 1157 // Create a mapping of the input triangles to the output triangles.1158 vtkIdType mark, numOutTris = 0;1159 for (vtkIdType triId = 0; triId < numTris; ++triId)1160 {1161 mark = triMap[triId];1162 triMap[triId] = numOutTris;1163 numOutTris += mark;1164 }1165 triMap[numTris] = numOutTris;1166 1167 // Generate the cell output (decimated list of triangles), with the1168 // triangle connectivity based on bin ids (not point ids). We'll directly1169 // create the offet and connectivity arrays for the output polydata.1170 vtkCellArray* outTrisArray = output->GetPolys();1171 vtkNew<vtkIdTypeArray> outConn;1172 vtkIdType* outTris = outConn->WritePointer(0, numOutTris * 3);1173 vtkNew<vtkIdTypeArray> outOffsets;1174 vtkIdType* outTriOffsets = outOffsets->WritePointer(0, numOutTris + 1);1175 outTriOffsets[numOutTris] = 3 * numOutTris;1176 1177 ArrayList arrays;1178 if (outCD) // copy cell data if requested1179 {1180 outCD->CopyAllocate(inCD, numOutTris);1181 arrays.AddArrays(numOutTris, inCD, outCD);1182 }1183 1184 BinAveTriangles<TIds> binTris(1185 binTuples, tris, triMap, outTris, outTriOffsets, (outCD != nullptr ? (&arrays) : nullptr));1186 vtkSMPTools::For(0, numTris, binTris);1187 outTrisArray->SetData(outOffsets, outConn);1188 1189 // Now sort the bin tuples by bin id. This is the first step in1190 // transforming the input points to new points. The sort operation1191 // will organize points into lists per bin.1192 vtkSMPTools::Sort(binTuples, binTuples + numPts);1193 1194 // To rapidly (random) access the bins, we have to build an offset array1195 // into the sorted bin tuples (i.e., runs of points in each bin).1196 TIds* offsets = new TIds[numBins + 1];1197 TIds numBatches = (numPts < 10000 ? 1 : 100); // totally arbitrary1198 MapOffsets<TIds> offMapper(binTuples, offsets, numPts, numBins, numBatches);1199 vtkSMPTools::For(0, numBatches, offMapper);1200 offsets[numBins] = numPts;