1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkSynchronizedTemplatesCutter3D.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 "vtkSynchronizedTemplatesCutter3D.h"
16
17 #include "vtkCellArray.h"
18 #include "vtkCellData.h"
19 #include "vtkCharArray.h"
20 #include "vtkDoubleArray.h"
21 #include "vtkFloatArray.h"
22 #include "vtkImplicitFunction.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 "vtkShortArray.h"
32 #include "vtkStreamingDemandDrivenPipeline.h"
33 #include "vtkStructuredPoints.h"
34 #include "vtkUnsignedCharArray.h"
35 #include "vtkUnsignedIntArray.h"
36 #include "vtkUnsignedLongArray.h"
37 #include "vtkUnsignedShortArray.h"
38 #include "vtkPolygonBuilder.h"
39 #include "vtkIdList.h"
40 #include "vtkIdListCollection.h"
41 #include "vtkSmartPointer.h"
42
43 #include <cmath>
44
45 vtkStandardNewMacro(vtkSynchronizedTemplatesCutter3D);
46 vtkCxxSetObjectMacro(vtkSynchronizedTemplatesCutter3D,CutFunction,vtkImplicitFunction);
47
48 //----------------------------------------------------------------------------
49 // Description:
50 // Construct object with initial scalar range (0,1) and single contour value
51 // of 0.0. The ImageRange are set to extract the first k-plane.
vtkSynchronizedTemplatesCutter3D()52 vtkSynchronizedTemplatesCutter3D::vtkSynchronizedTemplatesCutter3D()
53 {
54 this->CutFunction = nullptr;
55 this->OutputPointsPrecision = vtkAlgorithm::DEFAULT_PRECISION;
56 }
57
58 //----------------------------------------------------------------------------
~vtkSynchronizedTemplatesCutter3D()59 vtkSynchronizedTemplatesCutter3D::~vtkSynchronizedTemplatesCutter3D()
60 {
61 this->SetCutFunction(nullptr);
62 }
63
64 //----------------------------------------------------------------------------
vtkSynchronizedTemplatesCutter3DInitializeOutput(int * ext,int precision,vtkImageData * input,vtkPolyData * o)65 static void vtkSynchronizedTemplatesCutter3DInitializeOutput(
66 int *ext, int precision, vtkImageData *input, vtkPolyData *o)
67 {
68 vtkPoints *newPts;
69 vtkCellArray *newPolys;
70 long estimatedSize;
71
72 estimatedSize = (int) pow ((double)
73 ((ext[1]-ext[0]+1)*(ext[3]-ext[2]+1)*(ext[5]-ext[4]+1)), .75);
74 if (estimatedSize < 1024)
75 {
76 estimatedSize = 1024;
77 }
78 newPts = vtkPoints::New();
79
80 // set precision for the points in the output
81 if(precision == vtkAlgorithm::DEFAULT_PRECISION)
82 {
83 vtkPointSet *inputPointSet = vtkPointSet::SafeDownCast(input);
84 if(inputPointSet)
85 {
86 newPts->SetDataType(inputPointSet->GetPoints()->GetDataType());
87 }
88 else
89 {
90 newPts->SetDataType(VTK_FLOAT);
91 }
92 }
93 else if(precision == vtkAlgorithm::SINGLE_PRECISION)
94 {
95 newPts->SetDataType(VTK_FLOAT);
96 }
97 else if(precision == vtkAlgorithm::DOUBLE_PRECISION)
98 {
99 newPts->SetDataType(VTK_DOUBLE);
100 }
101
102 newPts->Allocate(estimatedSize,estimatedSize);
103 newPolys = vtkCellArray::New();
104 newPolys->Allocate(newPolys->EstimateSize(estimatedSize,3));
105
106 o->GetPointData()->CopyAllOn();
107
108 o->GetPointData()->InterpolateAllocate(input->GetPointData(),
109 estimatedSize,estimatedSize/2);
110 o->GetCellData()->CopyAllocate(input->GetCellData(),
111 estimatedSize,estimatedSize/2);
112
113 o->SetPoints(newPts);
114 newPts->Delete();
115
116 o->SetPolys(newPolys);
117 newPolys->Delete();
118 }
119
120 //----------------------------------------------------------------------------
121 //
122 // Contouring filter specialized for images
123 //
124 template <class T>
ContourImage(vtkSynchronizedTemplatesCutter3D * self,int * exExt,vtkImageData * data,vtkPolyData * output,T * ptr,bool outputTriangles)125 void ContourImage(vtkSynchronizedTemplatesCutter3D *self, int *exExt,
126 vtkImageData *data, vtkPolyData *output, T *ptr, bool outputTriangles)
127 {
128 int *inExt = data->GetExtent();
129 vtkIdType xdim = exExt[1] - exExt[0] + 1;
130 vtkIdType ydim = exExt[3] - exExt[2] + 1;
131 double *values = self->GetValues();
132 int numContours = self->GetNumberOfContours();
133 T *inPtrX, *inPtrY, *inPtrZ;
134 T *s0, *s1, *s2, *s3;
135 int xMin, xMax, yMin, yMax, zMin, zMax;
136 vtkIdType xInc, yInc, zInc;
137 int xIncFunc, yIncFunc, zIncFunc, scalarZIncFunc;
138 double *origin = data->GetOrigin();
139 double *spacing = data->GetSpacing();
140 vtkIdType *isect1Ptr, *isect2Ptr;
141 double y, z, t;
142 int i, j, k;
143 vtkIdType zstep, yisectstep;
144 vtkIdType offsets[12];
145 int *tablePtr;
146 vtkIdType idx;
147 int vidx;
148 double x[3], xz[3];
149 vtkIdType v0, v1, v2, v3;
150 vtkIdType ptIds[3];
151 double value;
152 // We need to know the edgePointId's for interpolating attributes.
153 vtkIdType edgePtId, inCellId, outCellId;
154 vtkPointData *inPD = data->GetPointData();
155 vtkCellData *inCD = data->GetCellData();
156 vtkPointData *outPD = output->GetPointData();
157 vtkCellData *outCD = output->GetCellData();
158 // Use to be arguments
159 vtkPoints *newPts;
160 vtkCellArray *newPolys;
161 ptr += self->GetArrayComponent();
162 vtkPolygonBuilder polyBuilder;
163 vtkSmartPointer<vtkIdListCollection> polys =
164 vtkSmartPointer<vtkIdListCollection>::New();
165
166 vtkSynchronizedTemplatesCutter3DInitializeOutput(exExt, self->GetOutputPointsPrecision(), data, output);
167 newPts = output->GetPoints();
168 newPolys = output->GetPolys();
169
170 xMin = exExt[0];
171 xMax = exExt[1];
172 yMin = exExt[2];
173 yMax = exExt[3];
174 zMin = exExt[4];
175 zMax = exExt[5];
176
177 vtkImplicitFunction *func = self->GetCutFunction();
178 if (!func)
179 {
180 return;
181 }
182
183 xInc = 1;
184 yInc = xInc*(inExt[1]-inExt[0]+1);
185 zInc = yInc*(inExt[3]-inExt[2]+1);
186
187 // Note that the implicit functions are specified
188 //over exExt so we need to compute the steps differently
189 xIncFunc = 1;
190 yIncFunc = xIncFunc*xdim;
191 zIncFunc = yIncFunc*ydim;
192 scalarZIncFunc = zIncFunc;
193
194 // Kens increments, probably to do with edge array
195 zstep = xdim*ydim;
196 yisectstep = xdim*3;
197 // compute offsets probably how to get to the edges in the edge array.
198 offsets[0] = -xdim*3;
199 offsets[1] = -xdim*3 + 1;
200 offsets[2] = -xdim*3 + 2;
201 offsets[3] = -xdim*3 + 4;
202 offsets[4] = -xdim*3 + 5;
203 offsets[5] = 0;
204 offsets[6] = 2;
205 offsets[7] = 5;
206 offsets[8] = (zstep - xdim)*3;
207 offsets[9] = (zstep - xdim)*3 + 1;
208 offsets[10] = (zstep - xdim)*3 + 4;
209 offsets[11] = zstep*3;
210
211 // allocate storage array
212 vtkIdType *isect1 = new vtkIdType [xdim*ydim*3*2];
213 // set impossible edges to -1
214 for (i = 0; i < ydim; i++)
215 {
216 isect1[(i+1)*xdim*3-3] = -1;
217 isect1[(i+1)*xdim*3*2-3] = -1;
218 }
219 for (i = 0; i < xdim; i++)
220 {
221 isect1[((ydim-1)*xdim + i)*3 + 1] = -1;
222 isect1[((ydim-1)*xdim + i)*3*2 + 1] = -1;
223 }
224
225 // allocate scalar storage for two slices
226 T *scalars = new T [xdim*ydim*2];
227 T *scalars1 = scalars;
228 T *scalars2 = scalars + xdim*ydim;
229
230 // for each contour
231 for (vidx = 0; vidx < numContours; vidx++)
232 {
233 value = values[vidx];
234 inPtrZ = ptr;
235
236 // fill the first slice
237 z = origin[2] + spacing[2]*zMin;
238 x[2] = z;
239 T *scalarsTmp = scalars1;
240 scalars1 = scalars2;
241 scalars2 = scalarsTmp;
242 for (j = yMin; j <= yMax; j++)
243 {
244 x[1] = origin[1] + spacing[1]*j;
245 for (i = xMin; i <= xMax; i++)
246 {
247 x[0] = origin[0] + spacing[0]*i;
248 *scalarsTmp = func->FunctionValue(x);
249 scalarsTmp++;
250 }
251 }
252 scalarZIncFunc = -scalarZIncFunc;
253
254 //==================================================================
255 for (k = zMin; k <= zMax; k++)
256 {
257 self->UpdateProgress((double)vidx/numContours +
258 (k-zMin)/((zMax - zMin+1.0)*numContours));
259 inPtrY = inPtrZ;
260
261 // for each slice compute the scalars
262 z = origin[2] + spacing[2]*(k+1);
263 x[2] = z;
264 scalarsTmp = scalars1;
265 scalars1 = scalars2;
266 scalars2 = scalarsTmp;
267 // if not the last slice then get more scalars
268 if (k < zMax)
269 {
270 for (j = yMin; j <= yMax; j++)
271 {
272 x[1] = origin[1] + spacing[1]*j;
273 for (i = xMin; i <= xMax; i++)
274 {
275 x[0] = origin[0] + spacing[0]*i;
276 *scalarsTmp = func->FunctionValue(x);
277 scalarsTmp++;
278 }
279 }
280 }
281 inPtrY = scalars1;
282 scalarZIncFunc = -scalarZIncFunc;
283
284 z = origin[2] + spacing[2]*k;
285 x[2] = z;
286
287 // swap the buffers
288 if (k%2)
289 {
290 offsets[8] = (zstep - xdim)*3;
291 offsets[9] = (zstep - xdim)*3 + 1;
292 offsets[10] = (zstep - xdim)*3 + 4;
293 offsets[11] = zstep*3;
294 isect1Ptr = isect1;
295 isect2Ptr = isect1 + xdim*ydim*3;
296 }
297 else
298 {
299 offsets[8] = (-zstep - xdim)*3;
300 offsets[9] = (-zstep - xdim)*3 + 1;
301 offsets[10] = (-zstep - xdim)*3 + 4;
302 offsets[11] = -zstep*3;
303 isect1Ptr = isect1 + xdim*ydim*3;
304 isect2Ptr = isect1;
305 }
306
307 for (j = yMin; j <= yMax; j++)
308 {
309 // Should not impact performance here/
310 edgePtId = (xMin-inExt[0])*xInc + (j-inExt[2])*yInc + (k-inExt[4])*zInc;
311
312 // Increments are different for cells. Since the cells are not
313 // contoured until the second row of templates, subtract 1 from
314 // i,j,and k. Note: first cube is formed when i=0, j=1, and k=1.
315 inCellId =
316 (xMin-inExt[0]) + (inExt[1]-inExt[0])*
317 ( (j-inExt[2]-1) + (k-inExt[4]-1)*(inExt[3]-inExt[2]) );
318
319 y = origin[1] + j*spacing[1];
320 xz[1] = y;
321
322 s1 = inPtrY;
323 v1 = (*s1 < value ? 0 : 1);
324
325 inPtrX = inPtrY;
326 for (i = xMin; i <= xMax; i++)
327 {
328 s0 = s1;
329 v0 = v1;
330 *isect2Ptr = -1;
331 *(isect2Ptr + 1) = -1;
332 *(isect2Ptr + 2) = -1;
333 if (i < xMax)
334 {
335 s1 = (inPtrX + xIncFunc);
336 v1 = (*s1 < value ? 0 : 1);
337 if (v0 ^ v1)
338 {
339 // watch for degenerate points
340 if (*s0 == value)
341 {
342 if (i > xMin && *(isect2Ptr-3) > -1)
343 {
344 *isect2Ptr = *(isect2Ptr-3);
345 }
346 else if (j > yMin && *(isect2Ptr - yisectstep + 1) > -1)
347 {
348 *isect2Ptr = *(isect2Ptr - yisectstep + 1);
349 }
350 else if (k > zMin && *(isect1Ptr+2) > -1)
351 {
352 *isect2Ptr = *(isect1Ptr+2);
353 }
354 }
355 else if (*s1 == value)
356 {
357 if (j > yMin && *(isect2Ptr - yisectstep +4) > -1)
358 {
359 *isect2Ptr = *(isect2Ptr - yisectstep + 4);
360 }
361 else if (k > zMin && i < xMax && *(isect1Ptr + 5) > -1)
362 {
363 *isect2Ptr = *(isect1Ptr + 5);
364 }
365 }
366 // if the edge has not been set yet then it is a new point
367 if (*isect2Ptr == -1)
368 {
369 t = (value - (double)(*s0)) / ((double)(*s1) - (double)(*s0));
370 x[0] = origin[0] + spacing[0]*(i+t);
371 x[1] = y;
372 *isect2Ptr = newPts->InsertNextPoint(x);
373 outPD->InterpolateEdge(inPD, *isect2Ptr, edgePtId, edgePtId+1, t);
374 }
375 }
376 }
377 if (j < yMax)
378 {
379 s2 = (inPtrX + yIncFunc);
380 v2 = (*s2 < value ? 0 : 1);
381 if (v0 ^ v2)
382 {
383 if (*s0 == value)
384 {
385 if (*isect2Ptr > -1)
386 {
387 *(isect2Ptr + 1) = *isect2Ptr;
388 }
389 else if (i > xMin && *(isect2Ptr-3) > -1)
390 {
391 *(isect2Ptr + 1) = *(isect2Ptr-3);
392 }
393 else if (j > yMin && *(isect2Ptr - yisectstep + 1) > -1)
394 {
395 *(isect2Ptr + 1) = *(isect2Ptr - yisectstep + 1);
396 }
397 else if (k > zMin && *(isect1Ptr+2) > -1)
398 {
399 *(isect2Ptr + 1) = *(isect1Ptr+2);
400 }
401 }
402 else if (*s2 == value && k > zMin && *(isect1Ptr + yisectstep + 2) > -1)
403 {
404 *(isect2Ptr+1) = *(isect1Ptr + yisectstep + 2);
405 }
406 // if the edge has not been set yet then it is a new point
407 if (*(isect2Ptr + 1) == -1)
408 {
409 t = (value - (double)(*s0)) / ((double)(*s2) - (double)(*s0));
410 x[0] = origin[0] + spacing[0]*i;
411 x[1] = y + spacing[1]*t;
412 *(isect2Ptr + 1) = newPts->InsertNextPoint(x);
413 outPD->InterpolateEdge(inPD, *(isect2Ptr+1), edgePtId, edgePtId+yInc, t);
414 }
415 }
416 }
417 if (k < zMax)
418 {
419 s3 = (inPtrX + scalarZIncFunc);
420 v3 = (*s3 < value ? 0 : 1);
421 if (v0 ^ v3)
422 {
423 if (*s0 == value)
424 {
425 if (*isect2Ptr > -1)
426 {
427 *(isect2Ptr + 2) = *isect2Ptr;
428 }
429 else if (*(isect2Ptr+1) > -1)
430 {
431 *(isect2Ptr + 2) = *(isect2Ptr+1);
432 }
433 else if (i > xMin && *(isect2Ptr-3) > -1)
434 {
435 *(isect2Ptr + 2) = *(isect2Ptr-3);
436 }
437 else if (j > yMin && *(isect2Ptr - yisectstep + 1) > -1)
438 {
439 *(isect2Ptr + 2) = *(isect2Ptr - yisectstep + 1);
440 }
441 else if (k > zMin && *(isect1Ptr+2) > -1)
442 {
443 *(isect2Ptr + 2) = *(isect1Ptr+2);
444 }
445 }
446 if (*(isect2Ptr + 2) == -1)
447 {
448 t = (value - (double)(*s0)) / ((double)(*s3) - (double)(*s0));
449 xz[0] = origin[0] + spacing[0]*i;
450 xz[2] = z + spacing[2]*t;
451 *(isect2Ptr + 2) = newPts->InsertNextPoint(xz);
452 outPD->InterpolateEdge(inPD, *(isect2Ptr+2), edgePtId, edgePtId+zInc, t);
453 }
454 }
455 }
456 // To keep track of ids for interpolating attributes.
457 ++edgePtId;
458
459 // now add any polys that need to be added
460 // basically look at the isect values,
461 // form an index and lookup the polys
462 if (j > yMin && i < xMax && k > zMin)
463 {
464 idx = (v0 ? 4096 : 0);
465 idx = idx + (*(isect1Ptr - yisectstep) > -1 ? 2048 : 0);
466 idx = idx + (*(isect1Ptr -yisectstep +1) > -1 ? 1024 : 0);
467 idx = idx + (*(isect1Ptr -yisectstep +2) > -1 ? 512 : 0);
468 idx = idx + (*(isect1Ptr -yisectstep +4) > -1 ? 256 : 0);
469 idx = idx + (*(isect1Ptr -yisectstep +5) > -1 ? 128 : 0);
470 idx = idx + (*(isect1Ptr) > -1 ? 64 : 0);
471 idx = idx + (*(isect1Ptr + 2) > -1 ? 32 : 0);
472 idx = idx + (*(isect1Ptr + 5) > -1 ? 16 : 0);
473 idx = idx + (*(isect2Ptr -yisectstep) > -1 ? 8 : 0);
474 idx = idx + (*(isect2Ptr -yisectstep +1) > -1 ? 4 : 0);
475 idx = idx + (*(isect2Ptr -yisectstep +4) > -1 ? 2 : 0);
476 idx = idx + (*(isect2Ptr) > -1 ? 1 : 0);
477
478 tablePtr = VTK_SYNCHRONIZED_TEMPLATES_3D_TABLE_2
479 + VTK_SYNCHRONIZED_TEMPLATES_3D_TABLE_1[idx];
480
481 if (!outputTriangles)
482 {
483 polyBuilder.Reset();
484 }
485 while (*tablePtr != -1)
486 {
487 ptIds[0] = *(isect1Ptr + offsets[*tablePtr]);
488 tablePtr++;
489 ptIds[1] = *(isect1Ptr + offsets[*tablePtr]);
490 tablePtr++;
491 ptIds[2] = *(isect1Ptr + offsets[*tablePtr]);
492 tablePtr++;
493 if (ptIds[0] != ptIds[1] &&
494 ptIds[0] != ptIds[2] &&
495 ptIds[1] != ptIds[2])
496 {
497 if(outputTriangles)
498 {
499 outCellId = newPolys->InsertNextCell(3,ptIds);
500 outCD->CopyData(inCD, inCellId, outCellId);
501 }
502 else
503 {
504 polyBuilder.InsertTriangle(ptIds);
505 }
506 }
507 }
508 if(!outputTriangles)
509 {
510 polyBuilder.GetPolygons(polys);
511 int nPolys = polys->GetNumberOfItems();
512 for (int polyId = 0; polyId < nPolys; ++polyId)
513 {
514 vtkIdList* poly = polys->GetItem(polyId);
515 if(poly->GetNumberOfIds()!=0)
516 {
517 outCellId = newPolys->InsertNextCell(poly);
518 outCD->CopyData(inCD, inCellId, outCellId);
519 }
520 poly->Delete();
521 }
522 polys->RemoveAllItems();
523 }
524 }
525 inPtrX += xIncFunc;
526 isect2Ptr += 3;
527 isect1Ptr += 3;
528 // To keep track of ids for copying cell attributes..
529 ++inCellId;
530 }
531 inPtrY += yIncFunc;
532 }
533 inPtrZ += zIncFunc;
534 }
535 }
536 delete [] isect1;
537
538 delete [] scalars;
539 }
540
541 //----------------------------------------------------------------------------
542 //
543 // Contouring filter specialized for images (or slices from images)
544 //
ThreadedExecute(vtkImageData * data,vtkInformation * outInfo,int)545 void vtkSynchronizedTemplatesCutter3D::ThreadedExecute(vtkImageData *data,
546 vtkInformation *outInfo,
547 int)
548 {
549 vtkPolyData *output;
550
551 vtkDebugMacro(<< "Executing Cutter3D structured contour");
552
553 output = vtkPolyData::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));
554
555 int* exExt = data->GetExtent();
556 if ( exExt[0] >= exExt[1] || exExt[2] >= exExt[3] || exExt[4] >= exExt[5] )
557 {
558 vtkDebugMacro(<<"Cutter3D structured contours requires Cutter3D data");
559 return;
560 }
561
562
563 // Check data type and execute appropriate function
564 ContourImage(this, exExt, data, output, (double *)nullptr, this->GenerateTriangles!=0);
565 }
566
567 //----------------------------------------------------------------------------
RequestData(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)568 int vtkSynchronizedTemplatesCutter3D::RequestData(
569 vtkInformation *vtkNotUsed(request),
570 vtkInformationVector **inputVector,
571 vtkInformationVector *outputVector)
572 {
573 // get the info objects
574 vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
575 vtkInformation *outInfo = outputVector->GetInformationObject(0);
576
577 // get the input and output
578 vtkImageData *input = vtkImageData::SafeDownCast(
579 inInfo->Get(vtkDataObject::DATA_OBJECT()));
580 vtkPolyData *output = vtkPolyData::SafeDownCast(
581 outInfo->Get(vtkDataObject::DATA_OBJECT()));
582
583 // Just call the threaded execute directly.
584 this->ThreadedExecute(input, outInfo, 0);
585
586 output->Squeeze();
587
588 return 1;
589 }
590
591
592 //----------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)593 void vtkSynchronizedTemplatesCutter3D::PrintSelf(ostream& os, vtkIndent indent)
594 {
595 this->Superclass::PrintSelf(os,indent);
596 os << indent << "Cut Function: " << this->CutFunction << "\n";
597 os << indent << "Precision of the output points: "
598 << this->OutputPointsPrecision << "\n";
599 }
600
601