vtkThreshold.cxx 15.8 KB
Newer Older
Will Schroeder's avatar
Will Schroeder committed
1 2
/*=========================================================================

Ken Martin's avatar
Ken Martin committed
3
  Program:   Visualization Toolkit
Ken Martin's avatar
Ken Martin committed
4
  Module:    vtkThreshold.cxx
Will Schroeder's avatar
Will Schroeder committed
5

6
  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7 8
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
Ken Martin's avatar
Ken Martin committed
9

10 11
     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12
     PURPOSE.  See the above copyright notice for more information.
Will Schroeder's avatar
Will Schroeder committed
13 14

=========================================================================*/
Ken Martin's avatar
Ken Martin committed
15
#include "vtkThreshold.h"
16

17 18 19
#include "vtkCell.h"
#include "vtkCellData.h"
#include "vtkIdList.h"
20 21
#include "vtkInformation.h"
#include "vtkInformationVector.h"
22
#include "vtkObjectFactory.h"
23
#include "vtkPointData.h"
24
#include "vtkUnstructuredGrid.h"
25
#include "vtkStreamingDemandDrivenPipeline.h"
26
#include "vtkMath.h"
27

28 29
#include <algorithm>

Brad King's avatar
Brad King committed
30
vtkStandardNewMacro(vtkThreshold);
31

32
// Construct with lower threshold=0, upper threshold=1, and threshold
Ken Martin's avatar
Ken Martin committed
33
// function=upper AllScalars=1.
Ken Martin's avatar
Ken Martin committed
34
vtkThreshold::vtkThreshold()
Will Schroeder's avatar
Will Schroeder committed
35
{
36 37 38
  this->LowerThreshold         = 0.0;
  this->UpperThreshold         = 1.0;
  this->AllScalars             = 1;
39
  this->AttributeMode          = -1;
40 41 42
  this->ThresholdFunction      = &vtkThreshold::Upper;
  this->ComponentMode          = VTK_COMPONENT_MODE_USE_SELECTED;
  this->SelectedComponent      = 0;
43
  this->OutputPointsPrecision = DEFAULT_PRECISION;
44 45 46 47

  // by default process active point scalars
  this->SetInputArrayToProcess(0,0,0,vtkDataObject::FIELD_ASSOCIATION_POINTS_THEN_CELLS,
                               vtkDataSetAttributes::SCALARS);
48 49 50

  this->GetInformation()->Set(vtkAlgorithm::PRESERVES_RANGES(), 1);
  this->GetInformation()->Set(vtkAlgorithm::PRESERVES_BOUNDS(), 1);
51
  this->UseContinuousCellRange = 0;
Charles Law's avatar
Charles Law committed
52 53 54 55
}

vtkThreshold::~vtkThreshold()
{
Will Schroeder's avatar
Will Schroeder committed
56 57
}

Ken Martin's avatar
Ken Martin committed
58
// Criterion is cells whose scalars are less or equal to lower threshold.
59
void vtkThreshold::ThresholdByLower(double lower)
Will Schroeder's avatar
Will Schroeder committed
60
{
61
  if ( this->LowerThreshold != lower ||
62
       this->ThresholdFunction != &vtkThreshold::Lower)
Will Schroeder's avatar
Will Schroeder committed
63
    {
64
    this->LowerThreshold = lower;
Ken Martin's avatar
Ken Martin committed
65
    this->ThresholdFunction = &vtkThreshold::Lower;
Will Schroeder's avatar
Will Schroeder committed
66 67 68
    this->Modified();
    }
}
69

Ken Martin's avatar
Ken Martin committed
70
// Criterion is cells whose scalars are greater or equal to upper threshold.
Ken Martin's avatar
Ken Martin committed
71
void vtkThreshold::ThresholdByUpper(double upper)
Will Schroeder's avatar
Will Schroeder committed
72
{
73 74
  if ( this->UpperThreshold != upper ||
       this->ThresholdFunction != &vtkThreshold::Upper)
Will Schroeder's avatar
Will Schroeder committed
75
    {
76
    this->UpperThreshold = upper;
Ken Martin's avatar
Ken Martin committed
77
    this->ThresholdFunction = &vtkThreshold::Upper;
Will Schroeder's avatar
Will Schroeder committed
78 79 80
    this->Modified();
    }
}
81

