1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkBoxClipDataSet.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 /*----------------------------------------------------------------------------
16 Copyright (c) Sandia Corporation
17 See Copyright.txt or http://www.paraview.org/HTML/Copyright.html for details.
18 ----------------------------------------------------------------------------*/
19 #include "vtkBoxClipDataSet.h"
20
21 #include "vtkCellArray.h"
22 #include "vtkCellData.h"
23 #include "vtkExecutive.h"
24 #include "vtkFloatArray.h"
25 #include "vtkGenericCell.h"
26 #include "vtkImageData.h"
27 #include "vtkInformation.h"
28 #include "vtkInformationVector.h"
29 #include "vtkIntArray.h"
30 #include "vtkMergePoints.h"
31 #include "vtkObjectFactory.h"
32 #include "vtkPointData.h"
33 #include "vtkUnsignedCharArray.h"
34 #include "vtkUnstructuredGrid.h"
35 #include "vtkIdList.h"
36 #include "vtkIncrementalPointLocator.h"
37
38 #include <math.h>
39 #include <vector>
40
41 vtkStandardNewMacro(vtkBoxClipDataSet);
vtkCxxSetObjectMacro(vtkBoxClipDataSet,Locator,vtkIncrementalPointLocator)42 vtkCxxSetObjectMacro(vtkBoxClipDataSet, Locator, vtkIncrementalPointLocator)
43 //----------------------------------------------------------------------------
44 vtkBoxClipDataSet::vtkBoxClipDataSet()
45 {
46 this->Locator = NULL;
47 this->GenerateClipScalars = 0;
48
49 this->GenerateClippedOutput = 0;
50 //this->MergeTolerance = 0.01;
51
52 this->SetNumberOfOutputPorts(2);
53
54 this->Orientation = 1;
55
56 this->PlaneNormal[0][0] = -1.0;
57 this->PlaneNormal[0][1] = 0.0;
58 this->PlaneNormal[0][2] = 0.0;
59
60 this->PlaneNormal[1][0] = 1.0;
61 this->PlaneNormal[1][1] = 0.0;
62 this->PlaneNormal[1][2] = 0.0;
63
64 this->PlaneNormal[2][0] = 0.0;
65 this->PlaneNormal[2][1] = -1.0;
66 this->PlaneNormal[2][2] = 0.0;
67
68 this->PlaneNormal[3][0] = 0.0;
69 this->PlaneNormal[3][1] = 1.0;
70 this->PlaneNormal[3][2] = 0.0;
71
72 this->PlaneNormal[4][0] = 0.0;
73 this->PlaneNormal[4][1] = 0.0;
74 this->PlaneNormal[4][2] = -1.0;
75
76 this->PlaneNormal[5][0] = 0.0;
77 this->PlaneNormal[5][1] = 0.0;
78 this->PlaneNormal[5][2] = 1.0;
79
80 this->PlanePoint[0][0] = 0.0;
81 this->PlanePoint[0][1] = 0.0;
82 this->PlanePoint[0][2] = 0.0;
83
84 this->PlanePoint[1][0] = 1.0;
85 this->PlanePoint[1][1] = 0.0;
86 this->PlanePoint[1][2] = 0.0;
87
88 this->PlanePoint[2][0] = 0.0;
89 this->PlanePoint[2][1] = 0.0;
90 this->PlanePoint[2][2] = 0.0;
91
92 this->PlanePoint[3][0] = 0.0;
93 this->PlanePoint[3][1] = 1.0;
94 this->PlanePoint[3][2] = 0.0;
95
96 this->PlanePoint[4][0] = 0.0;
97 this->PlanePoint[4][1] = 0.0;
98 this->PlanePoint[4][2] = 0.0;
99
100 this->PlanePoint[5][0] = 0.0;
101 this->PlanePoint[5][1] = 0.0;
102 this->PlanePoint[5][2] = 1.0;
103
104 int i=0;
105 while(i<3)
106 {
107 this->BoundBoxClip[i][0]=0.0;
108 this->BoundBoxClip[i][1]=1.0;
109 ++i;
110 }
111
112 // by default process active point scalars
113 this->SetInputArrayToProcess(0,0,0,vtkDataObject::FIELD_ASSOCIATION_POINTS,
114 vtkDataSetAttributes::SCALARS);
115 }
116
117 //----------------------------------------------------------------------------
~vtkBoxClipDataSet()118 vtkBoxClipDataSet::~vtkBoxClipDataSet()
119 {
120 this->SetLocator(NULL);
121 }
122
123 //----------------------------------------------------------------------------
124 // Do not say we have two outputs unless we are generating the clipped output.
GetNumberOfOutputs()125 int vtkBoxClipDataSet::GetNumberOfOutputs()
126 {
127 if (this->GenerateClippedOutput)
128 {
129 return 2;
130 }
131 return 1;
132 }
133
134 //----------------------------------------------------------------------------
135 // Overload standard modified time function. If Clip functions is modified,
136 // then this object is modified as well.
GetMTime()137 unsigned long vtkBoxClipDataSet::GetMTime()
138 {
139 unsigned long mTime = this->Superclass::GetMTime();
140 unsigned long time;
141
142 if ( this->Locator != NULL )
143 {
144 time = this->Locator->GetMTime();
145 mTime = ( time > mTime ? time : mTime );
146 }
147
148 return mTime;
149 }
150
151 //----------------------------------------------------------------------------
GetClippedOutput()152 vtkUnstructuredGrid *vtkBoxClipDataSet::GetClippedOutput()
153 {
154 if (this->GetNumberOfOutputPorts() < 2)
155 {
156 return NULL;
157 }
158
159 return vtkUnstructuredGrid::SafeDownCast(
160 this->GetExecutive()->GetOutputData(1));
161 }
162
163 //----------------------------------------------------------------------------
164 //
165 // Clip by box
166 //
RequestData(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)167 int vtkBoxClipDataSet::RequestData(vtkInformation *vtkNotUsed(request),
168 vtkInformationVector **inputVector,
169 vtkInformationVector *outputVector)
170 {
171 // get the info objects
172 vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
173 vtkInformation *outInfo = outputVector->GetInformationObject(0);
174
175 // get the input and output
176 vtkDataSet *input = vtkDataSet::SafeDownCast(
177 inInfo->Get(vtkDataObject::DATA_OBJECT()));
178 vtkUnstructuredGrid *output = vtkUnstructuredGrid::SafeDownCast(
179 outInfo->Get(vtkDataObject::DATA_OBJECT()));
180
181 vtkUnstructuredGrid *clippedOutput = this->GetClippedOutput();
182
183 vtkIdType i;
184 vtkIdType npts;
185 vtkIdType *pts;
186 vtkIdType estimatedSize;
187 vtkIdType newCellId;
188 vtkIdType numPts = input->GetNumberOfPoints();
189 vtkIdType numCells = input->GetNumberOfCells();
190 vtkPointData *inPD=input->GetPointData(), *outPD = output->GetPointData();
191 vtkCellData *inCD=input->GetCellData();
192 vtkCellData *outCD[2];
193 vtkPoints *newPoints;
194 vtkPoints *cellPts;
195 vtkIdTypeArray *locs[2];
196 vtkDebugMacro( << "Clip by Box\n" );
197 vtkUnsignedCharArray *types[2];
198
199 int j;
200 int cellType = 0;
201 int numOutputs = 1;
202 int inputObjectType = input->GetDataObjectType();
203
204 // if we have volumes
205 if (inputObjectType == VTK_STRUCTURED_POINTS ||
206 inputObjectType == VTK_IMAGE_DATA)
207 {
208 int dimension;
209 int *dims = vtkImageData::SafeDownCast(input)->GetDimensions();
210 for (dimension=3, i=0; i<3; i++)
211 {
212 if ( dims[i] <= 1 )
213 {
214 dimension--;
215 }
216 }
217 }
218
219 // Initialize self; create output objects
220 //
221 if ( numPts < 1 )
222 {
223 vtkDebugMacro(<<"No data to clip");
224 return 1;
225 }
226
227 // allocate the output and associated helper classes
228 estimatedSize = numCells;
229 estimatedSize = estimatedSize / 1024 * 1024; //multiple of 1024
230 if (estimatedSize < 1024)
231 {
232 estimatedSize = 1024;
233 }
234 vtkCellArray *conn[2];
235 conn[0] = vtkCellArray::New();
236 conn[0]->Allocate(estimatedSize,estimatedSize/2);
237 conn[0]->InitTraversal();
238 types[0] = vtkUnsignedCharArray::New();
239 types[0]->Allocate(estimatedSize,estimatedSize/2);
240 locs[0] = vtkIdTypeArray::New();
241 locs[0]->Allocate(estimatedSize,estimatedSize/2);
242
243 if ( this->GenerateClippedOutput )
244 {
245 numOutputs = 2;
246 conn[1] = vtkCellArray::New();
247 conn[1]->Allocate(estimatedSize,estimatedSize/2);
248 conn[1]->InitTraversal();
249 types[1] = vtkUnsignedCharArray::New();
250 types[1]->Allocate(estimatedSize,estimatedSize/2);
251 locs[1] = vtkIdTypeArray::New();
252 locs[1]->Allocate(estimatedSize,estimatedSize/2);
253 }
254
255 newPoints = vtkPoints::New();
256 newPoints->Allocate(numPts,numPts/2);
257
258 // locator used to merge potentially duplicate points
259 if ( this->Locator == NULL )
260 {
261 this->CreateDefaultLocator();
262 }
263 this->Locator->InitPointInsertion (newPoints, input->GetBounds());
264
265
266 vtkDataArray *scalars = this->GetInputArrayToProcess(0,inputVector);
267 if ( !this->GenerateClipScalars && !scalars)
268 {
269 outPD->CopyScalarsOff();
270 }
271 else
272 {
273 outPD->CopyScalarsOn();
274 }
275 outPD->InterpolateAllocate(inPD,estimatedSize,estimatedSize/2);
276 outCD[0] = output->GetCellData();
277 outCD[0]->CopyAllocate(inCD,estimatedSize,estimatedSize/2);
278 if ( this->GenerateClippedOutput )
279 {
280 outCD[1] = clippedOutput->GetCellData();
281 outCD[1]->CopyAllocate(inCD,estimatedSize,estimatedSize/2);
282 }
283
284 //Process all cells and clip each in turn
285
286 vtkIdType updateTime = numCells/20 + 1; // update roughly every 5%
287 vtkGenericCell *cell = vtkGenericCell::New();
288 vtkIdType cellId;
289
290 int abort = 0;
291 int num[2];
292 int numNew[2];
293
294 num[0] = num[1] = 0;
295 numNew[0] = numNew[1] = 0;
296
297 unsigned int orientation = this->GetOrientation(); //Test if there is a transformation
298
299 //clock_t init_tmp = clock();
300 for (cellId=0; cellId < numCells && !abort; cellId++)
301 {
302 if ( !(cellId % updateTime) )
303 {
304 this->UpdateProgress(static_cast<float>(cellId) / numCells);
305 abort = this->GetAbortExecute();
306 }
307
308 input->GetCell(cellId,cell);
309 cellPts = cell->GetPoints();
310 npts = cellPts->GetNumberOfPoints();
311
312 if (this->GenerateClippedOutput)
313 {
314 if((cell->GetCellDimension())==3 )
315 {
316 if (orientation)
317 {
318 this->ClipHexahedronInOut(newPoints,cell,this->Locator, conn,
319 inPD, outPD, inCD, cellId, outCD);
320 }
321 else
322 {
323 this->ClipBoxInOut(newPoints,cell, this->Locator, conn,
324 inPD, outPD, inCD, cellId, outCD);
325 }
326
327 numNew[0] = conn[0]->GetNumberOfCells() - num[0];
328 numNew[1] = conn[1]->GetNumberOfCells() - num[1];
329 num[0] = conn[0]->GetNumberOfCells();
330 num[1] = conn[1]->GetNumberOfCells();
331 }
332 else if((cell->GetCellDimension())==2 )
333 {
334 if (orientation)
335 {
336 this->ClipHexahedronInOut2D(newPoints,cell,this->Locator, conn,
337 inPD, outPD, inCD, cellId, outCD);
338 }
339 else
340 {
341 this->ClipBoxInOut2D(newPoints,cell, this->Locator, conn,
342 inPD, outPD, inCD, cellId, outCD);
343 }
344 numNew[0] = conn[0]->GetNumberOfCells() - num[0];
345 numNew[1] = conn[1]->GetNumberOfCells() - num[1];
346 num[0] = conn[0]->GetNumberOfCells();
347 num[1] = conn[1]->GetNumberOfCells();
348 }
349 else if (cell->GetCellDimension() == 1)
350 {
351 if (orientation)
352 {
353 this->ClipHexahedronInOut1D(newPoints,cell, this->Locator, conn,
354 inPD, outPD, inCD, cellId, outCD);
355 }
356 else
357 {
358 this->ClipBoxInOut1D(newPoints,cell, this->Locator, conn,
359 inPD, outPD, inCD, cellId, outCD);
360 }
361 numNew[0] = conn[0]->GetNumberOfCells() - num[0];
362 numNew[1] = conn[1]->GetNumberOfCells() - num[1];
363 num[0] = conn[0]->GetNumberOfCells();
364 num[1] = conn[1]->GetNumberOfCells();
365 }
366 else if (cell->GetCellDimension() == 0)
367 {
368 if (orientation)
369 {
370 this->ClipHexahedronInOut0D(cell, this->Locator, conn,
371 inPD, outPD, inCD, cellId, outCD);
372 }
373 else
374 {
375 this->ClipBoxInOut0D(cell, this->Locator, conn,
376 inPD, outPD, inCD, cellId, outCD);
377 }
378 numNew[0] = conn[0]->GetNumberOfCells() - num[0];
379 numNew[1] = conn[1]->GetNumberOfCells() - num[1];
380 num[0] = conn[0]->GetNumberOfCells();
381 num[1] = conn[1]->GetNumberOfCells();
382 }
383 else
384 {
385 vtkErrorMacro(<< "Do not support cells of dimension "
386 << cell->GetCellDimension());
387 }
388 }
389 else
390 {
391 if((cell->GetCellDimension())==3 )
392 {
393 if (orientation)
394 {
395 this->ClipHexahedron(newPoints,cell,this->Locator, conn[0],
396 inPD, outPD, inCD, cellId, outCD[0]);
397 }
398 else
399 {
400 this->ClipBox(newPoints,cell, this->Locator, conn[0],
401 inPD, outPD, inCD, cellId, outCD[0]);
402 }
403
404 numNew[0] = conn[0]->GetNumberOfCells() - num[0];
405 num[0] = conn[0]->GetNumberOfCells();
406 }
407 else if((cell->GetCellDimension())==2 )
408 {
409 if (orientation)
410 {
411 this->ClipHexahedron2D(newPoints,cell,this->Locator, conn[0],
412 inPD, outPD, inCD, cellId, outCD[0]);
413 }
414 else
415 {
416 this->ClipBox2D(newPoints,cell, this->Locator, conn[0],
417 inPD, outPD, inCD, cellId, outCD[0]);
418 }
419 numNew[0] = conn[0]->GetNumberOfCells() - num[0];
420 num[0] = conn[0]->GetNumberOfCells();
421 }
422 else if (cell->GetCellDimension() == 1)
423 {
424 if (orientation)
425 {
426 this->ClipHexahedron1D(newPoints,cell, this->Locator, conn[0],
427 inPD, outPD, inCD, cellId, outCD[0]);
428 }
429 else
430 {
431 this->ClipBox1D(newPoints,cell, this->Locator, conn[0],
432 inPD, outPD, inCD, cellId, outCD[0]);
433 }
434 numNew[0] = conn[0]->GetNumberOfCells() - num[0];
435 num[0] = conn[0]->GetNumberOfCells();
436 }
437 else if (cell->GetCellDimension() == 0)
438 {
439 if (orientation)
440 {
441 this->ClipHexahedron0D(cell, this->Locator, conn[0],
442 inPD, outPD, inCD, cellId, outCD[0]);
443 }
444 else
445 {
446 this->ClipBox0D(cell, this->Locator, conn[0],
447 inPD, outPD, inCD, cellId, outCD[0]);
448 }
449 numNew[0] = conn[0]->GetNumberOfCells() - num[0];
450 num[0] = conn[0]->GetNumberOfCells();
451 }
452 else
453 {
454 vtkErrorMacro(<< "Do not support cells of dimension "
455 << cell->GetCellDimension());
456 }
457 }
458
459 for (i=0 ; i<numOutputs; i++) // for both outputs
460 {
461 for (j=0; j < numNew[i]; j++)
462 {
463 locs[i]->InsertNextValue(conn[i]->GetTraversalLocation());
464 conn[i]->GetNextCell(npts,pts);
465
466 //For each new cell added, got to set the type of the cell
467 switch ( cell->GetCellDimension() )
468 {
469 case 0: //points are generated-------------------------------
470 cellType = (npts > 1 ? VTK_POLY_VERTEX : VTK_VERTEX);
471 break;
472
473 case 1: //lines are generated----------------------------------
474 cellType = (npts > 2 ? VTK_POLY_LINE : VTK_LINE);
475 break;
476
477 case 2: //polygons are generated------------------------------
478 cellType = (npts == 3 ? VTK_TRIANGLE :
479 (npts == 4 ? VTK_QUAD : VTK_POLYGON));
480 break;
481
482 case 3: //tetrahedra are generated------------------------------
483 cellType = VTK_TETRA;
484 break;
485 } //switch
486
487 newCellId = types[i]->InsertNextValue(cellType);
488 outCD[i]->CopyData(inCD, cellId, newCellId);
489 } //for each new cell
490 } // for both outputs
491 } //for each cell
492
493 cell->Delete();
494
495 output->SetPoints(newPoints);
496 output->SetCells(types[0], locs[0], conn[0]);
497
498 conn[0]->Delete();
499 types[0]->Delete();
500 locs[0]->Delete();
501
502 if ( this->GenerateClippedOutput )
503 {
504 clippedOutput->SetPoints(newPoints);
505 clippedOutput->SetCells(types[1], locs[1], conn[1]);
506 conn[1]->Delete();
507 types[1]->Delete();
508 locs[1]->Delete();
509 }
510 newPoints->Delete();
511 this->Locator->Initialize();//release any extra memory
512 output->Squeeze();
513
514 return 1;
515 }
516
517 //----------------------------------------------------------------------------
518 // Specify a spatial locator for merging points. By default,
519 // an instance of vtkMergePoints is used.
CreateDefaultLocator()520 void vtkBoxClipDataSet::CreateDefaultLocator()
521 {
522 if ( this->Locator == NULL )
523 {
524 this->Locator = vtkMergePoints::New();
525 this->Locator->Register(this);
526 this->Locator->Delete();
527 }
528 }
529
530 //----------------------------------------------------------------------------
531 // Set the box for clipping
532 // for each plane, specify the normal and one vertex on the plane.
533 //
SetBoxClip(const double * n0,const double * o0,const double * n1,const double * o1,const double * n2,const double * o2,const double * n3,const double * o3,const double * n4,const double * o4,const double * n5,const double * o5)534 void vtkBoxClipDataSet::SetBoxClip(const double *n0, const double *o0,
535 const double *n1, const double *o1,
536 const double *n2, const double *o2,
537 const double *n3, const double *o3,
538 const double *n4, const double *o4,
539 const double *n5, const double *o5)
540 {
541 if ( (this->Orientation == 1)
542
543 && (this->PlaneNormal[0][0] == n0[0])
544 && (this->PlaneNormal[0][1] == n0[1])
545 && (this->PlaneNormal[0][2] == n0[2])
546
547 && (this->PlaneNormal[1][0] == n1[0])
548 && (this->PlaneNormal[1][1] == n1[1])
549 && (this->PlaneNormal[1][2] == n1[2])
550
551 && (this->PlaneNormal[2][0] == n2[0])
552 && (this->PlaneNormal[2][1] == n2[1])
553 && (this->PlaneNormal[2][2] == n2[2])
554
555 && (this->PlaneNormal[3][0] == n3[0])
556 && (this->PlaneNormal[3][1] == n3[1])
557 && (this->PlaneNormal[3][2] == n3[2])
558
559 && (this->PlaneNormal[4][0] == n4[0])
560 && (this->PlaneNormal[4][1] == n4[1])
561 && (this->PlaneNormal[4][2] == n4[2])
562
563 && (this->PlaneNormal[5][0] == n5[0])
564 && (this->PlaneNormal[5][1] == n5[1])
565 && (this->PlaneNormal[5][2] == n5[2])
566
567 && (this->PlanePoint[0][0] == o0[0])
568 && (this->PlanePoint[0][1] == o0[1])
569 && (this->PlanePoint[0][2] == o0[2])
570
571 && (this->PlanePoint[1][0] == o1[0])
572 && (this->PlanePoint[1][1] == o1[1])
573 && (this->PlanePoint[1][2] == o1[2])
574
575 && (this->PlanePoint[2][0] == o2[0])
576 && (this->PlanePoint[2][1] == o2[1])
577 && (this->PlanePoint[2][2] == o2[2])
578
579 && (this->PlanePoint[3][0] == o3[0])
580 && (this->PlanePoint[3][1] == o3[1])
581 && (this->PlanePoint[3][2] == o3[2])
582
583 && (this->PlanePoint[4][0] == o4[0])
584 && (this->PlanePoint[4][1] == o4[1])
585 && (this->PlanePoint[4][2] == o4[2])
586
587 && (this->PlanePoint[5][0] == o5[0])
588 && (this->PlanePoint[5][1] == o5[1])
589 && (this->PlanePoint[5][2] == o5[2]) )
590 {
591 return;
592 }
593
594 this->SetOrientation(1);
595
596 this->PlaneNormal[0][0] = n0[0];
597 this->PlaneNormal[0][1] = n0[1];
598 this->PlaneNormal[0][2] = n0[2];
599
600 this->PlaneNormal[1][0] = n1[0];
601 this->PlaneNormal[1][1] = n1[1];
602 this->PlaneNormal[1][2] = n1[2];
603
604 this->PlaneNormal[2][0] = n2[0];
605 this->PlaneNormal[2][1] = n2[1];
606 this->PlaneNormal[2][2] = n2[2];
607
608 this->PlaneNormal[3][0] = n3[0];
609 this->PlaneNormal[3][1] = n3[1];
610 this->PlaneNormal[3][2] = n3[2];
611
612 this->PlaneNormal[4][0] = n4[0];
613 this->PlaneNormal[4][1] = n4[1];
614 this->PlaneNormal[4][2] = n4[2];
615
616 this->PlaneNormal[5][0] = n5[0];
617 this->PlaneNormal[5][1] = n5[1];
618 this->PlaneNormal[5][2] = n5[2];
619
620 this->PlanePoint[0][0] = o0[0];
621 this->PlanePoint[0][1] = o0[1];
622 this->PlanePoint[0][2] = o0[2];
623
624 this->PlanePoint[1][0] = o1[0];
625 this->PlanePoint[1][1] = o1[1];
626 this->PlanePoint[1][2] = o1[2];
627
628 this->PlanePoint[2][0] = o2[0];
629 this->PlanePoint[2][1] = o2[1];
630 this->PlanePoint[2][2] = o2[2];
631
632 this->PlanePoint[3][0] = o3[0];
633 this->PlanePoint[3][1] = o3[1];
634 this->PlanePoint[3][2] = o3[2];
635
636 this->PlanePoint[4][0] = o4[0];
637 this->PlanePoint[4][1] = o4[1];
638 this->PlanePoint[4][2] = o4[2];
639
640 this->PlanePoint[5][0] = o5[0];
641 this->PlanePoint[5][1] = o5[1];
642 this->PlanePoint[5][2] = o5[2];
643
644 this->Modified();
645 }
646
647 //----------------------------------------------------------------------------
648 // Specify the bounding box for clipping
649
SetBoxClip(double xmin,double xmax,double ymin,double ymax,double zmin,double zmax)650 void vtkBoxClipDataSet::SetBoxClip(double xmin,double xmax,
651 double ymin,double ymax,
652 double zmin,double zmax)
653 {
654 if ( (this->Orientation == 0)
655 && (this->BoundBoxClip[0][0] == xmin)
656 && (this->BoundBoxClip[0][1] == xmax)
657 && (this->BoundBoxClip[1][0] == ymin)
658 && (this->BoundBoxClip[1][1] == ymax)
659 && (this->BoundBoxClip[2][0] == zmin)
660 && (this->BoundBoxClip[2][1] == zmax) )
661 {
662 return;
663 }
664
665 this->SetOrientation(0);
666 this->BoundBoxClip[0][0] = xmin;
667 this->BoundBoxClip[0][1] = xmax;
668 this->BoundBoxClip[1][0] = ymin;
669 this->BoundBoxClip[1][1] = ymax;
670 this->BoundBoxClip[2][0] = zmin;
671 this->BoundBoxClip[2][1] = zmax;
672
673 this->Modified();
674 }
675
676 //----------------------------------------------------------------------------
FillInputPortInformation(int,vtkInformation * info)677 int vtkBoxClipDataSet::FillInputPortInformation(int, vtkInformation *info)
678 {
679 info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataSet");
680 return 1;
681 }
682
683 //----------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)684 void vtkBoxClipDataSet::PrintSelf(ostream& os, vtkIndent indent)
685 {
686 this->Superclass::PrintSelf(os,indent);
687
688 //os << indent << "Merge Tolerance: " << this->MergeTolerance << "\n";
689 os << indent << "Orientation: " << this->Orientation << "\n";
690
691 if ( this->Locator )
692 {
693 os << indent << "Locator: " << this->Locator << "\n";
694 }
695 else
696 {
697 os << indent << "Locator: (none)\n";
698 }
699
700 os << indent << "Generate Clipped Output: "
701 << (this->GenerateClippedOutput ? "Yes\n" : "Off\n");
702 os << indent << "Generate Clip Scalars: "
703 << (this->GenerateClipScalars ? "On\n" : "Off\n");
704 }
705
706 //-----------------------------------------------------------------------------
707 // InterpolateEdge: Interpolate the data in a vtkDataSetAttributes along a line
708 //
709 // This method works very much like vtkDataSetAttributes::InterpolateEdge
710 // except that rather than take the interpolation information from an
711 // input and copy it to an output, the values to interpolate are already
712 // placed in the output arrays. This is necessary because vtkBoxClipDataSet
713 // will continually cut (and consequently interpolate) the input until it
714 // fits within the bounds.
InterpolateEdge(vtkDataSetAttributes * attributes,vtkIdType toId,vtkIdType fromId1,vtkIdType fromId2,double t)715 void vtkBoxClipDataSet::InterpolateEdge(vtkDataSetAttributes *attributes,
716 vtkIdType toId,
717 vtkIdType fromId1, vtkIdType fromId2,
718 double t)
719 {
720 int numArrays = attributes->GetNumberOfArrays();
721 for (int i = 0; i < numArrays; i++)
722 {
723 vtkAbstractArray *array = attributes->GetAbstractArray(i);
724
725 // We are ignoring any special interpolate flags (such as nearest neighbor
726 // interpolation). That kind of interpolation is not linear and could be
727 // incorrect when multiple cuts are performed on the same primitive (which
728 // can happen).
729
730 array->InterpolateTuple(toId, fromId1, array, fromId2, array, t);
731 }
732 }
733
734 //----------------------------------------------------------------------------
735 // CellGrid: Subdivide cells in consistent tetrahedra.
736 // Case : Voxel(11) or Hexahedron(12).
737 //
738 // MinEdgF search the smallest vertex index in linear order of a face(4 vertices)
739 //
MinEdgeF(const unsigned int * id_v,const vtkIdType * cellIds,unsigned int * edgF)740 void vtkBoxClipDataSet::MinEdgeF(const unsigned int *id_v,
741 const vtkIdType *cellIds,
742 unsigned int *edgF)
743 {
744
745 int i;
746 unsigned int id;
747 int ids;
748 int min_f;
749
750 ids = 0;
751 id = id_v[0]; //Face index
752 min_f = cellIds[id_v[0]];
753
754 for(i=1; i<4; i++)
755 {
756 if(min_f > cellIds[id_v[i]])
757 {
758 min_f = cellIds[id_v[i]];
759 id = id_v[i];
760 ids= i;
761 }
762 }
763
764 switch(ids)
765 {
766 case 0:
767 if(id < id_v[2])
768 {
769 edgF[0] = id;
770 edgF[1] = id_v[2];
771 }
772 else
773 {
774 edgF[0] = id_v[2];
775 edgF[1] = id;
776 }
777 break;
778 case 1:
779 if(id < id_v[3])
780 {
781 edgF[0] = id;
782 edgF[1] = id_v[3];
783 }
784 else
785 {
786 edgF[0] = id_v[3];
787 edgF[1] = id;
788 }
789 break;
790 case 2:
791 if(id < id_v[0])
792 {
793 edgF[0] = id;
794 edgF[1] = id_v[0];
795 }
796 else
797 {
798 edgF[0] = id_v[0];
799 edgF[1] = id;
800 }
801 break;
802 case 3:
803 if(id < id_v[1])
804 {
805 edgF[0] = id;
806 edgF[1] = id_v[1];
807 }
808 else
809 {
810 edgF[0] = id_v[1];
811 edgF[1] = id;
812 }
813 break;
814 }
815 }
816
817 //----------------------------------------------------------------------------
818 // CellGrid: Subdivide cells in consistent tetrahedra.
819 //
820 // Case : Voxel or Hexahedron:
821 // if ( subdivide voxel in 6 tetrahedra)
822 // voxel : 2 wedges (*)
823 // else
824 // voxel : 5 tetrahedra
825 //
826 // Case : Wedge (*)
827 //
828 // ------------------------------------------------
829 //
830 //(*) WedgeToTetra: subdivide one wedge in 3 tetrahedra
831 //
832 // wedge : 1 tetrahedron + 1 pyramid = 3 tetrahedra.
833 //
WedgeToTetra(const vtkIdType * wedgeId,const vtkIdType * cellIds,vtkCellArray * newCellArray)834 void vtkBoxClipDataSet::WedgeToTetra(const vtkIdType *wedgeId,
835 const vtkIdType *cellIds,
836 vtkCellArray *newCellArray)
837 {
838 int i;
839 int id;
840 vtkIdType xmin;
841 vtkIdType tab[4];
842 vtkIdType tabpyram[5];
843
844 const vtkIdType vwedge[6][4]={ {0, 4, 3, 5}, {1, 4, 3, 5}, {2, 4, 3, 5},
845 {3, 0, 1, 2}, {4, 0, 1, 2}, {5, 0, 1, 2} };
846
847 // the table 'vwedge' set 6 possibilities of the smallest index
848 //
849 // v5
850 // /\ .
851 // v3 /..\ v4
852 // / /
853 // v2/\ /
854 // v0/__\/v1
855 // if(v0 index ) is the smallest index:
856 // wedge is subdivided in:
857 // 1 tetrahedron-> vwedge[0]: {v0,v4,v3,v5}
858 // and 1 pyramid -> vert[0] : {v1,v2,v5,v4,v0}
859 //
860
861 id = 0;
862 xmin = cellIds[wedgeId[0]];
863 for(i=1;i<6;i++)
864 {
865 if(xmin > cellIds[wedgeId[i]])
866 {
867 xmin = cellIds[wedgeId[i]];// the smallest global index
868 id = i; // local index
869 }
870 }
871 for (i =0;i<4;i++)
872 {
873 tab[i] = wedgeId[vwedge[id][i]];
874 }
875 newCellArray->InsertNextCell(4,tab);
876
877 // Pyramid :create 2 tetrahedra
878 const vtkIdType vert[6][5]={ {1, 2, 5, 4, 0}, {2, 0, 3, 5, 1},
879 {3, 0, 1, 4, 2}, {1, 2, 5, 4, 3},
880 {2, 0, 3, 5, 4}, {3, 0, 1, 4, 5} };
881 for(i=0;i<5;i++)
882 {
883 tabpyram[i] = wedgeId[vert[id][i]];
884 }
885 this->PyramidToTetra(tabpyram,cellIds,newCellArray);
886 }
887
888 //----------------------------------------------------------------------------
889 // CellGrid: Subdivide cells in consistent tetrahedra.
890 //
891 // PyramidToTetra :Subdivide the pyramid in consistent tetrahedra.
892 // Pyramid : 2 tetrahedra.
893 //
PyramidToTetra(const vtkIdType * pyramId,const vtkIdType * cellIds,vtkCellArray * newCellArray)894 void vtkBoxClipDataSet::PyramidToTetra(const vtkIdType *pyramId,
895 const vtkIdType *cellIds,
896 vtkCellArray *newCellArray)
897 {
898 vtkIdType xmin;
899 unsigned int i,j,idpy;
900 vtkIdType tab[4];
901
902 // the table 'vpy' set 3 possibilities of the smallest index
903 // vertices{v0,v1,v2,v3,v4}. {v0,v1,v2,v3} is a square face of pyramid
904 //
905 // v4
906 // ^
907 //
908 //
909 // v3 _ _ __ _ v2
910 // / /
911 // v0/_ _ _ _ _/v1
912 //
913 // if(v0 index ) is the smallest index:
914 // the pyramid is subdivided in:
915 // 2 tetrahedra-> vpy[0]: {v0,v1,v2,v4}
916 // vpy[1]: {v0,v2,v3,v4}
917 //
918 const vtkIdType vpy[8][4] ={{0,1,2,4},{0,2,3,4},{1,2,3,4},{1,3,0,4},
919 {2,3,0,4},{2,0,1,4},{3,0,1,4},{3,1,2,4}};
920
921 xmin = cellIds[pyramId[0]];
922 idpy = 0;
923 for(i=1;i<4;i++)
924 {
925 if(xmin > cellIds[pyramId[i]])
926 {
927 xmin = cellIds[pyramId[i]]; // global index
928 idpy = i; // local index
929 }
930 }
931 for(j = 0; j < 4 ; j++)
932 {
933 tab[j] = pyramId[vpy[2*idpy][j]];
934 }
935 newCellArray->InsertNextCell(4,tab);
936
937 for(j = 0; j < 4 ; j++)
938 {
939 tab[j] = pyramId[vpy[2*idpy+1][j]];
940 }
941 newCellArray->InsertNextCell(4,tab);
942 }
943
944
945 //----------------------------------------------------------------------------
946 //Tetra Grid : Subdivide cells in consistent tetrahedra.
947 // For each cell, search the smallest global index.
948 //
949 // Case Tetrahedron(10): Just insert this cell in the newCellArray
950 //
951 // Case Voxel(11) or Hexahedron(12):
952 // - for each face: looking for the diagonal edge with the smallest index
953 // - 2 possibilities: subdivide a cell in 5 or 6 tetrahedra
954 //
955 // (I)Case 6 tetrahedra:
956 // - 6 possibilities: subdivide de cell in 2 wedges:
957 // (1) diagonal edges (v0,v5),(v2,v7)
958 // vwedge[0]={0,5,4,2,7,6},
959 // vwedge[1]={0,1,5,2,3,7},
960 // (2) diagonal edges (v4,v7),(v0,v3)
961 // vwedge[2]={4,7,6,0,3,2},
962 // vwedge[3]={4,5,7,0,1,3},
963 // (3) diagonal edges (v0,v6),(v1,v7)
964 // vwedge[4]= {1,7,5,0,6,4},
965 // vwedge[5]{1,3,7,0,2,6},
966 // subdivide de cell in 2 wedges:
967 // (4) diagonal edges (v1,v2),(v5,v6)
968 // vwedge[6]={4,5,6,0,1,2},
969 // vwedge[7]={6,5,7,2,1,3},
970 // (5) diagonal edges (v2,v4),(v3,v5)
971 // vwedge[8]={3,7,5,2,6,4},
972 // vwedge[9]={1,3,5,0,2,4},
973 // (6) diagonal edges (v1,v4),(v3,v6)
974 // vwedge[10]={0,1,4,2,3,6}
975 // vwedge[11]={1,5,4,3,7,6}
976 // v6 _ _ __ _ v7
977 // /| /| VOXEL
978 // v4/_|_ _ _ _/ | opposite vertex of v0 is v7 and vice-versa
979 // | | v5| | diagonal edges Edg_f[i]
980 // |v2 _ _ _ |_| v3
981 // |/ |/
982 // v0/_ _ _ _ _|/v1
983 //
984 //
985 // (II)Case 5 tetrahedra:
986 // - search the smallest vertex vi
987 // - verify if the opposite vertices of vi do not belong to any diagonal edges Edg_f
988 // - 2 possibilites: create 5 tetraedra
989 // - if vi is ( 0 or 3 or 5 or 6)
990 // vtetra[]={v0,v5,v3,v6},{v0,v4,v5,v6},{v0,v1,v3,v5},{v5,v3,v6,v7},{v0,v3,v2,v6}};
991 // - if vi is ( 1 or 2 or 4 or 7)
992 // vtetra[]={v1,v2,v4,v7},{v0,v1,v2,v4},{v1,v4,v5,v7},{v1,v3,v2,v7},{v2,v6,v4,v7}};
993 // Case Wedge (13):
994 // the table 'vwedge' set 6 possibilities of the smallest index
995 //
996 // v5
997 // /\ .
998 // v3 /..\ v4
999 // / /
1000 // v2/\ /
1001 // v0/__\/v1
1002 //
1003 // if(v0 index ) is the smallest index:
1004 // wedge is subdivided in:
1005 // 1 tetrahedron-> vwedge[0]: {v0,v4,v3,v5}
1006 // and 1 pyramid-> vert[0] : {v1,v2,v5,v4,v0}
1007 //
1008 // Case Pyramid (14):
1009 // the table 'vpy' set 3 possibilities of the smallest index
1010 // vertices{v0,v1,v2,v3,v4}. {v0,v1,v2,v3} is a square face of pyramid
1011 //
1012 // v4 (opposite vertex of face with 4 vertices)
1013 // ^
1014 //
1015 //
1016 // v3 _ _ __ _ v2
1017 // / /
1018 // v0/_ _ _ _ _/v1
1019 //
1020 // if(v0 index ) is the smallest index:
1021 // the pyramid is subdivided in:
1022 // 2 tetrahedra-> vpyram[0]: {v0,v1,v2,v4}
1023 // vpyram[1]: {v0,v2,v3,v4}
1024 //
1025 //
1026 //----------------------------------------------------------------------------
CellGrid(vtkIdType typeobj,vtkIdType npts,const vtkIdType * cellIds,vtkCellArray * newCellArray)1027 void vtkBoxClipDataSet::CellGrid(vtkIdType typeobj, vtkIdType npts,
1028 const vtkIdType *cellIds,
1029 vtkCellArray *newCellArray)
1030 {
1031 vtkIdType tab[4];
1032 vtkIdType tabp[5];
1033 vtkIdType ptstriangle = 3;
1034 vtkIdType ptstetra = 4;
1035 vtkIdType xmin;
1036 vtkIdType idt;
1037 int i,j;
1038 unsigned int id =0;
1039 unsigned int idpy =0;
1040
1041 unsigned int Edg_f[6][2]; //edge selected of each face
1042 unsigned int idv[4];
1043
1044 unsigned int idopos;
1045 unsigned int numbertetra;
1046
1047 const vtkIdType triPassThrough[3] = {0, 1, 2};
1048 vtkIdType tri[3];
1049 vtkIdType line[2];
1050
1051 switch(typeobj)
1052 {
1053 case VTK_VERTEX:
1054 case VTK_POLY_VERTEX:
1055 for (idt = 0; idt < npts; idt++)
1056 {
1057 newCellArray->InsertNextCell(1, &idt);
1058 }
1059 break;
1060
1061 case VTK_LINE:
1062 case VTK_POLY_LINE:
1063 for (idt = 0; idt < npts-1; idt++)
1064 {
1065 line[0] = idt;
1066 line[1] = idt+1;
1067 newCellArray->InsertNextCell(2, line);
1068 }
1069 break;
1070
1071 case VTK_TRIANGLE: // 5
1072 case VTK_QUADRATIC_TRIANGLE:
1073 case VTK_BIQUADRATIC_TRIANGLE:
1074 newCellArray->InsertNextCell(ptstriangle, triPassThrough);
1075 break;
1076
1077 case VTK_TRIANGLE_STRIP: // 6
1078 for (idt=0 ; idt < npts-2; idt++)
1079 {
1080 if (idt%2 == 0)
1081 {
1082 tri[0] = idt;
1083 tri[1] = idt+1;
1084 tri[2] = idt+2;
1085 }
1086 else
1087 {
1088 tri[0] = idt;
1089 tri[1] = idt+2;
1090 tri[2] = idt+1;
1091 }
1092 newCellArray->InsertNextCell(3,tri);
1093 }
1094 break;
1095
1096 case VTK_POLYGON: // 7 (Convex case)
1097 tri[0] = 0;
1098 for (idt=2 ; idt < npts; idt++)
1099 {
1100 tri[1] = idt-1;
1101 tri[2] = idt;
1102 newCellArray->InsertNextCell(3,tri);
1103 }
1104 break;
1105
1106 case VTK_PIXEL: // 8
1107 {
1108 const vtkIdType vtrip[2][3] = {{0,1,3},{0,3,2}};
1109 newCellArray->InsertNextCell(3,vtrip[0]);
1110 newCellArray->InsertNextCell(3,vtrip[1]);
1111 }
1112 break;
1113
1114 case VTK_QUAD: // 9
1115 case VTK_QUADRATIC_QUAD:
1116 case VTK_BIQUADRATIC_QUAD:
1117 case VTK_QUADRATIC_LINEAR_QUAD:
1118 {
1119 const vtkIdType vtriq[2][3] = {{0,1,2},{0,2,3}};
1120 newCellArray->InsertNextCell(3,vtriq[0]);
1121 newCellArray->InsertNextCell(3,vtriq[1]);
1122 }
1123 break;
1124
1125 case VTK_TETRA: // 10
1126 case VTK_QUADRATIC_TETRA:
1127 {
1128 const vtkIdType tetra[4]={0,1,2,3};
1129 newCellArray->InsertNextCell(ptstetra,tetra);
1130 }
1131 break;
1132
1133 case VTK_VOXEL: // 11
1134 // each face: search edge with smallest global index
1135 // face 0 (0,1,5,4)
1136 idv[0]= 0;
1137 idv[1]= 1;
1138 idv[2]= 5;
1139 idv[3]= 4;
1140 this->MinEdgeF(idv,cellIds,Edg_f[0]);
1141
1142 // face 1 (0,1,3,2)
1143 idv[0]= 0;
1144 idv[1]= 1;
1145 idv[2]= 3;
1146 idv[3]= 2;
1147 this->MinEdgeF(idv,cellIds,Edg_f[1]);
1148 // face 2 (0,2,6,4)
1149 idv[0]= 0;
1150 idv[1]= 2;
1151 idv[2]= 6;
1152 idv[3]= 4;
1153 this->MinEdgeF(idv,cellIds,Edg_f[2]);
1154 // face 3 (4,5,7,6)
1155 idv[0]= 4;
1156 idv[1]= 5;
1157 idv[2]= 7;
1158 idv[3]= 6;
1159 this->MinEdgeF(idv,cellIds,Edg_f[3]);
1160 // face 4 (2,3,7,6)
1161 idv[0]= 2;
1162 idv[1]= 3;
1163 idv[2]= 7;
1164 idv[3]= 6;
1165 this->MinEdgeF(idv,cellIds,Edg_f[4]);
1166 // face 5 (1,3,7,5)
1167 idv[0]= 1;
1168 idv[1]= 3;
1169 idv[2]= 7;
1170 idv[3]= 5;
1171 this->MinEdgeF(idv,cellIds,Edg_f[5]);
1172
1173 //search the smallest global index of voxel
1174 xmin = cellIds[0];
1175 id = 0;
1176 for(i=1;i<8;i++)
1177 {
1178 if(xmin > cellIds[i])
1179 {
1180 xmin = cellIds[i];// the smallest global index
1181 id = i; // local index
1182 }
1183 }
1184 //two cases:
1185 idopos = 7 - id;
1186 numbertetra = 5;
1187 for(i=0;i<6;i++)
1188 {
1189 j=0;
1190 if (idopos == Edg_f[i][j])
1191 {
1192 numbertetra = 6;
1193 break;
1194 }
1195 j=1;
1196 if (idopos == Edg_f[i][j])
1197 {
1198 numbertetra = 6;
1199 break;
1200 }
1201 }
1202
1203 if(numbertetra == 5)
1204 {
1205 // case 1: create 5 tetraedra
1206 if((id == 0)||(id == 3)||(id == 5)||(id == 6))
1207 {
1208 const vtkIdType vtetra[5][4]={{0,5,3,6},{0,4,5,6},
1209 {0,1,3,5},{5,3,6,7},{0,3,2,6}};
1210 for(i=0; i<5; i++)
1211 {
1212 newCellArray->InsertNextCell(4,vtetra[i]);
1213 }
1214 }
1215 else
1216 {
1217 const vtkIdType vtetra[5][4]={{1,2,4,7},{0,1,2,4},
1218 {1,4,5,7},{1,3,2,7},{2,6,4,7}};
1219
1220 for(i=0; i<5; i++)
1221 {
1222 newCellArray->InsertNextCell(4,vtetra[i]);
1223 }
1224 }
1225 }
1226 else
1227 {
1228 //case 2: create 2 wedges-> 6 tetrahedra
1229 const vtkIdType vwedge[12][6]={{0,5,4,2,7,6},{0,1,5,2,3,7},
1230 {4,7,6,0,3,2},{4,5,7,0,1,3},
1231 {1,7,5,0,6,4},{1,3,7,0,2,6},
1232 {4,5,6,0,1,2},{6,5,7,2,1,3},
1233 {3,7,5,2,6,4},{1,3,5,0,2,4},
1234 {0,1,4,2,3,6},{1,5,4,3,7,6}};
1235 unsigned int edgeId = 10*Edg_f[i][0]+ Edg_f[i][1];
1236 switch(edgeId)
1237 {
1238 case 5: // edge(v0,v5):10*0 + 5
1239 case 27: // edge(v2,v7):10*2 + 7
1240 this->WedgeToTetra(vwedge[0],cellIds,newCellArray);
1241 this->WedgeToTetra(vwedge[1],cellIds,newCellArray);
1242 break;
1243 case 3: // edge(v0,v2)
1244 case 47: // edge(v4,v6)
1245 this->WedgeToTetra(vwedge[2],cellIds,newCellArray);
1246 this->WedgeToTetra(vwedge[3],cellIds,newCellArray);
1247 break;
1248 case 6:
1249 case 17:
1250 this->WedgeToTetra(vwedge[4],cellIds,newCellArray);
1251 this->WedgeToTetra(vwedge[5],cellIds,newCellArray);
1252 break;
1253 case 12:
1254 case 56:
1255 this->WedgeToTetra(vwedge[6],cellIds,newCellArray);
1256 this->WedgeToTetra(vwedge[7],cellIds,newCellArray);
1257 break;
1258 case 24:
1259 case 35:
1260 this->WedgeToTetra(vwedge[8],cellIds,newCellArray);
1261 this->WedgeToTetra(vwedge[9],cellIds,newCellArray);
1262 break;
1263 case 14:
1264 case 36:
1265 this->WedgeToTetra(vwedge[10],cellIds,newCellArray);
1266 this->WedgeToTetra(vwedge[11],cellIds,newCellArray);
1267 break;
1268 }
1269 }
1270 break;
1271
1272 case VTK_HEXAHEDRON: // 12
1273 case VTK_QUADRATIC_HEXAHEDRON:
1274 case VTK_TRIQUADRATIC_HEXAHEDRON:
1275 case VTK_BIQUADRATIC_QUADRATIC_HEXAHEDRON:
1276 {
1277 // each face: search edge with smallest global index
1278 // face 0 (0,1,5,4)
1279 idv[0]= 0;
1280 idv[1]= 1;
1281 idv[2]= 5;
1282 idv[3]= 4;
1283 this->MinEdgeF(idv,cellIds,Edg_f[0]);
1284
1285 // face 1 (0,1,2,3)
1286 idv[0]= 0;
1287 idv[1]= 1;
1288 idv[2]= 2;
1289 idv[3]= 3;
1290 this->MinEdgeF(idv,cellIds,Edg_f[1]);
1291 // face 2 (0,3,7,4)
1292 idv[0]= 0;
1293 idv[1]= 3;
1294 idv[2]= 7;
1295 idv[3]= 4;
1296 this->MinEdgeF(idv,cellIds,Edg_f[2]);
1297 // face 3 (4,5,6,7)
1298 idv[0]= 4;
1299 idv[1]= 5;
1300 idv[2]= 6;
1301 idv[3]= 7;
1302 this->MinEdgeF(idv,cellIds,Edg_f[3]);
1303 // face 4 (3,2,6,7)
1304 idv[0]= 3;
1305 idv[1]= 2;
1306 idv[2]= 6;
1307 idv[3]= 7;
1308 this->MinEdgeF(idv,cellIds,Edg_f[4]);
1309 // face 5 (1,2,6,5)
1310 idv[0]= 1;
1311 idv[1]= 2;
1312 idv[2]= 6;
1313 idv[3]= 5;
1314 this->MinEdgeF(idv,cellIds,Edg_f[5]);
1315
1316 //search the smallest global index of voxel
1317 xmin = cellIds[0];
1318 id = 0;
1319 for(i=1;i<8;i++)
1320 {
1321 if(xmin > cellIds[i])
1322 {
1323 xmin = cellIds[i];// the smallest global index
1324 id = i; // local index
1325 }
1326 }
1327
1328 //two cases:
1329 const unsigned int tabopos[8] = {6,7,4,5,2,3,0,1};
1330 idopos = tabopos[id];
1331 numbertetra = 5;
1332 for(i=0;i<6;i++)
1333 {
1334 j = 0;
1335 if (idopos == Edg_f[i][j])
1336 {
1337 numbertetra = 6;
1338 break;
1339 }
1340 j=1;
1341 if (idopos == Edg_f[i][j])
1342 {
1343 numbertetra = 6;
1344 break;
1345 }
1346 }
1347
1348 if(numbertetra == 5)
1349 {
1350 // case 1: create 5 tetraedra
1351 if((id == 0)||(id == 2)||(id == 5)||(id == 7))
1352 {
1353 const vtkIdType vtetra[5][4]={{0,5,2,7},{0,4,5,7},
1354 {0,1,2,5},{5,2,7,6},{0,2,3,7}};
1355 for(i=0; i<5; i++)
1356 {
1357 newCellArray->InsertNextCell(4,vtetra[i]);
1358 }
1359 }
1360 else
1361 {
1362 const vtkIdType vtetra[5][4]={{1,3,4,6},{0,1,3,4},
1363 {1,4,5,6},{1,2,3,6},{3,7,4,6}};
1364
1365 for(i=0; i<5; i++)
1366 {
1367 newCellArray->InsertNextCell(4,vtetra[i]);
1368 }
1369 }
1370 }
1371 else
1372 {
1373 //case 2: create 2 wedges-> 6 tetrahedra
1374 const vtkIdType vwedge[12][6]={{0,5,4,3,6,7},{0,1,5,3,2,6},
1375 {4,6,7,0,2,3},{4,5,6,0,1,2},
1376 {1,6,5,0,7,4},{1,2,6,0,3,7},
1377 {4,5,7,0,1,3},{7,5,6,3,1,2},
1378 {2,6,5,3,7,4},{1,2,5,0,3,4},
1379 {0,1,4,3,2,7},{1,5,4,2,6,7}};
1380 unsigned int edgeId = 10*Edg_f[i][0]+ Edg_f[i][1];
1381
1382 switch(edgeId)
1383 {
1384 case 5: // edge(v0,v5):10*0 + 5
1385 case 36: // edge(v3,v6):10*3 + 6
1386 this->WedgeToTetra(vwedge[0],cellIds,newCellArray);
1387 this->WedgeToTetra(vwedge[1],cellIds,newCellArray);
1388 break;
1389 case 2: // edge(v0,v2)
1390 case 46: // edge(v4,v6)
1391 this->WedgeToTetra(vwedge[2],cellIds,newCellArray);
1392 this->WedgeToTetra(vwedge[3],cellIds,newCellArray);
1393 break;
1394 case 7:
1395 case 16:
1396 this->WedgeToTetra(vwedge[4],cellIds,newCellArray);
1397 this->WedgeToTetra(vwedge[5],cellIds,newCellArray);
1398 break;
1399 case 13:
1400 case 57:
1401 this->WedgeToTetra(vwedge[6],cellIds,newCellArray);
1402 this->WedgeToTetra(vwedge[7],cellIds,newCellArray);
1403 break;
1404 case 34:
1405 case 25:
1406 this->WedgeToTetra(vwedge[8],cellIds,newCellArray);
1407 this->WedgeToTetra(vwedge[9],cellIds,newCellArray);
1408 break;
1409 case 14:
1410 case 27:
1411 this->WedgeToTetra(vwedge[10],cellIds,newCellArray);
1412 this->WedgeToTetra(vwedge[11],cellIds,newCellArray);
1413 break;
1414 }
1415 }
1416 }
1417 break;
1418
1419 case VTK_WEDGE:
1420 case VTK_QUADRATIC_WEDGE:
1421 case VTK_QUADRATIC_LINEAR_WEDGE:
1422 case VTK_BIQUADRATIC_QUADRATIC_WEDGE:
1423 if(npts == 6) //create 3 tetrahedra
1424 {
1425 //first tetrahedron
1426 const vtkIdType vwedge[6][4]={{0,4,3,5},{1,4,3,5},{2,4,3,5},
1427 {3,0,1,2},{4,0,1,2},{5,0,1,2}};
1428 xmin = cellIds[0];
1429 id = 0;
1430 for(i=1;i<6;i++)
1431 {
1432 if(xmin > cellIds[i])
1433 {
1434 xmin = cellIds[i]; // the smallest global index
1435 id = i; // local index
1436 }
1437 }
1438 newCellArray->InsertNextCell(4, vwedge[id]);
1439
1440 //Pyramid :create 2 tetrahedra
1441
1442 const vtkIdType vert[6][5]={{1,2,5,4,0},{2,0,3,5,1},{3,0,1,4,2},
1443 {1,2,5,4,3},{2,0,3,5,4},{3,0,1,4,5}};
1444 const vtkIdType vpy[8][4] ={{0,1,2,4},{0,2,3,4},{1,2,3,4},{1,3,0,4},
1445 {2,3,0,4},{2,0,1,4},{3,0,1,4},{3,1,2,4}};
1446 xmin = cellIds[vert[id][0]];
1447 tabp[0] = vert[id][0];
1448 idpy = 0;
1449 for(i=1;i<4;i++)
1450 {
1451 tabp[i] = vert[id][i];
1452 if(xmin > cellIds[vert[id][i]])
1453 {
1454 xmin = cellIds[vert[id][i]]; // global index
1455 idpy = i; // local index
1456 }
1457 }
1458 tabp[4] = vert[id][4];
1459 for(j = 0; j < 4 ; j++)
1460 {
1461 tab[j] = tabp[vpy[2*idpy][j]];
1462 }
1463 newCellArray->InsertNextCell(4,tab);
1464
1465 for(j = 0; j < 4 ; j++)
1466 {
1467 tab[j] = tabp[vpy[2*idpy+1][j]];
1468 }
1469 newCellArray->InsertNextCell(4,tab);
1470
1471 }
1472 else
1473 {
1474 vtkErrorMacro( << " This cell is not a wedge\n" );
1475 return;
1476 }
1477 break;
1478
1479 case VTK_PYRAMID: //Create 2 tetrahedra
1480 case VTK_QUADRATIC_PYRAMID:
1481 if(npts == 5)
1482 {
1483 //note: the first element vpyram[][0] is the smallest index of pyramid
1484 const vtkIdType vpyram[8][4]={{0,1,2,4},{0,2,3,4},{1,2,3,4},{1,3,0,4},
1485 {2,3,0,4},{2,0,1,4},{3,0,1,4},{3,1,2,4}};
1486 xmin = cellIds[0];
1487 id = 0;
1488 for(i=1;i<4;i++)
1489 {
1490 if(xmin > cellIds[i])
1491 {
1492 xmin = cellIds[i]; // the smallest global index of square face
1493 id = i; // local index
1494 }
1495 }
1496 newCellArray->InsertNextCell(4,vpyram[2*id]);
1497 newCellArray->InsertNextCell(4,vpyram[2*id+1]);
1498 }
1499 else
1500 {
1501 vtkErrorMacro( << " This cell is not a pyramid\n" );
1502 return;
1503 }
1504 break;
1505 }
1506 }
1507
1508 //----------------------------------------------------------------------------
1509 // The new cell created in intersection between tetrahedron and plane
1510 // are tetrahedron or wedges or pyramides.
1511 //
1512 // CreateTetra is used to subdivide wedges and pyramids in tetrahedron. The
1513 // proper vertex ordering for wedges and pyramids can be found in "The
1514 // Visualization Toolkit." In the third edition, they are in Figure 5-2 on page
1515 // 115 in section 5.4 ("Cell Types") in the "Basic Data Representation" chapter.
1516 //
CreateTetra(vtkIdType npts,const vtkIdType * cellIds,vtkCellArray * newCellArray)1517 void vtkBoxClipDataSet::CreateTetra(vtkIdType npts, const vtkIdType *cellIds,
1518 vtkCellArray *newCellArray)
1519 {
1520 vtkIdType tabp[5];
1521 vtkIdType tab[3][4];
1522
1523 unsigned int i,j;
1524 unsigned int id =0;
1525 unsigned int idpy;
1526
1527 vtkIdType xmin;
1528
1529 if (npts == 6)
1530 {
1531 //VTK_WEDGE: Create 3 tetrahedra
1532 //first tetrahedron
1533 const vtkIdType vwedge[6][4]={{0,4,3,5},{1,4,3,5},{2,4,3,5},
1534 {3,0,1,2},{4,0,1,2},{5,0,1,2}};
1535 xmin = cellIds[0];
1536 id = 0;
1537 for(i=1;i<6;i++)
1538 {
1539 if(xmin > cellIds[i])
1540 {
1541 xmin = cellIds[i];// the smallest global index
1542 id = i; // local index
1543 }
1544 }
1545 for(j = 0; j < 4 ; j++)
1546 {
1547 tab[0][j] = cellIds[vwedge[id][j]];
1548 }
1549 newCellArray->InsertNextCell(4,tab[0]);
1550
1551 //Pyramid: create 2 tetrahedra
1552
1553 const vtkIdType vert[6][5]= {{1,2,5,4,0},{2,0,3,5,1},{3,0,1,4,2},
1554 {1,2,5,4,3},{2,0,3,5,4},{3,0,1,4,5}};
1555 const vtkIdType vpy[8][4] = {{0,1,2,4},{0,2,3,4},{1,2,3,4},{1,3,0,4},
1556 {2,3,0,4},{2,0,1,4},{3,0,1,4},{3,1,2,4}};
1557 xmin = cellIds[vert[id][0]];
1558 tabp[0] = vert[id][0];
1559 idpy = 0;
1560 for(i=1;i<4;i++)
1561 {
1562 tabp[i] = vert[id][i];
1563 if(xmin > cellIds[vert[id][i]])
1564 {
1565 xmin = cellIds[vert[id][i]]; // global index
1566 idpy = i; // local index
1567 }
1568 }
1569 tabp[4] = vert[id][4];
1570 for(j = 0; j < 4 ; j++)
1571 {
1572 tab[1][j] = cellIds[tabp[vpy[2*idpy][j]]];
1573 }
1574 newCellArray->InsertNextCell(4,tab[1]);
1575
1576 for(j = 0; j < 4 ; j++)
1577 {
1578 tab[2][j] = cellIds[tabp[vpy[2*idpy+1][j]]];
1579 }
1580 newCellArray->InsertNextCell(4,tab[2]);
1581 }
1582 else
1583 {
1584 //VTK_PYRAMID: Create 2 tetrahedra
1585 //The first element in each set is the smallest index of pyramid
1586 const vtkIdType vpyram[8][4]={{0,1,2,4},{0,2,3,4},{1,2,3,4},{1,3,0,4},
1587 {2,3,0,4},{2,0,1,4},{3,0,1,4},{3,1,2,4}};
1588 xmin = cellIds[0];
1589 id = 0;
1590 for(i=1;i<4;i++)
1591 {
1592 if(xmin > cellIds[i])
1593 {
1594 xmin = cellIds[i]; // the smallest global index of face with 4 vertices
1595 id = i; // local index
1596 }
1597 }
1598 for(j = 0; j < 4 ; j++)
1599 {
1600 tab[0][j] = cellIds[vpyram[2*id][j]];
1601 }
1602 newCellArray->InsertNextCell(4,tab[0]);
1603
1604 for(j = 0; j < 4 ; j++)
1605 {
1606 tab[1][j] = cellIds[vpyram[2*id+1][j]];
1607 }
1608 newCellArray->InsertNextCell(4,tab[1]);
1609 }
1610 }
1611
1612 //----------------------------------------------------------------------------
1613 // Clip each cell of an unstructured grid.
1614 //
1615 //----------------------------------------------------------------------------
1616 //(1) How decide when the cell is NOT outside
1617 //
1618 // Explaining with an example in 2D.
1619 // Look at 9 region in the picture and the triangle represented there.
1620 // v0,v1,v2 are vertices of triangle T.
1621 //
1622 // | |
1623 // 1 | 2 | 3
1624 // _ _ _ _ _ _ _ _ _ _ _ _ _ _ ymax
1625 // | | v1
1626 // 4 | 5 |/\ 6
1627 // | / \ .
1628 // _ _ _ _ _ _ _ _ _ /| _ \_ _ _ymin
1629 // | v2/_|_ _ \v0
1630 // 7 | 8 | 9
1631 // xmin xmax
1632 //
1633 // set test={1,1,1,1} (one test for each plane: test={xmin,xmax,ymin,ymax} )
1634 // for each vertex, if the test is true set 0 in the test table:
1635 // vO > xmin ?, v0 < xmax ?, v0 > ymin ?, v0 < ymax ?
1636 // In the example: test={0,1,1,0}
1637 // v1 > xmin ?, v1 < xmax ?, v1 > ymin ?, v1 < ymax ?
1638 // In the example: test={0,1,0,0}
1639 // v2 > xmin ?, v2 < xmax ?, v2 > ymin ?, v2 < ymax ?
1640 // In the Example: test={0,0,0,0} -> triangle is NOT outside.
1641 //
1642 //
1643 // In general, look at the possibilities of each region
1644 //
1645 // (1,0,0,1) | (0,0,0,1) | (0,1,0,1)
1646 // - - - - - - - - - - - - - - - - - -
1647 // (1,0,0,0) | (0,0,0,0) | (0,1,0,0)
1648 // - - - - - - - - - - - - - - - - - -
1649 // (1,0,1,0) | (0,0,1,0) | (0,1,1,0)
1650 //
1651 // In the case above we have:(v0,v1,v2)
1652 // (1,1,1,1)->(0,1,1,0)->(0,1,0,0)->(0,0,0,0)
1653 //
1654 // If you have one vertex in region 5, this triangle is NOT outside
1655 // The triangle IS outside if (v0,v1,v2) are in region like
1656 // (1,2,3): (1,0,0,1)->(0,0,0,1)->(0,0,0,1)
1657 // (1,4,7): (1,0,0,1)->(1,0,0,0)->(1,0,0,0)
1658 // (7,8,9): (1,0,1,0)->(0,0,1,0)->(0,0,1,0)
1659 // (9,6,3): (0,1,1,0)->(0,1,0,0)->(0,1,0,0)
1660 //
1661 // Note that if the triangle T is on the region 5, T is NOT outside
1662 // In fact, T is inside
1663 //
1664 // You can extend this idea to 3D.
1665 //
1666 // Note: xmin = this->BoundBoxClip[0][0], xmax= this->BoundBoxClip[0][1],...
1667 //
1668 //----------------------------------------------------------------------------
1669 // (2) Intersection between Tetrahedron and Plane:
1670 // Description:
1671 // vertices of tetrahedron {v0,v1,v2,v3}
1672 // edge e1 : (v0,v1), edge e2 : (v1,v2)
1673 // edge e3 : (v2,v0), edge e4 : (v0,v3)
1674 // edge e5 : (v1,v3), edge e6 : (v2,v3)
1675 //
1676 // intersecting points p0, pi, ...
1677 //
1678 // Note: The algorithm search intersection
1679 // points with these edge order.
1680 //
1681 // v3 v3
1682 // e6/|\ e5 |
1683 // / | \ e2 |
1684 // v2/_ |_ \ v1 v2 - - -v1 |e4
1685 // \ | / |
1686 // e3\ | /e1 |
1687 // \|/ |
1688 // v0 v0
1689 //
1690 // (a) Intersecting 4 edges: see tab4[]
1691 // -> create: 2 wedges
1692 // -> 3 cases
1693 // ------------------------------------------
1694 // case 1246:
1695 // if (plane intersecting edges {e1,e2,e4,e6})
1696 // then p0 belongs to e1, p1 belongs to e2,
1697 // p2 belongs to e4, p3 belongs to e6.
1698 //
1699 // v3
1700 // | v3
1701 // | /
1702 // v1 v2 - *- -v1 | * p3
1703 // / p1 * p2 /
1704 // *p0 | v2/
1705 // / |
1706 // v0 v0
1707 // (e1) (e2) (e4) (e6)
1708 //
1709 // So, there are two wedges:
1710 // {p0,v1,p1,p2,v3,p3} {p2,v0,p0,p3,v2,p1}
1711 //
1712 // Note: if e1 and e2 are intersected by plane,
1713 // and the plane intersects 4 edges,
1714 // the edge e5 could not be intersected
1715 // (skew lines, do not create a planar face)
1716 // neither the edge e3((e1,e2,e3) is a face)
1717 //
1718 // ------------------------------------------
1719 // case 2345:
1720 // if (plane intersecting edges {e2,e3,e4,e5})
1721 // The two wedges are:
1722 // {p3,v3,p2,p0,v2,p1}, {p1,v0,p2,p0,v1,p3}
1723 //
1724 // ------------------------------------------
1725 // case 1356:
1726 // if (plane intersecting edges {e1,e3,e5,e6})
1727 // The two wedges are:
1728 // {p0,v0,p1,p2,v3,p3}, {p0,v1,p2,p1,v2,p3}
1729 //
1730 // -----------------------------------------
1731 //
1732 // (b) Intersecting 3 edges: tab3[]
1733 // ->create: 1 tetrahedron + 1 wedge
1734 // -> 4 cases
1735 // ------------------------------------------
1736 // case 134:
1737 // if (plane intersecting edges {e1,e3,e4})
1738 // then p0 belongs to e1, p1 belongs to e3,
1739 // p2 belongs to e4.
1740 //
1741 // v3
1742 // |
1743 // |
1744 // v2 *p2 v1
1745 // \ | /
1746 // p1* | *p0
1747 // \|/
1748 // v0
1749 //
1750 //
1751 // So, there are:
1752 // 1 tetrahedron {v0,p0,p2,p1)
1753 // and 1 wedge: {p0,p2,p1,v1,v3,v2}: tab3[0]
1754 //
1755 // Note:(a)if e1 and e3 are intersected by plane,
1756 // and the plane intersects 3 edges,
1757 // the edges e5 could not be intersected,
1758 // if so, e6 be intersected too and the
1759 // plane intersect 4 edges.
1760 // (b) tetrahedron vertices:
1761 // Use the first three indices of tab3:
1762 // {v0,p(0),p(2),p(1)}
1763 //
1764 // ------------------------------------------
1765 // case 125:
1766 // if (plane intersecting edges {e1,e2,e5})
1767 // There are :
1768 // 1 tetrahedron: {v1,p0,p2,p1}
1769 // 1 wedge : {p0,p2,p1,v0,v3,v2},
1770 // ------------------------------------------
1771 // case 236:
1772 // if (plane intersecting edges {e2,e3,e6})
1773 // There are :
1774 // 1 tetrahedron: {v2,p0,p1,p2}
1775 // 1 wedge : {p0,p1,p2,v1,v0,v3},
1776 // ------------------------------------------
1777 // case 456:
1778 // if (plane intersecting edges {e4,e5,e6})
1779 // There are :
1780 // 1 tetrahedron: {v3,p0,p1,p2}
1781 // 1 wedge : {p0,p1,p2,v0,v1,v2},
1782 // ------------------------------------------
1783 //
1784 // (c) Intersecting 2 edges and 1 vertex: tab2[]
1785 // -> create: 1 tetrahedron + 1 pyramid
1786 // -> 12 cases
1787 // ------------------------------------------
1788 // case 12: v3 indexes:{0,0,1,2,3},
1789 // * 0 1 2 3 4
1790 // e6/ | \ e5
1791 // / |p1\ tetrahedron: indices(*,4,1,2)
1792 // v2/_ _|_*_\ v1
1793 // \ | /
1794 // e3\ |p0*
1795 // \ | /e1
1796 // v0
1797 //
1798 // if (plane intersecting edges {e1,e2}
1799 // and one vertex)
1800 // There are
1801 // 1 tetrahedron: {v1,v3,p0,p1}
1802 // 1 pyramid : {v0,p0,p1,v2,v3},
1803 //
1804 // Note: if e1 and e2 are intersected by plane,
1805 // and the plane intersects 2 edges,
1806 // then the vertex is v3.
1807 // ------------------------------------------
1808 //
1809 // case 13: indexes{2,1,0,1,3}
1810 // Intersecting edges {e1,e3}
1811 // Intersectinf vertes v3,
1812 //
1813 // 1 tetrahedron: {v0,v3,p1,p0}
1814 // 1 pyramid : {v2,p1,p0,v1,v3},
1815 // ------------------------------------------
1816 // all cases see tab2[]
1817 // ------------------------------------------
1818 // (d) Intersecting 1 edges and 2 vertices: tab1[]
1819 // -> create: 2 tetrahedra
1820 // -> 6 cases
1821 //
1822 // - edge e1, vertices {v2,v3}(e6)
1823 //
1824 // v3
1825 // * tetrahedra: {p0,v2,v3,v1},{p0,v3,v2,v0}
1826 // e6/ | \ e5
1827 // / | \ .
1828 // v2 *_ _|_ _\ v1
1829 // \ | /
1830 // e3\ |p0*
1831 // \ | /e1
1832 // v0
1833 // - other cases see tab1[]
1834 //----------------------------------------------------------------------------
1835 //
ClipBox(vtkPoints * newPoints,vtkGenericCell * cell,vtkIncrementalPointLocator * locator,vtkCellArray * tets,vtkPointData * inPD,vtkPointData * outPD,vtkCellData * inCD,vtkIdType cellId,vtkCellData * outCD)1836 void vtkBoxClipDataSet::ClipBox(vtkPoints *newPoints,
1837 vtkGenericCell *cell,
1838 vtkIncrementalPointLocator *locator,
1839 vtkCellArray *tets,
1840 vtkPointData *inPD,
1841 vtkPointData *outPD,
1842 vtkCellData *inCD,
1843 vtkIdType cellId,
1844 vtkCellData *outCD)
1845 {
1846 vtkIdType cellType = cell->GetCellType();
1847 vtkIdList *cellIds = cell->GetPointIds();
1848 vtkCellArray *arraytetra = vtkCellArray::New();
1849 vtkPoints *cellPts = cell->GetPoints();
1850 vtkIdType npts = cellPts->GetNumberOfPoints();
1851 std::vector<vtkIdType> cellptId(npts);
1852 vtkIdType iid[4];
1853 vtkIdType *v_id = NULL;
1854 vtkIdType *verts, v1, v2;
1855 vtkIdType ptId;
1856 vtkIdType tab_id[6];
1857 vtkIdType ptstetra = 4;
1858
1859 vtkIdType i,j;
1860 unsigned int allInside;
1861 unsigned int idcellnew;
1862 unsigned int cutInd;
1863
1864 vtkIdType edges[6][2] = { {0,1}, {1,2}, {2,0},
1865 {0,3}, {1,3}, {2,3} }; /* Edges Tetrahedron */
1866 double value,deltaScalar;
1867 double t;
1868 double *p1, *p2;
1869 double x[3], v[3];
1870 double v_tetra[4][3];
1871
1872 for (i=0; i<npts; i++)
1873 {
1874 cellptId[i] = cellIds->GetId(i);
1875 }
1876
1877 // Convert all volume cells to tetrahedra
1878 this->CellGrid(cellType,npts,&cellptId[0],arraytetra);
1879 unsigned int totalnewtetra = arraytetra->GetNumberOfCells();
1880 unsigned int idtetranew;
1881
1882 for (idtetranew = 0 ; idtetranew < totalnewtetra; idtetranew++)
1883 {
1884 arraytetra->GetNextCell(ptstetra,v_id);
1885
1886 for (allInside=1, i=0; i<4; i++)
1887 {
1888 cellPts->GetPoint(v_id[i],v);
1889
1890 if (!(((v[0] >= this->BoundBoxClip[0][0])&&(v[0] <= this->BoundBoxClip[0][1]) &&
1891 (v[1] >= this->BoundBoxClip[1][0])&&(v[1] <= this->BoundBoxClip[1][1])&&
1892 (v[2] >= this->BoundBoxClip[2][0])&&(v[2] <= this->BoundBoxClip[2][1]))))
1893 {
1894 //outside,its type might change later (nearby intersection)
1895 allInside = 0;
1896 }
1897 }//for all points of the tetrahedron.
1898
1899 // Test Outside: see(1)
1900 if (!allInside)
1901 {
1902 unsigned int test[6] = {1,1,1,1,1,1};
1903 for (i=0; i<4; i++)
1904 {
1905 cellPts->GetPoint(v_id[i],v);
1906
1907 if (v[0] > this->BoundBoxClip[0][0])
1908 {
1909 test[0] = 0;
1910 }
1911 if (v[0] < this->BoundBoxClip[0][1])
1912 {
1913 test[1] = 0;
1914 }
1915 if (v[1] > this->BoundBoxClip[1][0])
1916 {
1917 test[2] = 0;
1918 }
1919 if (v[1] < this->BoundBoxClip[1][1])
1920 {
1921 test[3] = 0;
1922 }
1923 if (v[2] > this->BoundBoxClip[2][0])
1924 {
1925 test[4] = 0;
1926 }
1927 if (v[2] < this->BoundBoxClip[2][1])
1928 {
1929 test[5] = 0;
1930 }
1931
1932 }//for all points of the cell.
1933
1934 if ((test[0] == 1)|| (test[1] == 1) ||
1935 (test[2] == 1)|| (test[3] == 1) ||
1936 (test[4] == 1)|| (test[5] == 1))
1937 {
1938 continue; // Tetrahedron is outside.
1939 }
1940 }//if not all inside.
1941
1942 for (i=0; i<4; i++)
1943 {
1944 ptId = cellIds->GetId(v_id[i]);
1945 cellPts->GetPoint(v_id[i],v);
1946 // Currently all points are injected because of the possibility
1947 // of intersection point merging.
1948 if ( locator->InsertUniquePoint(v, iid[i]) )
1949 {
1950 outPD->CopyData(inPD,ptId, iid[i]);
1951 }
1952
1953 }//for all points of the tetrahedron.
1954
1955 if ( allInside )
1956 {
1957 // Tetrahedron inside.
1958 vtkIdType newCellId = tets->InsertNextCell(4,iid);
1959 outCD->CopyData(inCD,cellId,newCellId);
1960 continue;
1961 }
1962
1963 double *pedg1,*pedg2;
1964
1965
1966 // Tetrahedron Intersection Cases
1967 const unsigned int tab4[6][6] = { {1,1,0,3,3,2},
1968 {2,0,0,3,2,1},
1969 {3,3,2,0,2,1},
1970 {1,0,2,0,1,3},
1971 {0,0,1,2,3,3},
1972 {0,1,2,1,2,3}};
1973 const unsigned int tab3[4][6] = { {0,2,1,1,3,2},
1974 {0,1,2,0,2,3},
1975 {0,1,2,1,0,3},
1976 {0,1,2,0,1,2}};
1977 const unsigned int tab2[12][5] = { {0,0,1,2,3},
1978 {2,1,0,1,3},
1979 {1,0,1,0,3},
1980 {2,0,1,3,0},
1981 {3,1,0,1,0},
1982 {1,0,1,2,0},
1983 {3,1,0,2,1},
1984 {2,1,0,0,1},
1985 {0,0,1,3,1},
1986 {1,0,1,3,2},
1987 {3,1,0,0,2},
1988 {0,0,1,1,2}};
1989 const unsigned int tab1[12][3] = { {2,3,1},
1990 {3,2,0},
1991 {3,0,1},
1992 {0,3,2},
1993 {1,3,0},
1994 {3,1,2},
1995 {2,1,0},
1996 {1,2,3},
1997 {2,0,3},
1998 {0,2,1},
1999 {0,1,3},
2000 {1,0,2}};
2001
2002 vtkCellArray *cellarray = vtkCellArray::New();
2003 vtkIdType newCellId = cellarray->InsertNextCell(4,iid);
2004 unsigned int planes;
2005
2006 // Test Cell intersection with each plane of box
2007 for (planes = 0; planes < 6; planes++)
2008 {
2009 // The index of the dimension of the cut plane (x == 0, y == 1, z == 2).
2010 cutInd = planes/2;
2011
2012 // The plane is always parallel to unitary cube.
2013 value = this->BoundBoxClip[cutInd][planes%2];
2014
2015 unsigned int totalnewcells = cellarray->GetNumberOfCells();
2016 vtkCellArray *newcellArray = vtkCellArray::New();
2017 int edgeNum;
2018
2019 for (idcellnew = 0 ; idcellnew < totalnewcells; idcellnew++)
2020 {
2021 unsigned int num_inter = 0;
2022 unsigned int edges_inter = 0;
2023 unsigned int i0,i1;
2024 vtkIdType p_id[4];
2025 cellarray->GetNextCell(npts,v_id);
2026
2027 newPoints->GetPoint(v_id[0],v_tetra[0]); //coord (x,y,z)
2028 newPoints->GetPoint(v_id[1],v_tetra[1]);
2029 newPoints->GetPoint(v_id[2],v_tetra[2]);
2030 newPoints->GetPoint(v_id[3],v_tetra[3]);
2031
2032 for (edgeNum=0; edgeNum < 6; edgeNum++)
2033 {
2034 verts = edges[edgeNum];
2035
2036 p1 = v_tetra[verts[0]];
2037 p2 = v_tetra[verts[1]];
2038
2039 if ( (p1[cutInd] < value && value < p2[cutInd]) ||
2040 (p2[cutInd] < value && value < p1[cutInd]) )
2041 {
2042 deltaScalar = p2[cutInd] - p1[cutInd];
2043
2044 if (deltaScalar > 0)
2045 {
2046 pedg1 = p1; pedg2 = p2;
2047 v1 = verts[0]; v2 = verts[1];
2048 }
2049 else
2050 {
2051 pedg1 = p2; pedg2 = p1;
2052 v1 = verts[1]; v2 = verts[0];
2053 deltaScalar = -deltaScalar;
2054 }
2055
2056 // linear interpolation
2057 t = ( deltaScalar == 0.0 ? 0.0 :
2058 (value - pedg1[cutInd]) / deltaScalar );
2059
2060 for (j=0; j<3; j++)
2061 {
2062 x[j] = pedg1[j] + t*(pedg2[j] - pedg1[j]);
2063 }
2064
2065 // Incorporate point into output and interpolate edge data as necessary
2066 edges_inter = edges_inter * 10 + (edgeNum+1);
2067 if ( locator->InsertUniquePoint(x, p_id[num_inter]) )
2068 {
2069 this->InterpolateEdge(outPD, p_id[num_inter],
2070 v_id[v1], v_id[v2], t);
2071 }
2072 num_inter++;
2073 }//if edge intersects value
2074 }//for all edges
2075
2076 if (num_inter == 0)
2077 {
2078 unsigned int outside = 0;
2079 for(i=0; i<4; i++)
2080 {
2081 if (((v_tetra[i][cutInd] < value) && ((planes % 2) == 0)) ||
2082 ((v_tetra[i][cutInd] > value) && ((planes % 2) == 1)))
2083
2084 // If only one vertex is ouside, so the tetrahedron is outside
2085 // because there is not intersection.
2086 {
2087 outside = 1;
2088 break;
2089 }
2090 }
2091 if(outside == 0)
2092 {
2093 // else it is possible intersection if other plane
2094 newCellId = newcellArray->InsertNextCell(4,v_id);
2095 }
2096 continue;
2097 }
2098 switch(num_inter)
2099 {
2100 case 4: // We have two wedges
2101 switch(edges_inter)
2102 {
2103 case 1246:
2104 i0 = 0;
2105 break;
2106 case 2345:
2107 i0 = 2;
2108 break;
2109 case 1356:
2110 i0 = 4;
2111 break;
2112 default:
2113 vtkErrorMacro( << "Intersection not found: Num_inter = " <<
2114 num_inter << " Edges_inter = " << edges_inter );
2115 continue;
2116 }
2117
2118 if (((v_tetra[3][cutInd] < value) && ((planes % 2) == 0)) ||
2119 ((v_tetra[3][cutInd] > value) && ((planes % 2) == 1)))
2120 {
2121
2122 // The v_tetra[3] is outside box, so the first wedge is outside
2123
2124 tab_id[0] = p_id[tab4[i0+1][0]];
2125 // ps: v_tetra[3] is always in first wedge (see tab)
2126
2127 tab_id[1] = v_id[tab4[i0+1][1]];
2128 tab_id[2] = p_id[tab4[i0+1][2]];
2129 tab_id[3] = p_id[tab4[i0+1][3]];
2130 tab_id[4] = v_id[tab4[i0+1][4]];
2131 tab_id[5] = p_id[tab4[i0+1][5]];
2132 this->CreateTetra(6,tab_id,newcellArray);
2133 }
2134 else
2135 {
2136 tab_id[0] = p_id[tab4[i0][0]];
2137 tab_id[1] = v_id[tab4[i0][1]];
2138 tab_id[2] = p_id[tab4[i0][2]];
2139 tab_id[3] = p_id[tab4[i0][3]];
2140 tab_id[4] = v_id[tab4[i0][4]];
2141 tab_id[5] = p_id[tab4[i0][5]];
2142 this->CreateTetra(6,tab_id,newcellArray);
2143 }
2144 break;
2145 case 3: // We have one tetrahedron and one wedge
2146 // i0 gets the vertex on the tetrahedron.
2147 switch(edges_inter)
2148 {
2149 case 134:
2150 i0 = 0;
2151 break;
2152 case 125:
2153 i0 = 1;
2154 break;
2155 case 236:
2156 i0 = 2;
2157 break;
2158 case 456:
2159 i0 = 3;
2160 break;
2161 default:
2162 vtkErrorMacro( << "Intersection not found: Num_inter = " <<
2163 num_inter << " Edges_inter = " << edges_inter );
2164 continue;
2165 }
2166
2167 if (((v_tetra[i0][cutInd] < value) && ((planes % 2) == 0)) ||
2168 ((v_tetra[i0][cutInd] > value) && ((planes % 2) == 1)))
2169 {
2170
2171 // Isolate vertex is outside box, so the tetrahedron is outside
2172 tab_id[0] = p_id[tab3[i0][0]];
2173 tab_id[1] = p_id[tab3[i0][1]];
2174 tab_id[2] = p_id[tab3[i0][2]];
2175 tab_id[3] = v_id[tab3[i0][3]];
2176 tab_id[4] = v_id[tab3[i0][4]];
2177 tab_id[5] = v_id[tab3[i0][5]];
2178 this->CreateTetra(6,tab_id,newcellArray);
2179 }
2180 else
2181 {
2182 tab_id[0] = p_id[tab3[i0][0]];
2183 tab_id[1] = p_id[tab3[i0][1]];
2184 tab_id[2] = p_id[tab3[i0][2]];
2185 tab_id[3] = v_id[i0];
2186 newCellId = newcellArray->InsertNextCell(4,tab_id);
2187 }
2188 break;
2189
2190 case 2: // We have one tetrahedron and one pyramid
2191 switch(edges_inter) // i1 = vertex of the tetrahedron
2192 {
2193 case 12:
2194 i0 = 0; i1 = 1;
2195 break;
2196 case 13:
2197 i0 = 1; i1 = 0;
2198 break;
2199 case 23:
2200 i0 = 2; i1 = 2;
2201 break;
2202 case 25:
2203 i0 = 3; i1 = 1;
2204 break;
2205 case 26:
2206 i0 = 4; i1 = 2;
2207 break;
2208 case 56:
2209 i0 = 5; i1 = 3;
2210 break;
2211 case 34:
2212 i0 = 6; i1 = 0;
2213 break;
2214 case 46:
2215 i0 = 7; i1 = 3;
2216 break;
2217 case 36:
2218 i0 = 8; i1 = 2;
2219 break;
2220 case 14:
2221 i0 = 9; i1 = 0;
2222 break;
2223 case 15:
2224 i0 = 10; i1 = 1;
2225 break;
2226 case 45:
2227 i0 = 11; i1 = 3;
2228 break;
2229 default:
2230 vtkErrorMacro( << "Intersection not found: Num_inter = " <<
2231 num_inter << " Edges_inter = " << edges_inter );
2232 continue;
2233 }
2234 if (((v_tetra[i1][cutInd] < value) && ((planes % 2) == 0)) ||
2235 ((v_tetra[i1][cutInd] > value) && ((planes % 2) == 1)))
2236 {
2237 // Isolate vertex is outside box, so the tetrahedron is outside
2238 tab_id[0] = v_id[tab2[i0][0]];
2239 tab_id[1] = p_id[tab2[i0][1]];
2240 tab_id[2] = p_id[tab2[i0][2]];
2241 tab_id[3] = v_id[tab2[i0][3]];
2242 tab_id[4] = v_id[tab2[i0][4]];
2243 this->CreateTetra(5,tab_id,newcellArray);
2244 }
2245 else
2246 {
2247 tab_id[0] = v_id[i1];
2248 tab_id[1] = v_id[tab2[i0][4]];
2249 tab_id[2] = p_id[tab2[i0][2]];
2250 tab_id[3] = p_id[tab2[i0][1]];
2251 newCellId = newcellArray->InsertNextCell(4,tab_id);
2252 }
2253 break;
2254
2255 case 1: // We have two tetrahedron.
2256 if ((edges_inter > 6) || (edges_inter < 1))
2257 {
2258 vtkErrorMacro( << "Intersection not found: Num_inter = "
2259 << num_inter << " Edges_inter = " << edges_inter );
2260 continue;
2261 }
2262 if (((v_tetra[tab1[2*edges_inter-1][2]][cutInd] < value) && ((planes % 2) == 0)) ||
2263 ((v_tetra[tab1[2*edges_inter-1][2]][cutInd] > value) && ((planes % 2) == 1)))
2264 {
2265 // Isolate vertex is outside box, so the tetrahedron is outside
2266 tab_id[0] = p_id[0];
2267 tab_id[1] = v_id[tab1[2*edges_inter-2][0]];
2268 tab_id[2] = v_id[tab1[2*edges_inter-2][1]];
2269 tab_id[3] = v_id[tab1[2*edges_inter-2][2]];
2270 newCellId = newcellArray->InsertNextCell(4,tab_id);
2271 }
2272 else
2273 {
2274 tab_id[0] = p_id[0];
2275 tab_id[1] = v_id[tab1[2*edges_inter-1][0]];
2276 tab_id[2] = v_id[tab1[2*edges_inter-1][1]];
2277 tab_id[3] = v_id[tab1[2*edges_inter-1][2]];
2278 newCellId = newcellArray->InsertNextCell(4,tab_id);
2279 }
2280 break;
2281 }
2282 } // for all new cells
2283
2284 cellarray->Delete();
2285 cellarray = newcellArray;
2286 } //for all planes
2287
2288 unsigned int totalnewcells = cellarray->GetNumberOfCells();
2289
2290 for (idcellnew = 0 ; idcellnew < totalnewcells; idcellnew++)
2291 {
2292 cellarray->GetNextCell(npts,v_id);
2293 newCellId = tets->InsertNextCell(npts,v_id);
2294 outCD->CopyData(inCD,cellId,newCellId);
2295 }
2296 cellarray->Delete();
2297 }
2298 arraytetra->Delete();
2299 }
2300
2301 //----------------------------------------------------------------------------
2302 // ClipHexahedron: Box is like hexahedron.
2303 //
2304 // The difference between ClipBox and ClipHexahedron is the outside test.
2305 // The ClipHexahedron use plane equation to decide who is outside.
2306 //
ClipHexahedron(vtkPoints * newPoints,vtkGenericCell * cell,vtkIncrementalPointLocator * locator,vtkCellArray * tets,vtkPointData * inPD,vtkPointData * outPD,vtkCellData * inCD,vtkIdType cellId,vtkCellData * outCD)2307 void vtkBoxClipDataSet::ClipHexahedron(vtkPoints *newPoints,
2308 vtkGenericCell *cell,
2309 vtkIncrementalPointLocator *locator,
2310 vtkCellArray *tets,
2311 vtkPointData *inPD,
2312 vtkPointData *outPD,
2313 vtkCellData *inCD,
2314 vtkIdType cellId,
2315 vtkCellData *outCD)
2316 {
2317 vtkIdType idcellnew;
2318 vtkIdType cellType = cell->GetCellType();
2319 vtkIdList *cellIds = cell->GetPointIds();
2320 vtkCellArray *arraytetra = vtkCellArray::New();
2321 vtkPoints *cellPts = cell->GetPoints();
2322 vtkIdType npts = cellPts->GetNumberOfPoints();
2323 std::vector<vtkIdType> cellptId(npts);
2324 vtkIdType iid[4];
2325 vtkIdType *v_id = NULL;
2326 vtkIdType *verts, v1, v2;
2327 vtkIdType ptId;
2328 vtkIdType tab_id[6];
2329 vtkIdType ptstetra = 4;
2330
2331 vtkIdType i,j;
2332 unsigned int allInside, k;
2333 unsigned int planes;
2334
2335 vtkIdType edges[6][2] = { {0,1}, {1,2}, {2,0},
2336 {0,3}, {1,3}, {2,3} }; /* Edges Tetrahedron */
2337 double deltaScalar;
2338 double t;
2339 double *p1, *p2;
2340 double v[3], x[3];
2341 double p[6];
2342 double v_tetra[4][3];
2343
2344 for (i=0; i<npts; i++)
2345 {
2346 cellptId[i] = cellIds->GetId(i);
2347 }
2348
2349 this->CellGrid(cellType,npts,&cellptId[0],arraytetra); // Convert all volume cells to tetrahedra
2350
2351 unsigned int totalnewtetra = arraytetra->GetNumberOfCells();
2352 unsigned int idtetranew;
2353
2354 for (idtetranew = 0 ; idtetranew < totalnewtetra; idtetranew++)
2355 {
2356 arraytetra->GetNextCell(ptstetra,v_id);
2357
2358 for (allInside=1, i=0; i<4; i++)
2359 {
2360 cellPts->GetPoint(v_id[i],v);
2361
2362 for(k=0;k<6;k++)
2363 {
2364 p[k] = this->PlaneNormal[k][0]*(v[0] - this->PlanePoint[k][0])+
2365 this->PlaneNormal[k][1]*(v[1] - this->PlanePoint[k][1]) +
2366 this->PlaneNormal[k][2]*(v[2] - this->PlanePoint[k][2]);
2367 }
2368
2369 if (!((p[0] <= 0) && (p[1] <= 0) &&
2370 (p[2] <= 0) && (p[3] <= 0) &&
2371 (p[4] <= 0) && (p[5] <= 0)))
2372 {
2373 allInside = 0;
2374 }
2375 }//for all points of the cell.
2376
2377 // Test Outside
2378 unsigned int test[6] = {1,1,1,1,1,1};
2379 for (i=0; i<4; i++)
2380 {
2381 cellPts->GetPoint(v_id[i],v);
2382
2383 // Use plane equation
2384 for(k=0;k<6;k++)
2385 {
2386 p[k] = this->PlaneNormal[k][0]*(v[0] - this->PlanePoint[k][0])+
2387 this->PlaneNormal[k][1]*(v[1] - this->PlanePoint[k][1]) +
2388 this->PlaneNormal[k][2]*(v[2] - this->PlanePoint[k][2]);
2389 }
2390
2391
2392 for(k=0;k<3;k++)
2393 {
2394 if (p[2*k] < 0)
2395 {
2396 test[2*k] = 0;
2397 }
2398 if (p[2*k+1] < 0)
2399 {
2400 test[2*k+1] = 0;
2401 }
2402 }
2403
2404 }//for all points of the cell.
2405
2406 if (!allInside && ((test[0] == 1)|| (test[1] == 1) ||
2407 (test[2] == 1)|| (test[3] == 1) ||
2408 (test[4] == 1)|| (test[5] == 1)))
2409 {
2410 continue; // Tetrahedron is outside.
2411 }
2412
2413 for (i=0; i<4; i++)
2414 {
2415 ptId = cellIds->GetId(v_id[i]);
2416 cellPts->GetPoint(v_id[i],v);
2417
2418 // Currently all points are injected because of the possibility
2419 // of intersection point merging.
2420 if ( locator->InsertUniquePoint(v, iid[i]) )
2421 {
2422 outPD->CopyData(inPD,ptId, iid[i]);
2423 }
2424 }//for all points of the tetrahedron.
2425
2426 if ( allInside )
2427 {
2428 vtkIdType newCellId = tets->InsertNextCell(4,iid); // Tetrahedron inside.
2429 outCD->CopyData(inCD,cellId,newCellId);
2430 continue;
2431 }
2432
2433 double *pedg1,*pedg2;
2434
2435 // Tetrahedron Intersection Cases
2436 const unsigned int tab4[6][6] = { {1,1,0,3,3,2},
2437 {2,0,0,3,2,1},
2438 {3,3,2,0,2,1},
2439 {1,0,2,0,1,3},
2440 {0,0,1,2,3,3},
2441 {0,1,2,1,2,3}};
2442 const unsigned int tab3[4][6] = { {0,2,1,1,3,2},
2443 {0,1,2,0,2,3},
2444 {0,1,2,1,0,3},
2445 {0,1,2,0,1,2}};
2446 const unsigned int tab2[12][5] = { {0,0,1,2,3},
2447 {2,1,0,1,3},
2448 {1,0,1,0,3},
2449 {2,0,1,3,0},
2450 {3,1,0,1,0},
2451 {1,0,1,2,0},
2452 {3,1,0,2,1},
2453 {2,1,0,0,1},
2454 {0,0,1,3,1},
2455 {1,0,1,3,2},
2456 {3,1,0,0,2},
2457 {0,0,1,1,2}};
2458 const unsigned int tab1[12][3] = { {2,3,1},
2459 {3,2,0},
2460 {3,0,1},
2461 {0,3,2},
2462 {1,3,0},
2463 {3,1,2},
2464 {2,1,0},
2465 {1,2,3},
2466 {2,0,3},
2467 {0,2,1},
2468 {0,1,3},
2469 {1,0,2}};
2470
2471 vtkCellArray *cellarray = vtkCellArray::New();
2472 vtkIdType newCellId = cellarray->InsertNextCell(4,iid);
2473
2474 // Test Cell intersection with each plane of box
2475 for (planes = 0; planes < 6; planes++)
2476 {
2477 vtkIdType totalnewcells = cellarray->GetNumberOfCells();
2478 vtkCellArray *newcellArray = vtkCellArray::New();
2479
2480 for (idcellnew = 0 ; idcellnew < totalnewcells; idcellnew++)
2481 {
2482 unsigned int i0,i1;
2483 unsigned int num_inter = 0;
2484 unsigned int edges_inter = 0;
2485 vtkIdType p_id[4];
2486
2487 cellarray->GetNextCell(npts,v_id);
2488
2489 newPoints->GetPoint(v_id[0],v_tetra[0]); //coord (x,y,z)
2490 newPoints->GetPoint(v_id[1],v_tetra[1]);
2491 newPoints->GetPoint(v_id[2],v_tetra[2]);
2492 newPoints->GetPoint(v_id[3],v_tetra[3]);
2493
2494 p[0] = this->PlaneNormal[planes][0]*(v_tetra[0][0] - this->PlanePoint[planes][0]) +
2495 this->PlaneNormal[planes][1]*(v_tetra[0][1] - this->PlanePoint[planes][1]) +
2496 this->PlaneNormal[planes][2]*(v_tetra[0][2] - this->PlanePoint[planes][2]);
2497 p[1] = this->PlaneNormal[planes][0]*(v_tetra[1][0] - this->PlanePoint[planes][0]) +
2498 this->PlaneNormal[planes][1]*(v_tetra[1][1] - this->PlanePoint[planes][1]) +
2499 this->PlaneNormal[planes][2]*(v_tetra[1][2] - this->PlanePoint[planes][2]);
2500 p[2] = this->PlaneNormal[planes][0]*(v_tetra[2][0] - this->PlanePoint[planes][0]) +
2501 this->PlaneNormal[planes][1]*(v_tetra[2][1] - this->PlanePoint[planes][1]) +
2502 this->PlaneNormal[planes][2]*(v_tetra[2][2] - this->PlanePoint[planes][2]);
2503 p[3] = this->PlaneNormal[planes][0]*(v_tetra[3][0] - this->PlanePoint[planes][0]) +
2504 this->PlaneNormal[planes][1]*(v_tetra[3][1] - this->PlanePoint[planes][1]) +
2505 this->PlaneNormal[planes][2]*(v_tetra[3][2] - this->PlanePoint[planes][2]);
2506
2507 for (int edgeNum=0; edgeNum < 6; edgeNum++)
2508 {
2509 verts = edges[edgeNum];
2510
2511 p1 = v_tetra[verts[0]];
2512 p2 = v_tetra[verts[1]];
2513 double s1 = p[verts[0]];
2514 double s2 = p[verts[1]];
2515 if ( (s1 * s2) < 0)
2516 {
2517 deltaScalar = s2 - s1;
2518
2519 if (deltaScalar > 0)
2520 {
2521 pedg1 = p1; pedg2 = p2;
2522 v1 = verts[0]; v2 = verts[1];
2523 }
2524 else
2525 {
2526 pedg1 = p2; pedg2 = p1;
2527 v1 = verts[1]; v2 = verts[0];
2528 deltaScalar = -deltaScalar;
2529 t = s1; s1 = s2; s2 = t;
2530 }
2531
2532 // linear interpolation
2533 t = ( deltaScalar == 0.0 ? 0.0 : ( - s1) / deltaScalar );
2534
2535 for (j=0; j<3; j++)
2536 {
2537 x[j] = pedg1[j] + t*(pedg2[j] - pedg1[j]);
2538 }
2539
2540 // Incorporate point into output and interpolate edge data as necessary
2541 edges_inter = edges_inter * 10 + (edgeNum+1);
2542
2543 if ( locator->InsertUniquePoint(x, p_id[num_inter]) )
2544 {
2545 this->InterpolateEdge(outPD, p_id[num_inter],
2546 v_id[v1], v_id[v2], t);
2547 }
2548
2549 num_inter++;
2550 }//if edge intersects value
2551 }//for all edges
2552 if (num_inter == 0)
2553 {
2554 unsigned int outside = 0;
2555 for(i=0;i<4;i++)
2556 {
2557 if (p[i] > 0)
2558 {
2559 // If only one vertex is ouside, so the tetrahedron is outside
2560 // because there is not intersection.
2561 // some vertex could be on plane, so you need to test all vertex
2562
2563 outside = 1;
2564 break;
2565 }
2566 }
2567 if (outside == 0)
2568 {
2569 // else it is possible intersection if other plane
2570
2571 newCellId = newcellArray->InsertNextCell(4,v_id);
2572 }
2573 continue;
2574 }
2575 switch(num_inter)
2576 {
2577 case 4: // We have two wedges
2578 switch(edges_inter)
2579 {
2580 case 1246:
2581 i0 = 0;
2582 break;
2583 case 2345:
2584 i0 = 2;
2585 break;
2586 case 1356:
2587 i0 = 4;
2588 break;
2589 default:
2590 vtkErrorMacro( << "Intersection not found: Num_inter = " <<
2591 num_inter << " Edges_inter = " << edges_inter );
2592 continue;
2593 }
2594
2595 if (p[3] > 0)
2596 {
2597 // The v_tetra[3] is outside box, so the first wedge is outside
2598 // ps: v_tetra[3] is always in first wedge (see tab)
2599
2600 tab_id[0] = p_id[tab4[i0+1][0]];
2601 tab_id[1] = v_id[tab4[i0+1][1]];
2602 tab_id[2] = p_id[tab4[i0+1][2]];
2603 tab_id[3] = p_id[tab4[i0+1][3]];
2604 tab_id[4] = v_id[tab4[i0+1][4]];
2605 tab_id[5] = p_id[tab4[i0+1][5]];
2606 this->CreateTetra(6,tab_id,newcellArray);
2607 }
2608 else
2609 {
2610 tab_id[0] = p_id[tab4[i0][0]];
2611 tab_id[1] = v_id[tab4[i0][1]];
2612 tab_id[2] = p_id[tab4[i0][2]];
2613 tab_id[3] = p_id[tab4[i0][3]];
2614 tab_id[4] = v_id[tab4[i0][4]];
2615 tab_id[5] = p_id[tab4[i0][5]];
2616 this->CreateTetra(6,tab_id,newcellArray);
2617 }
2618 break;
2619 case 3: // We have one tetrahedron and one wedge
2620 switch(edges_inter)
2621 {
2622 case 134:
2623 i0 = 0;
2624 break;
2625 case 125:
2626 i0 = 1;
2627 break;
2628 case 236:
2629 i0 = 2;
2630 break;
2631 case 456:
2632 i0 = 3;
2633 break;
2634 default:
2635 vtkErrorMacro( << "Intersection not found: Num_inter = " <<
2636 num_inter << " Edges_inter = " << edges_inter );
2637 continue;
2638 }
2639
2640 if (p[i0] > 0)
2641 {
2642 // Isolate vertex is outside box, so the tetrahedron is outside
2643
2644 tab_id[0] = p_id[tab3[i0][0]];
2645 tab_id[1] = p_id[tab3[i0][1]];
2646 tab_id[2] = p_id[tab3[i0][2]];
2647 tab_id[3] = v_id[tab3[i0][3]];
2648 tab_id[4] = v_id[tab3[i0][4]];
2649 tab_id[5] = v_id[tab3[i0][5]];
2650 this->CreateTetra(6,tab_id,newcellArray);
2651 }
2652 else
2653 {
2654 tab_id[0] = p_id[tab3[i0][0]];
2655 tab_id[1] = p_id[tab3[i0][1]];
2656 tab_id[2] = p_id[tab3[i0][2]];
2657 tab_id[3] = v_id[i0];
2658 newCellId = newcellArray->InsertNextCell(4,tab_id);
2659 }
2660 break;
2661 case 2: // We have one tetrahedron and one pyramid
2662 switch(edges_inter) // i1 = vertex of the tetrahedron
2663 {
2664 case 12:
2665 i0 = 0; i1 = 1;
2666 break;
2667 case 13:
2668 i0 = 1; i1 = 0;
2669 break;
2670 case 23:
2671 i0 = 2; i1 = 2;
2672 break;
2673 case 25:
2674 i0 = 3; i1 = 1;
2675 break;
2676 case 26:
2677 i0 = 4; i1 = 2;
2678 break;
2679 case 56:
2680 i0 = 5; i1 = 3;
2681 break;
2682 case 34:
2683 i0 = 6; i1 = 0;
2684 break;
2685 case 46:
2686 i0 = 7; i1 = 3;
2687 break;
2688 case 36:
2689 i0 = 8; i1 = 2;
2690 break;
2691 case 14:
2692 i0 = 9; i1 = 0;
2693 break;
2694 case 15:
2695 i0 = 10; i1 = 1;
2696 break;
2697 case 45:
2698 i0 = 11; i1 = 3;
2699 break;
2700 default:
2701 vtkErrorMacro( << "Intersection not found: Num_inter = " <<
2702 num_inter << " Edges_inter = " << edges_inter );
2703 continue;
2704 }
2705 if (p[i1] > 0)
2706 {
2707 // Isolate vertex is outside box, so the tetrahedron is outside
2708
2709 tab_id[0] = v_id[tab2[i0][0]];
2710 tab_id[1] = p_id[tab2[i0][1]];
2711 tab_id[2] = p_id[tab2[i0][2]];
2712 tab_id[3] = v_id[tab2[i0][3]];
2713 tab_id[4] = v_id[tab2[i0][4]];
2714 this->CreateTetra(5,tab_id,newcellArray);
2715 }
2716 else
2717 {
2718 tab_id[0] = v_id[i1];
2719 tab_id[1] = v_id[tab2[i0][4]];
2720 tab_id[2] = p_id[tab2[i0][2]];
2721 tab_id[3] = p_id[tab2[i0][1]];
2722 newCellId = newcellArray->InsertNextCell(4,tab_id);
2723 }
2724 break;
2725 case 1: // We have two tetrahedron.
2726 if ((edges_inter > 6) || (edges_inter < 1))
2727 {
2728 vtkErrorMacro( << "Intersection not found: Num_inter = " <<
2729 num_inter << " Edges_inter = " << edges_inter );
2730 continue;
2731 }
2732 if (p[tab1[2*edges_inter-1][2]] > 0)
2733 {
2734 // Isolate vertex is outside box, so the tetrahedron is outside
2735 tab_id[0] = p_id[0];
2736 tab_id[1] = v_id[tab1[2*edges_inter-2][0]];
2737 tab_id[2] = v_id[tab1[2*edges_inter-2][1]];
2738 tab_id[3] = v_id[tab1[2*edges_inter-2][2]];
2739 newCellId = newcellArray->InsertNextCell(4,tab_id);
2740 }
2741 else
2742 {
2743 tab_id[0] = p_id[0];
2744 tab_id[1] = v_id[tab1[2*edges_inter-1][0]];
2745 tab_id[2] = v_id[tab1[2*edges_inter-1][1]];
2746 tab_id[3] = v_id[tab1[2*edges_inter-1][2]];
2747 newCellId = newcellArray->InsertNextCell(4,tab_id);
2748 }
2749 break;
2750 }
2751 } // for all new cells
2752 cellarray->Delete();
2753 cellarray = newcellArray;
2754 } //for all planes
2755
2756 vtkIdType totalnewcells = cellarray->GetNumberOfCells();
2757
2758 for (idcellnew = 0 ; idcellnew < totalnewcells; idcellnew++)
2759 {
2760 cellarray->GetNextCell(npts,v_id);
2761 newCellId = tets->InsertNextCell(npts,v_id);
2762 outCD->CopyData(inCD,cellId,newCellId);
2763 }
2764 cellarray->Delete();
2765 }
2766 arraytetra->Delete();
2767 }
2768 //----------------------------------------------------------------------------
2769 // ClipBoxInOut
2770 //
2771 // The difference between ClipBox and ClipBoxInOut is the outputs.
2772 // The ClipBoxInOut generate both outputs: inside and outside the clip box.
2773 //
ClipBoxInOut(vtkPoints * newPoints,vtkGenericCell * cell,vtkIncrementalPointLocator * locator,vtkCellArray ** tets,vtkPointData * inPD,vtkPointData * outPD,vtkCellData * inCD,vtkIdType cellId,vtkCellData ** outCD)2774 void vtkBoxClipDataSet::ClipBoxInOut(vtkPoints *newPoints,
2775 vtkGenericCell *cell,
2776 vtkIncrementalPointLocator *locator,
2777 vtkCellArray **tets,
2778 vtkPointData *inPD,
2779 vtkPointData *outPD,
2780 vtkCellData *inCD,
2781 vtkIdType cellId,
2782 vtkCellData **outCD)
2783 {
2784 vtkIdType cellType = cell->GetCellType();
2785 vtkIdList *cellIds = cell->GetPointIds();
2786 vtkCellArray *arraytetra = vtkCellArray::New();
2787 vtkPoints *cellPts = cell->GetPoints();
2788 vtkIdType npts = cellPts->GetNumberOfPoints();
2789 std::vector<vtkIdType> cellptId(npts);
2790 vtkIdType iid[4];
2791 vtkIdType *v_id = NULL;
2792 vtkIdType *verts, v1, v2;
2793 vtkIdType ptId;
2794 vtkIdType ptIdout[4];
2795 vtkIdType tab_id[6];
2796 vtkIdType ptstetra = 4;
2797
2798 int i,j;
2799 int allInside;
2800 int cutInd;
2801
2802 unsigned int planes;
2803 unsigned int idcellnew;
2804 unsigned int idtetranew;
2805
2806 vtkIdType edges[6][2] = { {0,1}, {1,2}, {2,0},
2807 {0,3}, {1,3}, {2,3} }; /* Edges Tetrahedron */
2808 double value,deltaScalar;
2809 double t;
2810 double v[3], x[3];
2811 double v_tetra[4][3];
2812 double *p1, *p2;
2813
2814 for (i=0; i<npts; i++)
2815 {
2816 cellptId[i] = cellIds->GetId(i);
2817 }
2818
2819 // Convert all volume cells to tetrahedra
2820 this->CellGrid(cellType,npts,&cellptId[0],arraytetra);
2821 unsigned int totalnewtetra = arraytetra->GetNumberOfCells();
2822
2823 for (idtetranew = 0 ; idtetranew < totalnewtetra; idtetranew++)
2824 {
2825 arraytetra->GetNextCell(ptstetra,v_id);
2826
2827 for (allInside=1, i=0; i<4; i++)
2828 {
2829 cellPts->GetPoint(v_id[i],v);
2830
2831 if (!(((v[0] >= this->BoundBoxClip[0][0])&&(v[0] <= this->BoundBoxClip[0][1]) &&
2832 (v[1] >= this->BoundBoxClip[1][0])&&(v[1] <= this->BoundBoxClip[1][1])&&
2833 (v[2] >= this->BoundBoxClip[2][0])&&(v[2] <= this->BoundBoxClip[2][1]))))
2834 {
2835 //outside,its type might change later (nearby intersection)
2836 allInside = 0;
2837 }
2838 }//for all points of the cell.
2839
2840 // Test Outside: see(1)
2841 if (!allInside)
2842 {
2843 unsigned int test[6] = {1,1,1,1,1,1};
2844 for (i=0; i<4; i++)
2845 {
2846 ptIdout[i] = cellIds->GetId(v_id[i]);
2847 cellPts->GetPoint(v_id[i],v_tetra[i]);
2848
2849 if (v_tetra[i][0] > this->BoundBoxClip[0][0])
2850 {
2851 test[0] = 0;
2852 }
2853 if (v_tetra[i][0] < this->BoundBoxClip[0][1])
2854 {
2855 test[1] = 0;
2856 }
2857 if (v_tetra[i][1] > this->BoundBoxClip[1][0])
2858 {
2859 test[2] = 0;
2860 }
2861 if (v_tetra[i][1] < this->BoundBoxClip[1][1])
2862 {
2863 test[3] = 0;
2864 }
2865 if (v_tetra[i][2] > this->BoundBoxClip[2][0])
2866 {
2867 test[4] = 0;
2868 }
2869 if (v_tetra[i][2] < this->BoundBoxClip[2][1])
2870 {
2871 test[5] = 0;
2872 }
2873 }//for all points of the cell.
2874
2875 if ((test[0] == 1)|| (test[1] == 1) ||
2876 (test[2] == 1)|| (test[3] == 1) ||
2877 (test[4] == 1)|| (test[5] == 1))
2878 {
2879 for (i=0; i<4; i++)
2880 {
2881 if ( locator->InsertUniquePoint(v_tetra[i], iid[i]) )
2882 {
2883 outPD->CopyData(inPD,ptIdout[i], iid[i]);
2884 }
2885 }
2886 int newCellId = tets[1]->InsertNextCell(4,iid);
2887 outCD[1]->CopyData(inCD,cellId,newCellId);
2888 continue; // Tetrahedron is outside.
2889 }
2890 }//if not allinside.
2891
2892 for (i=0; i<4; i++)
2893 {
2894 ptId = cellIds->GetId(v_id[i]);
2895 cellPts->GetPoint(v_id[i],v);
2896
2897 // Currently all points are injected because of the possibility
2898 // of intersection point merging.
2899 if ( locator->InsertUniquePoint(v, iid[i]) )
2900 {
2901 outPD->CopyData(inPD,ptId, iid[i]);
2902 }
2903
2904 }//for all points of the tetrahedron.
2905
2906 if ( allInside )
2907 {
2908 // Tetrahedron inside.
2909 int newCellId = tets[0]->InsertNextCell(4,iid);
2910 outCD[0]->CopyData(inCD,cellId,newCellId);
2911 continue;
2912 }
2913
2914 double *pedg1,*pedg2;
2915
2916 // Tetrahedron Intersection Cases
2917 const unsigned int tab4[6][6] = { {1,1,0,3,3,2},
2918 {2,0,0,3,2,1},
2919 {3,3,2,0,2,1},
2920 {1,0,2,0,1,3},
2921 {0,0,1,2,3,3},
2922 {0,1,2,1,2,3}};
2923 const unsigned int tab3[4][6] = { {0,2,1,1,3,2},
2924 {0,1,2,0,2,3},
2925 {0,1,2,1,0,3},
2926 {0,1,2,0,1,2}};
2927 const unsigned int tab2[12][5] = { {0,0,1,2,3},
2928 {2,1,0,1,3},
2929 {1,0,1,0,3},
2930 {2,0,1,3,0},
2931 {3,1,0,1,0},
2932 {1,0,1,2,0},
2933 {3,1,0,2,1},
2934 {2,1,0,0,1},
2935 {0,0,1,3,1},
2936 {1,0,1,3,2},
2937 {3,1,0,0,2},
2938 {0,0,1,1,2}};
2939 const unsigned int tab1[12][3] = { {2,3,1},
2940 {3,2,0},
2941 {3,0,1},
2942 {0,3,2},
2943 {1,3,0},
2944 {3,1,2},
2945 {2,1,0},
2946 {1,2,3},
2947 {2,0,3},
2948 {0,2,1},
2949 {0,1,3},
2950 {1,0,2}};
2951
2952 vtkCellArray *cellarray = vtkCellArray::New();
2953 vtkCellArray *cellarrayout = vtkCellArray::New();
2954 int newCellId = cellarray->InsertNextCell(4,iid);
2955
2956 // Test Cell intersection with each plane of box
2957 for (planes = 0; planes < 6; planes++)
2958 {
2959 // The index of the dimension of the cut plane (x == 0, y == 1, z == 2).
2960 cutInd = planes/2;
2961
2962 value = this->BoundBoxClip[cutInd][planes%2]; // The plane is always parallel to unitary cube.
2963
2964 unsigned int totalnewcells = cellarray->GetNumberOfCells();
2965 vtkCellArray *newcellArray = vtkCellArray::New();
2966
2967 for (idcellnew = 0 ; idcellnew < totalnewcells; idcellnew++)
2968 {
2969 unsigned int num_inter = 0;
2970 unsigned int edges_inter = 0;
2971 unsigned int i0,i1;
2972 vtkIdType p_id[4];
2973 cellarray->GetNextCell(npts,v_id);
2974
2975 newPoints->GetPoint(v_id[0],v_tetra[0]); //coord (x,y,z)
2976 newPoints->GetPoint(v_id[1],v_tetra[1]);
2977 newPoints->GetPoint(v_id[2],v_tetra[2]);
2978 newPoints->GetPoint(v_id[3],v_tetra[3]);
2979 for (int edgeNum=0; edgeNum < 6; edgeNum++)
2980 {
2981 verts = edges[edgeNum];
2982
2983 p1 = v_tetra[verts[0]];
2984 p2 = v_tetra[verts[1]];
2985
2986 if ( (p1[cutInd] < value && value < p2[cutInd]) ||
2987 (p2[cutInd] < value && value < p1[cutInd]) )
2988 {
2989 deltaScalar = p2[cutInd] - p1[cutInd];
2990
2991 if (deltaScalar > 0)
2992 {
2993 pedg1 = p1; pedg2 = p2;
2994 v1 = verts[0]; v2 = verts[1];
2995 }
2996 else
2997 {
2998 pedg1 = p2; pedg2 = p1;
2999 v1 = verts[1]; v2 = verts[0];
3000 deltaScalar = -deltaScalar;
3001 }
3002
3003 // linear interpolation
3004 t = ( deltaScalar == 0.0 ? 0.0 :
3005 (value - pedg1[cutInd]) / deltaScalar );
3006
3007 for (j=0; j<3; j++)
3008 {
3009 x[j] = pedg1[j] + t*(pedg2[j] - pedg1[j]);
3010 }
3011
3012 // Incorporate point into output and interpolate edge data as necessary
3013 edges_inter = edges_inter * 10 + (edgeNum+1);
3014 if ( locator->InsertUniquePoint(x, p_id[num_inter]) )
3015 {
3016 this->InterpolateEdge(outPD, p_id[num_inter],
3017 v_id[v1], v_id[v2], t);
3018 }
3019
3020 num_inter++;
3021 }//if edge intersects value
3022 }//for all edges
3023
3024 if (num_inter == 0)
3025 {
3026 unsigned int outside = 0;
3027 for(i=0; i<4; i++)
3028 {
3029 if (((v_tetra[i][cutInd] < value) && ((planes % 2) == 0)) ||
3030 ((v_tetra[i][cutInd] > value) && ((planes % 2) == 1)))
3031 {
3032 // If only one vertex is ouside, so the tetrahedron is outside
3033 // because there is not intersection.
3034 outside = 1;
3035 break;
3036 }
3037 }
3038 if(outside == 0)
3039 {
3040 // else it is possible intersection if other plane
3041 newCellId = newcellArray->InsertNextCell(4,v_id);
3042 }
3043 else
3044 {
3045 newCellId = tets[1]->InsertNextCell(4,v_id);
3046 outCD[1]->CopyData(inCD,cellId,newCellId);
3047 }
3048 continue;
3049 }
3050 switch(num_inter)
3051 {
3052 case 4: // We have two wedges
3053 switch(edges_inter)
3054 {
3055 case 1246:
3056 i0 = 0;
3057 break;
3058 case 2345:
3059 i0 = 2;
3060 break;
3061 case 1356:
3062 i0 = 4;
3063 break;
3064 default:
3065 vtkErrorMacro( << "Intersection not found: Num_inter = " <<
3066 num_inter << " Edges_inter = " << edges_inter );
3067 continue;
3068 }
3069 if (((v_tetra[3][cutInd] < value) && ((planes % 2) == 0)) ||
3070 ((v_tetra[3][cutInd] > value) && ((planes % 2) == 1)))
3071 {
3072 // The v_tetra[3] is outside box, so
3073 // the first wedge is outside
3074 // ps: v_tetra[3] is always in first wedge (see tab)
3075
3076 tab_id[0] = p_id[tab4[i0+1][0]]; // Inside
3077 tab_id[1] = v_id[tab4[i0+1][1]];
3078 tab_id[2] = p_id[tab4[i0+1][2]];
3079 tab_id[3] = p_id[tab4[i0+1][3]];
3080 tab_id[4] = v_id[tab4[i0+1][4]];
3081 tab_id[5] = p_id[tab4[i0+1][5]];
3082 this->CreateTetra(6,tab_id,newcellArray);
3083
3084 tab_id[0] = p_id[tab4[i0][0]]; // Outside
3085 tab_id[1] = v_id[tab4[i0][1]];
3086 tab_id[2] = p_id[tab4[i0][2]];
3087 tab_id[3] = p_id[tab4[i0][3]];
3088 tab_id[4] = v_id[tab4[i0][4]];
3089 tab_id[5] = p_id[tab4[i0][5]];
3090 this->CreateTetra(6,tab_id,cellarrayout);
3091 }
3092 else
3093 {
3094 tab_id[0] = p_id[tab4[i0][0]]; // Inside
3095 tab_id[1] = v_id[tab4[i0][1]];
3096 tab_id[2] = p_id[tab4[i0][2]];
3097 tab_id[3] = p_id[tab4[i0][3]];
3098 tab_id[4] = v_id[tab4[i0][4]];
3099 tab_id[5] = p_id[tab4[i0][5]];
3100 this->CreateTetra(6,tab_id,newcellArray);
3101
3102 tab_id[0] = p_id[tab4[i0+1][0]]; // Outside
3103 tab_id[1] = v_id[tab4[i0+1][1]];
3104 tab_id[2] = p_id[tab4[i0+1][2]];
3105 tab_id[3] = p_id[tab4[i0+1][3]];
3106 tab_id[4] = v_id[tab4[i0+1][4]];
3107 tab_id[5] = p_id[tab4[i0+1][5]];
3108 this->CreateTetra(6,tab_id,cellarrayout);
3109 }
3110 break;
3111 case 3: // We have one tetrahedron and one wedge
3112 switch(edges_inter)
3113 {
3114 case 134:
3115 i0 = 0;
3116 break;
3117 case 125:
3118 i0 = 1;
3119 break;
3120 case 236:
3121 i0 = 2;
3122 break;
3123 case 456:
3124 i0 = 3;
3125 break;
3126 default:
3127 vtkErrorMacro( << "Intersection not found: Num_inter = " <<
3128 num_inter << " Edges_inter = " << edges_inter );
3129 continue;
3130 }
3131 if (((v_tetra[i0][cutInd] < value) && ((planes % 2) == 0)) ||
3132 ((v_tetra[i0][cutInd] > value) && ((planes % 2) == 1)))
3133 {
3134 // Isolate vertex is outside box, so/ the tetrahedron is outside
3135 tab_id[0] = p_id[tab3[i0][0]]; // Inside
3136 tab_id[1] = p_id[tab3[i0][1]];
3137 tab_id[2] = p_id[tab3[i0][2]];
3138 tab_id[3] = v_id[tab3[i0][3]];
3139 tab_id[4] = v_id[tab3[i0][4]];
3140 tab_id[5] = v_id[tab3[i0][5]];
3141 this->CreateTetra(6,tab_id,newcellArray);
3142
3143 tab_id[0] = p_id[tab3[i0][0]]; // Outside
3144 tab_id[1] = p_id[tab3[i0][1]];
3145 tab_id[2] = p_id[tab3[i0][2]];
3146 tab_id[3] = v_id[i0];
3147 newCellId = cellarrayout->InsertNextCell(4,tab_id);
3148 }
3149 else
3150 {
3151 tab_id[0] = p_id[tab3[i0][0]];
3152 tab_id[1] = p_id[tab3[i0][1]];
3153 tab_id[2] = p_id[tab3[i0][2]]; // Inside
3154 tab_id[3] = v_id[i0];
3155 newCellId = newcellArray->InsertNextCell(4,tab_id);
3156
3157 tab_id[0] = p_id[tab3[i0][0]]; // Outside
3158 tab_id[1] = p_id[tab3[i0][1]];
3159 tab_id[2] = p_id[tab3[i0][2]];
3160 tab_id[3] = v_id[tab3[i0][3]];
3161 tab_id[4] = v_id[tab3[i0][4]];
3162 tab_id[5] = v_id[tab3[i0][5]];
3163 this->CreateTetra(6,tab_id,cellarrayout);
3164 }
3165 break;
3166 case 2: // We have one tetrahedron and one pyramid
3167 switch(edges_inter) // i1 = vertex of the tetrahedron
3168 {
3169 case 12:
3170 i0 = 0; i1 = 1;
3171 break;
3172 case 13:
3173 i0 = 1; i1 = 0;
3174 break;
3175 case 23:
3176 i0 = 2; i1 = 2;
3177 break;
3178 case 25:
3179 i0 = 3; i1 = 1;
3180 break;
3181 case 26:
3182 i0 = 4; i1 = 2;
3183 break;
3184 case 56:
3185 i0 = 5; i1 = 3;
3186 break;
3187 case 34:
3188 i0 = 6; i1 = 0;
3189 break;
3190 case 46:
3191 i0 = 7; i1 = 3;
3192 break;
3193 case 36:
3194 i0 = 8; i1 = 2;
3195 break;
3196 case 14:
3197 i0 = 9; i1 = 0;
3198 break;
3199 case 15:
3200 i0 = 10; i1 = 1;
3201 break;
3202 case 45:
3203 i0 = 11; i1 = 3;
3204 break;
3205 default:
3206 vtkErrorMacro( << "Intersection not found: Num_inter = " <<
3207 num_inter << " Edges_inter = " << edges_inter );
3208 continue;
3209 }
3210 if (((v_tetra[i1][cutInd] < value) && ((planes % 2) == 0)) ||
3211 ((v_tetra[i1][cutInd] > value) && ((planes % 2) == 1)))
3212 {
3213 // Isolate vertex is outside box, so the tetrahedron is outside
3214 tab_id[0] = v_id[tab2[i0][0]]; // Inside
3215 tab_id[1] = p_id[tab2[i0][1]];
3216 tab_id[2] = p_id[tab2[i0][2]];
3217 tab_id[3] = v_id[tab2[i0][3]];
3218 tab_id[4] = v_id[tab2[i0][4]];
3219 this->CreateTetra(5,tab_id,newcellArray);
3220
3221 tab_id[0] = v_id[i1]; // Outside
3222 tab_id[1] = v_id[tab2[i0][4]];
3223 tab_id[2] = p_id[tab2[i0][2]];
3224 tab_id[3] = p_id[tab2[i0][1]];
3225 newCellId = cellarrayout->InsertNextCell(4,tab_id);
3226 }
3227 else
3228 {
3229 tab_id[0] = v_id[i1]; // Inside
3230 tab_id[1] = v_id[tab2[i0][4]];
3231 tab_id[2] = p_id[tab2[i0][2]];
3232 tab_id[3] = p_id[tab2[i0][1]];
3233 newCellId = newcellArray->InsertNextCell(4,tab_id);
3234
3235 tab_id[0] = v_id[tab2[i0][0]]; // Outside
3236 tab_id[1] = p_id[tab2[i0][1]];
3237 tab_id[2] = p_id[tab2[i0][2]];
3238 tab_id[3] = v_id[tab2[i0][3]];
3239 tab_id[4] = v_id[tab2[i0][4]];
3240 this->CreateTetra(5,tab_id,cellarrayout);
3241 }
3242 break;
3243 case 1: // We have two tetrahedron.
3244 if ((edges_inter > 6) || (edges_inter < 1))
3245 {
3246 vtkErrorMacro( << "Intersection not found: Num_inter = " <<
3247 num_inter << " Edges_inter = " << edges_inter );
3248 continue;
3249 }
3250 if (((v_tetra[tab1[2*edges_inter-1][2]][cutInd] < value) && ((planes % 2) == 0)) ||
3251 ((v_tetra[tab1[2*edges_inter-1][2]][cutInd] > value) && ((planes % 2) == 1)))
3252 {
3253 // Isolate vertex is outside box, so the tetrahedron is outside
3254
3255 tab_id[0] = p_id[0]; // Inside
3256 tab_id[1] = v_id[tab1[2*edges_inter-2][0]];
3257 tab_id[2] = v_id[tab1[2*edges_inter-2][1]];
3258 tab_id[3] = v_id[tab1[2*edges_inter-2][2]];
3259 newCellId = newcellArray->InsertNextCell(4,tab_id);
3260
3261 tab_id[0] = p_id[0]; // Outside
3262 tab_id[1] = v_id[tab1[2*edges_inter-1][0]];
3263 tab_id[2] = v_id[tab1[2*edges_inter-1][1]];
3264 tab_id[3] = v_id[tab1[2*edges_inter-1][2]];
3265 newCellId = cellarrayout->InsertNextCell(4,tab_id);
3266 }
3267 else
3268 {
3269 tab_id[0] = p_id[0]; // Inside
3270 tab_id[1] = v_id[tab1[2*edges_inter-1][0]];
3271 tab_id[2] = v_id[tab1[2*edges_inter-1][1]];
3272 tab_id[3] = v_id[tab1[2*edges_inter-1][2]];
3273 newCellId = newcellArray->InsertNextCell(4,tab_id);
3274
3275 tab_id[0] = p_id[0]; // Outside
3276 tab_id[1] = v_id[tab1[2*edges_inter-2][0]];
3277 tab_id[2] = v_id[tab1[2*edges_inter-2][1]];
3278 tab_id[3] = v_id[tab1[2*edges_inter-2][2]];
3279 newCellId = cellarrayout->InsertNextCell(4,tab_id);
3280 }
3281 break;
3282 }
3283 } // for all new cells
3284 cellarray->Delete();
3285 cellarray = newcellArray;
3286 } //for all planes
3287
3288 unsigned int totalnewcells = cellarray->GetNumberOfCells(); // Inside
3289 for (idcellnew = 0 ; idcellnew < totalnewcells; idcellnew++)
3290 {
3291 cellarray->GetNextCell(npts,v_id);
3292 newCellId = tets[0]->InsertNextCell(npts,v_id);
3293 outCD[0]->CopyData(inCD,cellId,newCellId);
3294 }
3295 cellarray->Delete();
3296
3297 totalnewcells = cellarrayout->GetNumberOfCells(); // Outside
3298
3299 for (idcellnew = 0 ; idcellnew < totalnewcells; idcellnew++)
3300 {
3301 cellarrayout->GetNextCell(npts,v_id);
3302 newCellId = tets[1]->InsertNextCell(npts,v_id);
3303 outCD[1]->CopyData(inCD,cellId,newCellId);
3304 }
3305 cellarrayout->Delete();
3306 }
3307 arraytetra->Delete();
3308 }
3309
3310
3311 //----------------------------------------------------------------------------
3312 // ClipHexahedronInOut
3313 //
3314 // The difference between ClipHexahedron and ClipHexahedronInOut is the outputs.
3315 // The ClipHexahedronInOut generate both outputs: inside and outside the clip hexahedron.
3316 //
ClipHexahedronInOut(vtkPoints * newPoints,vtkGenericCell * cell,vtkIncrementalPointLocator * locator,vtkCellArray ** tets,vtkPointData * inPD,vtkPointData * outPD,vtkCellData * inCD,vtkIdType cellId,vtkCellData ** outCD)3317 void vtkBoxClipDataSet::ClipHexahedronInOut(vtkPoints *newPoints,
3318 vtkGenericCell *cell,
3319 vtkIncrementalPointLocator *locator,
3320 vtkCellArray **tets,
3321 vtkPointData *inPD,
3322 vtkPointData *outPD,
3323 vtkCellData *inCD,
3324 vtkIdType cellId,
3325 vtkCellData **outCD)
3326 {
3327 vtkIdType cellType = cell->GetCellType();
3328 vtkIdList *cellIds = cell->GetPointIds();
3329 vtkCellArray *arraytetra = vtkCellArray::New();
3330 vtkPoints *cellPts = cell->GetPoints();
3331 vtkIdType npts = cellPts->GetNumberOfPoints();
3332 std::vector<vtkIdType> cellptId(npts);
3333 vtkIdType iid[4];
3334 vtkIdType *v_id = NULL;
3335 vtkIdType *verts, v1, v2;
3336 vtkIdType ptId;
3337 vtkIdType ptIdout[4];
3338 vtkIdType tab_id[6];
3339 vtkIdType ptstetra = 4;
3340
3341 int i,j,k;
3342 int allInside;
3343
3344 vtkIdType edges[6][2] = { {0,1}, {1,2}, {2,0},
3345 {0,3}, {1,3}, {2,3} }; /* Edges Tetrahedron */
3346 double deltaScalar;
3347 double p[6], t;
3348 double v_tetra[4][3];
3349 double v[3], x[3];
3350 double *p1, *p2;
3351
3352 unsigned int idtetranew;
3353 unsigned int idcellnew;
3354 unsigned int planes;
3355
3356 for (i=0; i<npts; i++)
3357 {
3358 cellptId[i] = cellIds->GetId(i);
3359 }
3360
3361 this->CellGrid(cellType,npts,&cellptId[0],arraytetra); // Convert all volume cells to tetrahedra
3362
3363 unsigned int totalnewtetra = arraytetra->GetNumberOfCells();
3364 for (idtetranew = 0 ; idtetranew < totalnewtetra; idtetranew++)
3365 {
3366 arraytetra->GetNextCell(ptstetra,v_id);
3367
3368 for (allInside=1, i=0; i<4; i++)
3369 {
3370 cellPts->GetPoint(v_id[i],v);
3371
3372 for(k=0;k<6;k++)
3373 {
3374 p[k] = this->PlaneNormal[k][0]*(v[0] - this->PlanePoint[k][0]) +
3375 this->PlaneNormal[k][1]*(v[1] - this->PlanePoint[k][1]) +
3376 this->PlaneNormal[k][2]*(v[2] - this->PlanePoint[k][2]);
3377 }
3378
3379 if (!((p[0] <= 0) && (p[1] <= 0) &&
3380 (p[2] <= 0) && (p[3] <= 0) &&
3381 (p[4] <= 0) && (p[5] <= 0)))
3382 {
3383 allInside = 0;
3384 }
3385 }//for all points of the tetrahedron.
3386
3387 // Test Outside
3388 unsigned int test[6] = {1,1,1,1,1,1};
3389 for (i=0; i<4; i++)
3390 {
3391 ptIdout[i] = cellIds->GetId(v_id[i]);
3392 cellPts->GetPoint(v_id[i],v_tetra[i]);
3393
3394 // Use plane equation
3395 for(k=0;k<6;k++)
3396 {
3397 p[k] = this->PlaneNormal[k][0]*(v_tetra[i][0] - this->PlanePoint[k][0]) +
3398 this->PlaneNormal[k][1]*(v_tetra[i][1] - this->PlanePoint[k][1]) +
3399 this->PlaneNormal[k][2]*(v_tetra[i][2] - this->PlanePoint[k][2]);
3400 }
3401
3402
3403 for(k=0;k<3;k++)
3404 {
3405 if (p[2*k] < 0)
3406 {
3407 test[2*k] = 0;
3408 }
3409 if (p[2*k+1] < 0)
3410 {
3411 test[2*k+1] = 0;
3412 }
3413 }
3414
3415 }//for all points of the cell.
3416
3417 if (!allInside && ((test[0] == 1)|| (test[1] == 1) ||
3418 (test[2] == 1)|| (test[3] == 1) ||
3419 (test[4] == 1)|| (test[5] == 1)))
3420 {
3421 for (i=0; i<4; i++)
3422 {
3423 if ( locator->InsertUniquePoint(v_tetra[i], iid[i]) )
3424 {
3425 outPD->CopyData(inPD,ptIdout[i], iid[i]);
3426 }
3427 }
3428 int newCellId = tets[1]->InsertNextCell(4,iid);
3429 outCD[1]->CopyData(inCD,cellId,newCellId);
3430 continue; // Tetrahedron is outside.
3431 }
3432
3433 for (i=0; i<4; i++)
3434 {
3435 ptId = cellIds->GetId(v_id[i]);
3436 cellPts->GetPoint(v_id[i],v);
3437
3438 // Currently all points are injected because of the possibility
3439 // of intersection point merging.
3440
3441 if ( locator->InsertUniquePoint(v, iid[i]) )
3442 {
3443 outPD->CopyData(inPD,ptId, iid[i]);
3444 }
3445
3446 }//for all points of the tetrahedron.
3447
3448 if ( allInside )
3449 {
3450 int newCellId = tets[0]->InsertNextCell(4,iid); // Tetrahedron inside.
3451 outCD[0]->CopyData(inCD,cellId,newCellId);
3452 continue;
3453 }
3454
3455 //float *pc1 , *pc2;
3456 double *pedg1,*pedg2;
3457
3458 // Tetrahedron Intersection Cases
3459 const unsigned int tab4[6][6] = { {1,1,0,3,3,2},
3460 {2,0,0,3,2,1},
3461 {3,3,2,0,2,1},
3462 {1,0,2,0,1,3},
3463 {0,0,1,2,3,3},
3464 {0,1,2,1,2,3}};
3465 const unsigned int tab3[4][6] = { {0,2,1,1,3,2},
3466 {0,1,2,0,2,3},
3467 {0,1,2,1,0,3},
3468 {0,1,2,0,1,2}};
3469 const unsigned int tab2[12][5] = { {0,0,1,2,3},
3470 {2,1,0,1,3},
3471 {1,0,1,0,3},
3472 {2,0,1,3,0},
3473 {3,1,0,1,0},
3474 {1,0,1,2,0},
3475 {3,1,0,2,1},
3476 {2,1,0,0,1},
3477 {0,0,1,3,1},
3478 {1,0,1,3,2},
3479 {3,1,0,0,2},
3480 {0,0,1,1,2}};
3481 const unsigned int tab1[12][3] = { {2,3,1},
3482 {3,2,0},
3483 {3,0,1},
3484 {0,3,2},
3485 {1,3,0},
3486 {3,1,2},
3487 {2,1,0},
3488 {1,2,3},
3489 {2,0,3},
3490 {0,2,1},
3491 {0,1,3},
3492 {1,0,2}};
3493
3494 vtkCellArray *cellarray = vtkCellArray::New();
3495 vtkCellArray *cellarrayout = vtkCellArray::New();
3496 int newCellId = cellarray->InsertNextCell(4,iid);
3497
3498 // Test Cell intersection with each plane of box
3499 // FIXME: there is no difference in the for loop (planes has no influence!)
3500 for (planes = 0; planes < 6; planes++)
3501 {
3502 unsigned int totalnewcells = cellarray->GetNumberOfCells();
3503 vtkCellArray *newcellArray = vtkCellArray::New();
3504
3505 for (idcellnew = 0 ; idcellnew < totalnewcells; idcellnew++)
3506 {
3507 unsigned int i0,i1;
3508 unsigned int num_inter = 0;
3509 unsigned int edges_inter = 0;
3510 vtkIdType p_id[4];
3511
3512 cellarray->GetNextCell(npts,v_id);
3513
3514 newPoints->GetPoint(v_id[0],v_tetra[0]); //coord (x,y,z)
3515 newPoints->GetPoint(v_id[1],v_tetra[1]);
3516 newPoints->GetPoint(v_id[2],v_tetra[2]);
3517 newPoints->GetPoint(v_id[3],v_tetra[3]);
3518
3519 p[0] = this->PlaneNormal[planes][0]*(v_tetra[0][0] - this->PlanePoint[planes][0]) +
3520 this->PlaneNormal[planes][1]*(v_tetra[0][1] - this->PlanePoint[planes][1]) +
3521 this->PlaneNormal[planes][2]*(v_tetra[0][2] - this->PlanePoint[planes][2]);
3522 p[1] = this->PlaneNormal[planes][0]*(v_tetra[1][0] - this->PlanePoint[planes][0]) +
3523 this->PlaneNormal[planes][1]*(v_tetra[1][1] - this->PlanePoint[planes][1]) +
3524 this->PlaneNormal[planes][2]*(v_tetra[1][2] - this->PlanePoint[planes][2]);
3525 p[2] = this->PlaneNormal[planes][0]*(v_tetra[2][0] - this->PlanePoint[planes][0]) +
3526 this->PlaneNormal[planes][1]*(v_tetra[2][1] - this->PlanePoint[planes][1]) +
3527 this->PlaneNormal[planes][2]*(v_tetra[2][2] - this->PlanePoint[planes][2]);
3528 p[3] = this->PlaneNormal[planes][0]*(v_tetra[3][0] - this->PlanePoint[planes][0]) +
3529 this->PlaneNormal[planes][1]*(v_tetra[3][1] - this->PlanePoint[planes][1]) +
3530 this->PlaneNormal[planes][2]*(v_tetra[3][2] - this->PlanePoint[planes][2]);
3531
3532 for (int edgeNum=0; edgeNum < 6; edgeNum++)
3533 {
3534 verts = edges[edgeNum];
3535
3536 p1 = v_tetra[verts[0]];
3537 p2 = v_tetra[verts[1]];
3538 double s1 = p[verts[0]];
3539 double s2 = p[verts[1]];
3540 if ( (s1 * s2) < 0)
3541 {
3542 deltaScalar = s2 - s1;
3543
3544 if (deltaScalar > 0)
3545 {
3546 pedg1 = p1; pedg2 = p2;
3547 v1 = verts[0]; v2 = verts[1];
3548 }
3549 else
3550 {
3551 pedg1 = p2; pedg2 = p1;
3552 v1 = verts[1]; v2 = verts[0];
3553 deltaScalar = -deltaScalar;
3554 t = s1; s1 = s2; s2 = t;
3555 }
3556
3557 // linear interpolation
3558 t = ( deltaScalar == 0.0 ? 0.0 :
3559 ( - s1) / deltaScalar );
3560
3561 for (j=0; j<3; j++)
3562 {
3563 x[j] = pedg1[j] + t*(pedg2[j] - pedg1[j]);
3564 }
3565
3566 // Incorporate point into output and interpolate edge data as necessary
3567 edges_inter = edges_inter * 10 + (edgeNum+1);
3568
3569 if ( locator->InsertUniquePoint(x, p_id[num_inter]) )
3570 {
3571 this->InterpolateEdge(outPD, p_id[num_inter],
3572 v_id[v1], v_id[v2], t);
3573 }
3574
3575 num_inter++;
3576 }//if edge intersects value
3577 }//for all edges
3578
3579 if (num_inter == 0)
3580 {
3581 unsigned int outside = 0;
3582 for(i=0;i<4;i++)
3583 {
3584 if (p[i] > 0)
3585 { // If only one vertex is ouside, so the tetrahedron is outside
3586 outside = 1; // because there is not intersection.
3587 break; // some vertex could be on plane, so you need to test all vertex
3588 }
3589 }
3590 if (outside == 0)
3591 {
3592 // else it is possible intersection if other plane
3593 newCellId = newcellArray->InsertNextCell(4,v_id);
3594 }
3595 else
3596 {
3597 newCellId = tets[1]->InsertNextCell(4,v_id);
3598 outCD[1]->CopyData(inCD,cellId,newCellId);
3599 }
3600 continue;
3601 }
3602
3603 switch(num_inter)
3604 {
3605 case 4: // We have two wedges
3606 switch(edges_inter)
3607 {
3608 case 1246:
3609 i0 = 0;
3610 break;
3611 case 2345:
3612 i0 = 2;
3613 break;
3614 case 1356:
3615 i0 = 4;
3616 break;
3617 default:
3618 vtkErrorMacro( << "Intersection not found: Num_inter = " <<
3619 num_inter << " Edges_inter = " << edges_inter );
3620 continue;
3621 }
3622 if (p[3] > 0)
3623 {
3624 // The v_tetra[3] is outside box, so the first wedge is outside
3625 // ps: v_tetra[3] is always in first wedge (see tab)
3626
3627 tab_id[0] = p_id[tab4[i0+1][0]]; // Inside
3628 tab_id[1] = v_id[tab4[i0+1][1]];
3629 tab_id[2] = p_id[tab4[i0+1][2]];
3630 tab_id[3] = p_id[tab4[i0+1][3]];
3631 tab_id[4] = v_id[tab4[i0+1][4]];
3632 tab_id[5] = p_id[tab4[i0+1][5]];
3633 this->CreateTetra(6,tab_id,newcellArray);
3634
3635 tab_id[0] = p_id[tab4[i0][0]]; // Outside
3636 tab_id[1] = v_id[tab4[i0][1]];
3637 tab_id[2] = p_id[tab4[i0][2]];
3638 tab_id[3] = p_id[tab4[i0][3]];
3639 tab_id[4] = v_id[tab4[i0][4]];
3640 tab_id[5] = p_id[tab4[i0][5]];
3641 this->CreateTetra(6,tab_id,cellarrayout);
3642 }
3643 else
3644 {
3645 tab_id[0] = p_id[tab4[i0][0]]; // Inside
3646 tab_id[1] = v_id[tab4[i0][1]];
3647 tab_id[2] = p_id[tab4[i0][2]];
3648 tab_id[3] = p_id[tab4[i0][3]];
3649 tab_id[4] = v_id[tab4[i0][4]];
3650 tab_id[5] = p_id[tab4[i0][5]];
3651 this->CreateTetra(6,tab_id,newcellArray);
3652
3653 tab_id[0] = p_id[tab4[i0+1][0]]; // Outside
3654 tab_id[1] = v_id[tab4[i0+1][1]];
3655 tab_id[2] = p_id[tab4[i0+1][2]];
3656 tab_id[3] = p_id[tab4[i0+1][3]];
3657 tab_id[4] = v_id[tab4[i0+1][4]];
3658 tab_id[5] = p_id[tab4[i0+1][5]];
3659 this->CreateTetra(6,tab_id,cellarrayout);
3660 }
3661
3662 break;
3663 case 3: // We have one tetrahedron and one wedge
3664 switch(edges_inter)
3665 {
3666 case 134:
3667 i0 = 0;
3668 break;
3669 case 125:
3670 i0 = 1;
3671 break;
3672 case 236:
3673 i0 = 2;
3674 break;
3675 case 456:
3676 i0 = 3;
3677 break;
3678 default:
3679 vtkErrorMacro( << "Intersection not found: Num_inter = " <<
3680 num_inter << " Edges_inter = " << edges_inter );
3681 continue;
3682 }
3683 if (p[i0] > 0)
3684 {
3685 // Isolate vertex is outside box, so the tetrahedron is outside
3686 tab_id[0] = p_id[tab3[i0][0]]; // Inside
3687 tab_id[1] = p_id[tab3[i0][1]];
3688 tab_id[2] = p_id[tab3[i0][2]];
3689 tab_id[3] = v_id[tab3[i0][3]];
3690 tab_id[4] = v_id[tab3[i0][4]];
3691 tab_id[5] = v_id[tab3[i0][5]];
3692 this->CreateTetra(6,tab_id,newcellArray);
3693
3694 tab_id[0] = p_id[tab3[i0][0]]; // Outside
3695 tab_id[1] = p_id[tab3[i0][1]];
3696 tab_id[2] = p_id[tab3[i0][2]];
3697 tab_id[3] = v_id[i0];
3698 newCellId = cellarrayout->InsertNextCell(4,tab_id);
3699 }
3700 else
3701 {
3702 tab_id[0] = p_id[tab3[i0][0]]; // Inside
3703 tab_id[1] = p_id[tab3[i0][1]];
3704 tab_id[2] = p_id[tab3[i0][2]];
3705 tab_id[3] = v_id[i0];
3706 newCellId = newcellArray->InsertNextCell(4,tab_id);
3707
3708 tab_id[0] = p_id[tab3[i0][0]]; // Outside
3709 tab_id[1] = p_id[tab3[i0][1]];
3710 tab_id[2] = p_id[tab3[i0][2]];
3711 tab_id[3] = v_id[tab3[i0][3]];
3712 tab_id[4] = v_id[tab3[i0][4]];
3713 tab_id[5] = v_id[tab3[i0][5]];
3714 this->CreateTetra(6,tab_id,cellarrayout);
3715 }
3716 break;
3717 case 2: // We have one tetrahedron and one pyramid
3718 switch(edges_inter) // i1 = vertex of the tetrahedron
3719 {
3720 case 12:
3721 i0 = 0; i1 = 1;
3722 break;
3723 case 13:
3724 i0 = 1; i1 = 0;
3725 break;
3726 case 23:
3727 i0 = 2; i1 = 2;
3728 break;
3729 case 25:
3730 i0 = 3; i1 = 1;
3731 break;
3732 case 26:
3733 i0 = 4; i1 = 2;
3734 break;
3735 case 56:
3736 i0 = 5; i1 = 3;
3737 break;
3738 case 34:
3739 i0 = 6; i1 = 0;
3740 break;
3741 case 46:
3742 i0 = 7; i1 = 3;
3743 break;
3744 case 36:
3745 i0 = 8; i1 = 2;
3746 break;
3747 case 14:
3748 i0 = 9; i1 = 0;
3749 break;
3750 case 15:
3751 i0 = 10; i1 = 1;
3752 break;
3753 case 45:
3754 i0 = 11; i1 = 3;
3755 break;
3756 default:
3757 vtkErrorMacro( << "Intersection not found: Num_inter = " <<
3758 num_inter << " Edges_inter = %" << edges_inter );
3759 continue;
3760 }
3761
3762 if (p[i1] > 0)
3763 {
3764 // Isolate vertex is outside box, so the tetrahedron is outside
3765 tab_id[0] = v_id[tab2[i0][0]]; // Inside
3766 tab_id[1] = p_id[tab2[i0][1]];
3767 tab_id[2] = p_id[tab2[i0][2]];
3768 tab_id[3] = v_id[tab2[i0][3]];
3769 tab_id[4] = v_id[tab2[i0][4]];
3770 this->CreateTetra(5,tab_id,newcellArray);
3771
3772 tab_id[0] = v_id[i1]; // Outside
3773 tab_id[1] = v_id[tab2[i0][4]];
3774 tab_id[2] = p_id[tab2[i0][2]];
3775 tab_id[3] = p_id[tab2[i0][1]];
3776 newCellId = cellarrayout->InsertNextCell(4,tab_id);
3777 }
3778 else
3779 {
3780 tab_id[0] = v_id[i1]; // Inside
3781 tab_id[1] = v_id[tab2[i0][4]];
3782 tab_id[2] = p_id[tab2[i0][2]];
3783 tab_id[3] = p_id[tab2[i0][1]];
3784 newCellId = newcellArray->InsertNextCell(4,tab_id);
3785
3786 tab_id[0] = v_id[tab2[i0][0]]; // Outside
3787 tab_id[1] = p_id[tab2[i0][1]];
3788 tab_id[2] = p_id[tab2[i0][2]];
3789 tab_id[3] = v_id[tab2[i0][3]];
3790 tab_id[4] = v_id[tab2[i0][4]];
3791 this->CreateTetra(5,tab_id,cellarrayout);
3792 }
3793 break;
3794 case 1: // We have two tetrahedron.
3795 if ((edges_inter > 6) || (edges_inter < 1))
3796 {
3797 vtkErrorMacro( << "Intersection not found: Num_inter = " <<
3798 num_inter << " Edges_inter = " << edges_inter );
3799 continue;
3800 }
3801 if (p[tab1[2*edges_inter-1][2]] > 0)
3802 {
3803 // Isolate vertex is outside box, so the tetrahedron is outside
3804 tab_id[0] = p_id[0]; // Inside
3805 tab_id[1] = v_id[tab1[2*edges_inter-2][0]];
3806 tab_id[2] = v_id[tab1[2*edges_inter-2][1]];
3807 tab_id[3] = v_id[tab1[2*edges_inter-2][2]];
3808 newCellId = newcellArray->InsertNextCell(4,tab_id);
3809
3810 tab_id[0] = p_id[0]; // Outside
3811 tab_id[1] = v_id[tab1[2*edges_inter-1][0]];
3812 tab_id[2] = v_id[tab1[2*edges_inter-1][1]];
3813 tab_id[3] = v_id[tab1[2*edges_inter-1][2]];
3814 newCellId = cellarrayout->InsertNextCell(4,tab_id);
3815 }
3816 else
3817 {
3818 tab_id[0] = p_id[0]; // Inside
3819 tab_id[1] = v_id[tab1[2*edges_inter-1][0]];
3820 tab_id[2] = v_id[tab1[2*edges_inter-1][1]];
3821 tab_id[3] = v_id[tab1[2*edges_inter-1][2]];
3822 newCellId = newcellArray->InsertNextCell(4,tab_id);
3823
3824 tab_id[0] = p_id[0]; // Outside
3825 tab_id[1] = v_id[tab1[2*edges_inter-2][0]];
3826 tab_id[2] = v_id[tab1[2*edges_inter-2][1]];
3827 tab_id[3] = v_id[tab1[2*edges_inter-2][2]];
3828 newCellId = cellarrayout->InsertNextCell(4,tab_id);
3829 }
3830 break;
3831 }
3832 } // for all new cells
3833 cellarray->Delete();
3834 cellarray = newcellArray;
3835 } //for all planes
3836
3837 unsigned int totalnewcells = cellarray->GetNumberOfCells(); // Inside
3838
3839 for (idcellnew = 0 ; idcellnew < totalnewcells; idcellnew++)
3840 {
3841 cellarray->GetNextCell(npts,v_id);
3842 newCellId = tets[0]->InsertNextCell(npts,v_id);
3843 outCD[0]->CopyData(inCD,cellId,newCellId);
3844 }
3845 cellarray->Delete();
3846
3847 totalnewcells = cellarrayout->GetNumberOfCells(); // Outside
3848
3849 for (idcellnew = 0 ; idcellnew < totalnewcells; idcellnew++)
3850 {
3851 cellarrayout->GetNextCell(npts,v_id);
3852 newCellId = tets[1]->InsertNextCell(npts,v_id);
3853 outCD[1]->CopyData(inCD,cellId,newCellId);
3854 }
3855 cellarrayout->Delete();
3856 }
3857 arraytetra->Delete();
3858 }
3859
3860 //-------------------------------------------------------
3861
3862 //-------------------------------------------------------
ClipBox2D(vtkPoints * newPoints,vtkGenericCell * cell,vtkIncrementalPointLocator * locator,vtkCellArray * tets,vtkPointData * inPD,vtkPointData * outPD,vtkCellData * inCD,vtkIdType cellId,vtkCellData * outCD)3863 void vtkBoxClipDataSet::ClipBox2D(vtkPoints *newPoints,
3864 vtkGenericCell *cell,
3865 vtkIncrementalPointLocator *locator,
3866 vtkCellArray *tets,
3867 vtkPointData *inPD,
3868 vtkPointData *outPD,
3869 vtkCellData *inCD,
3870 vtkIdType cellId,
3871 vtkCellData *outCD)
3872 {
3873 vtkIdType cellType = cell->GetCellType();
3874 vtkIdList *cellIds = cell->GetPointIds();
3875 vtkCellArray *arraytriangle = vtkCellArray::New();
3876 vtkPoints *cellPts = cell->GetPoints();
3877 vtkIdType npts = cellPts->GetNumberOfPoints();
3878 std::vector<vtkIdType> cellptId(npts);
3879 vtkIdType iid[3];
3880 vtkIdType *v_id = NULL;
3881 vtkIdType *verts, v1, v2;
3882 vtkIdType ptId;
3883 vtkIdType tab_id[3];
3884 vtkIdType ptstriangle = 3;
3885
3886 int i,j;
3887 unsigned int allInside;
3888 unsigned int planes;
3889 unsigned int cutInd;
3890
3891 vtkIdType edges[3][2] = { {0,1}, {1,2}, {2,0}}; /* Edges Triangle*/
3892 double value,deltaScalar;
3893 double t, *p1, *p2;
3894 double v[3],x[3];
3895 double v_triangle[3][3];
3896
3897 // To avoid false positive warning.
3898 tab_id[0]=0;
3899 tab_id[1]=0;
3900 tab_id[2]=0;
3901
3902 for (i=0; i<npts; i++)
3903 {
3904 cellptId[i] = cellIds->GetId(i);
3905 }
3906
3907 // Convert all 2d cells to triangle
3908 this->CellGrid(cellType,npts,&cellptId[0],arraytriangle);
3909
3910 unsigned int totalnewtriangle= arraytriangle->GetNumberOfCells();
3911 unsigned int idtrianglenew;
3912
3913 for (idtrianglenew = 0 ; idtrianglenew < totalnewtriangle; idtrianglenew++)
3914 {
3915 arraytriangle->GetNextCell(ptstriangle,v_id);
3916
3917 for (allInside=1, i=0; i<3; i++)
3918 {
3919 cellPts->GetPoint(v_id[i],v);
3920
3921 if (!(((v[0] >= this->BoundBoxClip[0][0])&&(v[0] <= this->BoundBoxClip[0][1]) &&
3922 (v[1] >= this->BoundBoxClip[1][0])&&(v[1] <= this->BoundBoxClip[1][1])&&
3923 (v[2] >= this->BoundBoxClip[2][0])&&(v[2] <= this->BoundBoxClip[2][1]))))
3924 {
3925 //outside,its type might change later (nearby intersection)
3926 allInside = 0;
3927 }
3928 }
3929
3930 // Test Outside:
3931 if (!allInside)
3932 {
3933 unsigned int test[6] = {1,1,1,1,1,1};
3934 for (i=0; i<3; i++)
3935 {
3936 cellPts->GetPoint(v_id[i],v);
3937
3938 if (v[0] >= this->BoundBoxClip[0][0])
3939 {
3940 test[0] = 0;
3941 }
3942 if (v[0] <= this->BoundBoxClip[0][1])
3943 {
3944 test[1] = 0;
3945 }
3946 if (v[1] >= this->BoundBoxClip[1][0])
3947 {
3948 test[2] = 0;
3949 }
3950 if (v[1] <= this->BoundBoxClip[1][1])
3951 {
3952 test[3] = 0;
3953 }
3954 if (v[2] >= this->BoundBoxClip[2][0])
3955 {
3956 test[4] = 0;
3957 }
3958 if (v[2] <= this->BoundBoxClip[2][1])
3959 {
3960 test[5] = 0;
3961 }
3962 }//for all points of the triangle.
3963
3964 if ((test[0] == 1)|| (test[1] == 1) ||
3965 (test[2] == 1)|| (test[3] == 1) ||
3966 (test[4] == 1)|| (test[5] == 1))
3967 {
3968 continue; // Triangle is outside.
3969 }
3970 }
3971
3972 for (i=0; i<3; i++)
3973 {
3974 // Currently all points are injected because of the possibility
3975 // of intersection point merging.
3976 ptId = cellIds->GetId(v_id[i]);
3977 cellPts->GetPoint(v_id[i],v);
3978 if ( locator->InsertUniquePoint(v, iid[i]) )
3979 {
3980 outPD->CopyData(inPD,ptId, iid[i]);
3981 }
3982
3983 }//for all points of the triangle.
3984
3985 if ( allInside )
3986 {
3987 // Triangle inside.
3988 int newCellId = tets->InsertNextCell(3,iid);
3989 outCD->CopyData(inCD,cellId,newCellId);
3990 continue;
3991 }
3992
3993 //float *pc1 , *pc2;
3994 double *pedg1,*pedg2;
3995
3996 // Triangle intersection cases
3997
3998 unsigned int tab2[3][4] = { {1,2,1,0},
3999 {2,0,0,1},
4000 {0,1,0,1}};
4001 unsigned int tab1[3][2] = { {2,1},
4002 {0,2},
4003 {1,0}};
4004
4005 vtkCellArray *cellarray = vtkCellArray::New();
4006 int newCellId = cellarray->InsertNextCell(3,iid);
4007 unsigned int idcellnew;
4008
4009 // Test triangle intersection with each plane of box
4010 for (planes = 0; planes < 6; planes++)
4011 {
4012 // The index of the dimension of the cut plane (x == 0, y == 1, z == 2).
4013 cutInd = planes/2;
4014
4015 // The plane is always parallel to unitary cube.
4016 value = this->BoundBoxClip[cutInd][planes%2];
4017
4018 unsigned int totalnewcells = cellarray->GetNumberOfCells();
4019 vtkCellArray *newcellArray = vtkCellArray::New();
4020
4021 for (idcellnew = 0 ; idcellnew < totalnewcells; idcellnew++)
4022 {
4023 unsigned int num_inter = 0;
4024 unsigned int edges_inter = 0;
4025 unsigned int i0;
4026 vtkIdType p_id[3];
4027 cellarray->GetNextCell(npts,v_id);
4028
4029 newPoints->GetPoint(v_id[0],v_triangle[0]); //coord (x,y,z)
4030 newPoints->GetPoint(v_id[1],v_triangle[1]);
4031 newPoints->GetPoint(v_id[2],v_triangle[2]);
4032 for (int edgeNum=0; edgeNum < 3; edgeNum++)
4033 {
4034 verts = edges[edgeNum];
4035
4036 p1 = v_triangle[verts[0]];
4037 p2 = v_triangle[verts[1]];
4038
4039 if ( (p1[cutInd] < value && value < p2[cutInd]) ||
4040 (p2[cutInd] < value && value < p1[cutInd]) )
4041 {
4042 deltaScalar = p2[cutInd] - p1[cutInd];
4043
4044 if (deltaScalar > 0)
4045 {
4046 pedg1 = p1; pedg2 = p2;
4047 v1 = verts[0]; v2 = verts[1];
4048 }
4049 else
4050 {
4051 pedg1 = p2; pedg2 = p1;
4052 v1 = verts[1]; v2 = verts[0];
4053 deltaScalar = -deltaScalar;
4054 }
4055
4056 // linear interpolation
4057 t = ( deltaScalar == 0.0 ? 0.0 : (value - pedg1[cutInd]) / deltaScalar );
4058
4059 for (j=0; j<3; j++)
4060 {
4061 x[j] = pedg1[j] + t*(pedg2[j] - pedg1[j]);
4062 }
4063
4064 // Incorporate point into output and interpolate edge data as necessary
4065 edges_inter = edges_inter * 10 + (edgeNum+1);
4066 if ( locator->InsertUniquePoint(x, p_id[num_inter]) )
4067 {
4068 this->InterpolateEdge(outPD, p_id[num_inter],
4069 v_id[v1], v_id[v2], t);
4070 }
4071
4072 num_inter++;
4073 }//if edge intersects value
4074 }//for all edges
4075
4076 if (num_inter == 0)
4077 {
4078 unsigned int outside = 0;
4079 for(i=0; i<3; i++)
4080 {
4081 if (((v_triangle[i][cutInd] < value) && ((planes % 2) == 0)) ||
4082 // If only one vertex is ouside, so the triangle is outside
4083 ((v_triangle[i][cutInd] > value) && ((planes % 2) == 1))) // because there is not intersection.
4084 {
4085 outside = 1;
4086 break;
4087 }
4088 }
4089 if(outside == 0)
4090 {
4091 // else it is possible intersection if other plane
4092 newCellId = newcellArray->InsertNextCell(3,v_id);
4093 }
4094 continue;
4095 }
4096 switch(num_inter)
4097 {
4098 case 2: // We have one quad and one triangle
4099 switch(edges_inter)
4100 {
4101 case 12:
4102 i0 = 1;
4103 break;
4104 case 23:
4105 i0 = 2;
4106 break;
4107 case 13:
4108 i0 = 0;
4109 break;
4110 default:
4111 vtkErrorMacro( << "Intersection not found: Num_inter = " <<
4112 num_inter << " Edges_inter = " << edges_inter );
4113 continue;
4114 }
4115 if (((v_triangle[i0][cutInd] < value) && ((planes % 2) == 0)) ||
4116 ((v_triangle[i0][cutInd] > value) && ((planes % 2) == 1)))
4117 {
4118 // The v_triangle[i0] is outside box, so
4119 // The Quad is inside: two triangles: (v0,v1,p0) and (p0,p1,v1)
4120 tab_id[0] = v_id[tab2[i0][0]];
4121 tab_id[1] = v_id[tab2[i0][1]];
4122 tab_id[2] = p_id[tab2[i0][2]];
4123 newCellId = newcellArray->InsertNextCell(3,tab_id);
4124 tab_id[0] = p_id[tab2[i0][2]];
4125 tab_id[1] = p_id[tab2[i0][3]];
4126 tab_id[2] = v_id[tab2[i0][0]];
4127 newCellId = newcellArray->InsertNextCell(3,tab_id);
4128 }
4129 else
4130 {
4131 // The Triangle is inside: (v0,p0,p1)
4132 // The correct winding of the new triangle depends on where the
4133 // plane intersected the original triangle.
4134 switch (edges_inter)
4135 {
4136 case 12:
4137 case 23:
4138 tab_id[0] = v_id[i0];
4139 tab_id[1] = p_id[1];
4140 tab_id[2] = p_id[0];
4141 break;
4142 case 13:
4143 tab_id[0] = v_id[i0];
4144 tab_id[1] = p_id[0];
4145 tab_id[2] = p_id[1];
4146 break;
4147 }
4148 newCellId = newcellArray->InsertNextCell(3,tab_id);
4149 }
4150 break;
4151
4152 case 1: // We have two triangles
4153 switch(edges_inter)
4154 {
4155 case 1:
4156 i0 = 0;
4157 break;
4158 case 2:
4159 i0 = 1;
4160 break;
4161 case 3:
4162 i0 = 2;
4163 break;
4164 default:
4165 vtkErrorMacro( << "Intersection not found: Num_inter = " <<
4166 num_inter << " Edges_inter = " << edges_inter );
4167 continue;
4168 }
4169 if (((v_triangle[i0][cutInd] < value) && ((planes % 2) == 0)) ||
4170 ((v_triangle[i0][cutInd] > value) && ((planes % 2) == 1)))
4171 {
4172 // Test one of the vertices vertex i0 is outside
4173 tab_id[0] = v_id[tab1[i0][1]];
4174 tab_id[1] = v_id[tab1[i0][0]];
4175 tab_id[2] = p_id[0];
4176 newCellId = newcellArray->InsertNextCell(3,tab_id);
4177 }
4178 else
4179 {
4180 tab_id[0] = v_id[tab1[i0][0]];
4181 tab_id[1] = v_id[i0];
4182 tab_id[2] = p_id[0];
4183 newCellId = newcellArray->InsertNextCell(3,tab_id);
4184 }
4185 break;
4186 }
4187 } // for all new cells
4188 cellarray->Delete();
4189 cellarray = newcellArray;
4190 } //for all planes
4191
4192 unsigned int totalnewcells = cellarray->GetNumberOfCells();
4193
4194 for (idcellnew = 0 ; idcellnew < totalnewcells; idcellnew++)
4195 {
4196 cellarray->GetNextCell(npts,v_id);
4197 newCellId = tets->InsertNextCell(npts,v_id);
4198 outCD->CopyData(inCD,cellId,newCellId);
4199 }
4200 cellarray->Delete();
4201 }
4202 arraytriangle->Delete();
4203 }
4204 //-------------------------------------------------------
4205
ClipBoxInOut2D(vtkPoints * newPoints,vtkGenericCell * cell,vtkIncrementalPointLocator * locator,vtkCellArray ** tets,vtkPointData * inPD,vtkPointData * outPD,vtkCellData * inCD,vtkIdType cellId,vtkCellData ** outCD)4206 void vtkBoxClipDataSet::ClipBoxInOut2D(vtkPoints *newPoints,
4207 vtkGenericCell *cell,
4208 vtkIncrementalPointLocator *locator,
4209 vtkCellArray **tets,
4210 vtkPointData *inPD,
4211 vtkPointData *outPD,
4212 vtkCellData *inCD,
4213 vtkIdType cellId,
4214 vtkCellData **outCD)
4215 {
4216 vtkIdType cellType = cell->GetCellType();
4217 vtkIdList *cellIds = cell->GetPointIds();
4218 vtkCellArray *arraytriangle = vtkCellArray::New();
4219 vtkPoints *cellPts = cell->GetPoints();
4220 vtkIdType npts = cellPts->GetNumberOfPoints();
4221 std::vector<vtkIdType> cellptId(npts);
4222 vtkIdType iid[3];
4223 vtkIdType *v_id = NULL;
4224 vtkIdType *verts, v1, v2;
4225 vtkIdType ptId;
4226 vtkIdType ptIdout[4];
4227 vtkIdType tab_id[3];
4228 vtkIdType ptstriangle = 3;
4229
4230 int i,j;
4231 unsigned int allInside;
4232 unsigned int cutInd;
4233 unsigned int planes;
4234 unsigned int idcellnew;
4235
4236 vtkIdType edges[3][2] = { {0,1}, {1,2}, {2,0}}; /* Edges Triangle */
4237 double value,deltaScalar;
4238 double t, *p1, *p2;
4239 double v[3],x[3];
4240 double v_triangle[3][3];
4241
4242 // To avoid false positive warning.
4243 tab_id[0]=0;
4244 tab_id[1]=0;
4245 tab_id[2]=0;
4246
4247 for (i=0; i<npts; i++)
4248 {
4249 cellptId[i] = cellIds->GetId(i);
4250 }
4251
4252 // Convert all 2D cells to triangle
4253 this->CellGrid(cellType,npts, &cellptId[0],arraytriangle);
4254 unsigned int totalnewtriangle = arraytriangle->GetNumberOfCells();
4255 unsigned int idtrianglenew;
4256
4257 for (idtrianglenew = 0 ; idtrianglenew < totalnewtriangle; idtrianglenew++)
4258 {
4259 arraytriangle->GetNextCell(ptstriangle,v_id);
4260
4261 for (allInside=1, i=0; i<3; i++)
4262 {
4263 cellPts->GetPoint(v_id[i],v);
4264
4265 if (!(((v[0] >= this->BoundBoxClip[0][0])&&(v[0] <= this->BoundBoxClip[0][1]) &&
4266 (v[1] >= this->BoundBoxClip[1][0])&&(v[1] <= this->BoundBoxClip[1][1])&&
4267 (v[2] >= this->BoundBoxClip[2][0])&&(v[2] <= this->BoundBoxClip[2][1]))))
4268 {
4269 //outside,its type might change later (nearby intersection)
4270 allInside = 0;
4271 }
4272 }//for all points of the cell.
4273
4274 // Test Outside: see(1)
4275 if (!allInside)
4276 {
4277 unsigned int test[6] = {1,1,1,1,1,1};
4278 for (i=0; i<3; i++)
4279 {
4280 ptIdout[i] = cellIds->GetId(v_id[i]);
4281 cellPts->GetPoint(v_id[i],v_triangle[i]);
4282
4283 if (v_triangle[i][0] >= this->BoundBoxClip[0][0])
4284 {
4285 test[0] = 0;
4286 }
4287 if (v_triangle[i][0] <= this->BoundBoxClip[0][1])
4288 {
4289 test[1] = 0;
4290 }
4291 if (v_triangle[i][1] >= this->BoundBoxClip[1][0])
4292 {
4293 test[2] = 0;
4294 }
4295 if (v_triangle[i][1] <= this->BoundBoxClip[1][1])
4296 {
4297 test[3] = 0;
4298 }
4299 if (v_triangle[i][2] >= this->BoundBoxClip[2][0])
4300 {
4301 test[4] = 0;
4302 }
4303 if (v_triangle[i][2] <= this->BoundBoxClip[2][1])
4304 {
4305 test[5] = 0;
4306 }
4307
4308 }//for all points of the cell.
4309
4310 if ((test[0] == 1)|| (test[1] == 1) ||
4311 (test[2] == 1)|| (test[3] == 1) ||
4312 (test[4] == 1)|| (test[5] == 1))
4313 {
4314 for (i=0; i<3; i++)
4315 {
4316 if ( locator->InsertUniquePoint(v_triangle[i], iid[i]) )
4317 {
4318 outPD->CopyData(inPD,ptIdout[i], iid[i]);
4319 }
4320 }
4321
4322 int newCellId = tets[1]->InsertNextCell(3,iid);
4323 outCD[1]->CopyData(inCD,cellId,newCellId);
4324 continue; // Triangle is outside.
4325 }
4326 }//if not allInside.
4327
4328 for (i=0; i<3; i++)
4329 {
4330 ptId = cellIds->GetId(v_id[i]);
4331 cellPts->GetPoint(v_id[i],v);
4332
4333 // Currently all points are injected because of the possibility
4334 // of intersection point merging.
4335 if ( locator->InsertUniquePoint(v, iid[i]) )
4336 {
4337 outPD->CopyData(inPD,ptId, iid[i]);
4338 }
4339 }//for all points of the triangle.
4340
4341 if ( allInside )
4342 {
4343 // Triangle inside.
4344 int newCellId = tets[0]->InsertNextCell(3,iid);
4345 outCD[0]->CopyData(inCD,cellId,newCellId);
4346 continue;
4347 }
4348
4349 //float *pc1 , *pc2;
4350 double *pedg1,*pedg2;
4351
4352 // Triangle intersection cases
4353
4354 unsigned int tab2[3][4] = { {1,2,1,0},
4355 {2,0,0,1},
4356 {0,1,0,1}};
4357 unsigned int tab1[3][2] = { {2,1},
4358 {0,2},
4359 {1,0}};
4360
4361 vtkCellArray *cellarray = vtkCellArray::New();
4362 vtkCellArray *cellarrayout = vtkCellArray::New();
4363 int newCellId = cellarray->InsertNextCell(3,iid);
4364
4365 // Test Cell intersection with each plane of box
4366 for (planes = 0; planes < 6; planes++)
4367 {
4368 // The index of the dimension of the cut plane (x == 0, y == 1, z == 2).
4369 cutInd = planes/2;
4370
4371 // The plane is always parallel to unitary cube.
4372 value = this->BoundBoxClip[cutInd][planes%2];
4373
4374 unsigned int totalnewcells = cellarray->GetNumberOfCells();
4375 vtkCellArray *newcellArray = vtkCellArray::New();
4376
4377 for (idcellnew = 0 ; idcellnew < totalnewcells; idcellnew++)
4378 {
4379 unsigned int num_inter = 0;
4380 unsigned int edges_inter = 0;
4381 unsigned int i0;
4382 vtkIdType p_id[3];
4383 cellarray->GetNextCell(npts,v_id);
4384
4385 newPoints->GetPoint(v_id[0],v_triangle[0]); //coord (x,y,z)
4386 newPoints->GetPoint(v_id[1],v_triangle[1]);
4387 newPoints->GetPoint(v_id[2],v_triangle[2]);
4388 for (int edgeNum=0; edgeNum < 3; edgeNum++)
4389 {
4390 verts = edges[edgeNum];
4391
4392 p1 = v_triangle[verts[0]];
4393 p2 = v_triangle[verts[1]];
4394
4395 if ( (p1[cutInd] < value && value < p2[cutInd]) ||
4396 (p2[cutInd] < value && value < p1[cutInd]) )
4397 {
4398 deltaScalar = p2[cutInd] - p1[cutInd];
4399
4400 if (deltaScalar > 0)
4401 {
4402 pedg1 = p1; pedg2 = p2;
4403 v1 = verts[0]; v2 = verts[1];
4404 }
4405 else
4406 {
4407 pedg1 = p2; pedg2 = p1;
4408 v1 = verts[1]; v2 = verts[0];
4409 deltaScalar = -deltaScalar;
4410 }
4411
4412 // linear interpolation
4413 t = ( deltaScalar == 0.0 ? 0.0 : (value - pedg1[cutInd]) / deltaScalar );
4414
4415 for (j=0; j<3; j++)
4416 {
4417 x[j] = pedg1[j] + t*(pedg2[j] - pedg1[j]);
4418 }
4419
4420 // Incorporate point into output and interpolate edge data as necessary
4421 edges_inter = edges_inter * 10 + (edgeNum+1);
4422 if ( locator->InsertUniquePoint(x, p_id[num_inter]) )
4423 {
4424 this->InterpolateEdge(outPD, p_id[num_inter],
4425 v_id[v1], v_id[v2], t);
4426 }
4427
4428 num_inter++;
4429 }//if edge intersects value
4430 }//for all edges
4431
4432 if (num_inter == 0)
4433 {
4434 unsigned int outside = 0;
4435 for(i=0; i<3; i++)
4436 {
4437 if (((v_triangle[i][cutInd] < value) && ((planes % 2) == 0)) ||
4438 ((v_triangle[i][cutInd] > value) && ((planes % 2) == 1)))
4439 {
4440 // If only one vertex is ouside, so the triangle is outside
4441 // because there is not intersection.
4442 outside = 1;
4443 break;
4444 }
4445 }
4446 if(outside == 0)
4447 {
4448 // else it is possible intersection if other plane
4449 newCellId = newcellArray->InsertNextCell(3,v_id);
4450 }
4451 else
4452 {
4453 newCellId = tets[1]->InsertNextCell(3,v_id);
4454 outCD[1]->CopyData(inCD,cellId,newCellId);
4455 }
4456 continue;
4457 }
4458 switch(num_inter)
4459 {
4460 case 2: // We have one quad and one triangle
4461 // i0 gets the index of the triangle point that lies alone on
4462 // one side of the plane.
4463 switch(edges_inter)
4464 {
4465 case 12:
4466 i0 = 1;
4467 break;
4468 case 23:
4469 i0 = 2;
4470 break;
4471 case 13:
4472 i0 = 0;
4473 break;
4474 default:
4475 vtkErrorMacro( << "Intersection not found: Num_inter = " <<
4476 num_inter << " Edges_inter = " << edges_inter );
4477 continue;
4478 }
4479 if (((v_triangle[i0][cutInd] < value) && ((planes % 2) == 0)) ||
4480 ((v_triangle[i0][cutInd] > value) && ((planes % 2) == 1)))
4481 {
4482 // The v_triangle[i0] is outside box, so
4483
4484 tab_id[0] = v_id[tab2[i0][0]]; // Quad Inside
4485 tab_id[1] = v_id[tab2[i0][1]];
4486 tab_id[2] = p_id[tab2[i0][2]];
4487 newCellId = newcellArray->InsertNextCell(3,tab_id);
4488 tab_id[0] = p_id[tab2[i0][2]];
4489 tab_id[1] = p_id[tab2[i0][3]];
4490 tab_id[2] = v_id[tab2[i0][0]];
4491 newCellId = newcellArray->InsertNextCell(3,tab_id);
4492
4493 switch (edges_inter) // Triangle Outside
4494 {
4495 case 12:
4496 case 23:
4497 tab_id[0] = v_id[i0];
4498 tab_id[1] = p_id[1];
4499 tab_id[2] = p_id[0];
4500 break;
4501 case 13:
4502 tab_id[0] = v_id[i0];
4503 tab_id[1] = p_id[0];
4504 tab_id[2] = p_id[1];
4505 break;
4506 }
4507 newCellId = cellarrayout->InsertNextCell(3,tab_id);
4508 }
4509 else
4510 {
4511 // The Triangle is inside: (v0,p0,p1)
4512 switch (edges_inter)
4513 {
4514 case 12:
4515 case 23:
4516 tab_id[0] = v_id[i0];
4517 tab_id[1] = p_id[1];
4518 tab_id[2] = p_id[0];
4519 break;
4520 case 13:
4521 tab_id[0] = v_id[i0];
4522 tab_id[1] = p_id[0];
4523 tab_id[2] = p_id[1];
4524 break;
4525 }
4526 newCellId = newcellArray->InsertNextCell(3,tab_id);
4527
4528 tab_id[0] = v_id[tab2[i0][0]]; // Quad is Outside
4529 tab_id[1] = v_id[tab2[i0][1]];
4530 tab_id[2] = p_id[tab2[i0][2]];
4531 newCellId = cellarrayout->InsertNextCell(3,tab_id);
4532 tab_id[0] = p_id[tab2[i0][2]];
4533 tab_id[1] = p_id[tab2[i0][3]];
4534 tab_id[2] = v_id[tab2[i0][0]];
4535 newCellId = cellarrayout->InsertNextCell(3,tab_id);
4536 }
4537 break;
4538
4539 case 1: // We have two triangles
4540 switch(edges_inter)
4541 {
4542 case 1:
4543 i0 = 0;
4544 break;
4545 case 2:
4546 i0 = 1;
4547 break;
4548 case 3:
4549 i0 = 2;
4550 break;
4551 default:
4552 vtkErrorMacro( << "Intersection not found: Num_inter = " <<
4553 num_inter << " Edges_inter = " << edges_inter );
4554 continue;
4555 }
4556 if (((v_triangle[i0][cutInd] < value) && ((planes % 2) == 0)) ||
4557 ((v_triangle[i0][cutInd] > value) && ((planes % 2) == 1)))
4558 {
4559 // Test one of the vertices vertex i0 is outside
4560
4561 tab_id[0] = v_id[tab1[i0][1]]; // Inside
4562 tab_id[1] = v_id[tab1[i0][0]];
4563 tab_id[2] = p_id[0];
4564 newCellId = newcellArray->InsertNextCell(3,tab_id);
4565
4566 tab_id[0] = v_id[tab1[i0][0]]; // Outside
4567 tab_id[1] = v_id[i0];
4568 tab_id[2] = p_id[0];
4569 newCellId = cellarrayout->InsertNextCell(3,tab_id);
4570 }
4571 else
4572 {
4573 tab_id[0] = v_id[tab1[i0][0]]; // Inside
4574 tab_id[1] = v_id[i0];
4575 tab_id[2] = p_id[0];
4576 newCellId = newcellArray->InsertNextCell(3,tab_id);
4577
4578 tab_id[0] = v_id[tab1[i0][1]]; // Outside
4579 tab_id[1] = v_id[tab1[i0][0]];
4580 tab_id[2] = p_id[0];
4581 newCellId = cellarrayout->InsertNextCell(3,tab_id);
4582 }
4583 break;
4584 }
4585 } // for all new cells
4586 cellarray->Delete();
4587 cellarray = newcellArray;
4588 } //for all planes
4589
4590 unsigned int totalnewcells = cellarray->GetNumberOfCells(); // Inside
4591
4592 for (idcellnew = 0 ; idcellnew < totalnewcells; idcellnew++)
4593 {
4594 cellarray->GetNextCell(npts,v_id);
4595 newCellId = tets[0]->InsertNextCell(npts,v_id);
4596 outCD[0]->CopyData(inCD,cellId,newCellId);
4597 }
4598 cellarray->Delete();
4599
4600 totalnewcells = cellarrayout->GetNumberOfCells(); // Outside
4601
4602 for (idcellnew = 0 ; idcellnew < totalnewcells; idcellnew++)
4603 {
4604 cellarrayout->GetNextCell(npts,v_id);
4605 newCellId = tets[1]->InsertNextCell(npts,v_id);
4606 outCD[1]->CopyData(inCD,cellId,newCellId);
4607 }
4608 cellarrayout->Delete();
4609 }
4610 arraytriangle->Delete();
4611 }
4612
4613 //-------------------------------------------------------
4614
ClipHexahedron2D(vtkPoints * newPoints,vtkGenericCell * cell,vtkIncrementalPointLocator * locator,vtkCellArray * tets,vtkPointData * inPD,vtkPointData * outPD,vtkCellData * inCD,vtkIdType cellId,vtkCellData * outCD)4615 void vtkBoxClipDataSet::ClipHexahedron2D(vtkPoints *newPoints,
4616 vtkGenericCell *cell,
4617 vtkIncrementalPointLocator *locator,
4618 vtkCellArray *tets,
4619 vtkPointData *inPD,
4620 vtkPointData *outPD,
4621 vtkCellData *inCD,
4622 vtkIdType cellId,
4623 vtkCellData *outCD)
4624 {
4625 vtkIdType cellType = cell->GetCellType();
4626 vtkIdList *cellIds = cell->GetPointIds();
4627 vtkCellArray *arraytriangle = vtkCellArray::New();
4628 vtkPoints *cellPts = cell->GetPoints();
4629 vtkIdType npts = cellPts->GetNumberOfPoints();
4630 std::vector<vtkIdType> cellptId(npts);
4631 vtkIdType iid[3];
4632 vtkIdType *v_id = NULL;
4633 vtkIdType *verts, v1, v2;
4634 vtkIdType ptId;
4635 vtkIdType tab_id[3];
4636 vtkIdType ptstriangle = 3;
4637
4638 int i,j,k;
4639 unsigned int allInside;
4640 unsigned int idtrianglenew;
4641 unsigned int idcellnew;
4642 unsigned int planes;
4643
4644 vtkIdType edges[3][2] = { {0,1}, {1,2}, {2,0}};
4645 double deltaScalar;
4646 double t, *p1, *p2;
4647 double v[3],x[3];
4648 double p[6];
4649 double v_triangle[3][3];
4650
4651 // To avoid false positive warning.
4652 tab_id[0]=0;
4653 tab_id[1]=0;
4654 tab_id[2]=0;
4655
4656 for (i=0; i<npts; i++)
4657 {
4658 cellptId[i] = cellIds->GetId(i);
4659 }
4660
4661 this->CellGrid(cellType,npts,&cellptId[0],arraytriangle); // Convert all volume cells to triangle
4662
4663 unsigned int totalnewtriangle = arraytriangle->GetNumberOfCells();
4664 for (idtrianglenew = 0 ; idtrianglenew < totalnewtriangle; idtrianglenew++)
4665 {
4666 arraytriangle->GetNextCell(ptstriangle,v_id);
4667
4668 for (allInside=1, i=0; i<3; i++)
4669 {
4670 cellPts->GetPoint(v_id[i],v);
4671
4672 for(k=0;k<6;k++)
4673 {
4674 p[k] = this->PlaneNormal[k][0]*(v[0] - this->PlanePoint[k][0])+
4675 this->PlaneNormal[k][1]*(v[1] - this->PlanePoint[k][1]) +
4676 this->PlaneNormal[k][2]*(v[2] - this->PlanePoint[k][2]);
4677 }
4678
4679 if (!((p[0] <= 0) && (p[1] <= 0) &&
4680 (p[2] <= 0) && (p[3] <= 0) &&
4681 (p[4] <= 0) && (p[5] <= 0)))
4682 {
4683 allInside = 0;
4684 }
4685 }//for all points of the triangle.
4686
4687 // Test Outside
4688 unsigned int test[6] = {1,1,1,1,1,1};
4689 for (i=0; i<3; i++)
4690 {
4691 cellPts->GetPoint(v_id[i],v);
4692
4693 // Use plane equation
4694 for(k=0;k<6;k++)
4695 {
4696 p[k] = this->PlaneNormal[k][0]*(v[0] - this->PlanePoint[k][0]) +
4697 this->PlaneNormal[k][1]*(v[1] - this->PlanePoint[k][1]) +
4698 this->PlaneNormal[k][2]*(v[2] - this->PlanePoint[k][2]);
4699 }
4700
4701 for(k=0;k<3;k++)
4702 {
4703 if (p[2*k] <= 0)
4704 {
4705 test[2*k] = 0;
4706 }
4707 if (p[2*k+1] <= 0)
4708 {
4709 test[2*k+1] = 0;
4710 }
4711 }
4712 }//for all points of the cell.
4713
4714 if (!allInside && ((test[0] == 1)|| (test[1] == 1) ||
4715 (test[2] == 1)|| (test[3] == 1) ||
4716 (test[4] == 1)|| (test[5] == 1)))
4717 {
4718 continue; // Triangle is outside.
4719 }
4720
4721 for (i=0; i<3; i++)
4722 {
4723 ptId = cellIds->GetId(v_id[i]);
4724 cellPts->GetPoint(v_id[i],v);
4725
4726 // Currently all points are injected because of the possibility
4727 // of intersection point merging.
4728
4729 if ( locator->InsertUniquePoint(v, iid[i]) )
4730 {
4731 outPD->CopyData(inPD,ptId, iid[i]);
4732 }
4733 }//for all points of the triangle.
4734
4735 if ( allInside )
4736 {
4737 int newCellId = tets->InsertNextCell(3,iid); // Triangle inside.
4738 outCD->CopyData(inCD,cellId,newCellId);
4739 continue;
4740 }
4741
4742 //float *pc1 , *pc2;
4743 double *pedg1,*pedg2;
4744
4745 unsigned int tab2[3][4] = { {1,2,1,0},
4746 {2,0,0,1},
4747 {0,1,0,1}};
4748 unsigned int tab1[3][2] = { {2,1},
4749 {0,2},
4750 {1,0}};
4751
4752 vtkCellArray *cellarray = vtkCellArray::New();
4753 int newCellId = cellarray->InsertNextCell(3,iid);
4754
4755 // Test Cell intersection with each plane of box
4756 for (planes = 0; planes < 6; planes++)
4757 {
4758 unsigned int totalnewcells = cellarray->GetNumberOfCells();
4759 vtkCellArray *newcellArray = vtkCellArray::New();
4760
4761 for (idcellnew = 0 ; idcellnew < totalnewcells; idcellnew++)
4762 {
4763 unsigned int i0;
4764 unsigned int num_inter = 0;
4765 unsigned int edges_inter = 0;
4766 vtkIdType p_id[3];
4767
4768 cellarray->GetNextCell(npts,v_id);
4769
4770 newPoints->GetPoint(v_id[0],v_triangle[0]); //coord (x,y,z)
4771 newPoints->GetPoint(v_id[1],v_triangle[1]);
4772 newPoints->GetPoint(v_id[2],v_triangle[2]);
4773
4774 p[0] = this->PlaneNormal[planes][0]*(v_triangle[0][0] - this->PlanePoint[planes][0]) +
4775 this->PlaneNormal[planes][1]*(v_triangle[0][1] - this->PlanePoint[planes][1]) +
4776 this->PlaneNormal[planes][2]*(v_triangle[0][2] - this->PlanePoint[planes][2]);
4777 p[1] = this->PlaneNormal[planes][0]*(v_triangle[1][0] - this->PlanePoint[planes][0]) +
4778 this->PlaneNormal[planes][1]*(v_triangle[1][1] - this->PlanePoint[planes][1]) +
4779 this->PlaneNormal[planes][2]*(v_triangle[1][2] - this->PlanePoint[planes][2]);
4780 p[2] = this->PlaneNormal[planes][0]*(v_triangle[2][0] - this->PlanePoint[planes][0]) +
4781 this->PlaneNormal[planes][1]*(v_triangle[2][1] - this->PlanePoint[planes][1]) +
4782 this->PlaneNormal[planes][2]*(v_triangle[2][2] - this->PlanePoint[planes][2]);
4783
4784 for (int edgeNum=0; edgeNum < 3; edgeNum++)
4785 {
4786 verts = edges[edgeNum];
4787
4788 p1 = v_triangle[verts[0]];
4789 p2 = v_triangle[verts[1]];
4790 double s1 = p[verts[0]];
4791 double s2 = p[verts[1]];
4792 if ( (s1 * s2) < 0)
4793 {
4794 deltaScalar = s2 - s1;
4795
4796 if (deltaScalar > 0)
4797 {
4798 pedg1 = p1; pedg2 = p2;
4799 v1 = verts[0]; v2 = verts[1];
4800 }
4801 else
4802 {
4803 pedg1 = p2; pedg2 = p1;
4804 v1 = verts[1]; v2 = verts[0];
4805 deltaScalar = -deltaScalar;
4806 t = s1; s1 = s2; s2 = t;
4807 }
4808
4809 // linear interpolation
4810 t = ( deltaScalar == 0.0 ? 0.0 : ( - s1) / deltaScalar );
4811
4812 for (j=0; j<3; j++)
4813 {
4814 x[j] = pedg1[j] + t*(pedg2[j] - pedg1[j]);
4815 }
4816
4817 // Incorporate point into output and interpolate edge data as necessary
4818 edges_inter = edges_inter * 10 + (edgeNum+1);
4819
4820 if ( locator->InsertUniquePoint(x, p_id[num_inter]) )
4821 {
4822 this->InterpolateEdge(outPD, p_id[num_inter],
4823 v_id[v1], v_id[v2], t);
4824 }
4825
4826 num_inter++;
4827 }//if edge intersects value
4828 }//for all edges
4829
4830 if (num_inter == 0)
4831 {
4832 unsigned int outside = 0;
4833 for(i=0;i<3;i++)
4834 {
4835 if (p[i] > 0) // If only one vertex is ouside, so the triangle is outside
4836 {
4837 outside = 1; // because there is not intersection.
4838 break; // some vertex could be on plane, so you need to test all vertex
4839 }
4840 }
4841 if (outside == 0)
4842 {
4843 // else it is possible intersection if other plane
4844 newCellId = newcellArray->InsertNextCell(3,v_id);
4845 }
4846 continue;
4847 }
4848 switch(num_inter)
4849 {
4850 case 2: // We have one quad and one triangle
4851 // i0 gets the index of the triangle point that lies alone on
4852 // one side of the plane.
4853 switch(edges_inter)
4854 {
4855 case 12:
4856 i0 = 1;
4857 break;
4858 case 23:
4859 i0 = 2;
4860 break;
4861 case 13:
4862 i0 = 0;
4863 break;
4864 default:
4865 vtkErrorMacro(<< "Intersection not found: Num_inter = " <<
4866 num_inter << " Edges_inter = " << edges_inter);
4867 continue;
4868 }
4869 if (p[i0] > 0)
4870 {
4871 // The v_triangle[3] is outside box, so
4872 // the first wedge is outside
4873
4874 // The v_triangle[3] is outside box, so the quad is outside
4875 // The Quad is inside: two triangles: (v0,v1,p0) and (p0,p1,v1)
4876 tab_id[0] = v_id[tab2[i0][0]];
4877 tab_id[1] = v_id[tab2[i0][1]];
4878 tab_id[2] = p_id[tab2[i0][2]];
4879 newCellId = newcellArray->InsertNextCell(3,tab_id);
4880 tab_id[0] = p_id[tab2[i0][2]];
4881 tab_id[1] = p_id[tab2[i0][3]];
4882 tab_id[2] = v_id[tab2[i0][0]];
4883 newCellId = newcellArray->InsertNextCell(3,tab_id);
4884 }
4885 else
4886 {
4887 // The Triangle is inside: (v0,p0,p1)
4888 // The correct winding of the new triangle depends on where the
4889 // plane intersected the original triangle.
4890 switch (edges_inter)
4891 {
4892 case 12:
4893 case 23:
4894 tab_id[0] = v_id[i0];
4895 tab_id[1] = p_id[1];
4896 tab_id[2] = p_id[0];
4897 break;
4898 case 13:
4899 tab_id[0] = v_id[i0];
4900 tab_id[1] = p_id[0];
4901 tab_id[2] = p_id[1];
4902 break;
4903 }
4904 newCellId = newcellArray->InsertNextCell(3,tab_id);
4905 }
4906 break;
4907
4908 case 1: // We have two triangles
4909 switch(edges_inter)
4910 {
4911 case 1:
4912 i0 = 0;
4913 break;
4914 case 2:
4915 i0 = 1;
4916 break;
4917 case 3:
4918 i0 = 2;
4919 break;
4920 default:
4921 vtkErrorMacro(<< "Intersection not found: Num_inter = " <<
4922 num_inter << " Edges_inter = " << edges_inter);
4923 continue;
4924 }
4925 if (p[i0] > 0)
4926 {
4927 // Isolate vertex is outside box, so the triangle is outside
4928 tab_id[0] = v_id[tab1[i0][1]];
4929 tab_id[1] = v_id[tab1[i0][0]];
4930 tab_id[2] = p_id[0];
4931 newCellId = newcellArray->InsertNextCell(3,tab_id);
4932 }
4933 else
4934 {
4935 tab_id[0] = v_id[tab1[i0][0]];
4936 tab_id[1] = v_id[i0];
4937 tab_id[2] = p_id[0];
4938 newCellId = newcellArray->InsertNextCell(3,tab_id);
4939 }
4940 break;
4941 }
4942 } // for all new cells
4943 cellarray->Delete();
4944 cellarray = newcellArray;
4945 } //for all planes
4946
4947 unsigned int totalnewcells = cellarray->GetNumberOfCells();
4948
4949 for (idcellnew = 0 ; idcellnew < totalnewcells; idcellnew++)
4950 {
4951 cellarray->GetNextCell(npts,v_id);
4952 newCellId = tets->InsertNextCell(npts,v_id);
4953 outCD->CopyData(inCD,cellId,newCellId);
4954 }
4955 cellarray->Delete();
4956 }
4957 arraytriangle->Delete();
4958 }
4959
4960 //-------------------------------------------------------
4961
ClipHexahedronInOut2D(vtkPoints * newPoints,vtkGenericCell * cell,vtkIncrementalPointLocator * locator,vtkCellArray ** tets,vtkPointData * inPD,vtkPointData * outPD,vtkCellData * inCD,vtkIdType cellId,vtkCellData ** outCD)4962 void vtkBoxClipDataSet::ClipHexahedronInOut2D(vtkPoints *newPoints,
4963 vtkGenericCell *cell,
4964 vtkIncrementalPointLocator *locator,
4965 vtkCellArray **tets,
4966 vtkPointData *inPD,
4967 vtkPointData *outPD,
4968 vtkCellData *inCD,
4969 vtkIdType cellId,
4970 vtkCellData **outCD)
4971 {
4972 vtkIdType cellType = cell->GetCellType();
4973 vtkIdList *cellIds = cell->GetPointIds();
4974 vtkCellArray *arraytriangle = vtkCellArray::New();
4975 vtkPoints *cellPts = cell->GetPoints();
4976 vtkIdType npts = cellPts->GetNumberOfPoints();
4977 std::vector<vtkIdType> cellptId(npts);
4978 vtkIdType iid[3];
4979 vtkIdType *v_id = NULL;
4980 vtkIdType *verts, v1, v2;
4981 vtkIdType ptId;
4982 vtkIdType ptIdout[3];
4983 vtkIdType tab_id[3];
4984 vtkIdType ptstriangle = 3;
4985
4986 int i,j, k;
4987 unsigned int allInside;
4988 unsigned int idtrianglenew;
4989 unsigned int idcellnew;
4990 unsigned int planes;
4991
4992 vtkIdType edges[3][2] = { {0,1}, {1,2}, {2,0}}; /* Edges Triangle */
4993 double deltaScalar;
4994 double t, *p1, *p2;
4995 double v[3],x[3];
4996 double p[6];
4997 double v_triangle[3][3];
4998
4999 // To avoid false positive warning.
5000 tab_id[0]=0;
5001 tab_id[1]=0;
5002 tab_id[2]=0;
5003
5004 for (i=0; i<npts; i++)
5005 {
5006 cellptId[i] = cellIds->GetId(i);
5007 }
5008
5009 // Convert all polygon cells to triangles
5010 this->CellGrid(cellType,npts,&cellptId[0],arraytriangle);
5011
5012 unsigned int totalnewtriangle = arraytriangle->GetNumberOfCells();
5013 for (idtrianglenew = 0 ; idtrianglenew < totalnewtriangle; idtrianglenew++)
5014 {
5015 arraytriangle->GetNextCell(ptstriangle,v_id);
5016
5017 for (allInside=1, i=0; i<3; i++)
5018 {
5019 cellPts->GetPoint(v_id[i],v);
5020
5021 for(k=0;k<6;k++)
5022 {
5023 p[k] = this->PlaneNormal[k][0]*(v[0] - this->PlanePoint[k][0])+
5024 this->PlaneNormal[k][1]*(v[1] - this->PlanePoint[k][1]) +
5025 this->PlaneNormal[k][2]*(v[2] - this->PlanePoint[k][2]);
5026 }
5027
5028 if (!((p[0] <= 0) && (p[1] <= 0) &&
5029 (p[2] <= 0) && (p[3] <= 0) &&
5030 (p[4] <= 0) && (p[5] <= 0)))
5031 {
5032 allInside = 0;
5033 }
5034 }//for all points of the trianglehedron.
5035
5036 // Test Outside
5037 unsigned int test[6] = {1,1,1,1,1,1};
5038 for (i=0; i<3; i++)
5039 {
5040 ptIdout[i] = cellIds->GetId(v_id[i]);
5041 cellPts->GetPoint(v_id[i],v_triangle[i]);
5042
5043 // Use plane equation
5044 for(k=0;k<6;k++)
5045 {
5046 p[k] = this->PlaneNormal[k][0]*(v_triangle[i][0] - this->PlanePoint[k][0])+
5047 this->PlaneNormal[k][1]*(v_triangle[i][1] - this->PlanePoint[k][1]) +
5048 this->PlaneNormal[k][2]*(v_triangle[i][2] - this->PlanePoint[k][2]);
5049 }
5050
5051 for(k=0;k<3;k++)
5052 {
5053 if (p[2*k] <= 0)
5054 {
5055 test[2*k] = 0;
5056 }
5057 if (p[2*k+1] <= 0)
5058 {
5059 test[2*k+1] = 0;
5060 }
5061 }
5062 }//for all points of the cell.
5063
5064 if (!allInside && ((test[0] == 1)|| (test[1] == 1) ||
5065 (test[2] == 1)|| (test[3] == 1) ||
5066 (test[4] == 1)|| (test[5] == 1)))
5067 {
5068 for (i=0; i<3; i++)
5069 {
5070 if ( locator->InsertUniquePoint(v_triangle[i], iid[i]) )
5071 {
5072 outPD->CopyData(inPD,ptIdout[i], iid[i]);
5073 }
5074 }
5075
5076 int newCellId = tets[1]->InsertNextCell(3,iid);
5077 outCD[1]->CopyData(inCD,cellId,newCellId);
5078 continue; // Triangle is outside.
5079 }
5080
5081 for (i=0; i<3; i++)
5082 {
5083 ptId = cellIds->GetId(v_id[i]);
5084 cellPts->GetPoint(v_id[i],v);
5085
5086 // Currently all points are injected because of the possibility
5087 // of intersection point merging.
5088 if ( locator->InsertUniquePoint(v, iid[i]) )
5089 {
5090 outPD->CopyData(inPD,ptId, iid[i]);
5091 }
5092
5093 }//for all points of the trianglehedron.
5094
5095 if ( allInside )
5096 {
5097 int newCellId = tets[0]->InsertNextCell(3,iid); // Tetrahedron inside.
5098 outCD[0]->CopyData(inCD,cellId,newCellId);
5099 continue;
5100 }
5101
5102 double *pedg1,*pedg2;
5103
5104 unsigned int tab2[3][4] = { {1, 2, 1, 0},
5105 {2, 0, 0, 1},
5106 {0, 1, 0, 1} };
5107 unsigned int tab1[3][2] = { {2, 1},
5108 {0, 2},
5109 {1, 0} };
5110
5111 vtkCellArray *cellarray = vtkCellArray::New();
5112 vtkCellArray *cellarrayout = vtkCellArray::New();
5113 int newCellId = cellarray->InsertNextCell(3,iid);
5114
5115 // Test Cell intersection with each plane of box
5116 for (planes = 0; planes < 6; planes++)
5117 {
5118 unsigned int totalnewcells = cellarray->GetNumberOfCells();
5119 vtkCellArray *newcellArray = vtkCellArray::New();
5120
5121 for (idcellnew = 0 ; idcellnew < totalnewcells; idcellnew++)
5122 {
5123 unsigned int i0;
5124 unsigned int num_inter = 0;
5125 unsigned int edges_inter = 0;
5126 vtkIdType p_id[3];
5127
5128 cellarray->GetNextCell(npts,v_id);
5129
5130 newPoints->GetPoint(v_id[0],v_triangle[0]); //coord (x,y,z)
5131 newPoints->GetPoint(v_id[1],v_triangle[1]);
5132 newPoints->GetPoint(v_id[2],v_triangle[2]);
5133
5134 p[0] = this->PlaneNormal[planes][0]*(v_triangle[0][0] - this->PlanePoint[planes][0]) +
5135 this->PlaneNormal[planes][1]*(v_triangle[0][1] - this->PlanePoint[planes][1]) +
5136 this->PlaneNormal[planes][2]*(v_triangle[0][2] - this->PlanePoint[planes][2]);
5137 p[1] = this->PlaneNormal[planes][0]*(v_triangle[1][0] - this->PlanePoint[planes][0]) +
5138 this->PlaneNormal[planes][1]*(v_triangle[1][1] - this->PlanePoint[planes][1]) +
5139 this->PlaneNormal[planes][2]*(v_triangle[1][2] - this->PlanePoint[planes][2]);
5140 p[2] = this->PlaneNormal[planes][0]*(v_triangle[2][0] - this->PlanePoint[planes][0]) +
5141 this->PlaneNormal[planes][1]*(v_triangle[2][1] - this->PlanePoint[planes][1]) +
5142 this->PlaneNormal[planes][2]*(v_triangle[2][2] - this->PlanePoint[planes][2]);
5143
5144 for (int edgeNum=0; edgeNum < 3; edgeNum++)
5145 {
5146 verts = edges[edgeNum];
5147
5148 p1 = v_triangle[verts[0]];
5149 p2 = v_triangle[verts[1]];
5150 double s1 = p[verts[0]];
5151 double s2 = p[verts[1]];
5152 if ( (s1 * s2) < 0)
5153 {
5154 deltaScalar = s2 - s1;
5155
5156 if (deltaScalar > 0)
5157 {
5158 pedg1 = p1; pedg2 = p2;
5159 v1 = verts[0]; v2 = verts[1];
5160 }
5161 else
5162 {
5163 pedg1 = p2; pedg2 = p1;
5164 v1 = verts[1]; v2 = verts[0];
5165 deltaScalar = -deltaScalar;
5166 t = s1; s1 = s2; s2 = t;
5167 }
5168
5169 // linear interpolation
5170 t = ( deltaScalar == 0.0 ? 0.0 : ( - s1) / deltaScalar );
5171
5172 for (j=0; j<3; j++)
5173 {
5174 x[j] = pedg1[j] + t*(pedg2[j] - pedg1[j]);
5175 }
5176
5177 // Incorporate point into output and interpolate edge data as necessary
5178 edges_inter = edges_inter * 10 + (edgeNum+1);
5179
5180 if ( locator->InsertUniquePoint(x, p_id[num_inter]) )
5181 {
5182 this->InterpolateEdge(outPD, p_id[num_inter],
5183 v_id[v1], v_id[v2], t);
5184 }
5185
5186 num_inter++;
5187 }//if edge intersects value
5188 }//for all edges
5189
5190 if (num_inter == 0)
5191 {
5192 unsigned int outside = 0;
5193 for(i=0;i<3;i++)
5194 {
5195 if (p[i] > 0) // If only one vertex is ouside, so the trianglehedron is outside
5196 {
5197 outside = 1; // because there is not intersection.
5198 break; // some vertex could be on plane, so you need to test all vertex
5199 }
5200 }
5201 if (outside == 0)
5202 {
5203 // else it is possible intersection if other plane
5204 newCellId = newcellArray->InsertNextCell(3,v_id);
5205 }
5206 else
5207 {
5208 newCellId = tets[1]->InsertNextCell(3,v_id);
5209 outCD[1]->CopyData(inCD,cellId,newCellId);
5210 }
5211 continue;
5212 }
5213 switch(num_inter)
5214 {
5215 case 2: // We have one quad and one triangle
5216 switch(edges_inter)
5217 {
5218 case 12:
5219 i0 = 1;
5220 break;
5221 case 23:
5222 i0 = 2;
5223 break;
5224 case 13:
5225 i0 = 0;
5226 break;
5227 default:
5228 vtkErrorMacro(<< "Intersection not found: Num_inter = " <<
5229 num_inter << " Edges_inter = " << edges_inter);
5230 continue;
5231 }
5232
5233 if (p[i0] > 0)
5234 {
5235 // The v_triangle[3] is outside box, so
5236 // the first wedge is outside
5237
5238 // The Quad is inside: two triangles: (v0,v1,p0) and (p0,p1,v1)
5239 tab_id[0] = v_id[tab2[i0][0]];
5240 tab_id[1] = v_id[tab2[i0][1]];
5241 tab_id[2] = p_id[tab2[i0][2]];
5242 newCellId = newcellArray->InsertNextCell(3,tab_id);
5243 tab_id[0] = p_id[tab2[i0][2]];
5244 tab_id[1] = p_id[tab2[i0][3]];
5245 tab_id[2] = v_id[tab2[i0][0]];
5246 newCellId = newcellArray->InsertNextCell(3,tab_id);
5247
5248
5249 switch (edges_inter) // Triangle Outside
5250 {
5251 case 12:
5252 case 23:
5253 tab_id[0] = v_id[i0];
5254 tab_id[1] = p_id[1];
5255 tab_id[2] = p_id[0];
5256 break;
5257 case 13:
5258 tab_id[0] = v_id[i0];
5259 tab_id[1] = p_id[0];
5260 tab_id[2] = p_id[1];
5261 break;
5262 }
5263 newCellId = cellarrayout->InsertNextCell(3,tab_id);
5264 }
5265 else
5266 {
5267 // The Triangle is inside: (v0,p0,p1)
5268 switch (edges_inter)
5269 {
5270 case 12:
5271 case 23:
5272 tab_id[0] = v_id[i0];
5273 tab_id[1] = p_id[1];
5274 tab_id[2] = p_id[0];
5275 break;
5276 case 13:
5277 tab_id[0] = v_id[i0];
5278 tab_id[1] = p_id[0];
5279 tab_id[2] = p_id[1];
5280 break;
5281 }
5282 newCellId = newcellArray->InsertNextCell(3,tab_id);
5283
5284 // The Quad is outside: two triangles: (v0,v1,p0) and (p0,p1,v1)
5285 tab_id[0] = v_id[tab2[i0][0]];
5286 tab_id[1] = v_id[tab2[i0][1]];
5287 tab_id[2] = p_id[tab2[i0][2]];
5288 newCellId = cellarrayout->InsertNextCell(3,tab_id);
5289 tab_id[0] = p_id[tab2[i0][2]];
5290 tab_id[1] = p_id[tab2[i0][3]];
5291 tab_id[2] = v_id[tab2[i0][0]];
5292 newCellId = cellarrayout->InsertNextCell(3,tab_id);
5293 }
5294 break;
5295
5296 case 1: // We have two triangles
5297 switch(edges_inter)
5298 {
5299 case 1:
5300 i0 = 0;
5301 break;
5302 case 2:
5303 i0 = 1;
5304 break;
5305 case 3:
5306 i0 = 2;
5307 break;
5308 default:
5309 vtkErrorMacro(<< "Intersection not found: Num_inter = " <<
5310 num_inter << " Edges_inter = " << edges_inter);
5311 continue;
5312 }
5313 if (p[i0] > 0)
5314 {
5315 // Isolate vertex is outside box, so
5316 // the triangle is outside
5317 tab_id[0] = v_id[tab1[i0][1]]; // Inside
5318 tab_id[1] = v_id[tab1[i0][0]];
5319 tab_id[2] = p_id[0];
5320 newCellId = newcellArray->InsertNextCell(3,tab_id);
5321
5322 tab_id[0] = v_id[tab1[i0][0]]; // Outside
5323 tab_id[1] = v_id[i0];
5324 tab_id[2] = p_id[0];
5325 newCellId = cellarrayout->InsertNextCell(3,tab_id);
5326 }
5327 else
5328 {
5329 tab_id[0] = v_id[tab1[i0][0]]; // Inside
5330 tab_id[1] = v_id[i0];
5331 tab_id[2] = p_id[0];
5332 newCellId = newcellArray->InsertNextCell(3,tab_id);
5333
5334 tab_id[0] = v_id[tab1[i0][1]]; // Outside
5335 tab_id[1] = v_id[tab1[i0][0]];
5336 tab_id[2] = p_id[0];
5337 newCellId = cellarrayout->InsertNextCell(3,tab_id);
5338 }
5339 break;
5340 }
5341 } // for all new cells
5342 cellarray->Delete();
5343 cellarray = newcellArray;
5344 } //for all planes
5345
5346 unsigned int totalnewcells = cellarray->GetNumberOfCells(); // Inside
5347
5348 for (idcellnew = 0 ; idcellnew < totalnewcells; idcellnew++)
5349 {
5350 cellarray->GetNextCell(npts,v_id);
5351 newCellId = tets[0]->InsertNextCell(npts,v_id);
5352 outCD[0]->CopyData(inCD,cellId,newCellId);
5353 }
5354 cellarray->Delete();
5355
5356 totalnewcells = cellarrayout->GetNumberOfCells(); // Outside
5357
5358 for (idcellnew = 0 ; idcellnew < totalnewcells; idcellnew++)
5359 {
5360 cellarrayout->GetNextCell(npts,v_id);
5361 newCellId = tets[1]->InsertNextCell(npts,v_id);
5362 outCD[1]->CopyData(inCD,cellId,newCellId);
5363 }
5364 cellarrayout->Delete();
5365 }
5366 arraytriangle->Delete();
5367 }
5368
5369 //-----------------------------------------------------------------------------
5370
ClipBox1D(vtkPoints * newPoints,vtkGenericCell * cell,vtkIncrementalPointLocator * locator,vtkCellArray * lines,vtkPointData * inPD,vtkPointData * outPD,vtkCellData * inCD,vtkIdType cellId,vtkCellData * outCD)5371 void vtkBoxClipDataSet::ClipBox1D(vtkPoints *newPoints,
5372 vtkGenericCell *cell,
5373 vtkIncrementalPointLocator *locator,
5374 vtkCellArray *lines,
5375 vtkPointData *inPD,
5376 vtkPointData *outPD,
5377 vtkCellData *inCD,
5378 vtkIdType cellId,
5379 vtkCellData *outCD)
5380 {
5381 vtkIdType cellType = cell->GetCellType();
5382 vtkIdList *cellIds = cell->GetPointIds();
5383 vtkCellArray *arrayline = vtkCellArray::New();
5384 vtkPoints *cellPts = cell->GetPoints();
5385 vtkIdType npts = cellPts->GetNumberOfPoints();
5386 std::vector<vtkIdType> cellptId(npts);
5387 vtkIdType iid[2];
5388 vtkIdType *v_id = NULL;
5389 vtkIdType ptId;
5390 vtkIdType tab_id[2];
5391 vtkIdType ptsline = 2;
5392
5393 int i,j;
5394 unsigned int allInside;
5395 unsigned int planes;
5396 unsigned int cutInd;
5397
5398 double value;
5399 double t;
5400 double v[3],x[3];
5401 double v_line[2][3];
5402
5403 for (i = 0; i < npts; i++)
5404 {
5405 cellptId[i] = cellIds->GetId(i);
5406 }
5407
5408 // Convert all 1d cells to single line.
5409 this->CellGrid(cellType, npts, &cellptId[0], arrayline);
5410
5411 unsigned int totalnewline = arrayline->GetNumberOfCells();
5412 for (unsigned int idlinenew = 0; idlinenew < totalnewline; idlinenew++)
5413 {
5414 arrayline->GetNextCell(ptsline, v_id);
5415
5416 for (allInside=1, i=0; i<2; i++)
5417 {
5418 cellPts->GetPoint(v_id[i],v);
5419
5420 if (!( (v[0] >= this->BoundBoxClip[0][0])
5421 && (v[0] <= this->BoundBoxClip[0][1])
5422 && (v[1] >= this->BoundBoxClip[1][0])
5423 && (v[1] <= this->BoundBoxClip[1][1])
5424 && (v[2] >= this->BoundBoxClip[2][0])
5425 && (v[2] <= this->BoundBoxClip[2][1]) ))
5426 {
5427 //outside
5428 allInside = 0;
5429 }
5430 }
5431
5432 // Test Outside:
5433 if (!allInside)
5434 {
5435 unsigned int test[6] = {1,1,1,1,1,1};
5436 for (i=0; i<2; i++)
5437 {
5438 cellPts->GetPoint(v_id[i],v);
5439
5440 if (v[0] >= this->BoundBoxClip[0][0])
5441 {
5442 test[0] = 0;
5443 }
5444 if (v[0] <= this->BoundBoxClip[0][1])
5445 {
5446 test[1] = 0;
5447 }
5448 if (v[1] >= this->BoundBoxClip[1][0])
5449 {
5450 test[2] = 0;
5451 }
5452 if (v[1] <= this->BoundBoxClip[1][1])
5453 {
5454 test[3] = 0;
5455 }
5456 if (v[2] >= this->BoundBoxClip[2][0])
5457 {
5458 test[4] = 0;
5459 }
5460 if (v[2] <= this->BoundBoxClip[2][1])
5461 {
5462 test[5] = 0;
5463 }
5464 }//for all points of the line.
5465
5466 if ((test[0] == 1)|| (test[1] == 1) ||
5467 (test[2] == 1)|| (test[3] == 1) ||
5468 (test[4] == 1)|| (test[5] == 1))
5469 {
5470 continue; // Line is outside.
5471 }
5472 }//if not allInside
5473
5474 for (i=0; i<2; i++)
5475 {
5476 // Currently all points are injected because of the possibility
5477 // of intersection point merging.
5478 ptId = cellIds->GetId(v_id[i]);
5479 cellPts->GetPoint(v_id[i],v);
5480 if ( locator->InsertUniquePoint(v, iid[i]) )
5481 {
5482 outPD->CopyData(inPD,ptId, iid[i]);
5483 }
5484 }//for all points of the triangle.
5485
5486 if ( allInside )
5487 {
5488 // Triangle inside.
5489 int newCellId = lines->InsertNextCell(2,iid);
5490 outCD->CopyData(inCD,cellId,newCellId);
5491 continue;
5492 }
5493
5494 vtkCellArray *cellarray = vtkCellArray::New();
5495 int newCellId = cellarray->InsertNextCell(2,iid);
5496 unsigned int idcellnew;
5497
5498 // Test triangle intersection with each plane of box
5499 for (planes = 0; planes < 6; planes++)
5500 {
5501 // The index of the dimension of the cut plane (x == 0, y == 1, z == 2).
5502 cutInd = planes/2;
5503
5504 // The plane is always parallel to unitary cube.
5505 value = this->BoundBoxClip[cutInd][planes%2];
5506
5507 unsigned int totalnewcells = cellarray->GetNumberOfCells();
5508 vtkCellArray *newcellArray = vtkCellArray::New();
5509
5510 for (idcellnew = 0 ; idcellnew < totalnewcells; idcellnew++)
5511 {
5512 vtkIdType p_id;
5513 cellarray->GetNextCell(npts, v_id);
5514
5515 newPoints->GetPoint(v_id[0], v_line[0]); //coord (x,y,z)
5516 newPoints->GetPoint(v_id[1], v_line[1]);
5517
5518 // Check to see if line is inside plane.
5519 if ( ( (planes%2 == 0)
5520 && (v_line[0][cutInd] >= value)
5521 && (v_line[1][cutInd] >= value) )
5522 || ( (planes%2 == 1)
5523 && (v_line[0][cutInd] <= value)
5524 && (v_line[1][cutInd] <= value) ) )
5525 {
5526 newCellId = newcellArray->InsertNextCell(2, v_id);
5527 continue;
5528 }
5529
5530 // Check to see if line is outside plane.
5531 if ( ( (planes%2 == 0)
5532 && (v_line[0][cutInd] <= value)
5533 && (v_line[1][cutInd] <= value) )
5534 || ( (planes%2 == 1)
5535 && (v_line[0][cutInd] >= value)
5536 && (v_line[1][cutInd] >= value) ) )
5537 {
5538 continue;
5539 }
5540
5541 // If we are here, the plane intersects the line segment.
5542 t = (value - v_line[0][cutInd])/(v_line[1][cutInd] - v_line[0][cutInd]);
5543 for (j = 0; j < 3; j++)
5544 {
5545 x[j] = (v_line[1][j] - v_line[0][j])*t + v_line[0][j];
5546 }
5547
5548 // Incorporate point into output and interpolate edge data as
5549 // necessary.
5550 if (locator->InsertUniquePoint(x, p_id))
5551 {
5552 this->InterpolateEdge(outPD, p_id, v_id[0], v_id[0], t);
5553 }
5554
5555 // Add the clipped line to the output.
5556 if ( ((planes%2 == 0) && (v_line[0][cutInd] >= value))
5557 || ((planes%2 == 1) && (v_line[0][cutInd] <= value)) )
5558 {
5559 // First point of line is inside.
5560 tab_id[0] = v_id[0];
5561 tab_id[1] = p_id;
5562 newcellArray->InsertNextCell(2, tab_id);
5563 }
5564 else
5565 {
5566 // Second point of line is inside.
5567 tab_id[0] = p_id;
5568 tab_id[1] = v_id[1];
5569 newcellArray->InsertNextCell(2, tab_id);
5570 }
5571 } // for all new cells
5572 cellarray->Delete();
5573 cellarray = newcellArray;
5574 } // for all planes
5575
5576 unsigned int totalnewcells = cellarray->GetNumberOfCells();
5577
5578 for (idcellnew = 0 ; idcellnew < totalnewcells; idcellnew++)
5579 {
5580 cellarray->GetNextCell(npts,v_id);
5581 newCellId = lines->InsertNextCell(npts,v_id);
5582 outCD->CopyData(inCD,cellId,newCellId);
5583 }
5584 cellarray->Delete();
5585 }
5586 arrayline->Delete();
5587 }
5588
5589 //-----------------------------------------------------------------------------
5590
ClipBoxInOut1D(vtkPoints * newPoints,vtkGenericCell * cell,vtkIncrementalPointLocator * locator,vtkCellArray ** lines,vtkPointData * inPD,vtkPointData * outPD,vtkCellData * inCD,vtkIdType cellId,vtkCellData ** outCD)5591 void vtkBoxClipDataSet::ClipBoxInOut1D(vtkPoints *newPoints,
5592 vtkGenericCell *cell,
5593 vtkIncrementalPointLocator *locator,
5594 vtkCellArray **lines,
5595 vtkPointData *inPD,
5596 vtkPointData *outPD,
5597 vtkCellData *inCD,
5598 vtkIdType cellId,
5599 vtkCellData **outCD)
5600 {
5601 vtkIdType cellType = cell->GetCellType();
5602 vtkIdList *cellIds = cell->GetPointIds();
5603 vtkCellArray *arrayline = vtkCellArray::New();
5604 vtkPoints *cellPts = cell->GetPoints();
5605 vtkIdType npts = cellPts->GetNumberOfPoints();
5606 std::vector<vtkIdType> cellptId(npts);
5607 vtkIdType iid[2];
5608 vtkIdType *v_id = NULL;
5609 vtkIdType ptId;
5610 vtkIdType tab_id[2];
5611 vtkIdType ptsline = 2;
5612
5613 int i,j;
5614 unsigned int allInside;
5615 unsigned int planes;
5616 unsigned int cutInd;
5617
5618 double value;
5619 double t;
5620 double v[3],x[3];
5621 double v_line[2][3];
5622
5623 for (i = 0; i < npts; i++)
5624 {
5625 cellptId[i] = cellIds->GetId(i);
5626 }
5627
5628 // Convert all 1d cells to single line.
5629 this->CellGrid(cellType, npts, &cellptId[0], arrayline);
5630
5631 unsigned int totalnewline = arrayline->GetNumberOfCells();
5632 for (unsigned int idlinenew = 0; idlinenew < totalnewline; idlinenew++)
5633 {
5634 arrayline->GetNextCell(ptsline, v_id);
5635
5636 for (allInside=1, i=0; i<2; i++)
5637 {
5638 ptId = cellIds->GetId(v_id[i]);
5639 cellPts->GetPoint(v_id[i],v);
5640
5641 if (!( (v[0] >= this->BoundBoxClip[0][0])
5642 && (v[0] <= this->BoundBoxClip[0][1])
5643 && (v[1] >= this->BoundBoxClip[1][0])
5644 && (v[1] <= this->BoundBoxClip[1][1])
5645 && (v[2] >= this->BoundBoxClip[2][0])
5646 && (v[2] <= this->BoundBoxClip[2][1]) ))
5647 {
5648 //outside
5649 allInside = 0;
5650 }
5651
5652 if ( locator->InsertUniquePoint(v, iid[i]) )
5653 {
5654 outPD->CopyData(inPD, ptId, iid[i]);
5655 }
5656 }
5657
5658 // Test Outside:
5659 if (!allInside)
5660 {
5661 unsigned int test[6] = {1,1,1,1,1,1};
5662 for (i=0; i<2; i++)
5663 {
5664 cellPts->GetPoint(v_id[i],v);
5665
5666 if (v[0] >= this->BoundBoxClip[0][0])
5667 {
5668 test[0] = 0;
5669 }
5670 if (v[0] <= this->BoundBoxClip[0][1])
5671 {
5672 test[1] = 0;
5673 }
5674 if (v[1] >= this->BoundBoxClip[1][0])
5675 {
5676 test[2] = 0;
5677 }
5678 if (v[1] <= this->BoundBoxClip[1][1])
5679 {
5680 test[3] = 0;
5681 }
5682 if (v[2] >= this->BoundBoxClip[2][0])
5683 {
5684 test[4] = 0;
5685 }
5686 if (v[2] <= this->BoundBoxClip[2][1])
5687 {
5688 test[5] = 0;
5689 }
5690 }//for all points of the line.
5691
5692 if ((test[0] == 1)|| (test[1] == 1) ||
5693 (test[2] == 1)|| (test[3] == 1) ||
5694 (test[4] == 1)|| (test[5] == 1))
5695 {
5696 int newCellId = lines[1]->InsertNextCell(2, iid);
5697 outCD[1]->CopyData(inCD, cellId, newCellId);
5698 continue; // Line is outside.
5699 }
5700 }//if not allInside
5701
5702 if ( allInside )
5703 {
5704 // Triangle inside.
5705 int newCellId = lines[0]->InsertNextCell(2,iid);
5706 outCD[0]->CopyData(inCD,cellId,newCellId);
5707 continue;
5708 }
5709
5710 vtkCellArray *cellarray = vtkCellArray::New();
5711 vtkCellArray *cellarrayout = vtkCellArray::New();
5712 int newCellId = cellarray->InsertNextCell(2,iid);
5713 unsigned int idcellnew;
5714
5715 // Test triangle intersection with each plane of box
5716 for (planes = 0; planes < 6; planes++)
5717 {
5718 // The index of the dimension of the cut plane (x == 0, y == 1, z == 2).
5719 cutInd = planes/2;
5720
5721 // The plane is always parallel to unitary cube.
5722 value = this->BoundBoxClip[cutInd][planes%2];
5723
5724 unsigned int totalnewcells = cellarray->GetNumberOfCells();
5725 vtkCellArray *newcellArray = vtkCellArray::New();
5726
5727 for (idcellnew = 0 ; idcellnew < totalnewcells; idcellnew++)
5728 {
5729 vtkIdType p_id;
5730 cellarray->GetNextCell(npts, v_id);
5731
5732 newPoints->GetPoint(v_id[0], v_line[0]); //coord (x,y,z)
5733 newPoints->GetPoint(v_id[1], v_line[1]);
5734
5735 // Check to see if line is inside plane.
5736 if ( ( (planes%2 == 0)
5737 && (v_line[0][cutInd] >= value)
5738 && (v_line[1][cutInd] >= value) )
5739 || ( (planes%2 == 1)
5740 && (v_line[0][cutInd] <= value)
5741 && (v_line[1][cutInd] <= value) ) )
5742 {
5743 newCellId = newcellArray->InsertNextCell(2, v_id);
5744 continue;
5745 }
5746
5747 // Check to see if line is outside plane.
5748 if ( ( (planes%2 == 0)
5749 && (v_line[0][cutInd] <= value)
5750 && (v_line[1][cutInd] <= value) )
5751 || ( (planes%2 == 1)
5752 && (v_line[0][cutInd] >= value)
5753 && (v_line[1][cutInd] >= value) ) )
5754 {
5755 newCellId = lines[1]->InsertNextCell(2, v_id);
5756 outCD[1]->CopyData(inCD, cellId, newCellId);
5757 continue;
5758 }
5759
5760 // If we are here, the plane intersects the line segment.
5761 t = (value - v_line[0][cutInd])/(v_line[1][cutInd] - v_line[0][cutInd]);
5762 for (j = 0; j < 3; j++)
5763 {
5764 x[j] = (v_line[1][j] - v_line[0][j])*t + v_line[0][j];
5765 }
5766
5767 // Incorporate point into output and interpolate edge data as
5768 // necessary.
5769 if (locator->InsertUniquePoint(x, p_id))
5770 {
5771 this->InterpolateEdge(outPD, p_id, v_id[0], v_id[0], t);
5772 }
5773
5774 // Add the clipped line to the output.
5775 if ( ((planes%2 == 0) && (v_line[0][cutInd] >= value))
5776 || ((planes%2 == 1) && (v_line[0][cutInd] <= value)) )
5777 {
5778 // First point of line is inside.
5779 tab_id[0] = v_id[0];
5780 tab_id[1] = p_id;
5781 newcellArray->InsertNextCell(2, tab_id);
5782
5783 // Second point of line is outside.
5784 tab_id[0] = p_id;
5785 tab_id[1] = v_id[1];
5786 cellarrayout->InsertNextCell(2, tab_id);
5787 }
5788 else
5789 {
5790 // Second point of line is inside.
5791 tab_id[0] = p_id;
5792 tab_id[1] = v_id[1];
5793 newcellArray->InsertNextCell(2, tab_id);
5794
5795 // First point of line is outside.
5796 tab_id[0] = v_id[0];
5797 tab_id[1] = p_id;
5798 cellarrayout->InsertNextCell(2, tab_id);
5799 }
5800 } // for all new cells
5801 cellarray->Delete();
5802 cellarray = newcellArray;
5803 } // for all planes
5804
5805 unsigned int totalnewcells = cellarray->GetNumberOfCells(); // Inside
5806 for (idcellnew = 0 ; idcellnew < totalnewcells; idcellnew++)
5807 {
5808 cellarray->GetNextCell(npts, v_id);
5809 newCellId = lines[0]->InsertNextCell(npts, v_id);
5810 outCD[0]->CopyData(inCD, cellId, newCellId);
5811 }
5812 cellarray->Delete();
5813
5814 totalnewcells = cellarrayout->GetNumberOfCells(); // Outside
5815 for (idcellnew = 0; idcellnew < totalnewcells; idcellnew++)
5816 {
5817 cellarrayout->GetNextCell(npts, v_id);
5818 newCellId = lines[1]->InsertNextCell(npts, v_id);
5819 outCD[1]->CopyData(inCD, cellId, newCellId);
5820 }
5821 cellarrayout->Delete();
5822 }
5823 arrayline->Delete();
5824 }
5825
5826 //-----------------------------------------------------------------------------
5827
ClipHexahedron1D(vtkPoints * newPoints,vtkGenericCell * cell,vtkIncrementalPointLocator * locator,vtkCellArray * lines,vtkPointData * inPD,vtkPointData * outPD,vtkCellData * inCD,vtkIdType cellId,vtkCellData * outCD)5828 void vtkBoxClipDataSet::ClipHexahedron1D(vtkPoints *newPoints,
5829 vtkGenericCell *cell,
5830 vtkIncrementalPointLocator *locator,
5831 vtkCellArray *lines,
5832 vtkPointData *inPD,
5833 vtkPointData *outPD,
5834 vtkCellData *inCD,
5835 vtkIdType cellId,
5836 vtkCellData *outCD)
5837 {
5838 vtkIdType cellType = cell->GetCellType();
5839 vtkIdList *cellIds = cell->GetPointIds();
5840 vtkCellArray *arrayline = vtkCellArray::New();
5841 vtkPoints *cellPts = cell->GetPoints();
5842 vtkIdType npts = cellPts->GetNumberOfPoints();
5843 std::vector<vtkIdType> cellptId(npts);
5844 vtkIdType iid[2];
5845 vtkIdType *v_id = NULL;
5846 vtkIdType ptId;
5847 vtkIdType tab_id[2];
5848 vtkIdType ptsline = 2;
5849
5850 int i,j;
5851 unsigned int allInside;
5852 unsigned int planes;
5853
5854 double values[2];
5855 double t;
5856 double v[3],x[3];
5857 double v_line[2][3];
5858
5859 for (i = 0; i < npts; i++)
5860 {
5861 cellptId[i] = cellIds->GetId(i);
5862 }
5863
5864 // Convert all 1d cells to single line.
5865 this->CellGrid(cellType, npts, &cellptId[0], arrayline);
5866
5867 unsigned int totalnewline = arrayline->GetNumberOfCells();
5868 for (unsigned int idlinenew = 0; idlinenew < totalnewline; idlinenew++)
5869 {
5870 arrayline->GetNextCell(ptsline, v_id);
5871
5872 for (allInside=1, i=0; i<2; i++)
5873 {
5874 cellPts->GetPoint(v_id[i],v);
5875
5876 for (int k = 0; k < 6; k++)
5877 {
5878 values[0] = ( this->PlaneNormal[k][0]*(v[0] - this->PlanePoint[k][0])
5879 + this->PlaneNormal[k][1]*(v[1] - this->PlanePoint[k][1])
5880 + this->PlaneNormal[k][2]*(v[2] - this->PlanePoint[k][2]));
5881 if (values[0] > 0)
5882 {
5883 //outside
5884 allInside = 0;
5885 }
5886 }
5887 }
5888
5889 // Test Outside:
5890 if (!allInside)
5891 {
5892 unsigned int test[6] = {1,1,1,1,1,1};
5893 for (i=0; i<2; i++)
5894 {
5895 cellPts->GetPoint(v_id[i],v);
5896
5897 // Use plane equation.
5898 for (int k = 0; k < 6; k++)
5899 {
5900 values[0]
5901 = ( this->PlaneNormal[k][0]*(v[0] - this->PlanePoint[k][0])
5902 + this->PlaneNormal[k][1]*(v[1] - this->PlanePoint[k][1])
5903 + this->PlaneNormal[k][2]*(v[2] - this->PlanePoint[k][2]) );
5904 if (values[0] <= 0)
5905 {
5906 test[k] = 0;
5907 }
5908 }//for all planes of the hexahedron.
5909 }//for all points of the line.
5910
5911 if ((test[0] == 1)|| (test[1] == 1) ||
5912 (test[2] == 1)|| (test[3] == 1) ||
5913 (test[4] == 1)|| (test[5] == 1))
5914 {
5915 continue; // Line is outside.
5916 }
5917 }//if not allInside
5918
5919 for (i=0; i<2; i++)
5920 {
5921 // Currently all points are injected because of the possibility
5922 // of intersection point merging.
5923 ptId = cellIds->GetId(v_id[i]);
5924 cellPts->GetPoint(v_id[i],v);
5925 if ( locator->InsertUniquePoint(v, iid[i]) )
5926 {
5927 outPD->CopyData(inPD,ptId, iid[i]);
5928 }
5929 }//for all points of the triangle.
5930
5931 if ( allInside )
5932 {
5933 // Triangle inside.
5934 int newCellId = lines->InsertNextCell(2,iid);
5935 outCD->CopyData(inCD,cellId,newCellId);
5936 continue;
5937 }
5938
5939 vtkCellArray *cellarray = vtkCellArray::New();
5940 int newCellId = cellarray->InsertNextCell(2,iid);
5941 unsigned int idcellnew;
5942
5943 // Test triangle intersection with each plane of box
5944 for (planes = 0; planes < 6; planes++)
5945 {
5946 unsigned int totalnewcells = cellarray->GetNumberOfCells();
5947 vtkCellArray *newcellArray = vtkCellArray::New();
5948
5949 for (idcellnew = 0 ; idcellnew < totalnewcells; idcellnew++)
5950 {
5951 vtkIdType p_id;
5952 cellarray->GetNextCell(npts, v_id);
5953
5954 newPoints->GetPoint(v_id[0], v_line[0]); //coord (x,y,z)
5955 newPoints->GetPoint(v_id[1], v_line[1]);
5956
5957 const double *planeNormal = this->PlaneNormal[planes];
5958 const double *planePoint = this->PlanePoint[planes];
5959 values[0] = ( planeNormal[0]*(v_line[0][0] - planePoint[0])
5960 + planeNormal[1]*(v_line[0][1] - planePoint[1])
5961 + planeNormal[2]*(v_line[0][2] - planePoint[2]));
5962 values[1] = ( planeNormal[0]*(v_line[1][0] - planePoint[0])
5963 + planeNormal[1]*(v_line[1][1] - planePoint[1])
5964 + planeNormal[2]*(v_line[1][2] - planePoint[2]));
5965
5966 // Check to see if line is inside plane.
5967 if ((values[0] <= 0) && (values[1] <= 0))
5968 {
5969 newCellId = newcellArray->InsertNextCell(2, v_id);
5970 continue;
5971 }
5972
5973 // Check to see if line is outside plane.
5974 if ((values[0] >= 0) && (values[1] >= 0))
5975 {
5976 continue;
5977 }
5978
5979 // If we are here, the plane intersects the line segment.
5980 t = values[0]/(values[0] - values[1]);
5981 for (j = 0; j < 3; j++)
5982 {
5983 x[j] = (v_line[1][j] - v_line[0][j])*t + v_line[0][j];
5984 }
5985
5986 // Incorporate point into output and interpolate edge data as
5987 // necessary.
5988 if (locator->InsertUniquePoint(x, p_id))
5989 {
5990 this->InterpolateEdge(outPD, p_id, v_id[0], v_id[0], t);
5991 }
5992
5993 // Add the clipped line to the output.
5994 if (values[0] <= 0)
5995 {
5996 // First point of line is inside.
5997 tab_id[0] = v_id[0];
5998 tab_id[1] = p_id;
5999 newcellArray->InsertNextCell(2, tab_id);
6000 }
6001 else
6002 {
6003 // Second point of line is inside.
6004 tab_id[0] = p_id;
6005 tab_id[1] = v_id[1];
6006 newcellArray->InsertNextCell(2, tab_id);
6007 }
6008 } // for all new cells
6009 cellarray->Delete();
6010 cellarray = newcellArray;
6011 } // for all planes
6012
6013 unsigned int totalnewcells = cellarray->GetNumberOfCells();
6014
6015 for (idcellnew = 0 ; idcellnew < totalnewcells; idcellnew++)
6016 {
6017 cellarray->GetNextCell(npts,v_id);
6018 newCellId = lines->InsertNextCell(npts,v_id);
6019 outCD->CopyData(inCD,cellId,newCellId);
6020 }
6021 cellarray->Delete();
6022 }
6023 arrayline->Delete();
6024 }
6025
6026 //-----------------------------------------------------------------------------
6027
ClipHexahedronInOut1D(vtkPoints * newPoints,vtkGenericCell * cell,vtkIncrementalPointLocator * locator,vtkCellArray ** lines,vtkPointData * inPD,vtkPointData * outPD,vtkCellData * inCD,vtkIdType cellId,vtkCellData ** outCD)6028 void vtkBoxClipDataSet::ClipHexahedronInOut1D(vtkPoints *newPoints,
6029 vtkGenericCell *cell,
6030 vtkIncrementalPointLocator *locator,
6031 vtkCellArray **lines,
6032 vtkPointData *inPD,
6033 vtkPointData *outPD,
6034 vtkCellData *inCD,
6035 vtkIdType cellId,
6036 vtkCellData **outCD)
6037 {
6038 vtkIdType cellType = cell->GetCellType();
6039 vtkIdList *cellIds = cell->GetPointIds();
6040 vtkCellArray *arrayline = vtkCellArray::New();
6041 vtkPoints *cellPts = cell->GetPoints();
6042 vtkIdType npts = cellPts->GetNumberOfPoints();
6043 std::vector<vtkIdType> cellptId(npts);
6044 vtkIdType iid[2];
6045 vtkIdType *v_id = NULL;
6046 vtkIdType ptId;
6047 vtkIdType tab_id[2];
6048 vtkIdType ptsline = 2;
6049
6050 int i,j;
6051 unsigned int allInside;
6052 unsigned int planes;
6053
6054 double values[2];
6055 double t;
6056 double v[3],x[3];
6057 double v_line[2][3];
6058
6059 for (i = 0; i < npts; i++)
6060 {
6061 cellptId[i] = cellIds->GetId(i);
6062 }
6063
6064 // Convert all 1d cells to single line.
6065 this->CellGrid(cellType, npts, &cellptId[0], arrayline);
6066
6067 unsigned int totalnewline = arrayline->GetNumberOfCells();
6068 for (unsigned int idlinenew = 0; idlinenew < totalnewline; idlinenew++)
6069 {
6070 arrayline->GetNextCell(ptsline, v_id);
6071
6072 for (allInside=1, i=0; i<2; i++)
6073 {
6074 ptId = cellIds->GetId(v_id[i]);
6075 cellPts->GetPoint(v_id[i],v);
6076
6077 for (int k = 0; k < 6; k++)
6078 {
6079 values[0] = ( this->PlaneNormal[k][0]*(v[0] - this->PlanePoint[k][0])
6080 + this->PlaneNormal[k][1]*(v[1] - this->PlanePoint[k][1])
6081 + this->PlaneNormal[k][2]*(v[2] - this->PlanePoint[k][2]));
6082 if (values[0] > 0)
6083 {
6084 //outside
6085 allInside = 0;
6086 }
6087 }
6088
6089 if ( locator->InsertUniquePoint(v, iid[i]) )
6090 {
6091 outPD->CopyData(inPD, ptId, iid[i]);
6092 }
6093 }
6094
6095 // Test Outside:
6096 if (!allInside)
6097 {
6098 unsigned int test[6] = {1,1,1,1,1,1};
6099 for (i=0; i<2; i++)
6100 {
6101 cellPts->GetPoint(v_id[i],v);
6102
6103 // Use plane equation.
6104 for (int k = 0; k < 6; k++)
6105 {
6106 values[0]
6107 = ( this->PlaneNormal[k][0]*(v[0] - this->PlanePoint[k][0])
6108 + this->PlaneNormal[k][1]*(v[1] - this->PlanePoint[k][1])
6109 + this->PlaneNormal[k][2]*(v[2] - this->PlanePoint[k][2]) );
6110 if (values[0] <= 0)
6111 {
6112 test[k] = 0;
6113 }
6114 }//for all planes of the hexahedron.
6115 }//for all points of the line.
6116
6117 if ((test[0] == 1)|| (test[1] == 1) ||
6118 (test[2] == 1)|| (test[3] == 1) ||
6119 (test[4] == 1)|| (test[5] == 1))
6120 {
6121 int newCellId = lines[1]->InsertNextCell(2, iid);
6122 outCD[1]->CopyData(inCD, cellId, newCellId);
6123 continue; // Line is outside.
6124 }
6125 }//if not allInside
6126
6127 if ( allInside )
6128 {
6129 // Triangle inside.
6130 int newCellId = lines[0]->InsertNextCell(2,iid);
6131 outCD[0]->CopyData(inCD,cellId,newCellId);
6132 continue;
6133 }
6134
6135 vtkCellArray *cellarray = vtkCellArray::New();
6136 vtkCellArray *cellarrayout = vtkCellArray::New();
6137 int newCellId = cellarray->InsertNextCell(2,iid);
6138 unsigned int idcellnew;
6139
6140 // Test triangle intersection with each plane of box
6141 for (planes = 0; planes < 6; planes++)
6142 {
6143 unsigned int totalnewcells = cellarray->GetNumberOfCells();
6144 vtkCellArray *newcellArray = vtkCellArray::New();
6145
6146 for (idcellnew = 0 ; idcellnew < totalnewcells; idcellnew++)
6147 {
6148 vtkIdType p_id;
6149 cellarray->GetNextCell(npts, v_id);
6150
6151 newPoints->GetPoint(v_id[0], v_line[0]); //coord (x,y,z)
6152 newPoints->GetPoint(v_id[1], v_line[1]);
6153
6154 const double *planeNormal = this->PlaneNormal[planes];
6155 const double *planePoint = this->PlanePoint[planes];
6156 values[0] = ( planeNormal[0]*(v_line[0][0] - planePoint[0])
6157 + planeNormal[1]*(v_line[0][1] - planePoint[1])
6158 + planeNormal[2]*(v_line[0][2] - planePoint[2]));
6159 values[1] = ( planeNormal[0]*(v_line[1][0] - planePoint[0])
6160 + planeNormal[1]*(v_line[1][1] - planePoint[1])
6161 + planeNormal[2]*(v_line[1][2] - planePoint[2]));
6162
6163 // Check to see if line is inside plane.
6164 if ((values[0] <= 0) && (values[1] <= 0))
6165 {
6166 newCellId = newcellArray->InsertNextCell(2, v_id);
6167 continue;
6168 }
6169
6170 // Check to see if line is outside plane.
6171 if ((values[0] >= 0) && (values[1] >= 0))
6172 {
6173 newCellId = lines[1]->InsertNextCell(2, v_id);
6174 outCD[1]->CopyData(inCD, cellId, newCellId);
6175 continue;
6176 }
6177
6178 // If we are here, the plane intersects the line segment.
6179 t = values[0]/(values[0] - values[1]);
6180 for (j = 0; j < 3; j++)
6181 {
6182 x[j] = (v_line[1][j] - v_line[0][j])*t + v_line[0][j];
6183 }
6184
6185 // Incorporate point into output and interpolate edge data as
6186 // necessary.
6187 if (locator->InsertUniquePoint(x, p_id))
6188 {
6189 this->InterpolateEdge(outPD, p_id, v_id[0], v_id[0], t);
6190 }
6191
6192 // Add the clipped line to the output.
6193 if (values[0] <= 0)
6194 {
6195 // First point of line is inside.
6196 tab_id[0] = v_id[0];
6197 tab_id[1] = p_id;
6198 newcellArray->InsertNextCell(2, tab_id);
6199
6200 // Second point of line is outside.
6201 tab_id[0] = p_id;
6202 tab_id[1] = v_id[1];
6203 cellarrayout->InsertNextCell(2, tab_id);
6204 }
6205 else
6206 {
6207 // Second point of line is inside.
6208 tab_id[0] = p_id;
6209 tab_id[1] = v_id[1];
6210 newcellArray->InsertNextCell(2, tab_id);
6211
6212 // First point of line is outside.
6213 tab_id[0] = v_id[0];
6214 tab_id[1] = p_id;
6215 cellarrayout->InsertNextCell(2, tab_id);
6216 }
6217 } // for all new cells
6218 cellarray->Delete();
6219 cellarray = newcellArray;
6220 } // for all planes
6221
6222 unsigned int totalnewcells = cellarray->GetNumberOfCells(); // Inside
6223 for (idcellnew = 0 ; idcellnew < totalnewcells; idcellnew++)
6224 {
6225 cellarray->GetNextCell(npts, v_id);
6226 newCellId = lines[0]->InsertNextCell(npts, v_id);
6227 outCD[0]->CopyData(inCD, cellId, newCellId);
6228 }
6229 cellarray->Delete();
6230
6231 totalnewcells = cellarrayout->GetNumberOfCells(); // Outside
6232 for (idcellnew = 0; idcellnew < totalnewcells; idcellnew++)
6233 {
6234 cellarrayout->GetNextCell(npts, v_id);
6235 newCellId = lines[1]->InsertNextCell(npts, v_id);
6236 outCD[1]->CopyData(inCD, cellId, newCellId);
6237 }
6238 cellarrayout->Delete();
6239 }
6240 arrayline->Delete();
6241 }
6242
6243 //-----------------------------------------------------------------------------
6244
ClipBox0D(vtkGenericCell * cell,vtkIncrementalPointLocator * locator,vtkCellArray * verts,vtkPointData * inPD,vtkPointData * outPD,vtkCellData * inCD,vtkIdType cellId,vtkCellData * outCD)6245 void vtkBoxClipDataSet::ClipBox0D(vtkGenericCell *cell,
6246 vtkIncrementalPointLocator *locator,
6247 vtkCellArray *verts,
6248 vtkPointData *inPD,
6249 vtkPointData *outPD,
6250 vtkCellData *inCD,
6251 vtkIdType cellId,
6252 vtkCellData *outCD)
6253 {
6254 vtkIdType cellType = cell->GetCellType();
6255 vtkIdList *cellIds = cell->GetPointIds();
6256 vtkCellArray *arrayvert = vtkCellArray::New();
6257 vtkPoints *cellPts = cell->GetPoints();
6258 vtkIdType npts = cellPts->GetNumberOfPoints();
6259 std::vector<vtkIdType> cellptId(npts);
6260 vtkIdType iid;
6261 vtkIdType *v_id = NULL;
6262 vtkIdType ptId;
6263 vtkIdType ptsvert = 1;
6264
6265 int i;
6266
6267 double v[3];
6268
6269 for (i = 0; i < npts; i++)
6270 {
6271 cellptId[i] = cellIds->GetId(i);
6272 }
6273
6274 // Convert all 0d cells to single vert.
6275 this->CellGrid(cellType, npts, &cellptId[0], arrayvert);
6276
6277 unsigned int totalnewvert = arrayvert->GetNumberOfCells();
6278 for (unsigned int idlinenew = 0; idlinenew < totalnewvert; idlinenew++)
6279 {
6280 arrayvert->GetNextCell(ptsvert, v_id);
6281
6282 // Clipping verts is easy. Either it is inside the box or it isn't.
6283 cellPts->GetPoint(v_id[0], v);
6284 if ( (v[0] >= this->BoundBoxClip[0][0])
6285 && (v[0] <= this->BoundBoxClip[0][1])
6286 && (v[1] >= this->BoundBoxClip[1][0])
6287 && (v[1] <= this->BoundBoxClip[1][1])
6288 && (v[2] >= this->BoundBoxClip[2][0])
6289 && (v[2] <= this->BoundBoxClip[2][1]) )
6290 {
6291 // Vert is inside.
6292 ptId = cellIds->GetId(v_id[0]);
6293 if (locator->InsertUniquePoint(v, iid))
6294 {
6295 outPD->CopyData(inPD, ptId, iid);
6296 }
6297
6298 int newCellId = verts->InsertNextCell(1, &iid);
6299 outCD->CopyData(inCD, cellId, newCellId);
6300 }
6301 }
6302 arrayvert->Delete();
6303 }
6304
6305 //-----------------------------------------------------------------------------
6306
ClipBoxInOut0D(vtkGenericCell * cell,vtkIncrementalPointLocator * locator,vtkCellArray ** verts,vtkPointData * inPD,vtkPointData * outPD,vtkCellData * inCD,vtkIdType cellId,vtkCellData ** outCD)6307 void vtkBoxClipDataSet::ClipBoxInOut0D(vtkGenericCell *cell,
6308 vtkIncrementalPointLocator *locator,
6309 vtkCellArray **verts,
6310 vtkPointData *inPD,
6311 vtkPointData *outPD,
6312 vtkCellData *inCD,
6313 vtkIdType cellId,
6314 vtkCellData **outCD)
6315 {
6316 vtkIdType cellType = cell->GetCellType();
6317 vtkIdList *cellIds = cell->GetPointIds();
6318 vtkCellArray *arrayvert = vtkCellArray::New();
6319 vtkPoints *cellPts = cell->GetPoints();
6320 vtkIdType npts = cellPts->GetNumberOfPoints();
6321 std::vector<vtkIdType> cellptId(npts);
6322 vtkIdType iid;
6323 vtkIdType *v_id = NULL;
6324 vtkIdType ptId;
6325 vtkIdType ptsvert = 1;
6326
6327 int i;
6328
6329 double v[3];
6330
6331 for (i = 0; i < npts; i++)
6332 {
6333 cellptId[i] = cellIds->GetId(i);
6334 }
6335
6336 // Convert all 0d cells to single vert.
6337 this->CellGrid(cellType, npts, &cellptId[0], arrayvert);
6338
6339 unsigned int totalnewvert = arrayvert->GetNumberOfCells();
6340 for (unsigned int idlinenew = 0; idlinenew < totalnewvert; idlinenew++)
6341 {
6342 arrayvert->GetNextCell(ptsvert, v_id);
6343
6344 // One way or another, we are adding the point.
6345 ptId = cellIds->GetId(v_id[0]);
6346 cellPts->GetPoint(v_id[0], v);
6347
6348 if (locator->InsertUniquePoint(v, iid))
6349 {
6350 outPD->CopyData(inPD, ptId, iid);
6351 }
6352
6353 // Clipping verts is easy. Either it is inside the box or it isn't.
6354 if ( (v[0] >= this->BoundBoxClip[0][0])
6355 && (v[0] <= this->BoundBoxClip[0][1])
6356 && (v[1] >= this->BoundBoxClip[1][0])
6357 && (v[1] <= this->BoundBoxClip[1][1])
6358 && (v[2] >= this->BoundBoxClip[2][0])
6359 && (v[2] <= this->BoundBoxClip[2][1]) )
6360 {
6361 // Vert is inside.
6362 int newCellId = verts[0]->InsertNextCell(1, &iid);
6363 outCD[0]->CopyData(inCD, cellId, newCellId);
6364 }
6365 else
6366 {
6367 // Vert is outside.
6368 int newCellId = verts[1]->InsertNextCell(1, &iid);
6369 outCD[1]->CopyData(inCD, cellId, newCellId);
6370 }
6371 }
6372 arrayvert->Delete();
6373 }
6374
6375 //-----------------------------------------------------------------------------
6376
ClipHexahedron0D(vtkGenericCell * cell,vtkIncrementalPointLocator * locator,vtkCellArray * verts,vtkPointData * inPD,vtkPointData * outPD,vtkCellData * inCD,vtkIdType cellId,vtkCellData * outCD)6377 void vtkBoxClipDataSet::ClipHexahedron0D(vtkGenericCell *cell,
6378 vtkIncrementalPointLocator *locator,
6379 vtkCellArray *verts,
6380 vtkPointData *inPD,
6381 vtkPointData *outPD,
6382 vtkCellData *inCD,
6383 vtkIdType cellId,
6384 vtkCellData *outCD)
6385 {
6386 vtkIdType cellType = cell->GetCellType();
6387 vtkIdList *cellIds = cell->GetPointIds();
6388 vtkCellArray *arrayvert = vtkCellArray::New();
6389 vtkPoints *cellPts = cell->GetPoints();
6390 vtkIdType npts = cellPts->GetNumberOfPoints();
6391 std::vector<vtkIdType> cellptId(npts);
6392 vtkIdType iid;
6393 vtkIdType *v_id = NULL;
6394 vtkIdType ptId;
6395 vtkIdType ptsvert = 1;
6396
6397 int i;
6398
6399 double v[3];
6400
6401 for (i = 0; i < npts; i++)
6402 {
6403 cellptId[i] = cellIds->GetId(i);
6404 }
6405
6406 // Convert all 0d cells to single vert.
6407 this->CellGrid(cellType, npts, &cellptId[0], arrayvert);
6408
6409 unsigned int totalnewvert = arrayvert->GetNumberOfCells();
6410 for (unsigned int idlinenew = 0; idlinenew < totalnewvert; idlinenew++)
6411 {
6412 arrayvert->GetNextCell(ptsvert, v_id);
6413
6414 // Clipping verts is easy. Either it is inside the hexahedron or it isn't.
6415 cellPts->GetPoint(v_id[0], v);
6416 int inside = 1;
6417
6418 for (int k = 0; k < 6; k++)
6419 {
6420 double value
6421 = ( this->PlaneNormal[k][0]*(v[0] - this->PlanePoint[k][0])
6422 + this->PlaneNormal[k][1]*(v[1] - this->PlanePoint[k][1])
6423 + this->PlaneNormal[k][2]*(v[2] - this->PlanePoint[k][2]) );
6424 if (value > 0)
6425 {
6426 inside = 0;
6427 }
6428 }
6429
6430 if (inside)
6431 {
6432 // Vert is inside.
6433 ptId = cellIds->GetId(v_id[0]);
6434 if (locator->InsertUniquePoint(v, iid))
6435 {
6436 outPD->CopyData(inPD, ptId, iid);
6437 }
6438
6439 int newCellId = verts->InsertNextCell(1, &iid);
6440 outCD->CopyData(inCD, cellId, newCellId);
6441 }
6442 }
6443 arrayvert->Delete();
6444 }
6445
6446 //-----------------------------------------------------------------------------
6447
ClipHexahedronInOut0D(vtkGenericCell * cell,vtkIncrementalPointLocator * locator,vtkCellArray ** verts,vtkPointData * inPD,vtkPointData * outPD,vtkCellData * inCD,vtkIdType cellId,vtkCellData ** outCD)6448 void vtkBoxClipDataSet::ClipHexahedronInOut0D(vtkGenericCell *cell,
6449 vtkIncrementalPointLocator *locator,
6450 vtkCellArray **verts,
6451 vtkPointData *inPD,
6452 vtkPointData *outPD,
6453 vtkCellData *inCD,
6454 vtkIdType cellId,
6455 vtkCellData **outCD)
6456 {
6457 vtkIdType cellType = cell->GetCellType();
6458 vtkIdList *cellIds = cell->GetPointIds();
6459 vtkCellArray *arrayvert = vtkCellArray::New();
6460 vtkPoints *cellPts = cell->GetPoints();
6461 vtkIdType npts = cellPts->GetNumberOfPoints();
6462 std::vector<vtkIdType> cellptId(npts);
6463 vtkIdType iid;
6464 vtkIdType *v_id = NULL;
6465 vtkIdType ptId;
6466 vtkIdType ptsvert = 1;
6467
6468 int i;
6469
6470 double v[3];
6471
6472 for (i = 0; i < npts; i++)
6473 {
6474 cellptId[i] = cellIds->GetId(i);
6475 }
6476
6477 // Convert all 0d cells to single vert.
6478 this->CellGrid(cellType, npts, &cellptId[0], arrayvert);
6479
6480 unsigned int totalnewvert = arrayvert->GetNumberOfCells();
6481 for (unsigned int idlinenew = 0; idlinenew < totalnewvert; idlinenew++)
6482 {
6483 arrayvert->GetNextCell(ptsvert, v_id);
6484
6485 // One way or another, we are adding the point.
6486 ptId = cellIds->GetId(v_id[0]);
6487 cellPts->GetPoint(v_id[0], v);
6488
6489 if (locator->InsertUniquePoint(v, iid))
6490 {
6491 outPD->CopyData(inPD, ptId, iid);
6492 }
6493
6494 int inside = 1;
6495 for (int k = 0; k < 6; k++)
6496 {
6497 double value
6498 = ( this->PlaneNormal[k][0]*(v[0] - this->PlanePoint[k][0])
6499 + this->PlaneNormal[k][1]*(v[1] - this->PlanePoint[k][1])
6500 + this->PlaneNormal[k][2]*(v[2] - this->PlanePoint[k][2]) );
6501 if (value > 0)
6502 {
6503 inside = 0;
6504 }
6505 }
6506
6507 // Clipping verts is easy. Either it is inside the box or it isn't.
6508 if (inside)
6509 {
6510 // Vert is inside.
6511 int newCellId = verts[0]->InsertNextCell(1, &iid);
6512 outCD[0]->CopyData(inCD, cellId, newCellId);
6513 }
6514 else
6515 {
6516 // Vert is outside.
6517 int newCellId = verts[1]->InsertNextCell(1, &iid);
6518 outCD[1]->CopyData(inCD, cellId, newCellId);
6519 }
6520 }
6521 arrayvert->Delete();
6522 }
6523