1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkClipClosedSurface.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 "vtkClipClosedSurface.h"
16
17 #include "vtkCellArray.h"
18 #include "vtkCellArrayIterator.h"
19 #include "vtkCellData.h"
20 #include "vtkContourTriangulator.h"
21 #include "vtkDataSet.h"
22 #include "vtkDoubleArray.h"
23 #include "vtkImageData.h"
24 #include "vtkInformation.h"
25 #include "vtkInformationVector.h"
26 #include "vtkLine.h"
27 #include "vtkMath.h"
28 #include "vtkMatrix4x4.h"
29 #include "vtkObjectFactory.h"
30 #include "vtkPlaneCollection.h"
31 #include "vtkPointData.h"
32 #include "vtkPoints.h"
33 #include "vtkPolyData.h"
34 #include "vtkPolygon.h"
35 #include "vtkSignedCharArray.h"
36 #include "vtkTriangleStrip.h"
37 #include "vtkUnsignedCharArray.h"
38
39 #include "vtkIncrementalOctreePointLocator.h"
40
41 #include <algorithm>
42 #include <map>
43 #include <utility>
44 #include <vector>
45
46 vtkStandardNewMacro(vtkClipClosedSurface);
47
48 vtkCxxSetObjectMacro(vtkClipClosedSurface, ClippingPlanes, vtkPlaneCollection);
49
50 //------------------------------------------------------------------------------
vtkClipClosedSurface()51 vtkClipClosedSurface::vtkClipClosedSurface()
52 {
53 this->ClippingPlanes = nullptr;
54 this->Tolerance = 1e-6;
55 this->PassPointData = 0;
56
57 this->ScalarMode = VTK_CCS_SCALAR_MODE_NONE;
58 this->GenerateOutline = 0;
59 this->GenerateFaces = 1;
60 this->ActivePlaneId = -1;
61
62 this->BaseColor[0] = 1.0;
63 this->BaseColor[1] = 0.0;
64 this->BaseColor[2] = 0.0;
65
66 this->ClipColor[0] = 1.0;
67 this->ClipColor[1] = 0.5;
68 this->ClipColor[2] = 0.0;
69
70 this->ActivePlaneColor[0] = 1.0;
71 this->ActivePlaneColor[1] = 1.0;
72 this->ActivePlaneColor[2] = 0.0;
73
74 this->TriangulationErrorDisplay = 0;
75
76 // A whole bunch of objects needed during execution
77 this->IdList = nullptr;
78 }
79
80 //------------------------------------------------------------------------------
~vtkClipClosedSurface()81 vtkClipClosedSurface::~vtkClipClosedSurface()
82 {
83 if (this->ClippingPlanes)
84 {
85 this->ClippingPlanes->Delete();
86 }
87
88 if (this->IdList)
89 {
90 this->IdList->Delete();
91 }
92 }
93
94 //------------------------------------------------------------------------------
GetScalarModeAsString()95 const char* vtkClipClosedSurface::GetScalarModeAsString()
96 {
97 switch (this->ScalarMode)
98 {
99 case VTK_CCS_SCALAR_MODE_NONE:
100 return "None";
101 case VTK_CCS_SCALAR_MODE_COLORS:
102 return "Colors";
103 case VTK_CCS_SCALAR_MODE_LABELS:
104 return "Labels";
105 }
106 return "";
107 }
108
109 //------------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)110 void vtkClipClosedSurface::PrintSelf(ostream& os, vtkIndent indent)
111 {
112 this->Superclass::PrintSelf(os, indent);
113
114 os << indent << "ClippingPlanes: ";
115 if (this->ClippingPlanes)
116 {
117 os << this->ClippingPlanes << "\n";
118 }
119 else
120 {
121 os << "(none)\n";
122 }
123
124 os << indent << "Tolerance: " << this->Tolerance << "\n";
125
126 os << indent << "PassPointData: " << (this->PassPointData ? "On\n" : "Off\n");
127
128 os << indent << "GenerateOutline: " << (this->GenerateOutline ? "On\n" : "Off\n");
129
130 os << indent << "GenerateFaces: " << (this->GenerateFaces ? "On\n" : "Off\n");
131
132 os << indent << "ScalarMode: " << this->GetScalarModeAsString() << "\n";
133
134 os << indent << "BaseColor: " << this->BaseColor[0] << ", " << this->BaseColor[1] << ", "
135 << this->BaseColor[2] << "\n";
136
137 os << indent << "ClipColor: " << this->ClipColor[0] << ", " << this->ClipColor[1] << ", "
138 << this->ClipColor[2] << "\n";
139
140 os << indent << "ActivePlaneId: " << this->ActivePlaneId << "\n";
141
142 os << indent << "ActivePlaneColor: " << this->ActivePlaneColor[0] << ", "
143 << this->ActivePlaneColor[1] << ", " << this->ActivePlaneColor[2] << "\n";
144
145 os << indent
146 << "TriangulationErrorDisplay: " << (this->TriangulationErrorDisplay ? "On\n" : "Off\n");
147 }
148
149 //------------------------------------------------------------------------------
ComputePipelineMTime(vtkInformation * vtkNotUsed (request),vtkInformationVector ** vtkNotUsed (inputVector),vtkInformationVector * vtkNotUsed (outputVector),int vtkNotUsed (requestFromOutputPort),vtkMTimeType * mtime)150 int vtkClipClosedSurface::ComputePipelineMTime(vtkInformation* vtkNotUsed(request),
151 vtkInformationVector** vtkNotUsed(inputVector), vtkInformationVector* vtkNotUsed(outputVector),
152 int vtkNotUsed(requestFromOutputPort), vtkMTimeType* mtime)
153 {
154 vtkMTimeType mTime = this->GetMTime();
155
156 vtkPlaneCollection* planes = this->ClippingPlanes;
157 vtkPlane* plane = nullptr;
158
159 if (planes)
160 {
161 vtkMTimeType planesMTime = planes->GetMTime();
162 if (planesMTime > mTime)
163 {
164 mTime = planesMTime;
165 }
166
167 vtkCollectionSimpleIterator iter;
168 planes->InitTraversal(iter);
169 while ((plane = planes->GetNextPlane(iter)))
170 {
171 vtkMTimeType planeMTime = plane->GetMTime();
172 if (planeMTime > mTime)
173 {
174 mTime = planeMTime;
175 }
176 }
177 }
178
179 *mtime = mTime;
180
181 return 1;
182 }
183
184 //------------------------------------------------------------------------------
185 // A helper class to quickly locate an edge, given the endpoint ids.
186 // It uses an stl map rather than a table partitioning scheme, since
187 // we have no idea how many entries there will be when we start. So
188 // the performance is approximately log(n).
189
190 class vtkCCSEdgeLocatorNode
191 {
192 public:
vtkCCSEdgeLocatorNode()193 vtkCCSEdgeLocatorNode()
194 : ptId0(-1)
195 , ptId1(-1)
196 , edgeId(-1)
197 , next(nullptr)
198 {
199 }
200
FreeList()201 void FreeList()
202 {
203 vtkCCSEdgeLocatorNode* ptr = this->next;
204 while (ptr)
205 {
206 vtkCCSEdgeLocatorNode* tmp = ptr;
207 ptr = ptr->next;
208 tmp->next = nullptr;
209 delete tmp;
210 }
211 }
212
213 vtkIdType ptId0;
214 vtkIdType ptId1;
215 vtkIdType edgeId;
216 vtkCCSEdgeLocatorNode* next;
217 };
218
219 class vtkCCSEdgeLocator
220 {
221 private:
222 typedef std::map<vtkIdType, vtkCCSEdgeLocatorNode> MapType;
223 MapType EdgeMap;
224
225 public:
New()226 static vtkCCSEdgeLocator* New() { return new vtkCCSEdgeLocator; }
227
Delete()228 void Delete()
229 {
230 this->Initialize();
231 delete this;
232 }
233
234 // Description:
235 // Initialize the locator.
236 void Initialize();
237
238 // Description:
239 // If edge (i0, i1) is not in the list, then it will be added and
240 // a pointer for storing the new edgeId will be returned.
241 // If edge (i0, i1) is in the list, then edgeId will be set to the
242 // stored value and a null pointer will be returned.
243 vtkIdType* InsertUniqueEdge(vtkIdType i0, vtkIdType i1, vtkIdType& edgeId);
244 };
245
Initialize()246 void vtkCCSEdgeLocator::Initialize()
247 {
248 for (MapType::iterator i = this->EdgeMap.begin(); i != this->EdgeMap.end(); ++i)
249 {
250 i->second.FreeList();
251 }
252 this->EdgeMap.clear();
253 }
254
InsertUniqueEdge(vtkIdType i0,vtkIdType i1,vtkIdType & edgeId)255 vtkIdType* vtkCCSEdgeLocator::InsertUniqueEdge(vtkIdType i0, vtkIdType i1, vtkIdType& edgeId)
256 {
257 // Ensure consistent ordering of edge
258 if (i1 < i0)
259 {
260 vtkIdType tmp = i0;
261 i0 = i1;
262 i1 = tmp;
263 }
264
265 // Generate a integer key, try to make it unique
266 #ifdef VTK_USE_64BIT_IDS
267 vtkIdType key = ((i1 << 32) ^ i0);
268 #else
269 vtkIdType key = ((i1 << 16) ^ i0);
270 #endif
271
272 vtkCCSEdgeLocatorNode* node = &this->EdgeMap[key];
273
274 if (node->ptId1 < 0)
275 {
276 // Didn't find key, so add a new edge entry
277 node->ptId0 = i0;
278 node->ptId1 = i1;
279 return &node->edgeId;
280 }
281
282 // Search through the list for i0 and i1
283 if (node->ptId0 == i0 && node->ptId1 == i1)
284 {
285 edgeId = node->edgeId;
286 return nullptr;
287 }
288
289 int i = 1;
290 while (node->next != nullptr)
291 {
292 i++;
293 node = node->next;
294
295 if (node->ptId0 == i0 && node->ptId1 == i1)
296 {
297 edgeId = node->edgeId;
298 return nullptr;
299 }
300 }
301
302 // No entry for i1, so make one and return
303 node->next = new vtkCCSEdgeLocatorNode;
304 node = node->next;
305 node->ptId0 = i0;
306 node->ptId1 = i1;
307 node->edgeId = static_cast<vtkIdType>(this->EdgeMap.size() - 1);
308 return &node->edgeId;
309 }
310
311 //------------------------------------------------------------------------------
RequestData(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)312 int vtkClipClosedSurface::RequestData(vtkInformation* vtkNotUsed(request),
313 vtkInformationVector** inputVector, vtkInformationVector* outputVector)
314 {
315 // Get the info objects
316 vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
317 vtkInformation* outInfo = outputVector->GetInformationObject(0);
318
319 // Get the input and output
320 vtkPolyData* input = vtkPolyData::SafeDownCast(inInfo->Get(vtkDataObject::DATA_OBJECT()));
321 vtkPolyData* output = vtkPolyData::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));
322
323 // Create objects needed for temporary storage
324 if (this->IdList == nullptr)
325 {
326 this->IdList = vtkIdList::New();
327 }
328
329 // Get the input points
330 vtkPoints* inputPoints = input->GetPoints();
331 vtkIdType numPts = 0;
332 int inputPointsType = VTK_FLOAT;
333 if (inputPoints)
334 {
335 numPts = inputPoints->GetNumberOfPoints();
336 inputPointsType = inputPoints->GetDataType();
337 }
338
339 // Force points to double precision, copy the point attributes
340 vtkPoints* points = vtkPoints::New();
341 points->SetDataTypeToDouble();
342 points->SetNumberOfPoints(numPts);
343
344 vtkPointData* pointData = vtkPointData::New();
345 vtkPointData* inPointData = nullptr;
346
347 if (this->PassPointData)
348 {
349 inPointData = input->GetPointData();
350 pointData->InterpolateAllocate(inPointData, numPts, 0);
351 }
352
353 for (vtkIdType ptId = 0; ptId < numPts; ptId++)
354 {
355 double point[3];
356 inputPoints->GetPoint(ptId, point);
357 points->SetPoint(ptId, point);
358 // Point data is not copied from input
359 if (inPointData)
360 {
361 pointData->CopyData(inPointData, ptId, ptId);
362 }
363 }
364
365 // An edge locator to avoid point duplication while clipping
366 vtkCCSEdgeLocator* edgeLocator = vtkCCSEdgeLocator::New();
367
368 // A temporary polydata for the contour lines that are triangulated
369 vtkPolyData* tmpContourData = vtkPolyData::New();
370
371 // The cell scalars
372 vtkUnsignedCharArray* lineScalars = nullptr;
373 vtkUnsignedCharArray* polyScalars = nullptr;
374 vtkUnsignedCharArray* inputScalars = nullptr;
375
376 // For input scalars: the offsets to the various cell types
377 vtkIdType firstLineScalar = 0;
378 vtkIdType firstPolyScalar = 0;
379 vtkIdType firstStripScalar = 0;
380
381 // Make the colors to be used on the data.
382 int numberOfScalarComponents = 1;
383 unsigned char colors[3][3];
384
385 if (this->ScalarMode == VTK_CCS_SCALAR_MODE_COLORS)
386 {
387 numberOfScalarComponents = 3;
388 vtkClipClosedSurface::CreateColorValues(
389 this->BaseColor, this->ClipColor, this->ActivePlaneColor, colors);
390 }
391 else if (this->ScalarMode == VTK_CCS_SCALAR_MODE_LABELS)
392 {
393 colors[0][0] = 0;
394 colors[1][0] = 1;
395 colors[2][0] = 2;
396 }
397
398 // This is set if we have to work with scalars. The input scalars
399 // will be copied if they are unsigned char with 3 components, otherwise
400 // new scalars will be generated.
401 if (this->ScalarMode)
402 {
403 // Make the scalars
404 lineScalars = vtkUnsignedCharArray::New();
405 lineScalars->SetNumberOfComponents(numberOfScalarComponents);
406
407 vtkDataArray* tryInputScalars = input->GetCellData()->GetScalars();
408 // Get input scalars if they are RGB color scalars
409 if (tryInputScalars && tryInputScalars->IsA("vtkUnsignedCharArray") &&
410 numberOfScalarComponents == 3 && tryInputScalars->GetNumberOfComponents() == 3)
411 {
412 inputScalars = static_cast<vtkUnsignedCharArray*>(input->GetCellData()->GetScalars());
413
414 vtkIdType numVerts = 0;
415 vtkIdType numLines = 0;
416 vtkIdType numPolys = 0;
417 vtkCellArray* tmpCellArray = nullptr;
418 if ((tmpCellArray = input->GetVerts()))
419 {
420 numVerts = tmpCellArray->GetNumberOfCells();
421 }
422 if ((tmpCellArray = input->GetLines()))
423 {
424 numLines = tmpCellArray->GetNumberOfCells();
425 }
426 if ((tmpCellArray = input->GetPolys()))
427 {
428 numPolys = tmpCellArray->GetNumberOfCells();
429 }
430 firstLineScalar = numVerts;
431 firstPolyScalar = numVerts + numLines;
432 firstStripScalar = numVerts + numLines + numPolys;
433 }
434 }
435
436 // Break the input lines into segments, generate scalars for lines
437 vtkCellArray* lines = vtkCellArray::New();
438 if (input->GetLines() && input->GetLines()->GetNumberOfCells() > 0)
439 {
440 vtkClipClosedSurface::BreakPolylines(
441 input->GetLines(), lines, inputScalars, firstLineScalar, lineScalars, colors[0]);
442 }
443
444 // Copy the polygons, convert strips to triangles
445 vtkCellArray* polys = nullptr;
446 int polyMax = 3;
447 if ((input->GetPolys() && input->GetPolys()->GetNumberOfCells() > 0) ||
448 (input->GetStrips() && input->GetStrips()->GetNumberOfCells() > 0))
449 {
450 // If there are line scalars, then poly scalars are needed too
451 if (lineScalars)
452 {
453 polyScalars = vtkUnsignedCharArray::New();
454 polyScalars->SetNumberOfComponents(numberOfScalarComponents);
455 }
456
457 polys = vtkCellArray::New();
458 vtkClipClosedSurface::CopyPolygons(
459 input->GetPolys(), polys, inputScalars, firstPolyScalar, polyScalars, colors[0]);
460 vtkClipClosedSurface::BreakTriangleStrips(
461 input->GetStrips(), polys, inputScalars, firstStripScalar, polyScalars, colors[0]);
462
463 // Check if the input has polys and quads or just triangles
464 vtkIdType npts = 0;
465 const vtkIdType* pts = nullptr;
466 vtkCellArray* inPolys = input->GetPolys();
467 inPolys->InitTraversal();
468 while (inPolys->GetNextCell(npts, pts))
469 {
470 if (npts > polyMax)
471 {
472 polyMax = npts;
473 }
474 }
475 }
476
477 // Get the clipping planes
478 vtkPlaneCollection* planes = this->ClippingPlanes;
479
480 // Arrays for storing the clipped lines and polys.
481 vtkCellArray* newLines = vtkCellArray::New();
482 vtkCellArray* newPolys = nullptr;
483 if (polys)
484 {
485 newPolys = vtkCellArray::New();
486 }
487
488 // The point scalars, needed for clipping (not for the output!)
489 vtkDoubleArray* pointScalars = vtkDoubleArray::New();
490
491 // The line scalars, for coloring the outline
492 vtkCellData* inLineData = vtkCellData::New();
493 inLineData->CopyScalarsOn();
494 inLineData->SetScalars(lineScalars);
495 if (lineScalars)
496 {
497 lineScalars->Delete();
498 lineScalars = nullptr;
499 }
500
501 // The poly scalars, for coloring the faces
502 vtkCellData* inPolyData = vtkCellData::New();
503 inPolyData->CopyScalarsOn();
504 inPolyData->SetScalars(polyScalars);
505 if (polyScalars)
506 {
507 polyScalars->Delete();
508 polyScalars = nullptr;
509 }
510
511 // Also create output attribute data
512 vtkCellData* outLineData = vtkCellData::New();
513 outLineData->CopyScalarsOn();
514
515 vtkCellData* outPolyData = vtkCellData::New();
516 outPolyData->CopyScalarsOn();
517
518 // Go through the clipping planes and clip the input with each plane
519 vtkCollectionSimpleIterator iter;
520 int numPlanes = 0;
521 if (planes)
522 {
523 planes->InitTraversal(iter);
524 numPlanes = planes->GetNumberOfItems();
525 }
526
527 vtkPlane* plane = nullptr;
528 for (int planeId = 0; planes && (plane = planes->GetNextPlane(iter)); planeId++)
529 {
530 this->UpdateProgress((planeId + 1.0) / (numPlanes + 1.0));
531 if (this->GetAbortExecute())
532 {
533 break;
534 }
535
536 // Is this the last cut plane? If so, generate triangles.
537 int triangulate = 5;
538 if (planeId == numPlanes - 1)
539 {
540 triangulate = polyMax;
541 }
542
543 // Is this the active plane?
544 int active = (planeId == this->ActivePlaneId);
545
546 // Convert the plane into an easy-to-evaluate function
547 double pc[4];
548 plane->GetNormal(pc);
549 pc[3] = -vtkMath::Dot(pc, plane->GetOrigin());
550
551 // Create the clip scalars by evaluating the plane at each point
552 vtkIdType numPoints = points->GetNumberOfPoints();
553 pointScalars->SetNumberOfValues(numPoints);
554 for (vtkIdType pointId = 0; pointId < numPoints; pointId++)
555 {
556 double p[3];
557 points->GetPoint(pointId, p);
558 double val = p[0] * pc[0] + p[1] * pc[1] + p[2] * pc[2] + pc[3];
559 pointScalars->SetValue(pointId, val);
560 }
561
562 // Prepare the output scalars
563 outLineData->CopyAllocate(inLineData, 0, 0);
564 outPolyData->CopyAllocate(inPolyData, 0, 0);
565
566 // Reset the locator
567 edgeLocator->Initialize();
568
569 // Clip the lines
570 this->ClipLines(
571 points, pointScalars, pointData, edgeLocator, lines, newLines, inLineData, outLineData);
572
573 // Clip the polys
574 if (polys)
575 {
576 // Get the number of lines remaining after the clipping
577 vtkIdType numClipLines = newLines->GetNumberOfCells();
578
579 // Cut the polys to generate more lines
580 this->ClipAndContourPolys(points, pointScalars, pointData, edgeLocator, triangulate, polys,
581 newPolys, newLines, inPolyData, outPolyData, outLineData);
582
583 // Add scalars for the newly-created contour lines
584 vtkUnsignedCharArray* scalars =
585 vtkArrayDownCast<vtkUnsignedCharArray>(outLineData->GetScalars());
586
587 if (scalars)
588 {
589 // Set the color to the active color if plane is active
590 unsigned char* color = colors[1 + active];
591 unsigned char* activeColor = colors[2];
592
593 vtkIdType numLines = newLines->GetNumberOfCells();
594 for (vtkIdType lineId = numClipLines; lineId < numLines; lineId++)
595 {
596 unsigned char oldColor[3];
597 scalars->GetTypedTuple(lineId, oldColor);
598 if (numberOfScalarComponents != 3 || oldColor[0] != activeColor[0] ||
599 oldColor[1] != activeColor[1] || oldColor[2] != activeColor[2])
600 {
601 scalars->SetTypedTuple(lineId, color);
602 }
603 }
604 }
605
606 // Generate new polys from the cut lines
607 vtkIdType cellId = newPolys->GetNumberOfCells();
608 vtkIdType numClipAndContourLines = newLines->GetNumberOfCells();
609
610 // Create a polydata for the lines
611 tmpContourData->SetPoints(points);
612 tmpContourData->SetLines(newLines);
613 tmpContourData->BuildCells();
614
615 this->TriangulateContours(
616 tmpContourData, numClipLines, numClipAndContourLines - numClipLines, newPolys, pc);
617
618 // Add scalars for the newly-created polys
619 scalars = vtkArrayDownCast<vtkUnsignedCharArray>(outPolyData->GetScalars());
620
621 if (scalars)
622 {
623 unsigned char* color = colors[1 + active];
624
625 vtkIdType numCells = newPolys->GetNumberOfCells();
626 if (numCells > cellId)
627 {
628 // The insert allocates space up to numCells-1
629 scalars->InsertTypedTuple(numCells - 1, color);
630 for (; cellId < numCells; cellId++)
631 {
632 scalars->SetTypedTuple(cellId, color);
633 }
634 }
635 }
636
637 // Add scalars to any diagnostic lines that added by
638 // TriangulateContours(). In usual operation, no lines are added.
639 scalars = vtkArrayDownCast<vtkUnsignedCharArray>(outLineData->GetScalars());
640
641 if (scalars)
642 {
643 unsigned char color[3] = { 0, 255, 255 };
644
645 vtkIdType numCells = newLines->GetNumberOfCells();
646 if (numCells > numClipAndContourLines)
647 {
648 // The insert allocates space up to numCells-1
649 scalars->InsertTypedTuple(numCells - 1, color);
650 for (vtkIdType lineCellId = numClipAndContourLines; lineCellId < numCells; lineCellId++)
651 {
652 scalars->SetTypedTuple(lineCellId, color);
653 }
654 }
655 }
656 }
657
658 // Swap the lines, points, etcetera: old output becomes new input
659 vtkCellArray* tmp1 = lines;
660 lines = newLines;
661 newLines = tmp1;
662 newLines->Initialize();
663
664 if (polys)
665 {
666 vtkCellArray* tmp2 = polys;
667 polys = newPolys;
668 newPolys = tmp2;
669 newPolys->Initialize();
670 }
671
672 vtkCellData* tmp4 = inLineData;
673 inLineData = outLineData;
674 outLineData = tmp4;
675 outLineData->Initialize();
676
677 vtkCellData* tmp5 = inPolyData;
678 inPolyData = outPolyData;
679 outPolyData = tmp5;
680 outPolyData->Initialize();
681 }
682
683 // Delete the locator
684 edgeLocator->Delete();
685
686 // Delete the contour data container
687 tmpContourData->Delete();
688
689 // Delete the clip scalars
690 pointScalars->Delete();
691
692 // Get the line scalars
693 vtkUnsignedCharArray* scalars = vtkArrayDownCast<vtkUnsignedCharArray>(inLineData->GetScalars());
694
695 if (this->GenerateOutline)
696 {
697 output->SetLines(lines);
698 }
699 else if (scalars)
700 {
701 // If not adding lines to output, clear the line scalars
702 scalars->Initialize();
703 }
704
705 if (this->GenerateFaces)
706 {
707 output->SetPolys(polys);
708
709 if (polys && scalars)
710 {
711 vtkUnsignedCharArray* pScalars =
712 vtkArrayDownCast<vtkUnsignedCharArray>(inPolyData->GetScalars());
713
714 vtkIdType m = scalars->GetNumberOfTuples();
715 vtkIdType n = pScalars->GetNumberOfTuples();
716
717 if (n > 0)
718 {
719 unsigned char color[3];
720 color[0] = color[1] = color[2] = 0;
721
722 // This is just to expand the array
723 scalars->InsertTypedTuple(n + m - 1, color);
724
725 // Fill in the poly scalars
726 for (vtkIdType i = 0; i < n; i++)
727 {
728 pScalars->GetTypedTuple(i, color);
729 scalars->SetTypedTuple(i + m, color);
730 }
731 }
732 }
733 }
734
735 lines->Delete();
736
737 if (polys)
738 {
739 polys->Delete();
740 }
741
742 if (this->ScalarMode == VTK_CCS_SCALAR_MODE_COLORS)
743 {
744 scalars->SetName("Colors");
745 output->GetCellData()->SetScalars(scalars);
746 }
747 else if (this->ScalarMode == VTK_CCS_SCALAR_MODE_LABELS)
748 {
749 // Don't use VTK_UNSIGNED_CHAR or they will look like color scalars
750 vtkSignedCharArray* categories = vtkSignedCharArray::New();
751 categories->DeepCopy(scalars);
752 categories->SetName("Labels");
753 output->GetCellData()->SetScalars(categories);
754 categories->Delete();
755 }
756 else
757 {
758 output->GetCellData()->SetScalars(nullptr);
759 }
760
761 newLines->Delete();
762 if (newPolys)
763 {
764 newPolys->Delete();
765 }
766
767 inLineData->Delete();
768 outLineData->Delete();
769 inPolyData->Delete();
770 outPolyData->Delete();
771
772 // Finally, store the points in the output
773 vtkClipClosedSurface::SqueezeOutputPoints(output, points, pointData, inputPointsType);
774 output->Squeeze();
775
776 points->Delete();
777 pointData->Delete();
778
779 return 1;
780 }
781
782 //------------------------------------------------------------------------------
SqueezeOutputPoints(vtkPolyData * output,vtkPoints * points,vtkPointData * pointData,int outputPointDataType)783 void vtkClipClosedSurface::SqueezeOutputPoints(
784 vtkPolyData* output, vtkPoints* points, vtkPointData* pointData, int outputPointDataType)
785 {
786 // Create a list of points used by cells
787 vtkIdType n = points->GetNumberOfPoints();
788 vtkIdType numNewPoints = 0;
789
790 // The point data
791 vtkPointData* outPointData = output->GetPointData();
792
793 // A mapping from old pointIds to new pointIds
794 vtkIdType* pointMap = new vtkIdType[n];
795 for (vtkIdType i = 0; i < n; i++)
796 {
797 pointMap[i] = -1;
798 }
799
800 vtkIdType npts;
801 const vtkIdType* pts;
802 vtkCellArray* cellArrays[4];
803 cellArrays[0] = output->GetVerts();
804 cellArrays[1] = output->GetLines();
805 cellArrays[2] = output->GetPolys();
806 cellArrays[3] = output->GetStrips();
807 int arrayId;
808
809 // Find all the newPoints that are used by cells
810 for (arrayId = 0; arrayId < 4; arrayId++)
811 {
812 vtkCellArray* cellArray = cellArrays[arrayId];
813 if (cellArray)
814 {
815 cellArray->InitTraversal();
816 while (cellArray->GetNextCell(npts, pts))
817 {
818 for (vtkIdType ii = 0; ii < npts; ii++)
819 {
820 vtkIdType pointId = pts[ii];
821 if (pointMap[pointId] < 0)
822 {
823 pointMap[pointId] = numNewPoints++;
824 }
825 }
826 }
827 }
828 }
829
830 // Create exactly the number of points that are required
831 vtkPoints* newPoints = vtkPoints::New();
832 newPoints->SetDataType(outputPointDataType);
833 newPoints->SetNumberOfPoints(numNewPoints);
834 outPointData->CopyAllocate(pointData, numNewPoints, 0);
835
836 for (vtkIdType pointId = 0; pointId < n; pointId++)
837 {
838 vtkIdType newPointId = pointMap[pointId];
839 if (newPointId >= 0)
840 {
841 double p[3];
842 points->GetPoint(pointId, p);
843 newPoints->SetPoint(newPointId, p);
844 outPointData->CopyData(pointData, pointId, newPointId);
845 }
846 }
847
848 // Change the cell pointIds to reflect the new point array
849 vtkNew<vtkIdList> repCell;
850 for (arrayId = 0; arrayId < 4; arrayId++)
851 {
852 vtkCellArray* cellArray = cellArrays[arrayId];
853 if (cellArray)
854 {
855 auto cellIter = vtkSmartPointer<vtkCellArrayIterator>::Take(cellArray->NewIterator());
856 for (cellIter->GoToFirstCell(); !cellIter->IsDoneWithTraversal(); cellIter->GoToNextCell())
857 {
858 cellIter->GetCurrentCell(repCell);
859 for (vtkIdType ii = 0; ii < repCell->GetNumberOfIds(); ii++)
860 {
861 const vtkIdType pointId = repCell->GetId(ii);
862 repCell->SetId(ii, pointMap[pointId]);
863 }
864 cellIter->ReplaceCurrentCell(repCell);
865 }
866 }
867 }
868
869 output->SetPoints(newPoints);
870 newPoints->Delete();
871
872 delete[] pointMap;
873 }
874
875 //------------------------------------------------------------------------------
CreateColorValues(const double color1[3],const double color2[3],const double color3[3],unsigned char colors[3][3])876 void vtkClipClosedSurface::CreateColorValues(const double color1[3], const double color2[3],
877 const double color3[3], unsigned char colors[3][3])
878 {
879 // Convert colors from "double" to "unsigned char"
880
881 const double* dcolors[3];
882 dcolors[0] = color1;
883 dcolors[1] = color2;
884 dcolors[2] = color3;
885
886 for (int i = 0; i < 3; i++)
887 {
888 for (int j = 0; j < 3; j++)
889 {
890 double val = dcolors[i][j];
891 if (val < 0)
892 {
893 val = 0;
894 }
895 if (val > 1)
896 {
897 val = 1;
898 }
899 colors[i][j] = static_cast<unsigned char>(val * 255);
900 }
901 }
902 }
903
904 //------------------------------------------------------------------------------
905 // Point interpolation for clipping and contouring, given the scalar
906 // values (v0, v1) for the two endpoints (p0, p1). The use of this
907 // function guarantees perfect consistency in the results.
InterpolateEdge(vtkPoints * points,vtkPointData * pointData,vtkCCSEdgeLocator * locator,double tol,vtkIdType i0,vtkIdType i1,double v0,double v1,vtkIdType & i)908 int vtkClipClosedSurface::InterpolateEdge(vtkPoints* points, vtkPointData* pointData,
909 vtkCCSEdgeLocator* locator, double tol, vtkIdType i0, vtkIdType i1, double v0, double v1,
910 vtkIdType& i)
911 {
912 // This swap guarantees that exactly the same point is computed
913 // for both line directions, as long as the endpoints are the same.
914 if (v1 > 0)
915 {
916 vtkIdType tmpi = i0;
917 i0 = i1;
918 i1 = tmpi;
919
920 double tmp = v0;
921 v0 = v1;
922 v1 = tmp;
923 }
924
925 // After the above swap, i0 will be kept, and i1 will be clipped
926
927 // Check to see if this point has already been computed
928 vtkIdType* iptr = locator->InsertUniqueEdge(i0, i1, i);
929 if (iptr == nullptr)
930 {
931 return 0;
932 }
933
934 // Get the edge and interpolate the new point
935 double p0[3], p1[3], p[3];
936 points->GetPoint(i0, p0);
937 points->GetPoint(i1, p1);
938
939 double f = v0 / (v0 - v1);
940 double s = 1.0 - f;
941 double t = 1.0 - s;
942
943 p[0] = s * p0[0] + t * p1[0];
944 p[1] = s * p0[1] + t * p1[1];
945 p[2] = s * p0[2] + t * p1[2];
946
947 double tol2 = tol * tol;
948
949 // Make sure that new point is far enough from kept point
950 if (vtkMath::Distance2BetweenPoints(p, p0) < tol2)
951 {
952 i = i0;
953 *iptr = i0;
954 return 0;
955 }
956
957 if (vtkMath::Distance2BetweenPoints(p, p1) < tol2)
958 {
959 i = i1;
960 *iptr = i1;
961 return 0;
962 }
963
964 i = points->InsertNextPoint(p);
965 pointData->InterpolateEdge(pointData, i, i0, i1, t);
966
967 // Store the new index in the locator
968 *iptr = i;
969
970 return 1;
971 }
972
973 //------------------------------------------------------------------------------
ClipLines(vtkPoints * points,vtkDoubleArray * pointScalars,vtkPointData * pointData,vtkCCSEdgeLocator * edgeLocator,vtkCellArray * inputCells,vtkCellArray * outputLines,vtkCellData * inCellData,vtkCellData * outLineData)974 void vtkClipClosedSurface::ClipLines(vtkPoints* points, vtkDoubleArray* pointScalars,
975 vtkPointData* pointData, vtkCCSEdgeLocator* edgeLocator, vtkCellArray* inputCells,
976 vtkCellArray* outputLines, vtkCellData* inCellData, vtkCellData* outLineData)
977 {
978 vtkIdType numCells = inputCells->GetNumberOfCells();
979
980 inputCells->InitTraversal();
981 for (vtkIdType cellId = 0; cellId < numCells; cellId++)
982 {
983 vtkIdType numPts = 0;
984 const vtkIdType* pts = nullptr;
985 inputCells->GetNextCell(numPts, pts);
986
987 vtkIdType i1 = pts[0];
988 double v1 = pointScalars->GetValue(i1);
989 int c1 = (v1 > 0);
990
991 for (vtkIdType i = 1; i < numPts; i++)
992 {
993 vtkIdType i0 = i1;
994 double v0 = v1;
995 int c0 = c1;
996
997 i1 = pts[i];
998 v1 = pointScalars->GetValue(i1);
999 c1 = (v1 > 0);
1000
1001 // If at least one point wasn't clipped
1002 if ((c0 | c1))
1003 {
1004 vtkIdType linePts[2];
1005 linePts[0] = i0;
1006 linePts[1] = i1;
1007
1008 // If only one end was clipped, interpolate new point
1009 if ((c0 ^ c1))
1010 {
1011 vtkClipClosedSurface::InterpolateEdge(
1012 points, pointData, edgeLocator, this->Tolerance, i0, i1, v0, v1, linePts[c0]);
1013 }
1014
1015 // If endpoints are different, insert the line segment
1016 if (linePts[0] != linePts[1])
1017 {
1018 vtkIdType newCellId = outputLines->InsertNextCell(2, linePts);
1019 outLineData->CopyData(inCellData, cellId, newCellId);
1020 }
1021 }
1022 }
1023 }
1024 }
1025
1026 //------------------------------------------------------------------------------
ClipAndContourPolys(vtkPoints * points,vtkDoubleArray * pointScalars,vtkPointData * pointData,vtkCCSEdgeLocator * edgeLocator,int triangulate,vtkCellArray * inputCells,vtkCellArray * outputPolys,vtkCellArray * outputLines,vtkCellData * inCellData,vtkCellData * outPolyData,vtkCellData * outLineData)1027 void vtkClipClosedSurface::ClipAndContourPolys(vtkPoints* points, vtkDoubleArray* pointScalars,
1028 vtkPointData* pointData, vtkCCSEdgeLocator* edgeLocator, int triangulate,
1029 vtkCellArray* inputCells, vtkCellArray* outputPolys, vtkCellArray* outputLines,
1030 vtkCellData* inCellData, vtkCellData* outPolyData, vtkCellData* outLineData)
1031 {
1032 vtkIdList* idList = this->IdList;
1033
1034 // How many sides for output polygons?
1035 int polyMax = VTK_INT_MAX;
1036 if (triangulate)
1037 {
1038 if (triangulate < 4)
1039 { // triangles only
1040 polyMax = 3;
1041 }
1042 else if (triangulate == 4)
1043 { // allow triangles and quads
1044 polyMax = 4;
1045 }
1046 }
1047
1048 int triangulationFailure = 0;
1049
1050 // Go through all cells and clip them.
1051 vtkIdType numCells = inputCells->GetNumberOfCells();
1052
1053 inputCells->InitTraversal();
1054 for (vtkIdType cellId = 0; cellId < numCells; cellId++)
1055 {
1056 vtkIdType numPts = 0;
1057 const vtkIdType* pts = nullptr;
1058 inputCells->GetNextCell(numPts, pts);
1059 idList->Reset();
1060
1061 vtkIdType i1 = pts[numPts - 1];
1062 double v1 = pointScalars->GetValue(i1);
1063 int c1 = (v1 > 0);
1064
1065 // The ids for the current edge: init j0 to -1 if i1 will be clipped
1066 vtkIdType j0 = (c1 ? i1 : -1);
1067 vtkIdType j1 = 0;
1068
1069 // To store the ids of the contour line
1070 vtkIdType linePts[2];
1071 linePts[0] = 0;
1072 linePts[1] = 0;
1073
1074 for (vtkIdType i = 0; i < numPts; i++)
1075 {
1076 // Save previous point info
1077 vtkIdType i0 = i1;
1078 double v0 = v1;
1079 int c0 = c1;
1080
1081 // Generate new point info
1082 i1 = pts[i];
1083 v1 = pointScalars->GetValue(i1);
1084 c1 = (v1 > 0);
1085
1086 // If at least one edge end point wasn't clipped
1087 if ((c0 | c1))
1088 {
1089 // If only one end was clipped, interpolate new point
1090 if ((c0 ^ c1))
1091 {
1092 vtkClipClosedSurface::InterpolateEdge(
1093 points, pointData, edgeLocator, this->Tolerance, i0, i1, v0, v1, j1);
1094
1095 if (j1 != j0)
1096 {
1097 idList->InsertNextId(j1);
1098 j0 = j1;
1099 }
1100
1101 // Save as one end of the contour line
1102 linePts[c0] = j1;
1103 }
1104
1105 if (c1)
1106 {
1107 j1 = i1;
1108
1109 if (j1 != j0)
1110 {
1111 idList->InsertNextId(j1);
1112 j0 = j1;
1113 }
1114 }
1115 }
1116 }
1117
1118 // Insert the clipped poly
1119 vtkIdType numPoints = idList->GetNumberOfIds();
1120
1121 if (numPoints > polyMax)
1122 {
1123 vtkIdType newCellId = outputPolys->GetNumberOfCells();
1124
1125 // Triangulate the poly and insert triangles into output.
1126 if (!this->TriangulatePolygon(idList, points, outputPolys))
1127 {
1128 triangulationFailure = 1;
1129 }
1130
1131 // Copy the attribute data to the triangle cells
1132 vtkIdType nCells = outputPolys->GetNumberOfCells();
1133 for (; newCellId < nCells; newCellId++)
1134 {
1135 outPolyData->CopyData(inCellData, cellId, newCellId);
1136 }
1137 }
1138 else if (numPoints > 2)
1139 {
1140 // Insert the polygon without triangulating it
1141 vtkIdType newCellId = outputPolys->InsertNextCell(idList);
1142 outPolyData->CopyData(inCellData, cellId, newCellId);
1143 }
1144
1145 // Insert the contour line if one was created
1146 if (linePts[0] != linePts[1])
1147 {
1148 vtkIdType newCellId = outputLines->InsertNextCell(2, linePts);
1149 outLineData->CopyData(inCellData, cellId, newCellId);
1150 }
1151 }
1152
1153 if (triangulationFailure && this->TriangulationErrorDisplay)
1154 {
1155 vtkErrorMacro("Triangulation failed, output may not be watertight");
1156 }
1157
1158 // Free up the idList memory
1159 idList->Initialize();
1160 }
1161
1162 //------------------------------------------------------------------------------
BreakPolylines(vtkCellArray * inputLines,vtkCellArray * lines,vtkUnsignedCharArray * inputScalars,vtkIdType firstLineScalar,vtkUnsignedCharArray * scalars,const unsigned char color[3])1163 void vtkClipClosedSurface::BreakPolylines(vtkCellArray* inputLines, vtkCellArray* lines,
1164 vtkUnsignedCharArray* inputScalars, vtkIdType firstLineScalar, vtkUnsignedCharArray* scalars,
1165 const unsigned char color[3])
1166 {
1167 // The color for the lines
1168 unsigned char cellColor[3];
1169 cellColor[0] = color[0];
1170 cellColor[1] = color[1];
1171 cellColor[2] = color[2];
1172
1173 // Break the input lines into segments
1174 inputLines->InitTraversal();
1175 vtkIdType cellId = 0;
1176 vtkIdType npts;
1177 const vtkIdType* pts;
1178 while (inputLines->GetNextCell(npts, pts))
1179 {
1180 if (inputScalars)
1181 {
1182 inputScalars->GetTypedTuple(firstLineScalar + cellId++, cellColor);
1183 }
1184
1185 for (vtkIdType i = 1; i < npts; i++)
1186 {
1187 lines->InsertNextCell(2);
1188 lines->InsertCellPoint(pts[i - 1]);
1189 lines->InsertCellPoint(pts[i]);
1190
1191 if (scalars)
1192 {
1193 scalars->InsertNextTypedTuple(cellColor);
1194 }
1195 }
1196 }
1197 }
1198
1199 //------------------------------------------------------------------------------
CopyPolygons(vtkCellArray * inputPolys,vtkCellArray * polys,vtkUnsignedCharArray * inputScalars,vtkIdType firstPolyScalar,vtkUnsignedCharArray * polyScalars,const unsigned char color[3])1200 void vtkClipClosedSurface::CopyPolygons(vtkCellArray* inputPolys, vtkCellArray* polys,
1201 vtkUnsignedCharArray* inputScalars, vtkIdType firstPolyScalar, vtkUnsignedCharArray* polyScalars,
1202 const unsigned char color[3])
1203 {
1204 if (!inputPolys)
1205 {
1206 return;
1207 }
1208
1209 polys->DeepCopy(inputPolys);
1210
1211 if (polyScalars)
1212 {
1213 unsigned char scalarValue[3];
1214 scalarValue[0] = color[0];
1215 scalarValue[1] = color[1];
1216 scalarValue[2] = color[2];
1217
1218 vtkIdType n = polys->GetNumberOfCells();
1219 polyScalars->SetNumberOfTuples(n);
1220
1221 if (inputScalars)
1222 {
1223 for (vtkIdType i = 0; i < n; i++)
1224 {
1225 inputScalars->GetTypedTuple(i + firstPolyScalar, scalarValue);
1226 polyScalars->SetTypedTuple(i, scalarValue);
1227 }
1228 }
1229 else
1230 {
1231 for (vtkIdType i = 0; i < n; i++)
1232 {
1233 polyScalars->SetTypedTuple(i, scalarValue);
1234 }
1235 }
1236 }
1237 }
1238
1239 //------------------------------------------------------------------------------
BreakTriangleStrips(vtkCellArray * inputStrips,vtkCellArray * polys,vtkUnsignedCharArray * inputScalars,vtkIdType firstStripScalar,vtkUnsignedCharArray * polyScalars,const unsigned char color[3])1240 void vtkClipClosedSurface::BreakTriangleStrips(vtkCellArray* inputStrips, vtkCellArray* polys,
1241 vtkUnsignedCharArray* inputScalars, vtkIdType firstStripScalar, vtkUnsignedCharArray* polyScalars,
1242 const unsigned char color[3])
1243 {
1244 if (!inputStrips)
1245 {
1246 return;
1247 }
1248
1249 vtkIdType npts = 0;
1250 const vtkIdType* pts = nullptr;
1251
1252 inputStrips->InitTraversal();
1253
1254 for (vtkIdType cellId = firstStripScalar; inputStrips->GetNextCell(npts, pts); cellId++)
1255 {
1256 vtkTriangleStrip::DecomposeStrip(npts, pts, polys);
1257
1258 if (polyScalars)
1259 {
1260 unsigned char scalarValue[3];
1261 scalarValue[0] = color[0];
1262 scalarValue[1] = color[1];
1263 scalarValue[2] = color[2];
1264
1265 if (inputScalars)
1266 {
1267 // If there are input scalars, use them instead of "color"
1268 inputScalars->GetTypedTuple(cellId, scalarValue);
1269 }
1270
1271 vtkIdType n = npts - 3;
1272 vtkIdType m = polyScalars->GetNumberOfTuples();
1273 if (n >= 0)
1274 {
1275 // First insert is just to allocate space
1276 polyScalars->InsertTypedTuple(m + n, scalarValue);
1277
1278 for (vtkIdType i = 0; i < n; i++)
1279 {
1280 polyScalars->SetTypedTuple(m + i, scalarValue);
1281 }
1282 }
1283 }
1284 }
1285 }
1286
1287 //------------------------------------------------------------------------------
TriangulateContours(vtkPolyData * data,vtkIdType firstLine,vtkIdType numLines,vtkCellArray * outputPolys,const double normal[3])1288 void vtkClipClosedSurface::TriangulateContours(vtkPolyData* data, vtkIdType firstLine,
1289 vtkIdType numLines, vtkCellArray* outputPolys, const double normal[3])
1290 {
1291 // If no cut lines were generated, there's nothing to do
1292 if (numLines <= 0)
1293 {
1294 return;
1295 }
1296
1297 double nnormal[3] = { -normal[0], -normal[1], -normal[2] };
1298 int rval =
1299 vtkContourTriangulator::TriangulateContours(data, firstLine, numLines, outputPolys, nnormal);
1300
1301 if (rval == 0 && this->TriangulationErrorDisplay)
1302 {
1303 vtkErrorMacro("Triangulation failed, data may not be watertight.");
1304 }
1305 }
1306
1307 // ---------------------------------------------------
TriangulatePolygon(vtkIdList * polygon,vtkPoints * points,vtkCellArray * triangles)1308 int vtkClipClosedSurface::TriangulatePolygon(
1309 vtkIdList* polygon, vtkPoints* points, vtkCellArray* triangles)
1310 {
1311 return vtkContourTriangulator::TriangulatePolygon(polygon, points, triangles);
1312 }
1313