Will Schroeder's avatar
Will Schroeder committed
82
// Criterion is cells whose scalars are between lower and upper thresholds.
Ken Martin's avatar
Ken Martin committed
83
void vtkThreshold::ThresholdBetween(double lower, double upper)
Will Schroeder's avatar
Will Schroeder committed
84
{
85 86
  if ( this->LowerThreshold != lower || this->UpperThreshold != upper ||
       this->ThresholdFunction != &vtkThreshold::Between)
Will Schroeder's avatar
Will Schroeder committed
87
    {
88
    this->LowerThreshold = lower;
Will Schroeder's avatar
Will Schroeder committed
89
    this->UpperThreshold = upper;
Ken Martin's avatar
Ken Martin committed
90
    this->ThresholdFunction = &vtkThreshold::Between;
Will Schroeder's avatar
Will Schroeder committed
91 92 93
    this->Modified();
    }
}
94

95 96 97 98
int vtkThreshold::RequestData(
  vtkInformation *vtkNotUsed(request),
  vtkInformationVector **inputVector,
  vtkInformationVector *outputVector)
Will Schroeder's avatar
Will Schroeder committed
99
{
100 101 102 103
  // get the info objects
  vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation *outInfo = outputVector->GetInformationObject(0);

104
  // get the input and output
105 106 107 108 109
  vtkDataSet *input = vtkDataSet::SafeDownCast(
    inInfo->Get(vtkDataObject::DATA_OBJECT()));
  vtkUnstructuredGrid *output = vtkUnstructuredGrid::SafeDownCast(
    outInfo->Get(vtkDataObject::DATA_OBJECT()));

Amy Squillacote's avatar
Amy Squillacote committed
110
  vtkIdType cellId, newCellId;
Ken Martin's avatar
Ken Martin committed
111
  vtkIdList *cellPts, *pointMap;
112
  vtkIdList *newCellPts;
Ken Martin's avatar
Ken Martin committed
113
  vtkCell *cell;
114
  vtkPoints *newPoints;
Amy Squillacote's avatar
Amy Squillacote committed
115 116
  int i, ptId, newId, numPts;
  int numCellPts;
Ken Martin's avatar
Ken Martin committed
117
  double x[3];
118 119 120
  vtkPointData *pd=input->GetPointData(), *outPD=output->GetPointData();
  vtkCellData *cd=input->GetCellData(), *outCD=output->GetCellData();
  int keepCell, usePointScalars;
Charles Law's avatar
Charles Law committed
121

Ken Martin's avatar
Ken Martin committed
122
  vtkDebugMacro(<< "Executing threshold filter");
123

124 125 126 127 128
  if (this->AttributeMode != -1)
    {
    vtkErrorMacro(<<"You have set the attribute mode on vtkThreshold. This method is deprecated, please use SetInputArrayToProcess instead.");
    return 1;
    }
129

130
  vtkDataArray *inScalars = this->GetInputArrayToProcess(0,inputVector);
131

132
  if (!inScalars)
133
    {
134
    vtkDebugMacro(<<"No scalar data to threshold");
135
    return 1;
136
    }
137

138
  outPD->CopyGlobalIdsOn();
139
  outPD->CopyAllocate(pd);
140
  outCD->CopyGlobalIdsOn();
141 142
  outCD->CopyAllocate(cd);

143 144
  numPts = input->GetNumberOfPoints();
  output->Allocate(input->GetNumberOfCells());
145

146
  newPoints = vtkPoints::New();
147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165

  // 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 if(this->OutputPointsPrecision == vtkAlgorithm::SINGLE_PRECISION)
    {
    newPoints->SetDataType(VTK_FLOAT);
    }
  else if(this->OutputPointsPrecision == vtkAlgorithm::DOUBLE_PRECISION)
    {
    newPoints->SetDataType(VTK_DOUBLE);
    }

Will Schroeder's avatar
Will Schroeder committed
166
  newPoints->Allocate(numPts);
Will Schroeder's avatar
Will Schroeder committed
167

Will Schroeder's avatar
Will Schroeder committed
168
  pointMap = vtkIdList::New(); //maps old point ids into new
169
  pointMap->SetNumberOfIds(numPts);
Bill Lorensen's avatar
Bill Lorensen committed
170 171 172 173
  for (i=0; i < numPts; i++)
    {
    pointMap->SetId(i,-1);
    }
Charles Law's avatar
Charles Law committed
174

175
  newCellPts = vtkIdList::New();
176

177 178
  // are we using pointScalars?
  usePointScalars = (inScalars->GetNumberOfTuples() == numPts);
179

Will Schroeder's avatar
Will Schroeder committed
180
  // Check that the scalars of each cell satisfy the threshold criterion
181
  for (cellId=0; cellId < input->GetNumberOfCells(); cellId++)
Will Schroeder's avatar
Will Schroeder committed
182
    {
183
    cell = input->GetCell(cellId);
Will Schroeder's avatar
Will Schroeder committed
184 185
    cellPts = cell->GetPointIds();
    numCellPts = cell->GetNumberOfPoints();
186

187
    if ( usePointScalars )
Will Schroeder's avatar
Will Schroeder committed
188
      {
189
      if (this->AllScalars)
190 191 192 193 194
        {
        keepCell = 1;
        for ( i=0; keepCell && (i < numCellPts); i++)
          {
          ptId = cellPts->GetId(i);
195
          keepCell = this->EvaluateComponents( inScalars, ptId );
196 197
          }
        }
198
      else
199
        {
200
        if(!this->UseContinuousCellRange)
201
          {
202 203 204 205 206 207 208 209 210 211
          keepCell = 0;
          for ( i=0; (!keepCell) && (i < numCellPts); i++)
            {
            ptId = cellPts->GetId(i);
            keepCell = this->EvaluateComponents( inScalars, ptId );
            }
          }
        else
          {
          keepCell = this->EvaluateCell(inScalars, cellPts, numCellPts);
212 213
          }
        }
Ken Martin's avatar
Ken Martin committed
214
      }
215 216
    else //use cell scalars
      {
217
      keepCell = this->EvaluateComponents( inScalars, cellId );
218
      }
219

Berk Geveci's avatar
Berk Geveci committed
220
    if (  numCellPts > 0 && keepCell )
Will Schroeder's avatar
Will Schroeder committed
221
      {
Berk Geveci's avatar
Berk Geveci committed
222
      // satisfied thresholding (also non-empty cell, i.e. not VTK_EMPTY_CELL)
223 224 225 226 227
      for (i=0; i < numCellPts; i++)
        {
        ptId = cellPts->GetId(i);
        if ( (newId = pointMap->GetId(ptId)) < 0 )
          {
228
          input->GetPoint(ptId, x);
229 230 231 232 233 234
          newId = newPoints->InsertNextPoint(x);
          pointMap->SetId(ptId,newId);
          outPD->CopyData(pd,ptId,newId);
          }
        newCellPts->InsertId(i,newId);
        }
235 236 237 238 239 240 241 242 243 244
      // special handling for polyhedron cells
      if (vtkUnstructuredGrid::SafeDownCast(input) &&
          input->GetCellType(cellId) == VTK_POLYHEDRON)
        {
        newCellPts->Reset();
        vtkUnstructuredGrid::SafeDownCast(input)->
          GetFaceStream(cellId, newCellPts);
        vtkUnstructuredGrid::ConvertFaceStreamPointIds(
          newCellPts, pointMap->GetPointer(0));
        }
245
      newCellId = output->InsertNextCell(cell->GetCellType(),newCellPts);
246
      outCD->CopyData(cd,cellId,newCellId);
247
      newCellPts->Reset();
Will Schroeder's avatar
Will Schroeder committed
248 249 250
      } // satisfied thresholding
    } // for all cells

251
  vtkDebugMacro(<< "Extracted " << output->GetNumberOfCells()
252
                << " number of cells.");
Will Schroeder's avatar
Will Schroeder committed
253

254 255 256
  // now clean up / update ourselves
  pointMap->Delete();
  newCellPts->Delete();
257

258 259
  output->SetPoints(newPoints);
  newPoints->Delete();
Charles Law's avatar
Charles Law committed
260

261
  output->Squeeze();
262 263

  return 1;
Will Schroeder's avatar
Will Schroeder committed
264 265
}

