1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkRectilinearSynchronizedTemplates.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 "vtkRectilinearSynchronizedTemplates.h"
16
17 #include "vtkCellArray.h"
18 #include "vtkCellData.h"
19 #include "vtkCharArray.h"
20 #include "vtkDoubleArray.h"
21 #include "vtkFloatArray.h"
22 #include "vtkIdListCollection.h"
23 #include "vtkInformation.h"
24 #include "vtkInformationVector.h"
25 #include "vtkIntArray.h"
26 #include "vtkLongArray.h"
27 #include "vtkMath.h"
28 #include "vtkObjectFactory.h"
29 #include "vtkPointData.h"
30 #include "vtkPolyData.h"
31 #include "vtkRectilinearGrid.h"
32 #include "vtkShortArray.h"
33 #include "vtkStreamingDemandDrivenPipeline.h"
34 #include "vtkStructuredPoints.h"
35 #include "vtkSynchronizedTemplates3D.h"
36 #include "vtkUnsignedCharArray.h"
37 #include "vtkUnsignedIntArray.h"
38 #include "vtkUnsignedLongArray.h"
39 #include "vtkUnsignedShortArray.h"
40 #include "vtkPolygonBuilder.h"
41 #include "vtkSmartPointer.h"
42
43 #include <cmath>
44
45 vtkStandardNewMacro(vtkRectilinearSynchronizedTemplates);
46
47 //----------------------------------------------------------------------------
48 // Description:
49 // Construct object with initial scalar range (0,1) and single contour value
50 // of 0.0. The ImageRange are set to extract the first k-plane.
vtkRectilinearSynchronizedTemplates()51 vtkRectilinearSynchronizedTemplates::vtkRectilinearSynchronizedTemplates()
52 {
53 this->ContourValues = vtkContourValues::New();
54 this->ComputeNormals = 1;
55 this->ComputeGradients = 0;
56 this->ComputeScalars = 1;
57 this->GenerateTriangles = 1;
58
59 this->ArrayComponent = 0;
60
61 // by default process active point scalars
62 this->SetInputArrayToProcess(0,0,0,vtkDataObject::FIELD_ASSOCIATION_POINTS,
63 vtkDataSetAttributes::SCALARS);
64 }
65
66 //----------------------------------------------------------------------------
~vtkRectilinearSynchronizedTemplates()67 vtkRectilinearSynchronizedTemplates::~vtkRectilinearSynchronizedTemplates()
68 {
69 this->ContourValues->Delete();
70 }
71
72 //----------------------------------------------------------------------------
73 // Overload standard modified time function. If contour values are modified,
74 // then this object is modified as well.
GetMTime()75 vtkMTimeType vtkRectilinearSynchronizedTemplates::GetMTime()
76 {
77 vtkMTimeType mTime=this->Superclass::GetMTime();
78 vtkMTimeType mTime2=this->ContourValues->GetMTime();
79
80 mTime = ( mTime2 > mTime ? mTime2 : mTime );
81 return mTime;
82 }
83
84 //----------------------------------------------------------------------------
vtkRectilinearSynchronizedTemplatesInitializeOutput(int * ext,vtkRectilinearGrid * input,vtkPolyData * o,vtkFloatArray * scalars,vtkFloatArray * normals,vtkFloatArray * gradients,vtkDataArray * inScalars)85 static void vtkRectilinearSynchronizedTemplatesInitializeOutput(
86 int *ext, vtkRectilinearGrid *input, vtkPolyData *o, vtkFloatArray *scalars,
87 vtkFloatArray *normals, vtkFloatArray *gradients, vtkDataArray *inScalars)
88 {
89 vtkPoints *newPts;
90 vtkCellArray *newPolys;
91 long estimatedSize;
92
93 estimatedSize = (int) pow ((double)
94 ((ext[1]-ext[0]+1)*(ext[3]-ext[2]+1)*(ext[5]-ext[4]+1)), .75);
95 if (estimatedSize < 1024)
96 {
97 estimatedSize = 1024;
98 }
99 newPts = vtkPoints::New();
100 newPts->Allocate(estimatedSize,estimatedSize);
101 newPolys = vtkCellArray::New();
102 newPolys->Allocate(newPolys->EstimateSize(estimatedSize,3));
103
104 o->GetPointData()->CopyAllOn();
105 // It is more efficient to just create the scalar array
106 // rather than redundantly interpolate the scalars.
107 if (input->GetPointData()->GetScalars() == inScalars)
108 {
109 o->GetPointData()->CopyScalarsOff();
110 }
111 else
112 {
113 o->GetPointData()->CopyFieldOff(inScalars->GetName());
114 }
115
116 if (normals)
117 {
118 normals->SetNumberOfComponents(3);
119 normals->Allocate(3*estimatedSize,3*estimatedSize/2);
120 normals->SetName("Normals");
121 }
122 if (gradients)
123 {
124 gradients->SetNumberOfComponents(3);
125 gradients->Allocate(3*estimatedSize,3*estimatedSize/2);
126 gradients->SetName("Gradients");
127 }
128 if (scalars)
129 {
130 // A temporary name.
131 scalars->SetName("Scalars");
132 }
133
134 o->GetPointData()->InterpolateAllocate(input->GetPointData(),
135 estimatedSize,estimatedSize/2);
136 o->GetCellData()->CopyAllocate(input->GetCellData(),
137 estimatedSize,estimatedSize/2);
138
139 o->SetPoints(newPts);
140 newPts->Delete();
141
142 o->SetPolys(newPolys);
143 newPolys->Delete();
144 }
145
146 //----------------------------------------------------------------------------
147 // Calculate the gradient using central difference.
148 template <class T>
vtkRSTComputePointGradient(int i,int j,int k,T * s,int * inExt,int xInc,int yInc,int zInc,double * spacing,double n[3])149 void vtkRSTComputePointGradient(int i, int j, int k, T *s, int *inExt,
150 int xInc, int yInc, int zInc,
151 double *spacing, double n[3])
152 {
153 double sp, sm;
154
155 // x-direction
156 if ( i == inExt[0] )
157 {
158 sp = *(s+xInc);
159 sm = *s;
160 n[0] = (sp - sm) / spacing[1];
161 }
162 else if ( i == inExt[1] )
163 {
164 sp = *s;
165 sm = *(s-xInc);
166 n[0] = (sp - sm) / spacing[0];
167 }
168 else
169 {
170 sp = *(s+xInc);
171 sm = *(s-xInc);
172 n[0] = (sp - sm) / (spacing[0]+spacing[1]);
173 }
174
175 // y-direction
176 if ( j == inExt[2] )
177 {
178 sp = *(s+yInc);
179 sm = *s;
180 n[1] = (sp - sm) / spacing[3];
181 }
182 else if ( j == inExt[3] )
183 {
184 sp = *s;
185 sm = *(s-yInc);
186 n[1] = (sp - sm) / spacing[2];
187 }
188 else
189 {
190 sp = *(s+yInc);
191 sm = *(s-yInc);
192 n[1] = (sp - sm) / (spacing[2]+spacing[3]);
193 }
194
195 // z-direction
196 if ( k == inExt[4] )
197 {
198 sp = *(s+zInc);
199 sm = *s;
200 n[2] = (sp - sm) / spacing[5];
201 }
202 else if ( k == inExt[5] )
203 {
204 sp = *s;
205 sm = *(s-zInc);
206 n[2] = (sp - sm) / spacing[4];
207 }
208 else
209 {
210 sp = *(s+zInc);
211 sm = *(s-zInc);
212 n[2] = (sp - sm) / (spacing[4]+spacing[5]);
213 }
214 }
215
216 //----------------------------------------------------------------------------
217 #define VTK_RECT_CSP3PA(i2,j2,k2,s) \
218 if (NeedGradients) \
219 { \
220 if (!g0) \
221 { \
222 self->ComputeSpacing(data, i, j, k, exExt, spacing); \
223 vtkRSTComputePointGradient(i, j, k, s0, inExt, xInc, yInc, zInc, spacing, n0); \
224 g0 = 1; \
225 } \
226 self->ComputeSpacing(data, i2, j2, k2, exExt, spacing); \
227 vtkRSTComputePointGradient(i2, j2, k2, s, inExt, xInc, yInc, zInc, spacing, n1); \
228 for (jj=0; jj<3; jj++) \
229 { \
230 n[jj] = n0[jj] + t * (n1[jj] - n0[jj]); \
231 } \
232 if (ComputeGradients) \
233 { \
234 newGradients->InsertNextTuple(n); \
235 } \
236 if (ComputeNormals) \
237 { \
238 vtkMath::Normalize(n); \
239 n[0] = -n[0]; n[1] = -n[1]; n[2] = -n[2]; \
240 newNormals->InsertNextTuple(n); \
241 } \
242 } \
243 if (ComputeScalars) \
244 { \
245 newScalars->InsertNextTuple(&value); \
246 }
247
248 //----------------------------------------------------------------------------
249 //
250 // Contouring filter specialized for images
251 //
252 template <class T>
ContourRectilinearGrid(vtkRectilinearSynchronizedTemplates * self,int * exExt,vtkRectilinearGrid * data,vtkPolyData * output,T * ptr,vtkDataArray * inScalars,bool outputTriangles)253 void ContourRectilinearGrid(vtkRectilinearSynchronizedTemplates *self, int *exExt,
254 vtkRectilinearGrid *data, vtkPolyData *output, T *ptr,
255 vtkDataArray *inScalars, bool outputTriangles)
256 {
257 int *inExt = data->GetExtent();
258 int xdim = exExt[1] - exExt[0] + 1;
259 int ydim = exExt[3] - exExt[2] + 1;
260 double *values = self->GetValues();
261 int numContours = self->GetNumberOfContours();
262 T *inPtrX, *inPtrY, *inPtrZ;
263 T *s0, *s1, *s2, *s3;
264 int xMin, xMax, yMin, yMax, zMin, zMax;
265 int xInc, yInc, zInc;
266 int *isect1Ptr, *isect2Ptr;
267 double y, z, t;
268 int i, j, k;
269 int zstep, yisectstep;
270 int offsets[12];
271 int ComputeNormals = self->GetComputeNormals();
272 int ComputeGradients = self->GetComputeGradients();
273 int ComputeScalars = self->GetComputeScalars();
274 int NeedGradients = ComputeGradients || ComputeNormals;
275 double n[3], n0[3], n1[3];
276 int jj, g0;
277 int *tablePtr;
278 int idx, vidx;
279 double x[3], xz[3];
280 int v0, v1, v2, v3;
281 vtkIdType ptIds[3];
282 double value;
283 // We need to know the edgePointId's for interpolating attributes.
284 int edgePtId, inCellId, outCellId;
285 vtkPointData *inPD = data->GetPointData();
286 vtkCellData *inCD = data->GetCellData();
287 vtkPointData *outPD = output->GetPointData();
288 vtkCellData *outCD = output->GetCellData();
289 // Use to be arguments
290 vtkFloatArray *newScalars = nullptr;
291 vtkFloatArray *newNormals = nullptr;
292 vtkFloatArray *newGradients = nullptr;
293 vtkPoints *newPts;
294 vtkCellArray *newPolys;
295 ptr += self->GetArrayComponent();
296 vtkDataArray *xCoords = data->GetXCoordinates();
297 vtkDataArray *yCoords = data->GetYCoordinates();
298 vtkDataArray *zCoords = data->GetZCoordinates();
299 double x1, x2, y2, z2;
300 double spacing[6];
301 vtkPolygonBuilder polyBuilder;
302 vtkSmartPointer<vtkIdListCollection> polys =
303 vtkSmartPointer<vtkIdListCollection>::New();
304
305 if (ComputeScalars)
306 {
307 newScalars = vtkFloatArray::New();
308 }
309 if (ComputeNormals)
310 {
311 newNormals = vtkFloatArray::New();
312 }
313 if (ComputeGradients)
314 {
315 newGradients = vtkFloatArray::New();
316 }
317 vtkRectilinearSynchronizedTemplatesInitializeOutput(exExt, data, output,
318 newScalars, newNormals, newGradients, inScalars);
319 newPts = output->GetPoints();
320 newPolys = output->GetPolys();
321
322 // this is an exploded execute extent.
323 xMin = exExt[0];
324 xMax = exExt[1];
325 yMin = exExt[2];
326 yMax = exExt[3];
327 zMin = exExt[4];
328 zMax = exExt[5];
329
330 // increments to move through scalars Compute these ourself because
331 // we may be contouring an array other than scalars.
332
333 xInc = inScalars->GetNumberOfComponents();
334 yInc = xInc*(inExt[1]-inExt[0]+1);
335 zInc = yInc*(inExt[3]-inExt[2]+1);
336
337 // Kens increments, probably to do with edge array
338 zstep = xdim*ydim;
339 yisectstep = xdim*3;
340 // compute offsets probably how to get to the edges in the edge array.
341 offsets[0] = -xdim*3;
342 offsets[1] = -xdim*3 + 1;
343 offsets[2] = -xdim*3 + 2;
344 offsets[3] = -xdim*3 + 4;
345 offsets[4] = -xdim*3 + 5;
346 offsets[5] = 0;
347 offsets[6] = 2;
348 offsets[7] = 5;
349 offsets[8] = (zstep - xdim)*3;
350 offsets[9] = (zstep - xdim)*3 + 1;
351 offsets[10] = (zstep - xdim)*3 + 4;
352 offsets[11] = zstep*3;
353
354 // allocate storage array
355 int *isect1 = new int [xdim*ydim*3*2];
356 // set impossible edges to -1
357 for (i = 0; i < ydim; i++)
358 {
359 isect1[(i+1)*xdim*3-3] = -1;
360 isect1[(i+1)*xdim*3*2-3] = -1;
361 }
362 for (i = 0; i < xdim; i++)
363 {
364 isect1[((ydim-1)*xdim + i)*3 + 1] = -1;
365 isect1[((ydim-1)*xdim + i)*3*2 + 1] = -1;
366 }
367
368 // for each contour
369 for (vidx = 0; vidx < numContours; vidx++)
370 {
371 value = values[vidx];
372 inPtrZ = ptr;
373 s2 = inPtrZ;
374
375 //==================================================================
376 for (k = zMin; k <= zMax; k++)
377 {
378 self->UpdateProgress((double)vidx/numContours +
379 (k-zMin)/((zMax - zMin+1.0)*numContours));
380
381 z = zCoords->GetComponent(k-inExt[4], 0);
382 x[2] = z;
383
384 // swap the buffers
385 if (k%2)
386 {
387 offsets[8] = (zstep - xdim)*3;
388 offsets[9] = (zstep - xdim)*3 + 1;
389 offsets[10] = (zstep - xdim)*3 + 4;
390 offsets[11] = zstep*3;
391 isect1Ptr = isect1;
392 isect2Ptr = isect1 + xdim*ydim*3;
393 }
394 else
395 {
396 offsets[8] = (-zstep - xdim)*3;
397 offsets[9] = (-zstep - xdim)*3 + 1;
398 offsets[10] = (-zstep - xdim)*3 + 4;
399 offsets[11] = -zstep*3;
400 isect1Ptr = isect1 + xdim*ydim*3;
401 isect2Ptr = isect1;
402 }
403
404 inPtrY = inPtrZ;
405 for (j = yMin; j <= yMax; j++)
406 {
407 // Should not impact performance here/
408 edgePtId = (j-inExt[2])*yInc + (k-inExt[4])*zInc;
409 // Increments are different for cells.
410 // Since the cells are not contoured until the second row of templates,
411 // subtract 1 from i,j,and k. Note: first cube is formed when i=0, j=1, and k=1.
412 inCellId = (xMin-inExt[0]) + (inExt[1]-inExt[0])*( (j-inExt[2]-1) + (k-inExt[4]-1)*(inExt[3]-inExt[2]) );
413
414 y = yCoords->GetComponent(j-inExt[2], 0);
415 xz[1] = y;
416
417 s1 = inPtrY;
418 v1 = (*s1 < value ? 0 : 1);
419
420 inPtrX = inPtrY;
421 for (i = xMin; i <= xMax; i++)
422 {
423 s0 = s1;
424 v0 = v1;
425 // this flag keeps up from computing gradient for grid point 0 twice.
426 g0 = 0;
427 *isect2Ptr = -1;
428 *(isect2Ptr + 1) = -1;
429 *(isect2Ptr + 2) = -1;
430 if (i < xMax)
431 {
432 s1 = (inPtrX + xInc);
433 v1 = (*s1 < value ? 0 : 1);
434 if (v0 ^ v1)
435 {
436 // watch for degenerate points
437 if (*s0 == value)
438 {
439 if (i > xMin && *(isect2Ptr-3) > -1)
440 {
441 *isect2Ptr = *(isect2Ptr-3);
442 }
443 else if (j > yMin && *(isect2Ptr - yisectstep + 1) > -1)
444 {
445 *isect2Ptr = *(isect2Ptr - yisectstep + 1);
446 }
447 else if (k > zMin && *(isect1Ptr+2) > -1)
448 {
449 *isect2Ptr = *(isect1Ptr+2);
450 }
451 }
452 else if (*s1 == value)
453 {
454 if (j > yMin && *(isect2Ptr - yisectstep +4) > -1)
455 {
456 *isect2Ptr = *(isect2Ptr - yisectstep + 4);
457 }
458 else if (k > zMin && i < xMax && *(isect1Ptr + 5) > -1)
459 {
460 *isect2Ptr = *(isect1Ptr + 5);
461 }
462 }
463 // if the edge has not been set yet then it is a new point
464 if (*isect2Ptr == -1)
465 {
466 t = (value - (double)(*s0)) / ((double)(*s1) - (double)(*s0));
467 x1 = xCoords->GetComponent(i-inExt[0], 0);
468 x2 = xCoords->GetComponent(i-inExt[0]+1, 0);
469 x[0] = x1 + t*(x2-x1);
470 x[1] = y;
471
472 *isect2Ptr = newPts->InsertNextPoint(x);
473 VTK_RECT_CSP3PA(i+1,j,k,s1);
474 outPD->InterpolateEdge(inPD, *isect2Ptr, edgePtId, edgePtId+1, t);
475 }
476 }
477 }
478 if (j < yMax)
479 {
480 s2 = (inPtrX + yInc);
481 v2 = (*s2 < value ? 0 : 1);
482 if (v0 ^ v2)
483 {
484 // watch for degen points
485 if (*s0 == value)
486 {
487 if (*isect2Ptr > -1)
488 {
489 *(isect2Ptr + 1) = *isect2Ptr;
490 }
491 else if (i > xMin && *(isect2Ptr-3) > -1)
492 {
493 *(isect2Ptr + 1) = *(isect2Ptr-3);
494 }
495 else if (j > yMin && *(isect2Ptr - yisectstep + 1) > -1)
496 {
497 *(isect2Ptr + 1) = *(isect2Ptr - yisectstep + 1);
498 }
499 else if (k > zMin && *(isect1Ptr+2) > -1)
500 {
501 *(isect2Ptr + 1) = *(isect1Ptr+2);
502 }
503 }
504 else if (*s2 == value && k > zMin && *(isect1Ptr + yisectstep + 2) > -1)
505 {
506 *(isect2Ptr+1) = *(isect1Ptr + yisectstep + 2);
507 }
508 // if the edge has not been set yet then it is a new point
509 if (*(isect2Ptr + 1) == -1)
510 {
511 t = (value - (double)(*s0)) / ((double)(*s2) - (double)(*s0));
512 x[0] = xCoords->GetComponent(i-inExt[0], 0);
513
514 y2 = yCoords->GetComponent(j-inExt[2]+1, 0);
515 x[1] = y + t*(y2-y);
516
517 *(isect2Ptr + 1) = newPts->InsertNextPoint(x);
518 VTK_RECT_CSP3PA(i,j+1,k,s2);
519 outPD->InterpolateEdge(inPD, *(isect2Ptr+1), edgePtId, edgePtId+yInc, t);
520 }
521 }
522 }
523 if (k < zMax)
524 {
525 s3 = (inPtrX + zInc);
526 v3 = (*s3 < value ? 0 : 1);
527 if (v0 ^ v3)
528 {
529 // watch for degen points
530 if (*s0 == value)
531 {
532 if (*isect2Ptr > -1)
533 {
534 *(isect2Ptr + 2) = *isect2Ptr;
535 }
536 else if (*(isect2Ptr+1) > -1)
537 {
538 *(isect2Ptr + 2) = *(isect2Ptr+1);
539 }
540 else if (i > xMin && *(isect2Ptr-3) > -1)
541 {
542 *(isect2Ptr + 2) = *(isect2Ptr-3);
543 }
544 else if (j > yMin && *(isect2Ptr - yisectstep + 1) > -1)
545 {
546 *(isect2Ptr + 2) = *(isect2Ptr - yisectstep + 1);
547 }
548 else if (k > zMin && *(isect1Ptr+2) > -1)
549 {
550 *(isect2Ptr + 2) = *(isect1Ptr+2);
551 }
552 }
553 if (*(isect2Ptr + 2) == -1)
554 {
555 t = (value - (double)(*s0)) / ((double)(*s3) - (double)(*s0));
556 xz[0] = xCoords->GetComponent(i-inExt[0], 0);
557
558 z2 = zCoords->GetComponent(k-inExt[4]+1, 0);
559 xz[2] = z + t*(z2-z);
560
561 *(isect2Ptr + 2) = newPts->InsertNextPoint(xz);
562 VTK_RECT_CSP3PA(i,j,k+1,s3);
563 outPD->InterpolateEdge(inPD, *(isect2Ptr+2), edgePtId, edgePtId+zInc, t);
564 }
565 }
566 }
567 // To keep track of ids for interpolating attributes.
568 ++edgePtId;
569
570 // now add any polys that need to be added
571 // basically look at the isect values,
572 // form an index and lookup the polys
573 if (j > yMin && i < xMax && k > zMin)
574 {
575 idx = (v0 ? 4096 : 0);
576 idx = idx + (*(isect1Ptr - yisectstep) > -1 ? 2048 : 0);
577 idx = idx + (*(isect1Ptr -yisectstep +1) > -1 ? 1024 : 0);
578 idx = idx + (*(isect1Ptr -yisectstep +2) > -1 ? 512 : 0);
579 idx = idx + (*(isect1Ptr -yisectstep +4) > -1 ? 256 : 0);
580 idx = idx + (*(isect1Ptr -yisectstep +5) > -1 ? 128 : 0);
581 idx = idx + (*(isect1Ptr) > -1 ? 64 : 0);
582 idx = idx + (*(isect1Ptr + 2) > -1 ? 32 : 0);
583 idx = idx + (*(isect1Ptr + 5) > -1 ? 16 : 0);
584 idx = idx + (*(isect2Ptr -yisectstep) > -1 ? 8 : 0);
585 idx = idx + (*(isect2Ptr -yisectstep +1) > -1 ? 4 : 0);
586 idx = idx + (*(isect2Ptr -yisectstep +4) > -1 ? 2 : 0);
587 idx = idx + (*(isect2Ptr) > -1 ? 1 : 0);
588
589 tablePtr = VTK_SYNCHRONIZED_TEMPLATES_3D_TABLE_2
590 + VTK_SYNCHRONIZED_TEMPLATES_3D_TABLE_1[idx];
591
592 if (!outputTriangles)
593 {
594 polyBuilder.Reset();
595 }
596 while (*tablePtr != -1)
597 {
598 ptIds[0] = *(isect1Ptr + offsets[*tablePtr]);
599 tablePtr++;
600 ptIds[1] = *(isect1Ptr + offsets[*tablePtr]);
601 tablePtr++;
602 ptIds[2] = *(isect1Ptr + offsets[*tablePtr]);
603 tablePtr++;
604 if (ptIds[0] != ptIds[1] &&
605 ptIds[0] != ptIds[2] &&
606 ptIds[1] != ptIds[2])
607 {
608 if(outputTriangles)
609 {
610 outCellId = newPolys->InsertNextCell(3,ptIds);
611 outCD->CopyData(inCD, inCellId, outCellId);
612 }
613 else
614 {
615 polyBuilder.InsertTriangle(ptIds);
616 }
617 }
618 }
619 if(!outputTriangles)
620 {
621 polyBuilder.GetPolygons(polys);
622 int nPolys = polys->GetNumberOfItems();
623 for (int polyId = 0; polyId < nPolys; ++polyId)
624 {
625 vtkIdList* poly = polys->GetItem(polyId);
626 if(poly->GetNumberOfIds()!=0)
627 {
628 outCellId = newPolys->InsertNextCell(poly);
629 outCD->CopyData(inCD, inCellId, outCellId);
630 }
631 poly->Delete();
632 }
633 polys->RemoveAllItems();
634 }
635 }
636 inPtrX += xInc;
637 isect2Ptr += 3;
638 isect1Ptr += 3;
639 // To keep track of ids for copying cell attributes.
640 ++inCellId;
641 }
642 inPtrY += yInc;
643 }
644 inPtrZ += zInc;
645 }
646 }
647 delete [] isect1;
648
649 if (newScalars)
650 {
651 // Lets set the name of the scalars here.
652 if (inScalars)
653 {
654 newScalars->SetName(inScalars->GetName());
655 }
656 idx = output->GetPointData()->AddArray(newScalars);
657 output->GetPointData()->SetActiveAttribute(idx, vtkDataSetAttributes::SCALARS);
658 newScalars->Delete();
659 newScalars = nullptr;
660 }
661 if (newGradients)
662 {
663 output->GetPointData()->SetVectors(newGradients);
664 newGradients->Delete();
665 newGradients = nullptr;
666 }
667 if (newNormals)
668 {
669 output->GetPointData()->SetNormals(newNormals);
670 newNormals->Delete();
671 newNormals = nullptr;
672 }
673 }
674
675 //----------------------------------------------------------------------------
676 //
677 // Contouring filter specialized for images (or slices from images)
678 //
RequestData(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)679 int vtkRectilinearSynchronizedTemplates::RequestData(
680 vtkInformation *vtkNotUsed(request),
681 vtkInformationVector **inputVector,
682 vtkInformationVector *outputVector)
683 {
684 // get the info objects
685 vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
686 vtkInformation *outInfo = outputVector->GetInformationObject(0);
687
688 // get the input and output
689 vtkRectilinearGrid *data = vtkRectilinearGrid::SafeDownCast(
690 inInfo->Get(vtkDataObject::DATA_OBJECT()));
691 vtkPolyData *output = vtkPolyData::SafeDownCast(
692 outInfo->Get(vtkDataObject::DATA_OBJECT()));
693
694 void *ptr;
695 vtkDataArray *inScalars;
696
697 vtkDebugMacro(<< "Executing 3D structured contour");
698
699 //
700 // Check data type and execute appropriate function
701 //
702 inScalars = this->GetInputArrayToProcess(0,inputVector);
703 if (inScalars == nullptr)
704 {
705 vtkErrorMacro("No scalars for contouring.");
706 return 1;
707 }
708 int numComps = inScalars->GetNumberOfComponents();
709
710 if (this->ArrayComponent >= numComps)
711 {
712 vtkErrorMacro("Scalars have " << numComps << " components. "
713 "ArrayComponent must be smaller than " << numComps);
714 return 1;
715 }
716
717 int* inExt = data->GetExtent();
718 ptr = this->GetScalarsForExtent(inScalars, inExt, data);
719
720 int exExt[6];
721 inInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), exExt);
722 for (int i=0; i<3; i++)
723 {
724 if (inExt[2*i] > exExt[2*i])
725 {
726 exExt[2*i] = inExt[2*i];
727 }
728 if (inExt[2*i+1] < exExt[2*i+1])
729 {
730 exExt[2*i+1] = inExt[2*i+1];
731 }
732 }
733
734 switch (inScalars->GetDataType())
735 {
736 vtkTemplateMacro(
737 ContourRectilinearGrid(this, exExt, data,
738 output, (VTK_TT *)ptr, inScalars,this->GenerateTriangles!=0));
739 }
740
741 return 1;
742 }
743
744 //----------------------------------------------------------------------------
RequestUpdateExtent(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)745 int vtkRectilinearSynchronizedTemplates::RequestUpdateExtent(
746 vtkInformation *vtkNotUsed(request),
747 vtkInformationVector **inputVector,
748 vtkInformationVector *outputVector)
749 {
750 // These require extra ghost levels
751 if (this->ComputeGradients || this->ComputeNormals)
752 {
753 vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
754 vtkInformation *outInfo = outputVector->GetInformationObject(0);
755
756 int ghostLevels;
757 ghostLevels =
758 outInfo->Get(
759 vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS());
760 inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS(),
761 ghostLevels + 1);
762 }
763
764 return 1;
765 }
766
767 //----------------------------------------------------------------------------
GetScalarsForExtent(vtkDataArray * array,int extent[6],vtkRectilinearGrid * input)768 void* vtkRectilinearSynchronizedTemplates::GetScalarsForExtent(
769 vtkDataArray *array, int extent[6], vtkRectilinearGrid *input)
770 {
771 if ( ! array )
772 {
773 return nullptr;
774 }
775
776 int increments[3], iExt[6], idx;
777
778 input->GetExtent(iExt);
779
780 for (idx = 0; idx < 3; idx++)
781 {
782 if (extent[idx*2] < iExt[idx*2] ||
783 extent[idx*2] > iExt[idx*2+1])
784 {
785 vtkErrorMacro("requested extent not in input's extent");
786 return nullptr;
787 }
788 }
789
790 increments[0] = array->GetNumberOfComponents();
791 increments[1] = increments[0] * (iExt[1]-iExt[0]+1);
792 increments[2] = increments[1] * (iExt[3]-iExt[2]+1);
793
794 idx = (extent[0] - iExt[0]) * increments[0] +
795 (extent[2] - iExt[2]) * increments[1] +
796 (extent[4] - iExt[4]) * increments[2];
797
798 if (idx < 0 || idx > array->GetMaxId())
799 {
800 vtkErrorMacro("computed coordinate outside of array bounds");
801 return nullptr;
802 }
803
804 return array->GetVoidPointer(idx);
805 }
806
807 //----------------------------------------------------------------------------
ComputeSpacing(vtkRectilinearGrid * data,int i,int j,int k,int extent[6],double spacing[6])808 void vtkRectilinearSynchronizedTemplates::ComputeSpacing(
809 vtkRectilinearGrid *data, int i, int j, int k, int extent[6],
810 double spacing[6])
811 {
812 vtkDataArray *xCoords = data->GetXCoordinates();
813 vtkDataArray *yCoords = data->GetYCoordinates();
814 vtkDataArray *zCoords = data->GetZCoordinates();
815
816 spacing[0] = 0;
817 spacing[1] = 0;
818 spacing[2] = 0;
819 spacing[3] = 0;
820 spacing[4] = 0;
821 spacing[5] = 0;
822
823 if (i > extent[0])
824 {
825 spacing[0] = xCoords->GetComponent(i-extent[0], 0) -
826 xCoords->GetComponent(i-extent[0]-1, 0);
827 }
828 if (i < extent[1])
829 {
830 spacing[1] = xCoords->GetComponent(i-extent[0]+1, 0) -
831 xCoords->GetComponent(i-extent[0], 0);
832 }
833 if (j > extent[2])
834 {
835 spacing[2] = yCoords->GetComponent(j-extent[2], 0) -
836 yCoords->GetComponent(j-extent[2]-1, 0);
837 }
838 if (j < extent[3])
839 {
840 spacing[3] = yCoords->GetComponent(j-extent[2]+1, 0) -
841 yCoords->GetComponent(j-extent[2], 0);
842 }
843 if (k > extent[4])
844 {
845 spacing[4] = zCoords->GetComponent(k-extent[4], 0) -
846 zCoords->GetComponent(k-extent[4]-1, 0);
847 }
848 if (k < extent[5])
849 {
850 spacing[5] = zCoords->GetComponent(k-extent[4]+1, 0) -
851 zCoords->GetComponent(k-extent[4], 0);
852 }
853 }
854
855 //----------------------------------------------------------------------------
FillInputPortInformation(int,vtkInformation * info)856 int vtkRectilinearSynchronizedTemplates::FillInputPortInformation(
857 int, vtkInformation *info)
858 {
859 info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkRectilinearGrid");
860 return 1;
861 }
862
863 //----------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)864 void vtkRectilinearSynchronizedTemplates::PrintSelf(ostream& os, vtkIndent indent)
865 {
866 this->Superclass::PrintSelf(os,indent);
867
868 this->ContourValues->PrintSelf(os,indent.GetNextIndent());
869
870 os << indent << "Compute Normals: " << (this->ComputeNormals ? "On\n" : "Off\n");
871 os << indent << "Compute Gradients: " << (this->ComputeGradients ? "On\n" : "Off\n");
872 os << indent << "Compute Scalars: " << (this->ComputeScalars ? "On\n" : "Off\n");
873 os << indent << "ArrayComponent: " << this->ArrayComponent << endl;
874 }
875
876