1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkContour3DLinearGrid.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 "vtkContour3DLinearGrid.h"
16
17 #include "vtkUnstructuredGrid.h"
18 #include "vtkCellArray.h"
19 #include "vtkCellData.h"
20 #include "vtkContourValues.h"
21 #include "vtkTetra.h"
22 #include "vtkHexahedron.h"
23 #include "vtkWedge.h"
24 #include "vtkPyramid.h"
25 #include "vtkVoxel.h"
26 #include "vtkTriangle.h"
27 #include "vtkPointData.h"
28 #include "vtkPolyData.h"
29 #include "vtkFloatArray.h"
30 #include "vtkStaticPointLocator.h"
31 #include "vtkStaticEdgeLocatorTemplate.h"
32 #include "vtkArrayListTemplate.h" // For processing attribute data
33 #include "vtkStaticCellLinksTemplate.h"
34 #include "vtkGarbageCollector.h"
35 #include "vtkInformation.h"
36 #include "vtkInformationVector.h"
37 #include "vtkObjectFactory.h"
38 #include "vtkStreamingDemandDrivenPipeline.h"
39 #include "vtkSMPTools.h"
40 #include "vtkSMPThreadLocalObject.h"
41
42 #include <set>
43
44 vtkStandardNewMacro(vtkContour3DLinearGrid);
45
46 //-----------------------------------------------------------------------------
47 // Classes to support threaded execution. Note that there are different
48 // strategies implemented here: 1) a fast path that just produces output
49 // triangles and points, and 2) more general approach that supports point
50 // merging, field interpolation, and/or normal generation. There is also some
51 // cell-related machinery supporting faster contouring.
52
53 // Macros immediately below are just used to make code easier to
54 // read. Invokes functor _op _num times depending on serial (_seq==1) or
55 // parallel processing mode. The _REDUCE_ version is used to called functors
56 // with a Reduce() method).
57 #define EXECUTE_SMPFOR(_seq,_num,_op) \
58 if ( !_seq) \
59 {\
60 vtkSMPTools::For(0,_num,_op); \
61 }\
62 else \
63 {\
64 _op(0,_num);\
65 }
66
67 #define EXECUTE_REDUCED_SMPFOR(_seq,_num,_op) \
68 if ( !_seq) \
69 {\
70 vtkSMPTools::For(0,_num,_op); \
71 }\
72 else \
73 {\
74 _op.Initialize();\
75 _op(0,_num);\
76 _op.Reduce();\
77 }\
78
79 namespace {
80
81 //========================= CELL MACHINERY ====================================
82
83 // Implementation note: this filter currently handles 3D linear cells. It
84 // could be extended to handle other 3D cell types if the contouring
85 // operation can be expressed as transformation of scalar values from cell
86 // vertices to a case table of triangles.
87
88 // The maximum number of verts per cell
89 #define MAX_CELL_VERTS 8
90
91 // Base class to represent cells
92 struct BaseCell
93 {
94 unsigned char CellType;
95 unsigned char NumVerts;
96 unsigned char NumEdges;
97 unsigned short *Cases;
98 static unsigned char Mask[MAX_CELL_VERTS];
99
BaseCell__anon9cebd3fe0111::BaseCell100 BaseCell(int cellType) : CellType(cellType), NumVerts(0), NumEdges(0), Cases(nullptr) {}
~BaseCell__anon9cebd3fe0111::BaseCell101 virtual ~BaseCell() {}
102
103 // Set up the case table. This is done by accessing standard VTK cells and
104 // repackaging the case table for efficiency. The format of the case table
105 // is as follows: a linear array, organized into two parts: 1) offets into
106 // the second part, and 2) the cases. The first 2^NumVerts entries are the
107 // offsets which refer to the 2^NumVerts cases in the second part. Each
108 // case is represented by the number of edges, followed by pairs of
109 // vertices (v0,v1) for each edge. Note that groups of three contiguous
110 // edges form a triangle.
111 virtual void BuildCases() = 0;
112 void BuildCases(int numCases, int **edges, int **cases, unsigned short *caseArray);
113 };
114 // Used to generate case mask
115 unsigned char BaseCell::Mask[MAX_CELL_VERTS] = {1,2,4,8,16,32,64,128};
116 // Build repackaged case table and place into cases array.
BuildCases(int numCases,int ** edges,int ** cases,unsigned short * caseArray)117 void BaseCell::BuildCases(int numCases, int **edges, int **cases,
118 unsigned short *caseArray)
119 {
120 int caseOffset = numCases;
121 for (int caseNum=0; caseNum < numCases; ++caseNum)
122 {
123 caseArray[caseNum] = caseOffset;
124 int *triCases = cases[caseNum];
125
126 // Count the number of edges
127 int count;
128 for (count=0; triCases[count] != (-1); ++count) {}
129 caseArray[caseOffset++] = count;
130
131 // Now populate the edges
132 int *edge;
133 for (count=0; triCases[count] != (-1); ++count)
134 {
135 edge = edges[triCases[count]];
136 caseArray[caseOffset++] = edge[0];
137 caseArray[caseOffset++] = edge[1];
138 }
139 }//for all cases
140 }
141
142 // Contour tetrahedral cell------------------------------------------------------
143 // Repackages case table for more efficient processing.
144 struct TetCell : public BaseCell
145 {
146 static unsigned short TetCases[152];
147
TetCell__anon9cebd3fe0111::TetCell148 TetCell() : BaseCell(VTK_TETRA)
149 {
150 this->NumVerts = 4;
151 this->NumEdges = 6;
152 this->BuildCases();
153 this->Cases = this->TetCases;
154 }
~TetCell__anon9cebd3fe0111::TetCell155 ~TetCell() override {}
156 void BuildCases() override;
157 };
158 // Dummy initialization filled in later at initialization. The lengtth of the
159 // array is determined from the equation length=(2*NumCases + 3*2*NumTris).
160 unsigned short TetCell::TetCases[152] = { 0 };
161 // Load and transform vtkTetra case table. The case tables are repackaged for
162 // efficiency (e.g., support the GetCase() method).
BuildCases()163 void TetCell::BuildCases()
164 {
165 int **edges = new int*[this->NumEdges];
166 int numCases = std::pow(2,this->NumVerts);
167 int **cases = new int*[numCases];
168 for ( int i=0; i<this->NumEdges; ++i)
169 {
170 edges[i] = vtkTetra::GetEdgeArray(i);
171 }
172 for ( int i=0; i < numCases; ++i)
173 {
174 cases[i] = vtkTetra::GetTriangleCases(i);
175 }
176
177 BaseCell::BuildCases(numCases, edges, cases, this->TetCases);
178
179 delete [] edges;
180 delete [] cases;
181 }
182
183 // Contour hexahedral cell------------------------------------------------------
184 struct HexCell : public BaseCell
185 {
186 static unsigned short HexCases[5432];
187
HexCell__anon9cebd3fe0111::HexCell188 HexCell() : BaseCell(VTK_HEXAHEDRON)
189 {
190 this->NumVerts = 8;
191 this->NumEdges = 12;
192 this->BuildCases();
193 this->Cases = this->HexCases;
194 }
~HexCell__anon9cebd3fe0111::HexCell195 ~HexCell() override {}
196 void BuildCases() override;
197 };
198 // Dummy initialization filled in later at instantiation
199 unsigned short HexCell::HexCases[5432] = { 0 };
200 // Load and transform marching cubes case table. The case tables are
201 // repackaged for efficiency (e.g., support the GetCase() method).
BuildCases()202 void HexCell::BuildCases()
203 {
204 int **edges = new int*[this->NumEdges];
205 int numCases = std::pow(2,this->NumVerts);
206 int **cases = new int*[numCases];
207 for ( int i=0; i<this->NumEdges; ++i)
208 {
209 edges[i] = vtkHexahedron::GetEdgeArray(i);
210 }
211 for ( int i=0; i < numCases; ++i)
212 {
213 cases[i] = vtkHexahedron::GetTriangleCases(i);
214 }
215
216 BaseCell::BuildCases(numCases, edges, cases, this->HexCases);
217
218 delete [] edges;
219 delete [] cases;
220 }
221
222 // Contour wedge cell ------------------------------------------------------
223 struct WedgeCell : public BaseCell
224 {
225 static unsigned short WedgeCases[968];
226
WedgeCell__anon9cebd3fe0111::WedgeCell227 WedgeCell() : BaseCell(VTK_WEDGE)
228 {
229 this->NumVerts = 6;
230 this->NumEdges = 9;
231 this->BuildCases();
232 this->Cases = this->WedgeCases;
233 }
~WedgeCell__anon9cebd3fe0111::WedgeCell234 ~WedgeCell() override {}
235 void BuildCases() override;
236 };
237 // Dummy initialization filled in later at instantiation
238 unsigned short WedgeCell::WedgeCases[968] = { 0 };
239 // Load and transform marching cubes case table. The case tables are
240 // repackaged for efficiency (e.g., support the GetCase() method).
BuildCases()241 void WedgeCell::BuildCases()
242 {
243 int **edges = new int*[this->NumEdges];
244 int numCases = std::pow(2,this->NumVerts);
245 int **cases = new int*[numCases];
246 for ( int i=0; i<this->NumEdges; ++i)
247 {
248 edges[i] = vtkWedge::GetEdgeArray(i);
249 }
250 for ( int i=0; i < numCases; ++i)
251 {
252 cases[i] = vtkWedge::GetTriangleCases(i);
253 }
254
255 BaseCell::BuildCases(numCases, edges, cases, this->WedgeCases);
256
257 delete [] edges;
258 delete [] cases;
259 }
260
261 // Contour pyramid cell------------------------------------------------------
262 struct PyrCell : public BaseCell
263 {
264 static unsigned short PyrCases[448];
265
PyrCell__anon9cebd3fe0111::PyrCell266 PyrCell() : BaseCell(VTK_PYRAMID)
267 {
268 this->NumVerts = 5;
269 this->NumEdges = 8;
270 this->BuildCases();
271 this->Cases = this->PyrCases;
272 }
~PyrCell__anon9cebd3fe0111::PyrCell273 ~PyrCell() override {}
274 void BuildCases() override;
275 };
276 // Dummy initialization filled in later at instantiation
277 unsigned short PyrCell::PyrCases[448] = { 0 };
278 // Load and transform marching cubes case table. The case tables are
279 // repackaged for efficiency (e.g., support the GetCase() method).
BuildCases()280 void PyrCell::BuildCases()
281 {
282 int **edges = new int*[this->NumEdges];
283 int numCases = std::pow(2,this->NumVerts);
284 int **cases = new int*[numCases];
285 for ( int i=0; i<this->NumEdges; ++i)
286 {
287 edges[i] = vtkPyramid::GetEdgeArray(i);
288 }
289 for ( int i=0; i < numCases; ++i)
290 {
291 cases[i] = vtkPyramid::GetTriangleCases(i);
292 }
293
294 BaseCell::BuildCases(numCases, edges, cases, this->PyrCases);
295
296 delete [] edges;
297 delete [] cases;
298 }
299
300 // Contour voxel cell------------------------------------------------------
301 struct VoxCell : public BaseCell
302 {
303 static unsigned short VoxCases[5432];
304
VoxCell__anon9cebd3fe0111::VoxCell305 VoxCell() : BaseCell(VTK_VOXEL)
306 {
307 this->NumVerts = 8;
308 this->NumEdges = 12;
309 this->BuildCases();
310 this->Cases = this->VoxCases;
311 }
~VoxCell__anon9cebd3fe0111::VoxCell312 ~VoxCell() override {};
313 void BuildCases() override;
314 };
315 // Dummy initialization filled in later at instantiation
316 unsigned short VoxCell::VoxCases[5432] = { 0 };
317 // Load and transform marching cubes case table. The case tables are
318 // repackaged for efficiency (e.g., support the GetCase() method). Note that
319 // the MC cases (vtkMarchingCubesTriangleCases) are specified for the
320 // hexahedron; voxels require a transformation to produce correct output.
BuildCases()321 void VoxCell::BuildCases()
322 {
323 // Map the voxel points consistent with the hex edges and cases, Basically
324 // the hex points (2,3,6,7) are ordered (3,2,7,6) on the voxel.
325 int **edges = new int*[this->NumEdges];
326 int voxEdges[12][2] = { {0,1}, {1,3}, {2,3}, {0,2}, {4,5}, {5,7},
327 {6,7}, {4,6}, {0,4}, {1,5}, {2,6}, {3,7} };
328
329 for ( int i=0; i<this->NumEdges; ++i)
330 {
331 edges[i] = voxEdges[i];
332 }
333
334 // Build the voxel cases. Have to shuffle them around due to different
335 // vertex ordering.
336 unsigned int numCases = std::pow(2,this->NumVerts);
337 int **cases = new int*[numCases];
338 unsigned int hexCase, voxCase;
339 for ( hexCase=0; hexCase < numCases; ++hexCase)
340 {
341 voxCase = ((hexCase & BaseCell::Mask[0]) ? 1 : 0) << 0;
342 voxCase |= ((hexCase & BaseCell::Mask[1]) ? 1 : 0) << 1;
343 voxCase |= ((hexCase & BaseCell::Mask[2]) ? 1 : 0) << 3;
344 voxCase |= ((hexCase & BaseCell::Mask[3]) ? 1 : 0) << 2;
345 voxCase |= ((hexCase & BaseCell::Mask[4]) ? 1 : 0) << 4;
346 voxCase |= ((hexCase & BaseCell::Mask[5]) ? 1 : 0) << 5;
347 voxCase |= ((hexCase & BaseCell::Mask[6]) ? 1 : 0) << 7;
348 voxCase |= ((hexCase & BaseCell::Mask[7]) ? 1 : 0) << 6;
349 cases[voxCase] = vtkHexahedron::GetTriangleCases(hexCase);
350 }
351
352 BaseCell::BuildCases(numCases, edges, cases, this->VoxCases);
353
354 delete [] edges;
355 delete [] cases;
356 }
357
358 // Contour empty cell. These cells are skipped.---------------------------------
359 struct EmptyCell : public BaseCell
360 {
361 static unsigned short EmptyCases[2];
362
EmptyCell__anon9cebd3fe0111::EmptyCell363 EmptyCell() : BaseCell(VTK_EMPTY_CELL)
364 {
365 this->NumVerts = 0;
366 this->NumEdges = 0;
367 this->Cases = this->EmptyCases;
368 }
~EmptyCell__anon9cebd3fe0111::EmptyCell369 ~EmptyCell() override {}
BuildCases__anon9cebd3fe0111::EmptyCell370 void BuildCases() override {}
371 };
372 // No triangles generated
373 unsigned short EmptyCell::EmptyCases[2] = { 0,0 };
374
375
376 // This is a general iterator which assumes that the unstructured grid has a
377 // mix of cells. Any cell that is not processed by this contouring algorithm
378 // (i.e., not one of tet, hex, pyr, wedge, voxel) is skipped.
379 struct CellIter
380 {
381 // Current active cell, and whether it is a copy (which controls
382 // the destruction process).
383 bool Copy;
384 BaseCell *Cell;
385
386 // The iteration state.
387 unsigned char NumVerts;
388 const unsigned short *Cases;
389 vtkIdType Incr;
390
391 // Refernces to unstructured grid for cell traversal.
392 vtkIdType NumCells;
393 const unsigned char *Types;
394 const vtkIdType *Conn;
395 const vtkIdType *Locs;
396
397 // All possible cell types. The iterator switches between them when
398 // processing. All unsupported cells are of type EmptyCell.
399 TetCell *Tet;
400 HexCell *Hex;
401 PyrCell *Pyr;
402 WedgeCell *Wedge;
403 VoxCell *Vox;
404 EmptyCell *Empty;
405
CellIter__anon9cebd3fe0111::CellIter406 CellIter() : Copy(true), Cell(nullptr), NumVerts(0), Cases(nullptr), Incr(0),
407 NumCells(0), Types(nullptr), Conn(nullptr), Locs(nullptr), Tet(nullptr),
408 Hex(nullptr), Pyr(nullptr), Wedge(nullptr), Vox(nullptr), Empty(nullptr)
409 {}
410
CellIter__anon9cebd3fe0111::CellIter411 CellIter(vtkIdType numCells, unsigned char *types, vtkIdType *conn, vtkIdType *locs) :
412 Copy(false), Cell(nullptr), NumVerts(0), Cases(nullptr), Incr(0),
413 NumCells(numCells), Types(types), Conn(conn), Locs(locs)
414 {
415 this->Tet = new TetCell;
416 this->Hex = new HexCell;
417 this->Pyr = new PyrCell;
418 this->Wedge = new WedgeCell;
419 this->Vox = new VoxCell;
420 this->Empty = new EmptyCell;
421 }
422
~CellIter__anon9cebd3fe0111::CellIter423 ~CellIter()
424 {
425 if ( ! this->Copy )
426 {
427 delete this->Tet;
428 delete this->Hex;
429 delete this->Pyr;
430 delete this->Wedge;
431 delete this->Vox;
432 delete this->Empty;
433 }
434 }
435
436 CellIter(const CellIter &) = default; //remove compiler warnings
437
438 // Shallow copy to avoid new/delete.
operator =__anon9cebd3fe0111::CellIter439 CellIter& operator=(const CellIter& cellIter)
440 {
441 this->Copy = true;
442 this->Cell = nullptr;
443
444 this->NumVerts = cellIter.NumVerts;
445 this->Cases = cellIter.Cases;
446 this->Incr = cellIter.Incr;
447
448 this->NumCells = cellIter.NumCells;
449 this->Types = cellIter.Types;
450 this->Conn = cellIter.Conn;
451 this->Locs = cellIter.Locs;
452
453 this->Tet = cellIter.Tet;
454 this->Hex = cellIter.Hex;
455 this->Pyr = cellIter.Pyr;
456 this->Wedge = cellIter.Wedge;
457 this->Vox = cellIter.Vox;
458 this->Empty = cellIter.Empty;
459
460 return *this;
461 }
462
463 // Decode the case table. (See previous documentation of case table
464 // organization.) Note that bounds/range chacking is not performed
465 // for efficiency.
GetCase__anon9cebd3fe0111::CellIter466 const unsigned short *GetCase(unsigned char caseNum)
467 {
468 return (this->Cases + this->Cases[caseNum]);
469 }
470
471 // Methods for caching traversal. Initialize() sets up the traversal
472 // process; Next() advances to the next cell. Note that the public data
473 // members representing the iteration state (NumVerts, Cases, Incr) are
474 // modified by these methods, and then subsequently read during iteration.
Initialize__anon9cebd3fe0111::CellIter475 const vtkIdType* Initialize(vtkIdType cellId)
476 {
477 this->Cell = this->GetCell(this->Types[cellId]);
478 this->NumVerts = this->Cell->NumVerts;
479 this->Cases = this->Cell->Cases;
480
481 if ( this->Cell->CellType != VTK_EMPTY_CELL )
482 {
483 this->Incr = this->NumVerts + 1;
484 }
485 else // Else have to update the increment differently
486 {
487 this->Incr = (cellId >= (this->NumCells-1) ? 0 :
488 this->Locs[cellId+1] - this->Locs[cellId]);
489 }
490
491 return (this->Conn + this->Locs[cellId] + 1);
492 }
493
Next__anon9cebd3fe0111::CellIter494 vtkIdType Next(vtkIdType cellId)
495 {
496 // Guard against end of array condition; only update information if the
497 // cell type changes. Note however that empty cells may have to be
498 // treated specially.
499 if ( cellId >= (this->NumCells-1) ||
500 (this->Cell->CellType != VTK_EMPTY_CELL &&
501 this->Cell->CellType == this->Types[cellId+1]) )
502 {
503 return this->Incr;
504 }
505
506 // Need to look up new information as cell type has changed.
507 vtkIdType incr = this->Incr;
508 this->Cell = this->GetCell(this->Types[cellId+1]);
509 this->NumVerts = this->Cell->NumVerts;
510 this->Cases = this->Cell->Cases;
511
512 if ( this->Cell->CellType != VTK_EMPTY_CELL )
513 {
514 this->Incr = this->NumVerts + 1;
515 }
516 else // Else have to update the increment differently
517 {
518 this->Incr = this->Locs[cellId+2] - this->Locs[cellId+1];
519 }
520
521 return incr;
522 }
523
524 // Switch to the appropriate cell type.
GetCell__anon9cebd3fe0111::CellIter525 BaseCell *GetCell(int cellType)
526 {
527 switch ( cellType )
528 {
529 case VTK_TETRA:
530 return this->Tet;
531 case VTK_HEXAHEDRON:
532 return this->Hex;
533 case VTK_WEDGE:
534 return this->Wedge;
535 case VTK_PYRAMID:
536 return this->Pyr;
537 case VTK_VOXEL:
538 return this->Vox;
539 default:
540 return this->Empty;
541 }
542 }
543 };
544
545
546 //========================= FAST PATH =========================================
547 // Perform the contouring operation without merging coincident points.
548 template <typename TIP, typename TOP, typename TS>
549 struct ContourCells
550 {
551 typedef std::vector<TOP> LocalPtsType;
552
553 // Track local data on a per-thread basis. In the Reduce() method this
554 // information will be used to composite the data from each thread into a
555 // single vtkPolyData output.
556 struct LocalDataType
557 {
558 LocalPtsType LocalPts;
559 CellIter LocalCellIter;
560
LocalDataType__anon9cebd3fe0111::ContourCells::LocalDataType561 LocalDataType()
562 {
563 this->LocalPts.reserve(2048);
564 }
565 };
566
567 CellIter *Iter;
568 const TIP *InPts;
569 const TS *Scalars;
570 double Value;
571 vtkPoints *NewPts;
572 vtkCellArray *NewPolys;
573
574 // Keep track of generated points and triangles on a per thread basis
575 vtkSMPThreadLocal<LocalDataType> LocalData;
576
577 // Results from the compositing Reduce() method
578 vtkIdType NumPts;
579 vtkIdType NumTris;
580 int NumThreadsUsed;
581 vtkIdType TotalPts; //the total points thus far (support multiple contours)
582 vtkIdType TotalTris; //the total triangles thus far (support multiple contours)
583
ContourCells__anon9cebd3fe0111::ContourCells584 ContourCells(TIP *inPts, CellIter *iter, TS *s, double value,
585 vtkPoints* outPts, vtkCellArray *tris, vtkIdType totalPts,
586 vtkIdType totalTris) :
587 Iter(iter), InPts(inPts), Scalars(s), Value(value), NewPts(outPts),
588 NewPolys(tris), NumPts(0), NumTris(0), NumThreadsUsed(0),
589 TotalPts(totalPts), TotalTris(totalTris)
590 {
591 }
592
593 // Set up the iteration process.
Initialize__anon9cebd3fe0111::ContourCells594 void Initialize()
595 {
596 typename ContourCells<TIP,TOP,TS>::LocalDataType& localData = this->LocalData.Local();
597 localData.LocalCellIter = *(this->Iter);
598 }
599
600 // operator() method extracts points from cells (points taken three at a
601 // time form a triangle)
operator ()__anon9cebd3fe0111::ContourCells602 void operator()(vtkIdType cellId, vtkIdType endCellId)
603 {
604 typename ContourCells<TIP,TOP,TS>::LocalDataType& localData = this->LocalData.Local();
605 typename ContourCells<TIP,TOP,TS>::LocalPtsType& lPts = localData.LocalPts;
606 CellIter *cellIter = &localData.LocalCellIter;
607 const vtkIdType *c = cellIter->Initialize(cellId);
608 unsigned short isoCase, numEdges, i;
609 const unsigned short *edges;
610 double s[MAX_CELL_VERTS], value=this->Value, deltaScalar;
611 float t;
612 unsigned char v0, v1;
613 const TIP *x[MAX_CELL_VERTS];
614
615 for ( ; cellId < endCellId; ++cellId )
616 {
617 // Compute case by repeated masking of scalar value
618 for ( isoCase=0, i=0; i < cellIter->NumVerts; ++i )
619 {
620 s[i] = static_cast<double>(*(this->Scalars + c[i]));
621 isoCase |= ( s[i] >= value ? BaseCell::Mask[i] : 0 );
622 }
623 edges = cellIter->GetCase(isoCase);
624
625 if ( *edges > 0 )
626 {
627 numEdges = *edges++;
628 for ( i=0; i < cellIter->NumVerts; ++i )
629 {
630 x[i] = this->InPts + 3*c[i];
631 }
632
633 for (i=0; i<numEdges; ++i, edges+=2)
634 {
635 v0 = edges[0];
636 v1 = edges[1];
637 deltaScalar = s[v1] - s[v0];
638 t = (deltaScalar == 0.0 ? 0.0 : (value - s[v0]) / deltaScalar);
639 lPts.emplace_back(x[v0][0] + t*(x[v1][0] - x[v0][0]));
640 lPts.emplace_back(x[v0][1] + t*(x[v1][1] - x[v0][1]));
641 lPts.emplace_back(x[v0][2] + t*(x[v1][2] - x[v0][2]));
642 }//for all edges in this case
643 }//if contour passes through this cell
644 c += cellIter->Next(cellId); //move to the next cell
645 }//for all cells in this batch
646 }
647
648 // Composite results from each thread
Reduce__anon9cebd3fe0111::ContourCells649 void Reduce()
650 {
651 // Count the number of points. For fun keep track of the number of
652 // threads used.
653 vtkIdType numPts = 0;
654 this->NumThreadsUsed = 0;
655 typename vtkSMPThreadLocal<LocalDataType>::iterator ldItr;
656 typename vtkSMPThreadLocal<LocalDataType>::iterator ldEnd = this->LocalData.end();
657 for ( ldItr=this->LocalData.begin(); ldItr != ldEnd; ++ldItr )
658 {
659 numPts += static_cast<vtkIdType>(((*ldItr).LocalPts.size() / 3)); //x-y-z components
660 this->NumThreadsUsed++;
661 }
662
663 // (Re)Allocate space for output. Multiple contours require writing into
664 // the middle of arrays.
665 this->NumPts = numPts;
666 this->NumTris = numPts / 3;
667 this->NewPts->GetData()->WriteVoidPointer(0,3*(this->NumPts+this->TotalPts));
668 TOP *pts = static_cast<TOP*>(this->NewPts->GetVoidPointer(this->TotalPts*3));
669 this->NewPolys->WritePointer((this->NumTris+this->TotalTris),
670 4*(this->NumTris+this->TotalTris));
671
672 // Copy points output to VTK structures. Only point coordinates are
673 // copied for now; later we'll define the triangle topology.
674 typename LocalPtsType::iterator pItr, pEnd;
675 for ( ldItr=this->LocalData.begin(); ldItr != ldEnd; ++ldItr )
676 {
677 pEnd = (*ldItr).LocalPts.end();
678 for ( pItr = (*ldItr).LocalPts.begin(); pItr != pEnd; )
679 {
680 *pts++ = *pItr++;
681 }
682 }//For all threads
683 }//Reduce
684 };//ContourCells
685
686
687 // Fast path processing. Handles template dispatching and such.
688 template <typename TS>
ProcessFastPath(vtkIdType numTets,vtkPoints * inPts,CellIter * cellIter,TS * s,double isoValue,vtkPoints * outPts,vtkCellArray * tris,vtkTypeBool seqProcessing,int & numThreads,vtkIdType totalPts,vtkIdType totalTris)689 void ProcessFastPath(vtkIdType numTets, vtkPoints *inPts, CellIter *cellIter,
690 TS *s, double isoValue, vtkPoints *outPts,
691 vtkCellArray *tris, vtkTypeBool seqProcessing,
692 int &numThreads, vtkIdType totalPts, vtkIdType totalTris)
693 {
694 double val = static_cast<double>(isoValue);
695 int inPtsType = inPts->GetDataType();
696 void *inPtsPtr = inPts->GetVoidPointer(0);
697 int outPtsType = outPts->GetDataType();
698 if ( inPtsType == VTK_FLOAT && outPtsType == VTK_FLOAT )
699 {
700 ContourCells<float,float,TS> contour((float*)inPtsPtr,cellIter,
701 (TS*)s,val,outPts,tris,totalPts,totalTris);
702 EXECUTE_REDUCED_SMPFOR(seqProcessing,numTets,contour);
703 numThreads = contour.NumThreadsUsed;
704 }
705 else if ( inPtsType == VTK_DOUBLE && outPtsType == VTK_DOUBLE )
706 {
707 ContourCells<double,double,TS> contour((double*)inPtsPtr,cellIter,
708 (TS*)s,val,outPts,tris,totalPts,totalTris);
709 EXECUTE_REDUCED_SMPFOR(seqProcessing,numTets,contour);
710 numThreads = contour.NumThreadsUsed;
711 }
712 else if ( inPtsType == VTK_FLOAT && outPtsType == VTK_DOUBLE )
713 {
714 ContourCells<float,double,TS> contour((float*)inPtsPtr,cellIter,
715 (TS*)s,val,outPts,tris,totalPts,totalTris);
716 EXECUTE_REDUCED_SMPFOR(seqProcessing,numTets,contour);
717 numThreads = contour.NumThreadsUsed;
718 }
719 else //if ( inPtsType == VTK_DOUBLE && outPtsType == VTK_FLOAT )
720 {
721 ContourCells<double,float,TS> contour((double*)inPtsPtr,cellIter,
722 (TS*)s,val,outPts,tris,totalPts,totalTris);
723 EXECUTE_REDUCED_SMPFOR(seqProcessing,numTets,contour);
724 numThreads = contour.NumThreadsUsed;
725 }
726 };
727
728 // Produce triangles for non-merged points-------------------------
729 struct ProduceTriangles
730 {
731 vtkIdType *Tris;
732
ProduceTriangles__anon9cebd3fe0111::ProduceTriangles733 ProduceTriangles(vtkIdType *tris) : Tris(tris)
734 {}
735
operator ()__anon9cebd3fe0111::ProduceTriangles736 void operator()(vtkIdType triId, vtkIdType endTriId)
737 {
738 vtkIdType *tris = this->Tris + 4*triId;
739 vtkIdType id = 3*triId;
740 for ( ; triId < endTriId; ++triId )
741 {
742 *tris++ = 3;
743 *tris++ = id++;
744 *tris++ = id++;
745 *tris++ = id++;
746 }
747 }
748 };
749
750
751 //========================= GENERAL PATH (POINT MERGING) =======================
752 // Use vtkStaticEdgeLocatorTemplate for edge-based point merging
753 template <typename IDType, typename TS>
754 struct ExtractEdges
755 {
756 typedef std::vector<EdgeTuple<IDType,float>> EdgeVectorType;
757 typedef std::vector<MergeTuple<IDType,float>> MergeVectorType;
758
759 // Track local data on a per-thread basis. In the Reduce() method this
760 // information will be used to composite the data from each thread.
761 struct LocalDataType
762 {
763 EdgeVectorType LocalEdges;
764 CellIter LocalCellIter;
765
LocalDataType__anon9cebd3fe0111::ExtractEdges::LocalDataType766 LocalDataType()
767 {
768 this->LocalEdges.reserve(2048);
769 }
770 };
771
772 CellIter *Iter;
773 const TS *Scalars;
774 double Value;
775 MergeTuple<IDType,float> *Edges;
776 vtkCellArray *Tris;
777 vtkIdType NumTris;
778 int NumThreadsUsed;
779 vtkIdType TotalTris; //the total triangles thus far (support multiple contours)
780
781 // Keep track of generated points and triangles on a per thread basis
782 vtkSMPThreadLocal<LocalDataType> LocalData;
783
ExtractEdges__anon9cebd3fe0111::ExtractEdges784 ExtractEdges(CellIter *c, TS *s, double value, vtkCellArray *tris,
785 vtkIdType totalTris) :
786 Iter(c), Scalars(s), Value(value), Edges(nullptr), Tris(tris),
787 NumTris(0), NumThreadsUsed(0), TotalTris(totalTris)
788 {}
789
790 // Set up the iteration process
Initialize__anon9cebd3fe0111::ExtractEdges791 void Initialize()
792 {
793 typename ExtractEdges<IDType,TS>::LocalDataType& localData = this->LocalData.Local();
794 localData.LocalCellIter = *(this->Iter);
795 }
796
797 // operator() method extracts edges from cells (edges taken three at a
798 // time form a triangle)
operator ()__anon9cebd3fe0111::ExtractEdges799 void operator()(vtkIdType cellId, vtkIdType endCellId)
800 {
801 typename ExtractEdges<IDType,TS>::LocalDataType& localData = this->LocalData.Local();
802 typename ExtractEdges<IDType,TS>::EdgeVectorType& lEdges = localData.LocalEdges;
803 CellIter *cellIter = &localData.LocalCellIter;
804 const vtkIdType *c = cellIter->Initialize(cellId); //connectivity array
805 unsigned short isoCase, numEdges, i;
806 const unsigned short *edges;
807 double s[MAX_CELL_VERTS], value=this->Value, deltaScalar;
808 float t;
809 unsigned char v0, v1;
810
811 for ( ; cellId < endCellId; ++cellId )
812 {
813 // Compute case by repeated masking of scalar value
814 for ( isoCase=0, i=0; i < cellIter->NumVerts; ++i )
815 {
816 s[i] = static_cast<double>(*(this->Scalars + c[i]));
817 isoCase |= ( s[i] >= value ? BaseCell::Mask[i] : 0 );
818 }
819 edges = cellIter->GetCase(isoCase);
820
821 if ( *edges > 0 )
822 {
823 numEdges = *edges++;
824 for (i=0; i<numEdges; ++i, edges+=2)
825 {
826 v0 = edges[0];
827 v1 = edges[1];
828 deltaScalar = s[v1] - s[v0];
829 t = (deltaScalar == 0.0 ? 0.0 : (value - s[v0]) / deltaScalar);
830 t = ( c[v0] < c[v1] ? t : (1.0-t) ); //edges (v0,v1) must have v0<v1
831 lEdges.emplace_back(c[v0],c[v1],t); //edge constructor may swap v0<->v1
832 }//for all edges in this case
833 }//if contour passes through this cell
834 c += cellIter->Next(cellId); //move to the next cell
835 }//for all cells in this batch
836 }
837
838 // Composite local thread data
Reduce__anon9cebd3fe0111::ExtractEdges839 void Reduce()
840 {
841 // Count the number of triangles, and number of threads used.
842 vtkIdType numTris = 0;
843 this->NumThreadsUsed = 0;
844 typename vtkSMPThreadLocal<LocalDataType>::iterator ldItr;
845 typename vtkSMPThreadLocal<LocalDataType>::iterator ldEnd = this->LocalData.end();
846 for ( ldItr=this->LocalData.begin(); ldItr != ldEnd; ++ldItr )
847 {
848 numTris += static_cast<vtkIdType>(((*ldItr).LocalEdges.size() / 3)); //three edges per triangle
849 this->NumThreadsUsed++;
850 }
851
852 // Allocate space for VTK triangle output. Take into account previous
853 // contours.
854 this->NumTris = numTris;
855 this->Tris->WritePointer((this->NumTris+this->TotalTris),
856 4*(this->NumTris+this->TotalTris));
857
858 // Copy local edges to global edge array. Add in the originating edge id
859 // used later when merging.
860 EdgeVectorType emptyVector;
861 typename EdgeVectorType::iterator eItr, eEnd;
862 this->Edges = new MergeTuple<IDType,float>[3*this->NumTris]; //three edges per triangle
863 vtkIdType edgeNum=0;
864 for ( ldItr=this->LocalData.begin(); ldItr != ldEnd; ++ldItr )
865 {
866 eEnd = (*ldItr).LocalEdges.end();
867 for ( eItr = (*ldItr).LocalEdges.begin(); eItr != eEnd; ++eItr )
868 {
869 this->Edges[edgeNum].V0 = eItr->V0;
870 this->Edges[edgeNum].V1 = eItr->V1;
871 this->Edges[edgeNum].T = eItr->T;
872 this->Edges[edgeNum].EId = edgeNum;
873 edgeNum++;
874 }
875 (*ldItr).LocalEdges.swap(emptyVector); //frees memory
876 }//For all threads
877 }//Reduce
878 };//ExtractEdges
879
880 // This method generates the output isosurface triangle connectivity list.
881 template <typename IDType>
882 struct ProduceMergedTriangles
883 {
884 typedef MergeTuple<IDType,float> MergeTupleType;
885
886 const MergeTupleType *MergeArray;
887 const IDType *Offsets;
888 vtkIdType NumTris;
889 vtkIdType *Tris;
890 vtkIdType TotalPts;
891
ProduceMergedTriangles__anon9cebd3fe0111::ProduceMergedTriangles892 ProduceMergedTriangles(const MergeTupleType *merge, const IDType *offsets,
893 vtkIdType numTris, vtkIdType *tris,
894 vtkIdType totalPts, vtkIdType totalTris) :
895 MergeArray(merge), Offsets(offsets), NumTris(numTris), Tris(tris),
896 TotalPts(totalPts)
897 {
898 this->Tris = tris + 4*totalTris;
899 }
900
Initialize__anon9cebd3fe0111::ProduceMergedTriangles901 void Initialize()
902 {
903 ;//without this method Reduce() is not called
904 }
905
906 // Loop over all merged points and update the ids of the triangle
907 // connectivity. Offsets point to the beginning of a group of equal edges:
908 // all edges in the group are updated to the current merged point id.
operator ()__anon9cebd3fe0111::ProduceMergedTriangles909 void operator()(vtkIdType ptId, vtkIdType endPtId)
910 {
911 const MergeTupleType *mergeArray = this->MergeArray;
912 const IDType *offsets = this->Offsets;
913 IDType i, numPtsInGroup, eid, triId;
914
915 for ( ; ptId < endPtId; ++ptId )
916 {
917 numPtsInGroup = offsets[ptId+1] - offsets[ptId];
918 for ( i=0; i<numPtsInGroup; ++i )
919 {
920 eid = mergeArray[offsets[ptId]+i].EId;
921 triId = eid / 3;
922 *(this->Tris + 4*triId + eid-(3*triId) + 1) = ptId + this->TotalPts;
923 }//for this group of coincident edges
924 }//for all merged points
925 }
926
927 // Update the triangle connectivity (numPts for each triangle. This could
928 // be done in parallel but it's probably not faster.
Reduce__anon9cebd3fe0111::ProduceMergedTriangles929 void Reduce()
930 {
931 vtkIdType *tris = this->Tris;
932 for ( IDType triId=0; triId < this->NumTris; ++triId, tris+=4 )
933 {
934 *tris = 3;
935 }
936 }
937 };
938
939 // This method generates the output isosurface points. One point per
940 // merged edge is generated.
941 template <typename TIP, typename TOP, typename IDType>
942 struct ProduceMergedPoints
943 {
944 typedef MergeTuple<IDType,float> MergeTupleType;
945
946 const MergeTupleType *MergeArray;
947 const IDType *Offsets;
948 const TIP *InPts;
949 TOP *OutPts;
950
ProduceMergedPoints__anon9cebd3fe0111::ProduceMergedPoints951 ProduceMergedPoints(const MergeTupleType *merge, const IDType *offsets,
952 TIP *inPts, TOP *outPts, vtkIdType totalPts) :
953 MergeArray(merge), Offsets(offsets), InPts(inPts)
954 {
955 this->OutPts = outPts + 3*totalPts;
956 }
957
operator ()__anon9cebd3fe0111::ProduceMergedPoints958 void operator()(vtkIdType ptId, vtkIdType endPtId)
959 {
960 const MergeTupleType *mergeTuple;
961 IDType v0, v1;
962 const TIP *x0, *x1, *inPts=this->InPts;
963 TOP *x, *outPts=this->OutPts;
964 float t;
965
966 for ( ; ptId < endPtId; ++ptId )
967 {
968 mergeTuple = this->MergeArray + this->Offsets[ptId];
969 v0 = mergeTuple->V0;
970 v1 = mergeTuple->V1;
971 t = mergeTuple->T;
972 x0 = inPts + 3*v0;
973 x1 = inPts + 3*v1;
974 x = outPts + 3*ptId;
975 x[0] = x0[0] + t*(x1[0]-x0[0]);
976 x[1] = x0[1] + t*(x1[1]-x0[1]);
977 x[2] = x0[2] + t*(x1[2]-x0[2]);
978 }
979 }
980 };
981
982 // If requested, interpolate point data attributes. The merge tuple contains an
983 // interpolation value t for the merged edge.
984 template <typename TIds>
985 struct ProduceAttributes
986 {
987 const MergeTuple<TIds,float> *Edges; //all edges, sorted into groups of merged edges
988 const TIds *Offsets; //refer to single, unique, merged edge
989 ArrayList *Arrays; //carry list of attributes to interpolate
990 vtkIdType TotalPts; //total points / multiple contours computed previously
991
ProduceAttributes__anon9cebd3fe0111::ProduceAttributes992 ProduceAttributes(const MergeTuple<TIds,float> *mt, const TIds *offsets,
993 ArrayList *arrays, vtkIdType totalPts) :
994 Edges(mt), Offsets(offsets), Arrays(arrays), TotalPts(totalPts)
995 {}
996
operator ()__anon9cebd3fe0111::ProduceAttributes997 void operator()(vtkIdType ptId, vtkIdType endPtId)
998 {
999 const MergeTuple<TIds,float> *mergeTuple;
1000 TIds v0, v1;
1001 float t;
1002
1003 for ( ; ptId < endPtId; ++ptId )
1004 {
1005 mergeTuple = this->Edges + this->Offsets[ptId];
1006 v0 = mergeTuple->V0;
1007 v1 = mergeTuple->V1;
1008 t = mergeTuple->T;
1009 this->Arrays->InterpolateEdge(v0,v1,t,ptId+this->TotalPts);
1010 }
1011 }
1012 };
1013
1014 // Make the source code a little more readable
1015 #define EXTRACT_MERGED(VTK_type,_type) \
1016 case VTK_type: \
1017 { \
1018 ExtractEdges<TIds,_type> extractEdges(cellIter,(_type*)s,isoValue,newPolys,totalTris); \
1019 EXECUTE_REDUCED_SMPFOR(seqProcessing,numTets,extractEdges); \
1020 numThreads = extractEdges.NumThreadsUsed; \
1021 numTris = extractEdges.NumTris; \
1022 tris = newPolys->GetPointer(); \
1023 mergeEdges = extractEdges.Edges; \
1024 } \
1025 break;
1026
1027 // Wrapper to handle multiple template types for merged processing
1028 template <typename TIds>
ProcessMerged(vtkIdType numTets,vtkPoints * inPts,CellIter * cellIter,int sType,void * s,double isoValue,vtkPoints * outPts,vtkCellArray * newPolys,vtkTypeBool intAttr,vtkDataArray * inScalars,vtkPointData * inPD,vtkPointData * outPD,vtkTypeBool seqProcessing,int & numThreads,vtkIdType totalPts,vtkIdType totalTris)1029 int ProcessMerged(vtkIdType numTets, vtkPoints *inPts, CellIter *cellIter,
1030 int sType, void *s, double isoValue, vtkPoints *outPts,
1031 vtkCellArray *newPolys, vtkTypeBool intAttr, vtkDataArray *inScalars,
1032 vtkPointData *inPD, vtkPointData *outPD, vtkTypeBool seqProcessing,
1033 int &numThreads, vtkIdType totalPts, vtkIdType totalTris)
1034 {
1035 // Extract edges that the contour intersects. Templated on type of scalars.
1036 // List below the explicit choice of scalars that can be processed.
1037 vtkIdType numTris=0, *tris=nullptr;
1038 MergeTuple<TIds,float> *mergeEdges=nullptr; //may need reference counting
1039 switch ( sType )
1040 {
1041 EXTRACT_MERGED(VTK_UNSIGNED_CHAR,unsigned char);
1042 EXTRACT_MERGED(VTK_CHAR,char);
1043 EXTRACT_MERGED(VTK_UNSIGNED_SHORT,unsigned short);
1044 EXTRACT_MERGED(VTK_SHORT,short);
1045 EXTRACT_MERGED(VTK_UNSIGNED_INT,unsigned int);
1046 EXTRACT_MERGED(VTK_INT,int);
1047 EXTRACT_MERGED(VTK_FLOAT,float);
1048 EXTRACT_MERGED(VTK_DOUBLE,double);
1049 default:
1050 vtkGenericWarningMacro(<<"Scalar type not supported");
1051 return 0;
1052 };
1053
1054 // Make sure data was produced
1055 if ( numTris <= 0 )
1056 {
1057 outPts->SetNumberOfPoints(0);
1058 delete [] mergeEdges;
1059 return 1;
1060 }
1061
1062 // Merge coincident edges. The Offsets refer to the single unique edge
1063 // from the sorted group of duplicate edges.
1064 vtkIdType numPts;
1065 vtkStaticEdgeLocatorTemplate<TIds,float> loc;
1066 const TIds *offsets = loc.MergeEdges(3*numTris,mergeEdges,numPts);
1067
1068 // Generate triangles.
1069 ProduceMergedTriangles<TIds> produceTris(mergeEdges,offsets,numTris,tris,
1070 totalPts,totalTris);
1071 EXECUTE_REDUCED_SMPFOR(seqProcessing,numPts,produceTris);
1072
1073 // Generate points (one per unique edge)
1074 outPts->GetData()->WriteVoidPointer(0,3*(numPts+totalPts));
1075 int inPtsType = inPts->GetDataType();
1076 void *inPtsPtr = inPts->GetVoidPointer(0);
1077 int outPtsType = outPts->GetDataType();
1078 void *outPtsPtr = outPts->GetVoidPointer(0);
1079
1080 // Only handle combinations of real types
1081 if ( inPtsType == VTK_FLOAT && outPtsType == VTK_FLOAT )
1082 {
1083 ProduceMergedPoints<float,float,TIds>
1084 producePts(mergeEdges, offsets, (float*)inPtsPtr, (float*)outPtsPtr, totalPts );
1085 EXECUTE_SMPFOR(seqProcessing,numPts,producePts);
1086 }
1087 else if ( inPtsType == VTK_DOUBLE && outPtsType == VTK_DOUBLE )
1088 {
1089 ProduceMergedPoints<double,double,TIds>
1090 producePts(mergeEdges, offsets, (double*)inPtsPtr, (double*)outPtsPtr, totalPts );
1091 EXECUTE_SMPFOR(seqProcessing,numPts,producePts);
1092 }
1093 else if ( inPtsType == VTK_FLOAT && outPtsType == VTK_DOUBLE )
1094 {
1095 ProduceMergedPoints<float,double,TIds>
1096 producePts(mergeEdges, offsets, (float*)inPtsPtr, (double*)outPtsPtr, totalPts );
1097 EXECUTE_SMPFOR(seqProcessing,numPts,producePts);
1098 }
1099 else //if ( inPtsType == VTK_DOUBLE && outPtsType == VTK_FLOAT )
1100 {
1101 ProduceMergedPoints<double,float,TIds>
1102 producePts(mergeEdges, offsets, (double*)inPtsPtr, (float*)outPtsPtr, totalPts );
1103 EXECUTE_SMPFOR(seqProcessing,numPts,producePts);
1104 }
1105
1106 // Now process point data attributes if requested
1107 if ( intAttr )
1108 {
1109 ArrayList arrays;
1110 if ( totalPts <= 0 ) //first contour value generating output
1111 {
1112 outPD->InterpolateAllocate(inPD,numPts);
1113 outPD->RemoveArray(inScalars->GetName());
1114 arrays.ExcludeArray(inScalars);
1115 arrays.AddArrays(numPts,inPD,outPD);
1116 }
1117 else
1118 {
1119 arrays.Realloc(totalPts+numPts);
1120 }
1121 ProduceAttributes<TIds> interpolate(mergeEdges,offsets,&arrays,totalPts);
1122 EXECUTE_SMPFOR(seqProcessing,numPts,interpolate);
1123 }
1124
1125 // Clean up
1126 delete [] mergeEdges;
1127 return 1;
1128 };
1129 #undef EXTRACT_MERGED
1130
1131 // Functor for computing cell normals. Could easily be templated on output
1132 // point type but we are trying to control object size.
1133 struct ComputeCellNormals
1134 {
1135 vtkPoints *Points;
1136 vtkIdType *Tris;
1137 float *CellNormals;
1138
ComputeCellNormals__anon9cebd3fe0111::ComputeCellNormals1139 ComputeCellNormals(vtkPoints *pts, vtkIdType *tris, float *cellNormals) :
1140 Points(pts), Tris(tris), CellNormals(cellNormals)
1141 {}
1142
operator ()__anon9cebd3fe0111::ComputeCellNormals1143 void operator()(vtkIdType triId, vtkIdType endTriId)
1144 {
1145 vtkIdType *tri = this->Tris + 4*triId;
1146 tri++; //move the pointer to the begging of triangle connectivity
1147 float *n = this->CellNormals + 3*triId;
1148 double nd[3];
1149
1150 for ( ; triId < endTriId; ++triId, tri+=4 )
1151 {
1152 vtkTriangle::ComputeNormal(this->Points, 3, tri, nd);
1153 *n++ = nd[0];
1154 *n++ = nd[1];
1155 *n++ = nd[2];
1156 }
1157 }
1158 };
1159
1160 // Generate normals on output triangles
GenerateTriNormals(vtkTypeBool seqProcessing,vtkPoints * pts,vtkCellArray * tris)1161 vtkFloatArray *GenerateTriNormals(vtkTypeBool seqProcessing, vtkPoints *pts, vtkCellArray *tris)
1162 {
1163 vtkIdType numTris = tris->GetNumberOfCells();
1164
1165 vtkFloatArray *cellNormals = vtkFloatArray::New();
1166 cellNormals->SetNumberOfComponents(3);
1167 cellNormals->SetNumberOfTuples(numTris);
1168 float *n = static_cast<float*>(cellNormals->GetVoidPointer(0));
1169
1170 // Execute functor over all triangles
1171 ComputeCellNormals computeNormals(pts,tris->GetPointer(),n);
1172 EXECUTE_SMPFOR(seqProcessing, numTris, computeNormals);
1173
1174 return cellNormals;
1175 }
1176
1177 // Functor for averaging normals at each merged point.
1178 template <typename TId>
1179 struct AverageNormals
1180 {
1181 vtkStaticCellLinksTemplate<TId> *Links;
1182 const float *CellNormals;
1183 float *PointNormals;
1184
AverageNormals__anon9cebd3fe0111::AverageNormals1185 AverageNormals(vtkStaticCellLinksTemplate<TId> *links, float *cellNormals, float *ptNormals) :
1186 Links(links), CellNormals(cellNormals), PointNormals(ptNormals)
1187 {}
1188
operator ()__anon9cebd3fe0111::AverageNormals1189 void operator()(vtkIdType ptId, vtkIdType endPtId)
1190 {
1191 TId i, numTris;
1192 const TId *tris;
1193 const float *nc;
1194 float *n = this->PointNormals + 3*ptId;
1195
1196 for ( ; ptId < endPtId; ++ptId, n+=3 )
1197 {
1198 numTris = this->Links->GetNumberOfCells(ptId);
1199 tris = this->Links->GetCells(ptId);
1200 n[0] = n[1] = n[2] = 0.0;
1201 for ( i=0; i < numTris; ++i)
1202 {
1203 nc = this->CellNormals + 3*tris[i];
1204 n[0] += nc[0];
1205 n[1] += nc[1];
1206 n[2] += nc[2];
1207 }
1208 vtkMath::Normalize(n);
1209 }
1210 }
1211 };
1212
1213 // Generate normals on merged points. Average cell normals at each point.
1214 template <typename TId>
GeneratePointNormals(vtkTypeBool seqProcessing,vtkPoints * pts,vtkCellArray * tris,vtkFloatArray * cellNormals,vtkPointData * pd)1215 void GeneratePointNormals(vtkTypeBool seqProcessing, vtkPoints *pts, vtkCellArray *tris,
1216 vtkFloatArray *cellNormals, vtkPointData *pd)
1217 {
1218 vtkIdType numPts = pts->GetNumberOfPoints();
1219
1220 vtkFloatArray *ptNormals = vtkFloatArray::New();
1221 ptNormals->SetName("Normals");
1222 ptNormals->SetNumberOfComponents(3);
1223 ptNormals->SetNumberOfTuples(numPts);
1224 float *ptN = static_cast<float*>(ptNormals->GetVoidPointer(0));
1225
1226 // Grab the computed triangle normals
1227 float *triN = static_cast<float*>(cellNormals->GetVoidPointer(0));
1228
1229 // Build cell links
1230 vtkPolyData *dummy = vtkPolyData::New();
1231 dummy->SetPoints(pts);
1232 dummy->SetPolys(tris);
1233 vtkStaticCellLinksTemplate<TId> links;
1234 links.BuildLinks(dummy);
1235
1236 // Process all points, averaging normals
1237 AverageNormals<TId> average(&links, triN, ptN);
1238 EXECUTE_SMPFOR(seqProcessing,numPts,average);
1239
1240 // Clean up and get out
1241 dummy->Delete();
1242 pd->SetNormals(ptNormals);
1243 cellNormals->Delete();
1244 ptNormals->Delete();
1245 };
1246
1247 }//anonymous namespace
1248
1249
1250 //-----------------------------------------------------------------------------
1251 // Construct an instance of the class.
vtkContour3DLinearGrid()1252 vtkContour3DLinearGrid::vtkContour3DLinearGrid()
1253 {
1254 this->ContourValues = vtkContourValues::New();
1255
1256 this->OutputPointsPrecision = DEFAULT_PRECISION;
1257
1258 // By default process active point scalars
1259 this->SetInputArrayToProcess(0,0,0,vtkDataObject::FIELD_ASSOCIATION_POINTS,
1260 vtkDataSetAttributes::SCALARS);
1261
1262 this->MergePoints = false;
1263 this->InterpolateAttributes = false;
1264 this->ComputeNormals = false;
1265 this->SequentialProcessing = false;
1266 this->NumberOfThreadsUsed = 0;
1267 this->LargeIds = false;
1268 }
1269
1270 //-----------------------------------------------------------------------------
~vtkContour3DLinearGrid()1271 vtkContour3DLinearGrid::~vtkContour3DLinearGrid()
1272 {
1273 this->ContourValues->Delete();
1274 }
1275
1276 //-----------------------------------------------------------------------------
1277 // Overload standard modified time function. If contour values are modified,
1278 // then this object is modified as well.
GetMTime()1279 vtkMTimeType vtkContour3DLinearGrid::GetMTime()
1280 {
1281 vtkMTimeType mTime=this->Superclass::GetMTime();
1282 vtkMTimeType time;
1283
1284 if (this->ContourValues)
1285 {
1286 time = this->ContourValues->GetMTime();
1287 mTime = ( time > mTime ? time : mTime );
1288 }
1289
1290 return mTime;
1291 }
1292
1293 // Make code more readable
1294 #define EXTRACT_FAST_PATH(VTK_SCALAR_type,_type) \
1295 case VTK_SCALAR_type: \
1296 ProcessFastPath<_type>(numCells, inPts, cellIter, \
1297 (_type*)sPtr, value, outPts, newPolys, \
1298 this->SequentialProcessing,this->NumberOfThreadsUsed, \
1299 totalPts, totalTris); \
1300 break;
1301
1302 //-----------------------------------------------------------------------------
1303 // Specialized contouring filter to handle unstructured grids with 3D linear
1304 // cells (tetrahedras, hexes, wedges, pyradmids, voxels)
1305 //
1306 int vtkContour3DLinearGrid::
RequestData(vtkInformation *,vtkInformationVector ** inputVector,vtkInformationVector * outputVector)1307 RequestData(vtkInformation*, vtkInformationVector** inputVector,
1308 vtkInformationVector* outputVector)
1309 {
1310 // Get the input and output
1311 vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
1312 vtkInformation* outInfo = outputVector->GetInformationObject(0);
1313 vtkUnstructuredGrid *input =
1314 vtkUnstructuredGrid::SafeDownCast(inInfo->Get(vtkDataObject::DATA_OBJECT()));
1315 vtkPolyData *output =
1316 vtkPolyData::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));
1317 if (!input || !output)
1318 {
1319 return 0;
1320 }
1321
1322 // Make sure there is data to process
1323 vtkCellArray *cells = input->GetCells();
1324 vtkIdType numPts, numTris, numCells = cells->GetNumberOfCells();
1325 vtkDataArray *inScalars = this->GetInputArrayToProcess(0, inputVector);
1326 if (!inScalars)
1327 {
1328 return 1;
1329 }
1330 int sType = inScalars->GetDataType();
1331 void *sPtr = inScalars->GetVoidPointer(0);
1332
1333 // Get the contour values.
1334 int numContours = this->ContourValues->GetNumberOfContours();
1335 double value, *values=this->ContourValues->GetValues();
1336 if ( numContours < 1 )
1337 {
1338 vtkDebugMacro(<<"No contour values defined");
1339 return 1;
1340 }
1341
1342 // Check the input point type. Only real types are supported.
1343 vtkPoints *inPts = input->GetPoints();
1344 numPts = inPts->GetNumberOfPoints();
1345 int inPtsType = inPts->GetDataType();
1346 if ( (inPtsType != VTK_FLOAT && inPtsType != VTK_DOUBLE) )
1347 {
1348 vtkErrorMacro(<<"Input point type not supported");
1349 return 0;
1350 }
1351
1352 // Create the output points. Only real types are supported.
1353 vtkPoints *outPts = vtkPoints::New();
1354 if ( this->OutputPointsPrecision == vtkAlgorithm::DEFAULT_PRECISION )
1355 {
1356 outPts->SetDataType(inPts->GetDataType());
1357 }
1358 else if(this->OutputPointsPrecision == vtkAlgorithm::SINGLE_PRECISION)
1359 {
1360 outPts->SetDataType(VTK_FLOAT);
1361 }
1362 else if(this->OutputPointsPrecision == vtkAlgorithm::DOUBLE_PRECISION)
1363 {
1364 outPts->SetDataType(VTK_DOUBLE);
1365 }
1366
1367 // Output triangles go here.
1368 vtkCellArray *newPolys = vtkCellArray::New();
1369
1370 // Process all contour values
1371 vtkIdType totalPts = 0;
1372 vtkIdType totalTris = 0;
1373
1374 // Set up the cells for processing. A specialized iterator is used to traverse the cells.
1375 vtkIdType *conn = cells->GetPointer();
1376 unsigned char *cellTypes = static_cast<unsigned char*>(input->GetCellTypesArray()->GetVoidPointer(0));
1377 vtkIdType *locs = static_cast<vtkIdType*>(input->GetCellLocationsArray()->GetVoidPointer(0));
1378 CellIter *cellIter = new CellIter(numCells,cellTypes,conn,locs);
1379
1380 // Now produce the output: fast path or general path
1381 int mergePoints = this->MergePoints | this->ComputeNormals | this->InterpolateAttributes;
1382 if ( ! mergePoints )
1383 {//fast path
1384 // Generate all of the points at once (for multiple contours) and then produce the triangles.
1385 for (int vidx = 0; vidx < numContours; vidx++)
1386 {
1387 value = values[vidx];
1388 switch ( sType ) //process these scalar types
1389 {
1390 EXTRACT_FAST_PATH(VTK_UNSIGNED_CHAR,unsigned char);
1391 EXTRACT_FAST_PATH(VTK_CHAR,char);
1392 EXTRACT_FAST_PATH(VTK_UNSIGNED_SHORT,unsigned short);
1393 EXTRACT_FAST_PATH(VTK_SHORT,short);
1394 EXTRACT_FAST_PATH(VTK_UNSIGNED_INT,unsigned int);
1395 EXTRACT_FAST_PATH(VTK_INT,int);
1396 EXTRACT_FAST_PATH(VTK_FLOAT,float);
1397 EXTRACT_FAST_PATH(VTK_DOUBLE,double);
1398 default:
1399 vtkGenericWarningMacro(<<"Scalar type not supported");
1400 return 0;
1401 };
1402
1403 // Multiple contour values require accumulating points & triangles
1404 totalPts = outPts->GetNumberOfPoints();
1405 totalTris = newPolys->GetNumberOfCells();
1406 }//for all contours
1407
1408 // From the points create the output triangles. In the fast path there
1409 // are three points for every triangle. Many points are typically
1410 // duplicates but point merging is a significant cost so ignored in the
1411 // fast path.
1412 numTris = newPolys->GetNumberOfCells();
1413 vtkIdType *tris = newPolys->GetPointer();
1414 ProduceTriangles produceTris(tris);
1415 EXECUTE_SMPFOR(this->SequentialProcessing,numTris,produceTris)
1416 }
1417
1418 else //Need to merge points, and possibly perform attribute interpolation
1419 //and generate normals. Hence use the slower path.
1420 {
1421 vtkPointData *inPD = input->GetPointData();
1422 vtkPointData *outPD = output->GetPointData();
1423
1424 // Determine the size/type of point and cell ids needed to index points
1425 // and cells. Using smaller ids results in a greatly reduced memory footprint
1426 // and faster processing.
1427 this->LargeIds = ( numPts >= VTK_INT_MAX || numCells >= VTK_INT_MAX ? true : false );
1428
1429 // Generate all of the merged points and triangles at once (for multiple
1430 // contours) and then produce the normals if requested.
1431 for (int vidx = 0; vidx < numContours; vidx++)
1432 {
1433 value = values[vidx];
1434 if ( this->LargeIds == false )
1435 {
1436 if ( ! ProcessMerged<int>(numCells, inPts, cellIter, sType, sPtr, value,
1437 outPts, newPolys, this->InterpolateAttributes,
1438 inScalars, inPD, outPD, this->SequentialProcessing,
1439 this->NumberOfThreadsUsed, totalPts, totalTris) )
1440 {
1441 return 0;
1442 }
1443 }
1444 else
1445 {
1446 if ( ! ProcessMerged<vtkIdType>(numCells, inPts, cellIter, sType, sPtr, value,
1447 outPts, newPolys, this->InterpolateAttributes,
1448 inScalars, inPD, outPD, this->SequentialProcessing,
1449 this->NumberOfThreadsUsed, totalPts, totalTris) )
1450 {
1451 return 0;
1452 }
1453 }
1454
1455 // Multiple contour values require accumulating points & triangles
1456 totalPts = outPts->GetNumberOfPoints();
1457 totalTris = newPolys->GetNumberOfCells();
1458 }//for all contour values
1459
1460 // If requested, compute normals. Basically triangle normals are averaged
1461 // on each merged point. Requires building static CellLinks so it is a
1462 // relatively expensive operation. (This block of code is separate to
1463 // control .obj object bloat.)
1464 if ( this->ComputeNormals )
1465 {
1466 vtkFloatArray *triNormals =
1467 GenerateTriNormals(this->SequentialProcessing, outPts,newPolys);
1468 if ( this->LargeIds )
1469 {
1470 GeneratePointNormals<vtkIdType>(this->SequentialProcessing,
1471 outPts,newPolys,triNormals,outPD);
1472 }
1473 else
1474 {
1475 GeneratePointNormals<int>(this->SequentialProcessing,
1476 outPts,newPolys,triNormals,outPD);
1477 }
1478 }
1479 }// slower path requires point merging
1480
1481 // Report the results of execution
1482 vtkDebugMacro(<<"Created: " << outPts->GetNumberOfPoints() << " points, "
1483 << newPolys->GetNumberOfCells() << " triangles");
1484
1485 // Clean up
1486 delete cellIter;
1487 output->SetPoints(outPts);
1488 outPts->Delete();
1489 output->SetPolys(newPolys);
1490 newPolys->Delete();
1491
1492 return 1;
1493 }
1494
1495 //-----------------------------------------------------------------------------
SetOutputPointsPrecision(int precision)1496 void vtkContour3DLinearGrid::SetOutputPointsPrecision(int precision)
1497 {
1498 this->OutputPointsPrecision = precision;
1499 this->Modified();
1500 }
1501
1502 //-----------------------------------------------------------------------------
GetOutputPointsPrecision() const1503 int vtkContour3DLinearGrid::GetOutputPointsPrecision() const
1504 {
1505 return this->OutputPointsPrecision;
1506 }
1507
1508 //-----------------------------------------------------------------------------
FillInputPortInformation(int,vtkInformation * info)1509 int vtkContour3DLinearGrid::FillInputPortInformation(int, vtkInformation *info)
1510 {
1511 info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkUnstructuredGrid");
1512 return 1;
1513 }
1514
1515 //-----------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)1516 void vtkContour3DLinearGrid::PrintSelf(ostream& os, vtkIndent indent)
1517 {
1518 this->Superclass::PrintSelf(os,indent);
1519
1520 this->ContourValues->PrintSelf(os,indent.GetNextIndent());
1521
1522 os << indent << "Precision of the output points: "
1523 << this->OutputPointsPrecision << "\n";
1524
1525 os << indent << "Merge Points: "
1526 << (this->MergePoints ? "true\n" : "false\n");
1527 os << indent << "Interpolate Attributes: "
1528 << (this->InterpolateAttributes ? "true\n" : "false\n");
1529 os << indent << "Compute Normals: "
1530 << (this->ComputeNormals ? "true\n" : "false\n");
1531
1532 os << indent << "Sequential Processing: "
1533 << (this->SequentialProcessing ? "true\n" : "false\n");
1534 os << indent << "Large Ids: "
1535 << (this->LargeIds ? "true\n" : "false\n");
1536 }
1537
1538 #undef EXECUTE_SMPFOR
1539 #undef EXECUTE_REDUCED_SMPFOR
1540 #undef MAX_CELL_VERTS
1541 #undef EXTRACT_MERGED
1542 #undef EXTRACT_FAST_PATH
1543