266 267 268 269 270 271 272 273 274 275
int vtkThreshold::EvaluateCell( vtkDataArray *scalars,vtkIdList* cellPts, int numCellPts )
{
  int c(0);
  int numComp = scalars->GetNumberOfComponents();
  int keepCell(0);
  switch (this->ComponentMode)
    {
    case VTK_COMPONENT_MODE_USE_SELECTED:
      c  =   (this->SelectedComponent < numComp)?(this->SelectedComponent):(0);
      keepCell = EvaluateCell(scalars,c,cellPts,numCellPts);
276
      break;
277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309
    case VTK_COMPONENT_MODE_USE_ANY:
      keepCell = 0;
      for ( c = 0; (!keepCell) && (c < numComp); c++ )
        {
        keepCell =EvaluateCell(scalars,c,cellPts,numCellPts);
        }
      break;
    case VTK_COMPONENT_MODE_USE_ALL:
      keepCell = 1;
      for ( c = 0; keepCell && (c < numComp); c++ )
        {
        keepCell =EvaluateCell(scalars,c,cellPts,numCellPts);
        }
      break;
    }
  return keepCell;
}

int vtkThreshold::EvaluateCell( vtkDataArray *scalars, int c, vtkIdList* cellPts, int numCellPts )
{
  double minScalar=DBL_MAX, maxScalar=DBL_MIN;
  for (int i=0; i < numCellPts; i++)
    {
    int ptId = cellPts->GetId(i);
    double s = scalars->GetComponent(ptId,c);
    minScalar = std::min(s,minScalar);
    maxScalar = std::max(s,maxScalar);
    }

  int keepCell =  !(this->LowerThreshold > maxScalar || this->UpperThreshold < minScalar);
  return keepCell;
}

