1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkContourGrid.cxx
5
6 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7 All rights reserved.
8 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9
10 This software is distributed WITHOUT ANY WARRANTY; without even
11 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12 PURPOSE. See the above copyright notice for more information.
13
14 =========================================================================*/
15 #include "vtkContourGrid.h"
16
17 #include "vtkCell.h"
18 #include "vtkCellArray.h"
19 #include "vtkCellData.h"
20 #include "vtkCellIterator.h"
21 #include "vtkContourValues.h"
22 #include "vtkFloatArray.h"
23 #include "vtkGenericCell.h"
24 #include "vtkInformation.h"
25 #include "vtkInformationVector.h"
26 #include "vtkNew.h"
27 #include "vtkObjectFactory.h"
28 #include "vtkPointData.h"
29 #include "vtkPolyData.h"
30 #include "vtkPolyDataNormals.h"
31 #include "vtkSimpleScalarTree.h"
32 #include "vtkSmartPointer.h"
33 #include "vtkStreamingDemandDrivenPipeline.h"
34 #include "vtkUnstructuredGridBase.h"
35 #include "vtkCutter.h"
36 #include "vtkMergePoints.h"
37 #include "vtkPointLocator.h"
38 #include "vtkIncrementalPointLocator.h"
39 #include "vtkContourHelper.h"
40 #include <cmath>
41
42 vtkStandardNewMacro(vtkContourGrid);
43
44 //-----------------------------------------------------------------------------
45 // Construct object with initial range (0,1) and single contour value
46 // of 0.0.
vtkContourGrid()47 vtkContourGrid::vtkContourGrid()
48 {
49 this->ContourValues = vtkContourValues::New();
50
51 this->ComputeNormals = 0;
52 this->ComputeGradients = 0;
53 this->ComputeScalars = 1;
54 this->GenerateTriangles = 1;
55
56 this->Locator = nullptr;
57
58 this->UseScalarTree = 0;
59 this->ScalarTree = nullptr;
60
61 this->OutputPointsPrecision = DEFAULT_PRECISION;
62
63 // by default process active point scalars
64 this->SetInputArrayToProcess(0,0,0,vtkDataObject::FIELD_ASSOCIATION_POINTS,
65 vtkDataSetAttributes::SCALARS);
66
67 this->EdgeTable = nullptr;
68 }
69
70 //-----------------------------------------------------------------------------
~vtkContourGrid()71 vtkContourGrid::~vtkContourGrid()
72 {
73 this->ContourValues->Delete();
74 if ( this->Locator )
75 {
76 this->Locator->UnRegister(this);
77 this->Locator = nullptr;
78 }
79 if ( this->ScalarTree )
80 {
81 this->ScalarTree->Delete();
82 }
83 }
84
85 //-----------------------------------------------------------------------------
86 // Overload standard modified time function. If contour values are modified,
87 // then this object is modified as well.
GetMTime()88 vtkMTimeType vtkContourGrid::GetMTime()
89 {
90 vtkMTimeType mTime=this->Superclass::GetMTime();
91 vtkMTimeType time;
92
93 if (this->ContourValues)
94 {
95 time = this->ContourValues->GetMTime();
96 mTime = ( time > mTime ? time : mTime );
97 }
98 if (this->Locator)
99 {
100 time = this->Locator->GetMTime();
101 mTime = ( time > mTime ? time : mTime );
102 }
103
104 return mTime;
105 }
106
107 //-----------------------------------------------------------------------------
108 template <class Scalar>
vtkContourGridExecute(vtkContourGrid * self,vtkDataSet * input,vtkPolyData * output,vtkDataArray * inScalars,int numContours,double * values,int computeScalars,int useScalarTree,vtkScalarTree * scalarTree,bool generateTriangles)109 void vtkContourGridExecute(vtkContourGrid *self, vtkDataSet *input,
110 vtkPolyData *output,
111 vtkDataArray *inScalars,
112 int numContours, double *values,
113 int computeScalars,
114 int useScalarTree, vtkScalarTree *scalarTree,
115 bool generateTriangles)
116 {
117 vtkIdType i;
118 int abortExecute=0;
119 vtkIncrementalPointLocator *locator = self->GetLocator();
120 vtkNew<vtkGenericCell> cell;
121 Scalar range[2];
122 vtkCellArray *newVerts, *newLines, *newPolys;
123 vtkPoints *newPts;
124 vtkIdType numCells, estimatedSize;
125 vtkDataArray *cellScalars;
126 Scalar *cellScalarPtr;
127 vtkIdType numCellScalars;
128
129 vtkPointData *inPdOriginal = input->GetPointData();
130
131 // We don't want to change the active scalars in the input, but we
132 // need to set the active scalars to match the input array to
133 // process so that the point data copying works as expected. Create
134 // a shallow copy of point data so that we can do this without
135 // changing the input.
136 vtkSmartPointer<vtkPointData> inPd = vtkSmartPointer<vtkPointData>::New();
137 inPd->ShallowCopy(inPdOriginal);
138
139 // Keep track of the old active scalars because when we set the new
140 // scalars, the old scalars are removed from the point data entirely
141 // and we have to add them back.
142 vtkAbstractArray* oldScalars = inPd->GetScalars();
143 inPd->SetScalars(inScalars);
144 if (oldScalars)
145 {
146 inPd->AddArray(oldScalars);
147 }
148 vtkPointData *outPd = output->GetPointData();
149
150 vtkCellData *inCd = input->GetCellData();
151 vtkCellData *outCd = output->GetCellData();
152
153 //In this case, we know that the input is an unstructured grid.
154 vtkUnstructuredGridBase *grid = static_cast<vtkUnstructuredGridBase *>(input);
155 int needCell = 0;
156 vtkSmartPointer<vtkCellIterator> cellIter =
157 vtkSmartPointer<vtkCellIterator>::Take(input->NewCellIterator());
158
159 numCells = input->GetNumberOfCells();
160
161 //
162 // Create objects to hold output of contour operation. First estimate
163 // allocation size.
164 //
165 estimatedSize=static_cast<vtkIdType>(pow(static_cast<double>(numCells),.75));
166 estimatedSize *= numContours;
167 estimatedSize = estimatedSize / 1024 * 1024; //multiple of 1024
168 if (estimatedSize < 1024)
169 {
170 estimatedSize = 1024;
171 }
172
173 newPts = vtkPoints::New();
174
175 // set precision for the points in the output
176 if(self->GetOutputPointsPrecision() == vtkAlgorithm::DEFAULT_PRECISION)
177 {
178 newPts->SetDataType(grid->GetPoints()->GetDataType());
179 }
180 else if(self->GetOutputPointsPrecision() == vtkAlgorithm::SINGLE_PRECISION)
181 {
182 newPts->SetDataType(VTK_FLOAT);
183 }
184 else if(self->GetOutputPointsPrecision() == vtkAlgorithm::DOUBLE_PRECISION)
185 {
186 newPts->SetDataType(VTK_DOUBLE);
187 }
188
189 newPts->Allocate(estimatedSize,estimatedSize);
190 newVerts = vtkCellArray::New();
191 newVerts->Allocate(estimatedSize,estimatedSize);
192 newLines = vtkCellArray::New();
193 newLines->Allocate(estimatedSize,estimatedSize);
194 newPolys = vtkCellArray::New();
195 newPolys->Allocate(estimatedSize,estimatedSize);
196 cellScalars = inScalars->NewInstance();
197 cellScalars->SetNumberOfComponents(inScalars->GetNumberOfComponents());
198 cellScalars->Allocate(VTK_CELL_SIZE*inScalars->GetNumberOfComponents());
199
200 // locator used to merge potentially duplicate points
201 locator->InitPointInsertion (newPts, input->GetBounds(),
202 input->GetNumberOfPoints());
203
204 // interpolate data along edge
205 // if we did not ask for scalars to be computed, don't copy them
206 if (!computeScalars)
207 {
208 outPd->CopyScalarsOff();
209 }
210 outPd->InterpolateAllocate(inPd,estimatedSize,estimatedSize);
211 outCd->CopyAllocate(inCd,estimatedSize,estimatedSize);
212
213 vtkContourHelper helper(locator, newVerts, newLines, newPolys, inPd, inCd,
214 outPd, outCd, estimatedSize, generateTriangles);
215 // If enabled, build a scalar tree to accelerate search
216 //
217 vtkIdType numCellsContoured = 0;
218 if ( !useScalarTree )
219 {
220 // Three passes over the cells to process lower dimensional cells first.
221 // For poly data output cells need to be added in the order:
222 // verts, lines and then polys, or cell data gets mixed up.
223 // A better solution is to have an unstructured grid output.
224 // I create a table that maps cell type to cell dimensionality,
225 // because I need a fast way to get cell dimensionality.
226 // This assumes GetCell is slow and GetCellType is fast.
227 // I do not like hard coding a list of cell types here,
228 // but I do not want to add GetCellDimension(vtkIdType cellId)
229 // to the vtkDataSet API. Since I anticipate that the output
230 // will change to vtkUnstructuredGrid. This temporary solution
231 // is acceptable.
232 //
233 int cellType;
234 unsigned char cellTypeDimensions[VTK_NUMBER_OF_CELL_TYPES];
235 vtkCutter::GetCellTypeDimensions(cellTypeDimensions);
236 int dimensionality;
237 // We skip 0d cells (points), because they cannot be cut (generate no data).
238 for (dimensionality = 1; dimensionality <= 3; ++dimensionality)
239 {
240 // Loop over all cells; get scalar values for all cell points
241 // and process each cell.
242 //
243 for (cellIter->InitTraversal(); !cellIter->IsDoneWithTraversal();
244 cellIter->GoToNextCell())
245 {
246 if (abortExecute)
247 {
248 break;
249 }
250
251 cellType = cellIter->GetCellType();
252 if (cellType >= VTK_NUMBER_OF_CELL_TYPES)
253 { // Protect against new cell types added.
254 vtkGenericWarningMacro("Unknown cell type " << cellType);
255 continue;
256 }
257 if (cellTypeDimensions[cellType] != dimensionality)
258 {
259 continue;
260 }
261
262 cellScalars->SetNumberOfTuples(cellIter->GetNumberOfPoints());
263 inScalars->GetTuples(cellIter->GetPointIds(), cellScalars);
264 numCellScalars = cellScalars->GetNumberOfComponents()
265 * cellScalars->GetNumberOfTuples();
266 cellScalarPtr = static_cast<Scalar*>(cellScalars->GetVoidPointer(0));
267
268 //find min and max values in scalar data
269 range[0] = range[1] = cellScalarPtr[0];
270
271 for (Scalar *it = cellScalarPtr + 1,
272 *itEnd = cellScalarPtr + numCellScalars; it != itEnd; ++it)
273 {
274 if (*it <= range[0])
275 {
276 range[0] = *it;
277 } //if scalar <= min range value
278 if (*it >= range[1])
279 {
280 range[1] = *it;
281 } //if scalar >= max range value
282 } // for all cellScalars
283
284 if (dimensionality == 3 && ! (cellIter->GetCellId() % 5000) )
285 {
286 self->UpdateProgress(static_cast<double>(cellIter->GetCellId())
287 / numCells);
288 if (self->GetAbortExecute())
289 {
290 abortExecute = 1;
291 break;
292 }
293 }
294
295 for (i = 0; i < numContours; i++)
296 {
297 if ((values[i] >= range[0]) && (values[i] <= range[1]))
298 {
299 needCell = 1;
300 } // if contour value in range for this cell
301 } // end for numContours
302
303 if (needCell)
304 {
305 cellIter->GetCell(cell);
306
307 for (i=0; i < numContours; i++)
308 {
309 if ((values[i] >= range[0]) && (values[i] <= range[1]))
310 {
311 helper.Contour(cell, values[i], cellScalars,
312 cellIter->GetCellId());
313 } // if contour value in range of values for this cell
314 } // for all contour values
315 } // if contour goes through this cell
316 needCell = 0;
317 } // for all cells
318 } // For all dimensions.
319 } //if using scalar tree
320 else
321 {
322 // Note: This will have problems when input contains 2D and 3D cells.
323 // CellData will get scrambled because of the implicit ordering of
324 // verts, lines and polys in vtkPolyData. The solution
325 // is to convert this filter to create unstructured grid.
326 //
327
328 //
329 // Loop over all contour values. Then for each contour value,
330 // loop over all cells.
331 //
332 vtkCell *tmpCell;
333 vtkIdList *dummyIdList = nullptr;
334 vtkIdType cellId = cellIter->GetCellId();
335 for (i=0; i < numContours; i++)
336 {
337 for (scalarTree->InitTraversal(values[i]);
338 (tmpCell = scalarTree->GetNextCell(cellId, dummyIdList,
339 cellScalars)); )
340 {
341 helper.Contour(tmpCell, values[i], cellScalars, cellId);
342 numCellsContoured++;
343
344 //don't want to call Contour any more than necessary
345 } //for all cells
346 } //for all contour values
347 } //using scalar tree
348
349 //
350 // Update ourselves. Because we don't know up front how many verts, lines,
351 // polys we've created, take care to reclaim memory.
352 //
353 output->SetPoints(newPts);
354 newPts->Delete();
355 cellScalars->Delete();
356
357 if (newVerts->GetNumberOfCells())
358 {
359 output->SetVerts(newVerts);
360 }
361 newVerts->Delete();
362
363 if (newLines->GetNumberOfCells())
364 {
365 output->SetLines(newLines);
366 }
367 newLines->Delete();
368
369 if (newPolys->GetNumberOfCells())
370 {
371 output->SetPolys(newPolys);
372 }
373 newPolys->Delete();
374
375 locator->Initialize();//releases leftover memory
376 output->Squeeze();
377 }
378
379 //-----------------------------------------------------------------------------
380 // Contouring filter for unstructured grids.
381 //
RequestData(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)382 int vtkContourGrid::RequestData(
383 vtkInformation *vtkNotUsed(request),
384 vtkInformationVector **inputVector,
385 vtkInformationVector *outputVector)
386 {
387 // get the info objects
388 vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
389 vtkInformation *outInfo = outputVector->GetInformationObject(0);
390
391 // get the input and output
392 vtkUnstructuredGridBase *input = vtkUnstructuredGridBase::SafeDownCast(
393 inInfo->Get(vtkDataObject::DATA_OBJECT()));
394 vtkPolyData *output = vtkPolyData::SafeDownCast(
395 outInfo->Get(vtkDataObject::DATA_OBJECT()));
396
397 vtkDataArray *inScalars;
398 vtkIdType numCells;
399 int numContours = this->ContourValues->GetNumberOfContours();
400 double *values = this->ContourValues->GetValues();
401 int computeScalars = this->ComputeScalars;
402
403 vtkDebugMacro(<< "Executing contour filter");
404
405 if ( this->Locator == nullptr )
406 {
407 this->CreateDefaultLocator();
408 }
409
410 numCells = input->GetNumberOfCells();
411 inScalars = this->GetInputArrayToProcess(0,inputVector);
412 if ( ! inScalars || numCells < 1 )
413 {
414 vtkDebugMacro(<<"No data to contour");
415 return 1;
416 }
417
418 // Create scalar tree if necessary and if requested
419 int useScalarTree = this->GetUseScalarTree();
420 vtkScalarTree *scalarTree = this->ScalarTree;
421 if ( useScalarTree )
422 {
423 if ( scalarTree == nullptr )
424 {
425 this->ScalarTree = scalarTree = vtkSimpleScalarTree::New();
426 }
427 scalarTree->SetDataSet(input);
428 scalarTree->SetScalars(inScalars);
429 }
430
431 switch (inScalars->GetDataType())
432 {
433 vtkTemplateMacro(vtkContourGridExecute<VTK_TT>(
434 this, input, output, inScalars, numContours, values,
435 computeScalars, useScalarTree, scalarTree,
436 this->GenerateTriangles != 0));
437 default:
438 vtkErrorMacro(<< "Execute: Unknown ScalarType");
439 return 1;
440 }
441
442 if(this->ComputeNormals)
443 {
444 vtkInformation* info = outputVector->GetInformationObject(0);
445 vtkNew<vtkPolyDataNormals> normalsFilter;
446 normalsFilter->SetOutputPointsPrecision(this->OutputPointsPrecision);
447 vtkNew<vtkPolyData> tempInput;
448 tempInput->ShallowCopy(output);
449 normalsFilter->SetInputData(tempInput);
450 normalsFilter->SetFeatureAngle(180.);
451 normalsFilter->UpdatePiece(
452 info->Get(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER()),
453 info->Get(vtkStreamingDemandDrivenPipeline:: UPDATE_NUMBER_OF_PIECES()),
454 info->Get(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS()));
455 output->ShallowCopy(normalsFilter->GetOutput());
456 }
457
458 return 1;
459 }
460
461 //-----------------------------------------------------------------------------
462 // Specify a spatial locator for merging points. By default,
463 // an instance of vtkMergePoints is used.
SetScalarTree(vtkScalarTree * sTree)464 void vtkContourGrid::SetScalarTree(vtkScalarTree *sTree)
465 {
466 if ( this->ScalarTree == sTree )
467 {
468 return;
469 }
470 if ( this->ScalarTree )
471 {
472 this->ScalarTree->UnRegister(this);
473 this->ScalarTree = nullptr;
474 }
475 if ( sTree )
476 {
477 sTree->Register(this);
478 }
479 this->ScalarTree = sTree;
480 this->Modified();
481 }
482
483
484 //-----------------------------------------------------------------------------
485 // Specify a spatial locator for merging points. By default,
486 // an instance of vtkMergePoints is used.
SetLocator(vtkIncrementalPointLocator * locator)487 void vtkContourGrid::SetLocator(vtkIncrementalPointLocator *locator)
488 {
489 if ( this->Locator == locator )
490 {
491 return;
492 }
493 if ( this->Locator )
494 {
495 this->Locator->UnRegister(this);
496 this->Locator = nullptr;
497 }
498 if ( locator )
499 {
500 locator->Register(this);
501 }
502 this->Locator = locator;
503 this->Modified();
504 }
505
506 //-----------------------------------------------------------------------------
CreateDefaultLocator()507 void vtkContourGrid::CreateDefaultLocator()
508 {
509 if ( this->Locator == nullptr )
510 {
511 this->Locator = vtkMergePoints::New();
512 this->Locator->Register(this);
513 this->Locator->Delete();
514 }
515 }
516
517 //-----------------------------------------------------------------------------
SetOutputPointsPrecision(int precision)518 void vtkContourGrid::SetOutputPointsPrecision(int precision)
519 {
520 this->OutputPointsPrecision = precision;
521 this->Modified();
522 }
523
524 //-----------------------------------------------------------------------------
GetOutputPointsPrecision() const525 int vtkContourGrid::GetOutputPointsPrecision() const
526 {
527 return this->OutputPointsPrecision;
528 }
529
530 //-----------------------------------------------------------------------------
FillInputPortInformation(int,vtkInformation * info)531 int vtkContourGrid::FillInputPortInformation(int, vtkInformation *info)
532 {
533 info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(),
534 "vtkUnstructuredGridBase");
535 return 1;
536 }
537
538 //-----------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)539 void vtkContourGrid::PrintSelf(ostream& os, vtkIndent indent)
540 {
541 this->Superclass::PrintSelf(os,indent);
542
543 os << indent << "Compute Gradients: "
544 << (this->ComputeGradients ? "On\n" : "Off\n");
545 os << indent << "Compute Normals: "
546 << (this->ComputeNormals ? "On\n" : "Off\n");
547 os << indent << "Compute Scalars: "
548 << (this->ComputeScalars ? "On\n" : "Off\n");
549 os << indent << "Use Scalar Tree: "
550 << (this->UseScalarTree ? "On\n" : "Off\n");
551
552 this->ContourValues->PrintSelf(os,indent.GetNextIndent());
553
554 if ( this->ScalarTree )
555 {
556 os << indent << "Scalar Tree: " << this->ScalarTree << "\n";
557 }
558 else
559 {
560 os << indent << "Scalar Tree: (none)\n";
561 }
562
563 if ( this->Locator )
564 {
565 os << indent << "Locator: " << this->Locator << "\n";
566 }
567 else
568 {
569 os << indent << "Locator: (none)\n";
570 }
571
572 os << indent << "Precision of the output points: "
573 << this->OutputPointsPrecision << "\n";
574 }
575