/*========================================================================= Program: Visualization Toolkit Module: vtkExtractSelectedThresholds.cxx Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen All rights reserved. See Copyright.txt or http://www.kitware.com/Copyright.htm for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notice for more information. =========================================================================*/ #include "vtkExtractSelectedThresholds.h" #include "vtkDataSet.h" #include "vtkThreshold.h" #include "vtkIdList.h" #include "vtkIdTypeArray.h" #include "vtkInformation.h" #include "vtkInformationVector.h" #include "vtkObjectFactory.h" #include "vtkSelection.h" #include "vtkSelectionNode.h" #include "vtkUnstructuredGrid.h" #include "vtkCell.h" #include "vtkCellData.h" #include "vtkPointData.h" #include "vtkCellData.h" #include "vtkDoubleArray.h" #include "vtkSignedCharArray.h" vtkStandardNewMacro(vtkExtractSelectedThresholds); //---------------------------------------------------------------------------- vtkExtractSelectedThresholds::vtkExtractSelectedThresholds() { this->SetNumberOfInputPorts(2); } //---------------------------------------------------------------------------- vtkExtractSelectedThresholds::~vtkExtractSelectedThresholds() { } //---------------------------------------------------------------------------- int vtkExtractSelectedThresholds::RequestData( vtkInformation *vtkNotUsed(request), vtkInformationVector **inputVector, vtkInformationVector *outputVector) { // get the info objects vtkInformation *selInfo = inputVector[1]->GetInformationObject(0); vtkInformation *inInfo = inputVector[0]->GetInformationObject(0); vtkInformation *outInfo = outputVector->GetInformationObject(0); // verify the input, selection and ouptut vtkDataSet *input = vtkDataSet::SafeDownCast( inInfo->Get(vtkDataObject::DATA_OBJECT())); if (!input) { vtkErrorMacro(<<"No input specified"); return 0; } if ( ! selInfo ) { //When not given a selection, quietly select nothing. return 1; } if (input->GetNumberOfCells() == 0 && input->GetNumberOfPoints() == 0) { // empty input, nothing to do.. return 1; } vtkSelection *sel = vtkSelection::SafeDownCast( selInfo->Get(vtkDataObject::DATA_OBJECT())); vtkSelectionNode *node = 0; if (sel->GetNumberOfNodes() == 1) { node = sel->GetNode(0); } if (!node) { vtkErrorMacro("Selection must have a single node."); return 1; } if (!node->GetProperties()->Has(vtkSelectionNode::CONTENT_TYPE()) || node->GetProperties()->Get(vtkSelectionNode::CONTENT_TYPE()) != vtkSelectionNode::THRESHOLDS) { vtkErrorMacro("Missing or invalid CONTENT_TYPE."); return 1; } vtkDataSet *output = vtkDataSet::SafeDownCast( outInfo->Get(vtkDataObject::DATA_OBJECT())); vtkDebugMacro(<< "Extracting from dataset"); int thresholdByPointVals = 0; int fieldType = vtkSelectionNode::CELL; if (node->GetProperties()->Has(vtkSelectionNode::FIELD_TYPE())) { fieldType = node->GetProperties()->Get(vtkSelectionNode::FIELD_TYPE()); if (fieldType == vtkSelectionNode::POINT) { if (node->GetProperties()->Has(vtkSelectionNode::CONTAINING_CELLS())) { thresholdByPointVals = node->GetProperties()->Get(vtkSelectionNode::CONTAINING_CELLS()); } } } if (thresholdByPointVals || fieldType == vtkSelectionNode::CELL) { return this->ExtractCells(node, input, output, thresholdByPointVals); } if (fieldType == vtkSelectionNode::POINT) { return this->ExtractPoints(node, input, output); } return 1; } //---------------------------------------------------------------------------- int vtkExtractSelectedThresholds::ExtractCells( vtkSelectionNode *sel, vtkDataSet *input, vtkDataSet *output, int usePointScalars) { //find the values to threshold within vtkDataArray *lims = vtkDataArray::SafeDownCast(sel->GetSelectionList()); if (lims == NULL) { vtkErrorMacro(<<"No values to threshold with"); return 1; } //find out what array we are supposed to threshold in vtkDataArray *inScalars = NULL; bool use_ids = false; if (usePointScalars) { if (sel->GetSelectionList()->GetName()) { if (strcmp(sel->GetSelectionList()->GetName(), "vtkGlobalIds") == 0) { inScalars = input->GetPointData()->GetGlobalIds(); } else if (strcmp(sel->GetSelectionList()->GetName(), "vtkIndices") == 0) { use_ids = true; } else { inScalars = input->GetPointData()->GetArray( sel->GetSelectionList()->GetName()); } } else { inScalars = input->GetPointData()->GetScalars(); } } else { if (sel->GetSelectionList()->GetName()) { if (strcmp(sel->GetSelectionList()->GetName(), "vtkGlobalIds") == 0) { inScalars = input->GetCellData()->GetGlobalIds(); } else if (strcmp(sel->GetSelectionList()->GetName(), "vtkIndices") == 0) { use_ids = true; } else { inScalars = input->GetCellData()->GetArray( sel->GetSelectionList()->GetName()); } } else { inScalars = input->GetCellData()->GetScalars(); } } if (inScalars == NULL && !use_ids) { vtkErrorMacro("Could not figure out what array to threshold in."); return 1; } int inverse = 0; if (sel->GetProperties()->Has(vtkSelectionNode::INVERSE())) { inverse = sel->GetProperties()->Get(vtkSelectionNode::INVERSE()); } int passThrough = 0; if (this->PreserveTopology) { passThrough = 1; } int comp_no = 0; if (sel->GetProperties()->Has(vtkSelectionNode::COMPONENT_NUMBER())) { comp_no = sel->GetProperties()->Get(vtkSelectionNode::COMPONENT_NUMBER()); } vtkIdType cellId, newCellId; vtkIdList *cellPts, *pointMap = NULL; vtkIdList *newCellPts = NULL; vtkCell *cell = 0; vtkPoints *newPoints = 0; vtkIdType i, ptId, newId, numPts, numCells; int numCellPts; double x[3]; vtkPointData *pd=input->GetPointData(), *outPD=output->GetPointData(); vtkCellData *cd=input->GetCellData(), *outCD=output->GetCellData(); int keepCell; outPD->CopyGlobalIdsOn(); outPD->CopyAllocate(pd); outCD->CopyGlobalIdsOn(); outCD->CopyAllocate(cd); numPts = input->GetNumberOfPoints(); numCells = input->GetNumberOfCells(); vtkDataSet *outputDS = output; vtkSignedCharArray *pointInArray = NULL; vtkSignedCharArray *cellInArray = NULL; vtkUnstructuredGrid *outputUG = NULL; vtkIdTypeArray *originalCellIds = NULL; vtkIdTypeArray *originalPointIds = NULL; signed char flag = inverse ? 1 : -1; if (passThrough) { outputDS->ShallowCopy(input); pointInArray = vtkSignedCharArray::New(); pointInArray->SetNumberOfComponents(1); pointInArray->SetNumberOfTuples(numPts); for (i=0; i < numPts; i++) { pointInArray->SetValue(i, flag); } pointInArray->SetName("vtkInsidedness"); outPD->AddArray(pointInArray); outPD->SetScalars(pointInArray); cellInArray = vtkSignedCharArray::New(); cellInArray->SetNumberOfComponents(1); cellInArray->SetNumberOfTuples(numCells); for (i=0; i < numCells; i++) { cellInArray->SetValue(i, flag); } cellInArray->SetName("vtkInsidedness"); outCD->AddArray(cellInArray); outCD->SetScalars(cellInArray); } else { outputUG = vtkUnstructuredGrid::SafeDownCast(output); outputUG->Allocate(input->GetNumberOfCells()); newPoints = vtkPoints::New(); newPoints->Allocate(numPts); pointMap = vtkIdList::New(); //maps old point ids into new pointMap->SetNumberOfIds(numPts); for (i=0; i < numPts; i++) { pointMap->SetId(i,-1); } newCellPts = vtkIdList::New(); originalCellIds = vtkIdTypeArray::New(); originalCellIds->SetName("vtkOriginalCellIds"); originalCellIds->SetNumberOfComponents(1); outCD->AddArray(originalCellIds); originalPointIds = vtkIdTypeArray::New(); originalPointIds->SetName("vtkOriginalPointIds"); originalPointIds->SetNumberOfComponents(1); outPD->AddArray(originalPointIds); originalPointIds->Delete(); } flag = -flag; // Check that the scalars of each cell satisfy the threshold criterion for (cellId=0; cellId < input->GetNumberOfCells(); cellId++) { cell = input->GetCell(cellId); cellPts = cell->GetPointIds(); numCellPts = cell->GetNumberOfPoints(); // BUG: This code misses the case where the threshold is contained // completely within the cell but none of its points are inside // the range. Consider as an example the threshold range [1, 2] // with a cell [0, 3]. if ( usePointScalars ) { keepCell = 0; int totalAbove = 0; int totalBelow = 0; for ( i=0; (i < numCellPts) && (passThrough || !keepCell); i++) { int above = 0; int below = 0; ptId = cellPts->GetId(i); int inside = this->EvaluateValue( inScalars, comp_no, ptId, lims, &above, &below, NULL); totalAbove += above; totalBelow += below; // Have we detected a cell that straddles the threshold? if ((!inside) && (totalAbove && totalBelow)) { inside = 1; } if (passThrough && (inside ^ inverse)) { pointInArray->SetValue(ptId, flag); cellInArray->SetValue(cellId, flag); } keepCell |= inside; } } else //use cell scalars { keepCell = this->EvaluateValue(inScalars, comp_no, cellId, lims); if (passThrough && (keepCell ^ inverse)) { cellInArray->SetValue(cellId, flag); } } if ( !passThrough && (numCellPts > 0) && (keepCell + inverse == 1) ) // Poor man's XOR { // satisfied thresholding (also non-empty cell, i.e. not VTK_EMPTY_CELL) originalCellIds->InsertNextValue(cellId); for (i=0; i < numCellPts; i++) { ptId = cellPts->GetId(i); if ( (newId = pointMap->GetId(ptId)) < 0 ) { input->GetPoint(ptId, x); newId = newPoints->InsertNextPoint(x); pointMap->SetId(ptId,newId); outPD->CopyData(pd,ptId,newId); originalPointIds->InsertNextValue(ptId); } newCellPts->InsertId(i,newId); } newCellId = outputUG->InsertNextCell(cell->GetCellType(),newCellPts); outCD->CopyData(cd,cellId,newCellId); newCellPts->Reset(); } // satisfied thresholding } // for all cells // now clean up / update ourselves if (passThrough) { pointInArray->Delete(); cellInArray->Delete(); } else { outputUG->SetPoints(newPoints); newPoints->Delete(); pointMap->Delete(); newCellPts->Delete(); originalCellIds->Delete(); } output->Squeeze(); return 1; } //---------------------------------------------------------------------------- int vtkExtractSelectedThresholds::ExtractPoints( vtkSelectionNode *sel, vtkDataSet *input, vtkDataSet *output) { //find the values to threshold within vtkDataArray *lims = vtkDataArray::SafeDownCast(sel->GetSelectionList()); if (lims == NULL) { vtkErrorMacro(<<"No values to threshold with"); return 1; } //find out what array we are supposed to threshold in vtkDataArray *inScalars = NULL; bool use_ids = false; if (sel->GetSelectionList()->GetName()) { if (strcmp(sel->GetSelectionList()->GetName(), "vtkGlobalIds") == 0) { inScalars = input->GetPointData()->GetGlobalIds(); } else if (strcmp(sel->GetSelectionList()->GetName(), "vtkIndices") == 0) { use_ids = true; } else { inScalars = input->GetPointData()->GetArray( sel->GetSelectionList()->GetName()); } } else { inScalars = input->GetPointData()->GetScalars(); } if (inScalars == NULL && !use_ids) { vtkErrorMacro("Could not figure out what array to threshold in."); return 1; } int inverse = 0; if (sel->GetProperties()->Has(vtkSelectionNode::INVERSE())) { inverse = sel->GetProperties()->Get(vtkSelectionNode::INVERSE()); } int passThrough = 0; if (this->PreserveTopology) { passThrough = 1; } int comp_no = 0; if (sel->GetProperties()->Has(vtkSelectionNode::COMPONENT_NUMBER())) { comp_no = sel->GetProperties()->Get(vtkSelectionNode::COMPONENT_NUMBER()); } vtkIdType numPts = input->GetNumberOfPoints(); vtkPointData *inputPD = input->GetPointData(); vtkPointData *outPD = output->GetPointData(); vtkDataSet *outputDS = output; vtkSignedCharArray *pointInArray = NULL; vtkUnstructuredGrid * outputUG = NULL; vtkPoints *newPts = vtkPoints::New(); vtkIdTypeArray* originalPointIds = 0; signed char flag = inverse ? 1 : -1; if (passThrough) { outputDS->ShallowCopy(input); pointInArray = vtkSignedCharArray::New(); pointInArray->SetNumberOfComponents(1); pointInArray->SetNumberOfTuples(numPts); for (vtkIdType i=0; i < numPts; i++) { pointInArray->SetValue(i, flag); } pointInArray->SetName("vtkInsidedness"); outPD->AddArray(pointInArray); outPD->SetScalars(pointInArray); } else { outputUG = vtkUnstructuredGrid::SafeDownCast(output); outputUG->Allocate(numPts); newPts->Allocate(numPts); outputUG->SetPoints(newPts); outPD->CopyGlobalIdsOn(); outPD->CopyAllocate(inputPD); originalPointIds = vtkIdTypeArray::New(); originalPointIds->SetNumberOfComponents(1); originalPointIds->SetName("vtkOriginalPointIds"); outPD->AddArray(originalPointIds); originalPointIds->Delete(); } flag = -flag; vtkIdType outPtCnt = 0; for (vtkIdType ptId = 0; ptId < numPts; ptId++) { int keepPoint = this->EvaluateValue( inScalars, comp_no, ptId, lims ); if (keepPoint ^ inverse) { if (passThrough) { pointInArray->SetValue(ptId, flag); } else { double X[4]; input->GetPoint(ptId, X); newPts->InsertNextPoint(X); outPD->CopyData(inputPD, ptId, outPtCnt); originalPointIds->InsertNextValue(ptId); outputUG->InsertNextCell(VTK_VERTEX, 1, &outPtCnt); outPtCnt++; } } } if (passThrough) { pointInArray->Delete(); } newPts->Delete(); output->Squeeze(); return 1; } //---------------------------------------------------------------------------- void vtkExtractSelectedThresholds::PrintSelf(ostream& os, vtkIndent indent) { this->Superclass::PrintSelf(os,indent); } namespace { template <class daT> bool TestItem(vtkIdType numLims, daT* limsPtr, double value) { for (int i = 0; i < numLims; i+=2) { if (value >= limsPtr[i] && value <= limsPtr[i+1]) { return true; } } return false; } template <class daT> bool TestItem(vtkIdType numLims, daT* limsPtr, double value, int &above, int &below, int& inside) { bool keepCell = false; for (vtkIdType i = 0; i < numLims; i+=2) { daT low = limsPtr[i]; daT high = limsPtr[i+1]; if (value >= low && value <= high) { keepCell = true; ++inside; } else if (value < low) { ++below; } else if (value > high) { ++above; } } return keepCell; } }; //---------------------------------------------------------------------------- int vtkExtractSelectedThresholds::EvaluateValue( vtkDataArray *scalars, int comp_no, vtkIdType id, vtkDataArray *lims) { int keepCell = 0; //check the value in the array against all of the thresholds in lims //if it is inside any, return true double value = 0.0; if (comp_no < 0 && scalars) { // use magnitude. int numComps = scalars->GetNumberOfComponents(); const double *tuple = scalars->GetTuple(id); for (int cc=0; cc < numComps; cc++) { value += tuple[cc]*tuple[cc]; } value = sqrt(value); } else { value = scalars? scalars->GetComponent(id, comp_no) : static_cast<double>(id); /// <=== precision loss when using id. } void* rawLimsPtr = lims->GetVoidPointer(0); vtkIdType numLims = lims->GetNumberOfComponents() * lims->GetNumberOfTuples(); switch (lims->GetDataType()) { vtkTemplateMacro( keepCell = TestItem<VTK_TT>(numLims, static_cast<VTK_TT*>(rawLimsPtr), value)); } return keepCell; } //---------------------------------------------------------------------------- int vtkExtractSelectedThresholds::EvaluateValue( vtkDataArray *scalars, int comp_no, vtkIdType id, vtkDataArray *lims, int *AboveCount, int *BelowCount, int *InsideCount) { double value = 0.0; if (comp_no < 0 && scalars) { // use magnitude. int numComps = scalars->GetNumberOfComponents(); const double *tuple = scalars->GetTuple(id); for (int cc=0; cc < numComps; cc++) { value += tuple[cc]*tuple[cc]; } value = sqrt(value); } else { value = scalars? scalars->GetComponent(id, comp_no) : static_cast<double>(id); /// <=== precision loss when using id. } int keepCell = 0; //check the value in the array against all of the thresholds in lims //if it is inside any, return true int above = 0; int below = 0; int inside = 0; void* rawLimsPtr = lims->GetVoidPointer(0); vtkIdType numLims = lims->GetNumberOfComponents() * lims->GetNumberOfTuples(); switch (lims->GetDataType()) { vtkTemplateMacro( keepCell = TestItem<VTK_TT>(numLims, static_cast<VTK_TT*>(rawLimsPtr), value, above, below, inside)); } if (AboveCount) *AboveCount = above; if (BelowCount) *BelowCount = below; if (InsideCount) *InsideCount = inside; return keepCell; }