310 311 312 313 314
int vtkThreshold::EvaluateComponents( vtkDataArray *scalars, vtkIdType id )
{
  int keepCell = 0;
  int numComp = scalars->GetNumberOfComponents();
  int c;
315

316 317 318 319
  switch ( this->ComponentMode )
    {
    case VTK_COMPONENT_MODE_USE_SELECTED:
      c = (this->SelectedComponent < numComp)?(this->SelectedComponent):(0);
320
      keepCell =
321 322 323 324 325 326
        (this->*(this->ThresholdFunction))(scalars->GetComponent(id,c));
      break;
    case VTK_COMPONENT_MODE_USE_ANY:
      keepCell = 0;
      for ( c = 0; (!keepCell) && (c < numComp); c++ )
        {
327
        keepCell =
328 329 330 331 332 333 334
          (this->*(this->ThresholdFunction))(scalars->GetComponent(id,c));
        }
      break;
    case VTK_COMPONENT_MODE_USE_ALL:
      keepCell = 1;
      for ( c = 0; keepCell && (c < numComp); c++ )
        {
335
        keepCell =
336 337 338 339 340 341 342
          (this->*(this->ThresholdFunction))(scalars->GetComponent(id,c));
        }
      break;
    }
  return keepCell;
}

343

344
// Return the method for manipulating scalar data as a string.
345
const char *vtkThreshold::GetAttributeModeAsString(void)
346 347 348 349 350 351 352 353 354
{
  if ( this->AttributeMode == VTK_ATTRIBUTE_MODE_DEFAULT )
    {
    return "Default";
    }
  else if ( this->AttributeMode == VTK_ATTRIBUTE_MODE_USE_POINT_DATA )
    {
    return "UsePointData";
    }
355
  else
356 357 358 359 360
    {
    return "UseCellData";
    }
}

