// SPDX-FileCopyrightText: Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen // SPDX-License-Identifier: BSD-3-Clause #include "vtkClipDataSet.h" #include "vtkCallbackCommand.h" #include "vtkCellArray.h" #include "vtkCellData.h" #include "vtkClipVolume.h" #include "vtkExecutive.h" #include "vtkFloatArray.h" #include "vtkGenericCell.h" #include "vtkImageData.h" #include "vtkImplicitFunction.h" #include "vtkIncrementalPointLocator.h" #include "vtkInformation.h" #include "vtkInformationVector.h" #include "vtkMergePoints.h" #include "vtkNonLinearCell.h" #include "vtkObjectFactory.h" #include "vtkPointData.h" #include "vtkPolyhedron.h" #include "vtkSmartPointer.h" #include "vtkUnsignedCharArray.h" #include "vtkUnstructuredGrid.h" VTK_ABI_NAMESPACE_BEGIN vtkStandardNewMacro(vtkClipDataSet); vtkCxxSetObjectMacro(vtkClipDataSet, ClipFunction, vtkImplicitFunction); //------------------------------------------------------------------------------ // Construct with user-specified implicit function; InsideOut turned off; value // set to 0.0; and generate clip scalars turned off. vtkClipDataSet::vtkClipDataSet(vtkImplicitFunction* cf) { this->ClipFunction = cf; this->InsideOut = 0; this->Locator = nullptr; this->Value = 0.0; this->UseValueAsOffset = true; this->GenerateClipScalars = 0; this->OutputPointsPrecision = DEFAULT_PRECISION; this->GenerateClippedOutput = 0; this->MergeTolerance = 0.01; this->SetNumberOfOutputPorts(2); vtkUnstructuredGrid* output2 = vtkUnstructuredGrid::New(); this->GetExecutive()->SetOutputData(1, output2); output2->Delete(); // by default process active point scalars this->SetInputArrayToProcess( 0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS, vtkDataSetAttributes::SCALARS); // Setup a callback for the internal readers to report progress. this->InternalProgressObserver = vtkCallbackCommand::New(); this->InternalProgressObserver->SetCallback(&vtkClipDataSet::InternalProgressCallbackFunction); this->InternalProgressObserver->SetClientData(this); } //------------------------------------------------------------------------------ vtkClipDataSet::~vtkClipDataSet() { if (this->Locator) { this->Locator->UnRegister(this); this->Locator = nullptr; } this->SetClipFunction(nullptr); this->InternalProgressObserver->Delete(); } //------------------------------------------------------------------------------ void vtkClipDataSet::InternalProgressCallbackFunction( vtkObject* arg, unsigned long, void* clientdata, void*) { reinterpret_cast(clientdata) ->InternalProgressCallback(static_cast(arg)); } //------------------------------------------------------------------------------ void vtkClipDataSet::InternalProgressCallback(vtkAlgorithm* algorithm) { float progress = algorithm->GetProgress(); this->UpdateProgress(progress); } //------------------------------------------------------------------------------ // Overload standard modified time function. If Clip functions is modified, // then this object is modified as well. vtkMTimeType vtkClipDataSet::GetMTime() { vtkMTimeType mTime = this->Superclass::GetMTime(); vtkMTimeType time; if (this->ClipFunction != nullptr) { time = this->ClipFunction->GetMTime(); mTime = (time > mTime ? time : mTime); } if (this->Locator != nullptr) { time = this->Locator->GetMTime(); mTime = (time > mTime ? time : mTime); } return mTime; } //------------------------------------------------------------------------------ vtkUnstructuredGrid* vtkClipDataSet::GetClippedOutput() { if (!this->GenerateClippedOutput) { return nullptr; } return vtkUnstructuredGrid::SafeDownCast(this->GetExecutive()->GetOutputData(1)); } //------------------------------------------------------------------------------ // // Clip through data generating surface. // int vtkClipDataSet::RequestData(vtkInformation* vtkNotUsed(request), vtkInformationVector** inputVector, vtkInformationVector* outputVector) { // get the info objects vtkInformation* inInfo = inputVector[0]->GetInformationObject(0); vtkInformation* outInfo = outputVector->GetInformationObject(0); // get the input and output vtkDataSet* realInput = vtkDataSet::SafeDownCast(inInfo->Get(vtkDataObject::DATA_OBJECT())); // We have to create a copy of the input because clip requires being // able to InterpolateAllocate point data from the input that is // exactly the same as output. If the input arrays and output arrays // are different vtkCell3D's Clip will fail. By calling InterpolateAllocate // here, we make sure that the output will look exactly like the output // (unwanted arrays will be eliminated in InterpolateAllocate). The // last argument of InterpolateAllocate makes sure that arrays are shallow // copied from realInput to input. vtkSmartPointer input; input.TakeReference(realInput->NewInstance()); input->CopyStructure(realInput); input->GetCellData()->PassData(realInput->GetCellData()); input->GetPointData()->InterpolateAllocate(realInput->GetPointData(), 0, 0, 1); vtkUnstructuredGrid* output = vtkUnstructuredGrid::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT())); vtkUnstructuredGrid* clippedOutput = this->GetClippedOutput(); vtkIdType numPts = input->GetNumberOfPoints(); vtkIdType numCells = input->GetNumberOfCells(); vtkPointData *inPD = input->GetPointData(), *outPD = output->GetPointData(); vtkCellData* inCD = input->GetCellData(); vtkCellData* outCD[2]; vtkSmartPointer newPoints; vtkSmartPointer cellScalars; vtkDataArray* clipScalars; vtkPoints* cellPts; vtkIdList* cellIds; double s; vtkIdType npts; const vtkIdType* pts; vtkIdType i; int j; vtkIdType estimatedSize; vtkSmartPointer types[2]; types[0] = types[1] = nullptr; int numOutputs = 1; outCD[0] = nullptr; outCD[1] = nullptr; vtkDebugMacro(<< "Clipping dataset"); int inputObjectType = input->GetDataObjectType(); // if we have volumes if (inputObjectType == VTK_STRUCTURED_POINTS || inputObjectType == VTK_IMAGE_DATA) { int dimension; int* dims = vtkImageData::SafeDownCast(input)->GetDimensions(); for (dimension = 3, i = 0; i < 3; i++) { if (dims[i] <= 1) { dimension--; } } if (dimension >= 3) { this->ClipVolume(input, output); return 1; } } // Initialize self; create output objects // if (numPts < 1) { vtkDebugMacro(<< "No data to clip"); return 1; } if (!this->ClipFunction && this->GenerateClipScalars) { vtkErrorMacro(<< "Cannot generate clip scalars if no clip function defined"); return 1; } if (numCells < 1) { return this->ClipPoints(input, output, inputVector); } // allocate the output and associated helper classes estimatedSize = numCells; estimatedSize = estimatedSize / 1024 * 1024; // multiple of 1024 estimatedSize = std::max(estimatedSize, 1024); cellScalars = vtkSmartPointer::New(); cellScalars->Allocate(VTK_CELL_SIZE); vtkSmartPointer conn[2]; conn[0] = conn[1] = nullptr; conn[0] = vtkSmartPointer::New(); conn[0]->AllocateEstimate(estimatedSize, 1); conn[0]->InitTraversal(); types[0] = vtkSmartPointer::New(); types[0]->Allocate(estimatedSize, estimatedSize / 2); if (this->GenerateClippedOutput) { numOutputs = 2; conn[1] = vtkSmartPointer::New(); conn[1]->AllocateEstimate(estimatedSize, 1); conn[1]->InitTraversal(); types[1] = vtkSmartPointer::New(); types[1]->Allocate(estimatedSize, estimatedSize / 2); } newPoints = vtkSmartPointer::New(); // set precision for the points in the output if (this->OutputPointsPrecision == vtkAlgorithm::DEFAULT_PRECISION) { vtkPointSet* inputPointSet = vtkPointSet::SafeDownCast(input); if (inputPointSet) { newPoints->SetDataType(inputPointSet->GetPoints()->GetDataType()); } else { newPoints->SetDataType(VTK_FLOAT); } } else if (this->OutputPointsPrecision == vtkAlgorithm::SINGLE_PRECISION) { newPoints->SetDataType(VTK_FLOAT); } else if (this->OutputPointsPrecision == vtkAlgorithm::DOUBLE_PRECISION) { newPoints->SetDataType(VTK_DOUBLE); } newPoints->Allocate(numPts, numPts / 2); // locator used to merge potentially duplicate points if (this->Locator == nullptr) { this->CreateDefaultLocator(); } // Determine whether we're clipping with input scalars or a clip function // and do necessary setup. if (this->ClipFunction) { vtkFloatArray* tmpScalars = vtkFloatArray::New(); tmpScalars->SetNumberOfTuples(numPts); tmpScalars->SetName("ClipDataSetScalars"); inPD = vtkPointData::New(); inPD->ShallowCopy(input->GetPointData()); // copies original if (this->GenerateClipScalars) { inPD->SetScalars(tmpScalars); } double pt[3]; for (i = 0; i < numPts; i++) { input->GetPoint(i, pt); tmpScalars->SetValue(i, this->ClipFunction->FunctionValue(pt)); } clipScalars = tmpScalars; } else // using input scalars { clipScalars = this->GetInputArrayToProcess(0, inputVector); if (!clipScalars) { // When processing composite datasets with partial arrays, this warning is // not applicable, hence disabling it. // vtkErrorMacro(<<"Cannot clip without clip function or input scalars"); return 1; } } // Refer to BUG #8494 and BUG #11016. I cannot see any reason why one would // want to turn CopyScalars Off. My understanding is that this was done to // avoid copying of "ClipDataSetScalars" to the output when // this->GenerateClipScalars is false. But, if GenerateClipScalars is false, // then "ClipDataSetScalars" is not added as scalars to the input at all // (refer to code above) so it's a non-issue. Leaving CopyScalars untouched // i.e. ON avoids dropping of arrays (#8484) as well as segfaults (#11016). // if ( !this->GenerateClipScalars && // !this->GetInputArrayToProcess(0,inputVector)) // { // outPD->CopyScalarsOff(); // } // else // { // outPD->CopyScalarsOn(); // } vtkDataSetAttributes* tempDSA = vtkDataSetAttributes::New(); tempDSA->InterpolateAllocate(inPD, 1, 2); outPD->InterpolateAllocate(inPD, estimatedSize, estimatedSize / 2); tempDSA->Delete(); outCD[0] = output->GetCellData(); outCD[0]->CopyAllocate(inCD, estimatedSize, estimatedSize / 2); if (this->GenerateClippedOutput) { outCD[1] = clippedOutput->GetCellData(); outCD[1]->CopyAllocate(inCD, estimatedSize, estimatedSize / 2); } std::vector pointMap(numPts, -1); std::vector clippingCellIds; std::array, 2> intactCellIds; intactCellIds[0].reserve(numPts / 2); if (this->GenerateClippedOutput) { intactCellIds[1].reserve(numPts / 2); } double clipValue = 0.0; if (this->UseValueAsOffset || !this->ClipFunction) { clipValue = this->Value; } // Identify inside/outside cells and cells at the border vtkIdType newPointsNumber = this->SeparateCellIds( numCells, input, pointMap, clippingCellIds, clipValue, clipScalars, intactCellIds); this->UpdateProgress(0.15); // Fill points and copy point data newPoints->SetNumberOfPoints(newPointsNumber); for (vtkIdType pointId = 0; pointId < numPts; ++pointId) { vtkIdType newId = pointMap[pointId]; if (newId >= 0) { double coords[3]; input->GetPoint(pointId, coords); newPoints->SetPoint(newId, coords); outPD->CopyData(inPD, pointId, newId); } } this->UpdateProgress(0.2); this->FillIntactCells(input, numOutputs, intactCellIds, pointMap, conn, outCD, inCD, types); this->UpdateProgress(0.4); this->Locator->InitPointInsertion(newPoints, input->GetBounds()); int traversalId[2]; traversalId[0] = conn[0]->GetNumberOfCells(); if (this->GenerateClippedOutput) { traversalId[1] = conn[1]->GetNumberOfCells(); } auto getGeneratedCellType = [](vtkGenericCell* gCell, vtkIdType nPts, bool isSameCell) { if (isSameCell) { return static_cast(gCell->GetCellType()); } else if (gCell->GetCellType() == VTK_POLYHEDRON) { return VTK_POLYHEDRON; } else { switch (gCell->GetCellDimension()) { case 0: // points are generated-------------------------------- return (nPts > 1 ? VTK_POLY_VERTEX : VTK_VERTEX); case 1: // lines are generated--------------------------------- return (nPts > 2 ? VTK_POLY_LINE : VTK_LINE); case 2: // polygons are generated------------------------------ // NOLINTNEXTLINE(readability-avoid-nested-conditional-operator) return (nPts == 3 ? VTK_TRIANGLE : (nPts == 4 ? VTK_QUAD : VTK_POLYGON)); case 3: // tetrahedra or wedges are generated------------------ return (nPts == 4 ? VTK_TETRA : VTK_WEDGE); default: vtkErrorWithObjectMacro(nullptr, "Dimension cannot be lower than 0 or higher than 3"); break; } } return VTK_EMPTY_CELL; }; // Process all cells and clip each in turn // bool abort = false; vtkIdType updateTime = static_cast(clippingCellIds.size()) / 20 + 1; // update roughly every 5% int num[2]; num[0] = conn[0]->GetNumberOfCells(); num[1] = 0; if (this->GenerateClippedOutput) { num[1] = conn[1]->GetNumberOfCells(); } int numNew[2]; numNew[0] = numNew[1] = 0; bool sameCell[2] = { false, false }; vtkSmartPointer cell = vtkSmartPointer::New(); for (vtkIdType cellId = 0; cellId < static_cast(clippingCellIds.size()) && !abort; ++cellId) { if (!(cellId % updateTime)) { this->UpdateProgress(0.4 + (static_cast(cellId) / clippingCellIds.size()) * 0.45); abort = this->CheckAbort(); } vtkIdType cid = clippingCellIds[cellId]; input->GetCell(cid, cell); cellPts = cell->GetPoints(); cellIds = cell->GetPointIds(); npts = cellPts->GetNumberOfPoints(); vtkNonLinearCell* nonLinearCell = vtkNonLinearCell::SafeDownCast(cell->GetRepresentativeCell()); // evaluate implicit cutting function for (i = 0; i < npts; i++) { s = clipScalars->GetComponent(cellIds->GetId(i), 0); cellScalars->InsertTuple(i, &s); } // perform the clipping for (i = 0; i < numOutputs; ++i) { if (this->StableClipNonLinear && nonLinearCell != nullptr) { sameCell[i] = nonLinearCell->StableClip(clipValue, cellScalars, this->Locator, conn[i], inPD, outPD, inCD, cid, outCD[i], this->InsideOut); numNew[i] = conn[i]->GetNumberOfCells() - num[i]; num[i] = conn[i]->GetNumberOfCells(); } else { cell->Clip(clipValue, cellScalars, this->Locator, conn[i], inPD, outPD, inCD, cid, outCD[i], this->InsideOut); numNew[i] = conn[i]->GetNumberOfCells() - num[i]; num[i] = conn[i]->GetNumberOfCells(); sameCell[i] = false; } } for (i = 0; i < numOutputs; i++) { for (j = 0; j < numNew[i]; j++) { conn[i]->GetCellAtId(traversalId[i], npts, pts); traversalId[i]++; types[i]->InsertNextValue(getGeneratedCellType(cell, npts, sameCell[i])); } } } if (this->ClipFunction) { clipScalars->Delete(); inPD->Delete(); } output->SetPoints(newPoints); output->SetCells(types[0], conn[0]); if (this->GenerateClippedOutput) { clippedOutput->SetPoints(newPoints); clippedOutput->SetCells(types[1], conn[1]); } this->UpdateProgress(1.0); this->Locator->Initialize(); // release any extra memory output->Squeeze(); return 1; } //------------------------------------------------------------------------------ vtkIdType vtkClipDataSet::SeparateCellIds(vtkIdType numCells, vtkDataSet* input, std::vector& pointMap, std::vector& clippingCellIds, double clipValue, vtkDataArray* clipScalars, std::array, 2>& intactCellIds) { vtkSmartPointer cellPointIds = vtkSmartPointer::New(); vtkIdType outputPointCounter = 0; auto registerCell = [&outputPointCounter, &pointMap, &cellPointIds]( vtkIdType cId, std::vector& cellList) { cellList.push_back(cId); for (vtkIdType idx = 0; idx < cellPointIds->GetNumberOfIds(); ++idx) { vtkIdType pid = cellPointIds->GetId(idx); if (pointMap[pid] < 0) { pointMap[pid] = outputPointCounter++; } } }; for (vtkIdType cellId = 0; cellId < numCells; ++cellId) { input->GetCellPoints(cellId, cellPointIds); int type = input->GetCellType(cellId); // Only polyhedrons are optimized for now // Other cell types should probably be considered if (type != VTK_POLYHEDRON) { clippingCellIds.push_back(cellId); continue; } bool isInside = false; bool isOutside = false; for (vtkIdType idx = 0; idx < cellPointIds->GetNumberOfIds(); ++idx) { double val = clipScalars->GetComponent(cellPointIds->GetId(idx), 0); bool inDomain = this->InsideOut ? (val <= clipValue) : (val >= clipValue); isInside |= inDomain; isOutside |= !inDomain; } if (isInside && !isOutside) { registerCell(cellId, intactCellIds[0]); } if (!isInside && isOutside && this->GenerateClippedOutput) { registerCell(cellId, intactCellIds[1]); } if (isInside && isOutside) { clippingCellIds.push_back(cellId); } } return outputPointCounter; } //------------------------------------------------------------------------------ void vtkClipDataSet::FillIntactCells(vtkDataSet* input, int numOutputs, std::array, 2>& intactCellIds, std::vector& pointMap, vtkSmartPointer conn[2], vtkCellData* outCD[2], vtkCellData* inCD, vtkSmartPointer types[2]) { auto ug = vtkUnstructuredGrid::SafeDownCast(input); vtkSmartPointer cell = vtkSmartPointer::New(); for (int i = 0; i < numOutputs; ++i) { for (vtkIdType idx = 0; idx < static_cast(intactCellIds[i].size()); ++idx) { vtkIdType cellId = intactCellIds[i][idx]; vtkNew faceStream; // This is safe because this loop is entered only if there are polyhedrons // which are necessary attached to vtkUnstructuredGrid ug->GetFaceStream(cellId, faceStream); // Update the new faceStream ids int faceIdCursor = 1; for (int _ = 0; _ < faceStream->GetId(0); ++_) { vtkIdType nFaces = faceStream->GetId(faceIdCursor); for (int j = 1; j < faceStream->GetId(faceIdCursor) + 1; ++j) { vtkIdType newId = pointMap[faceStream->GetId(faceIdCursor + j)]; faceStream->SetId(faceIdCursor + j, newId); } faceIdCursor += nFaces + 1; } vtkIdType newCellId = conn[i]->InsertNextCell(faceStream); outCD[i]->CopyData(inCD, cellId, newCellId); types[i]->InsertNextValue(VTK_POLYHEDRON); } } } //------------------------------------------------------------------------------ int vtkClipDataSet::ClipPoints( vtkDataSet* input, vtkUnstructuredGrid* output, vtkInformationVector** inputVector) { vtkPoints* outPoints = vtkPoints::New(); vtkPointData* inPD = input->GetPointData(); vtkPointData* outPD = output->GetPointData(); vtkIdType numPts = input->GetNumberOfPoints(); outPD->CopyAllocate(inPD, numPts / 2, numPts / 4); double value = 0.0; if (this->UseValueAsOffset || !this->ClipFunction) { value = this->Value; } if (this->ClipFunction) { double pt[3]; for (vtkIdType i = 0; i < numPts; i++) { input->GetPoint(i, pt); double fv = this->ClipFunction->FunctionValue(pt); int addPoint = 0; if (this->InsideOut) { if (fv <= value) { addPoint = 1; } } else { if (fv > value) { addPoint = 1; } } if (addPoint) { vtkIdType id = outPoints->InsertNextPoint(pt); outPD->CopyData(inPD, i, id); } } } else { vtkDataArray* clipScalars = this->GetInputArrayToProcess(0, inputVector); if (clipScalars) { double pt[3]; for (vtkIdType i = 0; i < numPts; i++) { int addPoint = 0; double fv; clipScalars->GetTuple(i, &fv); if (this->InsideOut) { if (fv <= value) { addPoint = 1; } } else { if (fv > value) { addPoint = 1; } } if (addPoint) { input->GetPoint(i, pt); vtkIdType id = outPoints->InsertNextPoint(pt); outPD->CopyData(inPD, i, id); } } } } output->SetPoints(outPoints); outPoints->Delete(); return 1; } //------------------------------------------------------------------------------ // Specify a spatial locator for merging points. By default, // an instance of vtkMergePoints is used. void vtkClipDataSet::SetLocator(vtkIncrementalPointLocator* locator) { if (this->Locator == locator) { return; } if (this->Locator) { this->Locator->UnRegister(this); this->Locator = nullptr; } if (locator) { locator->Register(this); } this->Locator = locator; this->Modified(); } //------------------------------------------------------------------------------ void vtkClipDataSet::CreateDefaultLocator() { if (this->Locator == nullptr) { this->Locator = vtkMergePoints::New(); this->Locator->Register(this); this->Locator->Delete(); } } //------------------------------------------------------------------------------ void vtkClipDataSet::ClipVolume(vtkDataSet* input, vtkUnstructuredGrid* output) { vtkClipVolume* clipVolume = vtkClipVolume::New(); clipVolume->AddObserver(vtkCommand::ProgressEvent, this->InternalProgressObserver); // We cannot set the input directly. This messes up the partitioning. // output->UpdateNumberOfPieces gets set to 1. vtkImageData* tmp = vtkImageData::New(); tmp->ShallowCopy(vtkImageData::SafeDownCast(input)); clipVolume->SetInputData(tmp); double value = 0.0; if (this->UseValueAsOffset || !this->ClipFunction) { value = this->Value; } clipVolume->SetValue(value); clipVolume->SetInsideOut(this->InsideOut); clipVolume->SetClipFunction(this->ClipFunction); clipVolume->SetGenerateClipScalars(this->GenerateClipScalars); clipVolume->SetGenerateClippedOutput(this->GenerateClippedOutput); clipVolume->SetMergeTolerance(this->MergeTolerance); clipVolume->SetDebug(this->Debug); clipVolume->SetInputArrayToProcess(0, this->GetInputArrayInformation(0)); clipVolume->SetContainerAlgorithm(this); clipVolume->Update(); clipVolume->RemoveObserver(this->InternalProgressObserver); vtkUnstructuredGrid* clipOutput = clipVolume->GetOutput(); output->CopyStructure(clipOutput); output->GetPointData()->ShallowCopy(clipOutput->GetPointData()); output->GetCellData()->ShallowCopy(clipOutput->GetCellData()); clipVolume->Delete(); tmp->Delete(); } //------------------------------------------------------------------------------ int vtkClipDataSet::FillInputPortInformation(int, vtkInformation* info) { info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataSet"); return 1; } //------------------------------------------------------------------------------ void vtkClipDataSet::PrintSelf(ostream& os, vtkIndent indent) { this->Superclass::PrintSelf(os, indent); os << indent << "Merge Tolerance: " << this->MergeTolerance << "\n"; if (this->ClipFunction) { os << indent << "Clip Function: " << this->ClipFunction << "\n"; } else { os << indent << "Clip Function: (none)\n"; } os << indent << "InsideOut: " << (this->InsideOut ? "On\n" : "Off\n"); os << indent << "Value: " << this->Value << "\n"; if (this->Locator) { os << indent << "Locator: " << this->Locator << "\n"; } else { os << indent << "Locator: (none)\n"; } os << indent << "Generate Clip Scalars: " << (this->GenerateClipScalars ? "On\n" : "Off\n"); os << indent << "Generate Clipped Output: " << (this->GenerateClippedOutput ? "On\n" : "Off\n"); os << indent << "UseValueAsOffset: " << (this->UseValueAsOffset ? "On\n" : "Off\n"); os << indent << "Precision of the output points: " << this->OutputPointsPrecision << "\n"; } VTK_ABI_NAMESPACE_END