361 362 363 364 365 366 367 368 369 370 371
// Return a string representation of the component mode
const char *vtkThreshold::GetComponentModeAsString(void)
{
  if ( this->ComponentMode == VTK_COMPONENT_MODE_USE_SELECTED )
    {
    return "UseSelected";
    }
  else if ( this->ComponentMode == VTK_COMPONENT_MODE_USE_ANY )
    {
    return "UseAny";
    }
372
  else
373 374 375 376 377
    {
    return "UseAll";
    }
}

378 379 380 381
void vtkThreshold::SetPointsDataType(int type)
{
  if(type == VTK_FLOAT)
    {
382
    this->SetOutputPointsPrecision(SINGLE_PRECISION);
383 384 385
    }
  else if(type == VTK_DOUBLE)
    {
386
    this->SetOutputPointsPrecision(DOUBLE_PRECISION);
387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403
    }
}

int vtkThreshold::GetPointsDataType()
{
  if(this->OutputPointsPrecision == SINGLE_PRECISION)
    {
    return VTK_FLOAT;
    }
  else if(this->OutputPointsPrecision == DOUBLE_PRECISION)
    {
    return VTK_DOUBLE;
    }

  return 0;
}

404
void vtkThreshold::SetOutputPointsPrecision(int precision)
405 406 407 408 409 410 411 412 413 414
{
  this->OutputPointsPrecision = precision;
  this->Modified();
}

int vtkThreshold::GetOutputPointsPrecision() const
{
  return this->OutputPointsPrecision;
}

415 416 417 418 419 420
int vtkThreshold::FillInputPortInformation(int, vtkInformation *info)
{
  info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataSet");
  return 1;
}

421
void vtkThreshold::PrintSelf(ostream& os, vtkIndent indent)
Will Schroeder's avatar
Will Schroeder committed
422
{
Brad King's avatar
Brad King committed
423
  this->Superclass::PrintSelf(os,indent);
Will Schroeder's avatar
Will Schroeder committed
424

425
  os << indent << "Attribute Mode: " << this->GetAttributeModeAsString() << endl;
426 427
  os << indent << "Component Mode: " << this->GetComponentModeAsString() << endl;
  os << indent << "Selected Component: " << this->SelectedComponent << endl;
428

Charles Law's avatar
Charles Law committed
429
  os << indent << "All Scalars: " << this->AllScalars << "\n";
Ken Martin's avatar
Ken Martin committed
430
  if ( this->ThresholdFunction == &vtkThreshold::Upper )
Bill Lorensen's avatar
Bill Lorensen committed
431
    {
Will Schroeder's avatar
Will Schroeder committed
432
    os << indent << "Threshold By Upper\n";
Bill Lorensen's avatar
Bill Lorensen committed
433
    }
Will Schroeder's avatar
Will Schroeder committed
434

Ken Martin's avatar
Ken Martin committed
435
  else if ( this->ThresholdFunction == &vtkThreshold::Lower )
Bill Lorensen's avatar
Bill Lorensen committed
436
    {
Will Schroeder's avatar
Will Schroeder committed
437
    os << indent << "Threshold By Lower\n";
Bill Lorensen's avatar
Bill Lorensen committed
438
    }
Will Schroeder's avatar
Will Schroeder committed
439

Ken Martin's avatar
Ken Martin committed
440
  else if ( this->ThresholdFunction == &vtkThreshold::Between )
Bill Lorensen's avatar
Bill Lorensen committed
441
    {
Will Schroeder's avatar
Will Schroeder committed
442
    os << indent << "Threshold Between\n";
Bill Lorensen's avatar
Bill Lorensen committed
443
    }
Will Schroeder's avatar
Will Schroeder committed
444

Charles Law's avatar
Charles Law committed
445 446
  os << indent << "Lower Threshold: " << this->LowerThreshold << "\n";
  os << indent << "Upper Threshold: " << this->UpperThreshold << "\n";
447 448
  os << indent << "Precision of the output points: "
     << this->OutputPointsPrecision << "\n";
449
  os << indent << "Use Continuous Cell Range: "<<this->UseContinuousCellRange<<endl;
Will Schroeder's avatar
Will Schroeder committed
450
}
451 452 453 454 455 456 457 458 459 460 461 462 463 464

//----------------------------------------------------------------------------
int vtkThreshold::ProcessRequest(vtkInformation* request,
                                     vtkInformationVector** inputVector,
                                     vtkInformationVector* outputVector)
{
  // generate the data
  if(request->Has(vtkStreamingDemandDrivenPipeline::REQUEST_UPDATE_EXTENT_INFORMATION()))
    {
    // compute the priority for this UE
    vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
    vtkInformation *outInfo = outputVector->GetInformationObject(0);

    // get the range of the input if available
465 466
    vtkInformation *fInfo = NULL;
    vtkDataArray *inscalars = this->GetInputArrayToProcess(0, inputVector);
467
    if (inscalars)
468
      {
469 470
      vtkInformationVector *miv = inInfo->Get(vtkDataObject::POINT_DATA_VECTOR());
      for (int index = 0; index < miv->GetNumberOfInformationObjects(); index++)
471
        {
472 473 474 475 476 477 478 479
        vtkInformation *mInfo = miv->GetInformationObject(index);
        const char *minfo_arrayname =
          mInfo->Get(vtkDataObject::FIELD_ARRAY_NAME());
        if (minfo_arrayname && !strcmp(minfo_arrayname, inscalars->GetName()))
          {
          fInfo = mInfo;
          break;
          }
480 481
        }
      }
482 483 484 485 486 487
    else
      {
      fInfo = vtkDataObject::GetActiveFieldInformation
        (inInfo, vtkDataObject::FIELD_ASSOCIATION_POINTS,
         vtkDataSetAttributes::SCALARS);
      }
488 489 490 491
    if (!fInfo)
      {
      return 1;
      }
492

493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510
    double *range = fInfo->Get(vtkDataObject::PIECE_FIELD_RANGE());
    if (range)
      {
      // compute the priority
      // get the incoming priority if any
      double inPriority = 1;
      if (inInfo->Has(vtkStreamingDemandDrivenPipeline::PRIORITY()))
        {
        inPriority = inInfo->Get(vtkStreamingDemandDrivenPipeline::PRIORITY());
        }
      outInfo->Set(vtkStreamingDemandDrivenPipeline::PRIORITY(),inPriority);
      if (!inPriority)
        {
        return 1;
        }

      // do any contours intersect the range?
      if (this->ThresholdFunction == &vtkThreshold::Upper)
511
        {
512 513 514 515 516 517 518 519 520 521
        if ((this->*(this->ThresholdFunction))(range[0]))
          {
          return 1;
          }
        }
      if (this->ThresholdFunction == &vtkThreshold::Between)
        {
        if (
            (this->*(this->ThresholdFunction))(range[0]) ||
            (this->*(this->ThresholdFunction))(range[1]) ||
522 523
            (range[0] < this->LowerThreshold
             &&
524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556
             range[1] > this->UpperThreshold))
          {
          return 1;
          };
        }
      if (this->ThresholdFunction == &vtkThreshold::Lower)
        {
        if ((this->*(this->ThresholdFunction))(range[1]))
          {
          return 1;
          }
        }

      double inRes = 1.0;
      if (inInfo->Has(
                      vtkStreamingDemandDrivenPipeline::UPDATE_RESOLUTION()))
        {
        inRes = inInfo->Get(
                            vtkStreamingDemandDrivenPipeline::UPDATE_RESOLUTION());
        }
      if (inRes == 1.0)
        {
        outInfo->Set(vtkStreamingDemandDrivenPipeline::PRIORITY(),0.0);
        }
      else
        {
        outInfo->Set(vtkStreamingDemandDrivenPipeline::PRIORITY(),inPriority*0.1);
        }
      }
    return 1;
    }
  return this->Superclass::ProcessRequest(request, inputVector, outputVector);
}