vtkClipHyperOctree.cxx 45.6 KB
Newer Older
Francois Bertel's avatar
Francois Bertel committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkClipHyperOctree.cxx

  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notice for more information.

=========================================================================*/
#include "vtkClipHyperOctree.h"

#include "vtkHyperOctree.h"
#include "vtkCellArray.h"
#include "vtkCellData.h"
#include "vtkClipVolume.h"
#include "vtkExecutive.h"
#include "vtkDoubleArray.h"
#include "vtkGenericCell.h"
#include "vtkImageData.h"
#include "vtkImplicitFunction.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkIntArray.h"
#include "vtkMergePoints.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkUnsignedCharArray.h"
#include "vtkUnstructuredGrid.h"
#include "vtkHyperOctreeCursor.h"
#include "vtkVoxel.h"
#include "vtkPixel.h"
#include "vtkLine.h"
#include "vtkTetra.h"
#include "vtkPolygon.h"
#include "vtkOrderedTriangulator.h"
#include "vtkHyperOctreeClipCutPointsGrabber.h"
42
#include "vtkIncrementalPointLocator.h"
Francois Bertel's avatar
Francois Bertel committed
43

Sean McBride's avatar
Sean McBride committed
44
#include <cmath>
45
#include <cassert>
Francois Bertel's avatar
Francois Bertel committed
46
47
48
49
50
51
52
53
54

vtkStandardNewMacro(vtkClipHyperOctree);
vtkCxxSetObjectMacro(vtkClipHyperOctree,ClipFunction,vtkImplicitFunction);

//----------------------------------------------------------------------------
// Construct with user-specified implicit function; InsideOut turned off; value
// set to 0.0; and generate clip scalars turned off.
vtkClipHyperOctree::vtkClipHyperOctree(vtkImplicitFunction *cf)
{
David E. DeMarle's avatar
David E. DeMarle committed
55
  VTK_LEGACY_BODY(vtkClipHyperOctree, "VTK 8.1");
56

Francois Bertel's avatar
Francois Bertel committed
57
58
  this->ClipFunction = cf;
  this->InsideOut = 0;
59
60
  this->Locator = nullptr;
  this->Locator2=nullptr;
Francois Bertel's avatar
Francois Bertel committed
61
62
63
64
  this->Value = 0.0;
  this->GenerateClipScalars = 0;

  this->GenerateClippedOutput = 0;
65

Francois Bertel's avatar
Francois Bertel committed
66
67
68
69
70
71
72
73
  this->SetNumberOfOutputPorts(2);
  vtkUnstructuredGrid *output2 = vtkUnstructuredGrid::New();
  this->GetExecutive()->SetOutputData(1, output2);
  output2->Delete();

  // by default process active point scalars
  this->SetInputArrayToProcess(0,0,0,vtkDataObject::FIELD_ASSOCIATION_POINTS,
                               vtkDataSetAttributes::SCALARS);
74

75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
  this->Input=nullptr;
  this->Output=nullptr;
  this->ClippedOutput=nullptr;
  this->Conn[0]=nullptr;
  this->Conn[1]=nullptr;
  this->Types[0]=nullptr;
  this->Types[1]=nullptr;
  this->Locs[0]=nullptr;
  this->Locs[1]=nullptr;
  this->InCD=nullptr;
  this->OutCD[0]=nullptr;
  this->OutCD[1]=nullptr;
  this->OutPD[0]=nullptr;
  this->OutPD[1]=nullptr;
  this->Triangulator=nullptr;
  this->Sibling=nullptr;

  this->Tetra=nullptr;
  this->Polygon=nullptr;
  this->TetScalars=nullptr;
  this->CellScalars=nullptr;
  this->Pts=nullptr;
  this->Grabber=nullptr;
Francois Bertel's avatar
Francois Bertel committed
98
99
100
101
102
103
}

//----------------------------------------------------------------------------
vtkClipHyperOctree::~vtkClipHyperOctree()
{
  if ( this->Locator )
104
  {
Francois Bertel's avatar
Francois Bertel committed
105
    this->Locator->UnRegister(this);
106
    this->Locator = nullptr;
107
  }
108
  this->SetClipFunction(nullptr);
Francois Bertel's avatar
Francois Bertel committed
109
110
111
112
113
}

//----------------------------------------------------------------------------
// Overload standard modified time function. If Clip functions is modified,
// then this object is modified as well.
Bill Lorensen's avatar
Bill Lorensen committed
114
vtkMTimeType vtkClipHyperOctree::GetMTime()
Francois Bertel's avatar
Francois Bertel committed
115
{
Bill Lorensen's avatar
Bill Lorensen committed
116
117
  vtkMTimeType mTime=this->Superclass::GetMTime();
  vtkMTimeType time;
Francois Bertel's avatar
Francois Bertel committed
118

119
  if ( this->ClipFunction != nullptr )
120
  {
Francois Bertel's avatar
Francois Bertel committed
121
122
    time = this->ClipFunction->GetMTime();
    mTime = ( time > mTime ? time : mTime );
123
  }
124
  if ( this->Locator != nullptr )
125
  {
Francois Bertel's avatar
Francois Bertel committed
126
127
    time = this->Locator->GetMTime();
    mTime = ( time > mTime ? time : mTime );
128
  }
Francois Bertel's avatar
Francois Bertel committed
129
130
131
132
133
134
135

  return mTime;
}

vtkUnstructuredGrid *vtkClipHyperOctree::GetClippedOutput()
{
  if (!this->GenerateClippedOutput)
136
  {
137
    return nullptr;
138
  }
Francois Bertel's avatar
Francois Bertel committed
139
140
141
142
143
144
145
146
147
148
149
150
151
  return vtkUnstructuredGrid::SafeDownCast(
    this->GetExecutive()->GetOutputData(1));
}

//----------------------------------------------------------------------------
//
// Clip through data generating surface.
//
int vtkClipHyperOctree::RequestData(vtkInformation *vtkNotUsed(request),
                                    vtkInformationVector **inputVector,
                                    vtkInformationVector *outputVector)
{
  if ( !this->ClipFunction)
152
  {
Francois Bertel's avatar
Francois Bertel committed
153
154
    vtkErrorMacro(<<"As HyperOctree does not support point data yet, a clip function has to be provided.");
    return 1;
155
  }
156

Francois Bertel's avatar
Francois Bertel committed
157
  if ( !this->ClipFunction && this->GenerateClipScalars )
158
  {
Francois Bertel's avatar
Francois Bertel committed
159
160
    vtkErrorMacro(<<"Cannot generate clip scalars if no clip function defined");
    return 1;
161
  }
162

Francois Bertel's avatar
Francois Bertel committed
163
164
165
166
  // get the info objects
  vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation *outInfo = outputVector->GetInformationObject(0);

167
  // get the input and output
Francois Bertel's avatar
Francois Bertel committed
168
169
170
171
172
173
  this->Input = vtkHyperOctree::SafeDownCast(
    inInfo->Get(vtkDataObject::DATA_OBJECT()));
  this->Output=vtkUnstructuredGrid::SafeDownCast(
    outInfo->Get(vtkDataObject::DATA_OBJECT()));

  this->ClippedOutput = this->GetClippedOutput();
174

Francois Bertel's avatar
Francois Bertel committed
175
  vtkIdType numPts=this->Input->GetMaxNumberOfPoints(0);
176
  vtkIdType numCells = this->Input->GetNumberOfLeaves();
177

Francois Bertel's avatar
Francois Bertel committed
178
179
  vtkPoints *newPoints = vtkPoints::New();
  newPoints->Allocate(numPts,numPts/2);
180

Francois Bertel's avatar
Francois Bertel committed
181
182
183
184
  // allocate the output and associated helper classes
  vtkIdType estimatedSize = numCells;
  estimatedSize = estimatedSize / 1024 * 1024; //multiple of 1024
  if (estimatedSize < 1024)
185
  {
Francois Bertel's avatar
Francois Bertel committed
186
    estimatedSize = 1024;
187
  }
Francois Bertel's avatar
Francois Bertel committed
188
189
190
191
192
193
194
195
  this->Conn[0] = vtkCellArray::New();
  this->Conn[0]->Allocate(estimatedSize,estimatedSize/2);
  this->Conn[0]->InitTraversal();
  this->Types[0] = vtkUnsignedCharArray::New();
  this->Types[0]->Allocate(estimatedSize,estimatedSize/2);
  this->Locs[0] = vtkIdTypeArray::New();
  this->Locs[0]->Allocate(estimatedSize,estimatedSize/2);
  if ( this->GenerateClippedOutput )
196
  {
Francois Bertel's avatar
Francois Bertel committed
197
198
199
200
201
202
203
204
    // numOutputs = 2;
    this->Conn[1] = vtkCellArray::New();
    this->Conn[1]->Allocate(estimatedSize,estimatedSize/2);
    this->Conn[1]->InitTraversal();
    this->Types[1] = vtkUnsignedCharArray::New();
    this->Types[1]->Allocate(estimatedSize,estimatedSize/2);
    this->Locs[1] = vtkIdTypeArray::New();
    this->Locs[1]->Allocate(estimatedSize,estimatedSize/2);
205
  }
206

207
  vtkPoints *newPoints2=nullptr;
208

Francois Bertel's avatar
Francois Bertel committed
209
  // locator used to merge potentially duplicate points
210
  if ( this->Locator == nullptr )
211
  {
Francois Bertel's avatar
Francois Bertel committed
212
    this->CreateDefaultLocator();
213
  }
214

Francois Bertel's avatar
Francois Bertel committed
215
  if(this->GenerateClippedOutput)
216
  {
Francois Bertel's avatar
Francois Bertel committed
217
218
219
220
    this->Locator2=this->Locator->NewInstance();
    newPoints2 = vtkPoints::New();
    newPoints2->Allocate(numPts,numPts/2);
    this->Locator2->InitPointInsertion (newPoints2, this->Input->GetBounds());
221
  }
222

Francois Bertel's avatar
Francois Bertel committed
223
  this->Locator->InitPointInsertion (newPoints, this->Input->GetBounds());
224

Francois Bertel's avatar
Francois Bertel committed
225
  this->InCD=static_cast<vtkCellData *>(this->Input->GetLeafData());
Francois Bertel's avatar
Francois Bertel committed
226
227
228
  this->OutCD[0] = this->Output->GetCellData();
  this->OutCD[0]->CopyAllocate(this->InCD,estimatedSize,estimatedSize/2);
  if ( this->GenerateClippedOutput )
229
  {
Francois Bertel's avatar
Francois Bertel committed
230
231
    this->OutCD[1] = this->ClippedOutput->GetCellData();
    this->OutCD[1]->CopyAllocate(this->InCD,estimatedSize,estimatedSize/2);
232
  }
233

Francois Bertel's avatar
Francois Bertel committed
234
235
236
  this->OutPD[0]=this->Output->GetPointData();
  if ( !this->GenerateClipScalars &&
       !this->GetInputArrayToProcess(0,inputVector))
237
  {
Francois Bertel's avatar
Francois Bertel committed
238
    this->OutPD[0]->CopyScalarsOff();
239
  }
Francois Bertel's avatar
Francois Bertel committed
240
  else
241
  {
Francois Bertel's avatar
Francois Bertel committed
242
    this->OutPD[0]->CopyScalarsOn();
243
  }
244

Francois Bertel's avatar
Francois Bertel committed
245
  if(this->GenerateClippedOutput)
246
  {
Francois Bertel's avatar
Francois Bertel committed
247
248
249
    this->OutPD[1]=this->ClippedOutput->GetPointData();
    if ( !this->GenerateClipScalars &&
         !this->GetInputArrayToProcess(0,inputVector))
250
    {
Francois Bertel's avatar
Francois Bertel committed
251
      this->OutPD[1]->CopyScalarsOff();
252
    }
Francois Bertel's avatar
Francois Bertel committed
253
    else
254
    {
Francois Bertel's avatar
Francois Bertel committed
255
256
      this->OutPD[1]->CopyScalarsOn();
    }
257
  }
258

Francois Bertel's avatar
Francois Bertel committed
259
260
  vtkHyperOctreeCursor *cursor=this->Input->NewCellCursor();
  this->Sibling=cursor->Clone();
261

Francois Bertel's avatar
Francois Bertel committed
262
263
264
265
266
  cursor->ToRoot();
  double bounds[6];
  this->Input->GetBounds(bounds);

  switch(this->Input->GetDimension())
267
  {
Francois Bertel's avatar
Francois Bertel committed
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
    case 3:
      this->Tetra=vtkTetra::New();
      this->TetScalars=vtkDoubleArray::New();
      this->TetScalars->SetNumberOfComponents(1);
      this->TetScalars->SetNumberOfTuples(4);
      this->Grabber=vtkHyperOctreeClipCutPointsGrabber::New();
      this->Grabber->SetDimension(3);
      this->Triangulator=this->Grabber->GetTriangulator();
      break;
    case 2:
      this->Grabber=vtkHyperOctreeClipCutPointsGrabber::New();
      this->Grabber->SetDimension(2);
      this->Polygon=this->Grabber->GetPolygon();
      break;
    default:
      // do nothing
      break;
285
  }
Francois Bertel's avatar
Francois Bertel committed
286
287
  this->CellScalars=vtkDoubleArray::New();
  this->Pts=vtkPoints::New();
288

Francois Bertel's avatar
Francois Bertel committed
289
290
  this->TotalCounter=0;
  this->TemplateCounter=0;
291

Francois Bertel's avatar
Francois Bertel committed
292
293
  int j=0;
  while(j<65536)
294
  {
Francois Bertel's avatar
Francois Bertel committed
295
296
    this->CellTypeCounter[j]=0; // up-to-65536 points per octant
    ++j;
297
  }
Francois Bertel's avatar
Francois Bertel committed
298
  this->ClipNode(cursor,0,bounds);
299

Francois Bertel's avatar
Francois Bertel committed
300
//  cout<<"ClipHyperOctree: "<<this->TemplateCounter<<" templates over "<<this->TotalCounter<<" octants, ratio="<<(this->TemplateCounter/static_cast<double>(this->TotalCounter))<<endl;
301

Francois Bertel's avatar
Francois Bertel committed
302
303
  j=0;
  while(j<65536)
304
  {
Francois Bertel's avatar
Francois Bertel committed
305
    if(this->CellTypeCounter[j]>0)
306
    {
Francois Bertel's avatar
Francois Bertel committed
307
308
//      cout<<this->CellTypeCounter[j]<<" with "<<j+1<<"points"<<endl;
    }
309
310
    ++j;
  }
311

Francois Bertel's avatar
Francois Bertel committed
312
  switch(this->Input->GetDimension())
313
  {
Francois Bertel's avatar
Francois Bertel committed
314
315
    case 3:
      this->Tetra->UnRegister(this);
316
      this->Tetra=nullptr;
Francois Bertel's avatar
Francois Bertel committed
317
      this->TetScalars->UnRegister(this);
318
319
      this->TetScalars=nullptr;
      this->Triangulator=nullptr;
Francois Bertel's avatar
Francois Bertel committed
320
      this->Grabber->UnRegister(this);
321
      this->Grabber=nullptr;
Francois Bertel's avatar
Francois Bertel committed
322
323
      break;
    case 2:
324
      this->Polygon=nullptr;
Francois Bertel's avatar
Francois Bertel committed
325
      this->Grabber->UnRegister(this);
326
      this->Grabber=nullptr;
Francois Bertel's avatar
Francois Bertel committed
327
328
329
      break;
    default:
      break;
330
  }
331

Francois Bertel's avatar
Francois Bertel committed
332
  this->CellScalars->UnRegister(this);
333
  this->CellScalars=nullptr;
Francois Bertel's avatar
Francois Bertel committed
334
  this->Pts->UnRegister(this);
335
  this->Pts=nullptr;
336

Francois Bertel's avatar
Francois Bertel committed
337
338
  cursor->UnRegister(this);
  this->Sibling->UnRegister(this);
339
  this->Sibling=nullptr;
340

341
342
343
  this->OutPD[0]=nullptr;
  this->Input=nullptr;
  this->InCD=nullptr;
Francois Bertel's avatar
Francois Bertel committed
344
345
346
347
  this->Output->SetPoints(newPoints);
  newPoints->Delete();
  this->Output->SetCells(this->Types[0], this->Locs[0], this->Conn[0]);
  this->Conn[0]->Delete();
348
  this->Conn[0]=nullptr;
Francois Bertel's avatar
Francois Bertel committed
349
  this->Types[0]->Delete();
350
  this->Types[0]=nullptr;
Francois Bertel's avatar
Francois Bertel committed
351
  this->Locs[0]->Delete();
352
353
  this->Locs[0]=nullptr;
  this->OutCD[0]=nullptr;
354

Francois Bertel's avatar
Francois Bertel committed
355
  if(this->GenerateClippedOutput)
356
  {
Francois Bertel's avatar
Francois Bertel committed
357
358
359
    this->ClippedOutput->SetPoints(newPoints2);
    this->ClippedOutput->SetCells(this->Types[1], this->Locs[1],this->Conn[1]);
    this->Conn[1]->Delete();
360
    this->Conn[1]=nullptr;
Francois Bertel's avatar
Francois Bertel committed
361
    this->Types[1]->Delete();
362
    this->Types[1]=nullptr;
Francois Bertel's avatar
Francois Bertel committed
363
    this->Locs[1]->Delete();
364
    this->Locs[1]=nullptr;
Francois Bertel's avatar
Francois Bertel committed
365
366
    newPoints2->Delete();
    this->Locator2->Delete();
367
368
369
    this->Locator2=nullptr;
    this->OutCD[1]=nullptr;
    this->OutPD[1]=nullptr;
370
  }
371

Francois Bertel's avatar
Francois Bertel committed
372
373
  this->Locator->Initialize();//release any extra memory
  this->Output->Squeeze();
374
  this->Output=nullptr;
Francois Bertel's avatar
Francois Bertel committed
375
  if(this->GenerateClippedOutput)
376
  {
Francois Bertel's avatar
Francois Bertel committed
377
    this->ClippedOutput->Squeeze();
378
    this->ClippedOutput=nullptr;
379
  }
380

381
382
383
384
385
386
387
388
389
390
  assert("post: input_is_null" && this->Input==nullptr);
  assert("post: output_is_null" && this->Output==nullptr);
  assert("post: clipped_output_is_null" && this->ClippedOutput==nullptr);
  assert("post: locator2_is_null" && this->Locator2==nullptr);
  assert("post: types_are_null" && this->Types[0]==nullptr && this->Types[1]==nullptr);
  assert("post: conn_are_null" && this->Conn[0]==nullptr && this->Conn[1]==nullptr);
  assert("post: locs_are_null" && this->Locs[0]==nullptr && this->Locs[1]==nullptr);
  assert("post: incd_is_null" && this->InCD==nullptr);
  assert("post: outpd_are_null" && this->OutPD[0]==nullptr && this->OutPD[1]==nullptr);
  assert("post: outcd_are_null" && this->OutCD[0]==nullptr && this->OutCD[1]==nullptr);
391

Francois Bertel's avatar
Francois Bertel committed
392
393
394
395
396
397
398
  return 1;
}

//----------------------------------------------------------------------------
void vtkClipHyperOctree::ClipNode(vtkHyperOctreeCursor *cursor,
                                  int level,
                                  double bounds[6])
399
{
400
  assert("pre: cursor_exists" && cursor!=nullptr);
Francois Bertel's avatar
Francois Bertel committed
401
  assert("pre: positive_level" && level>=0);
402

Francois Bertel's avatar
Francois Bertel committed
403
  if(cursor->CurrentIsLeaf())
404
  {
Francois Bertel's avatar
Francois Bertel committed
405
    if(cursor->CurrentIsRoot() || (this->Input->GetDimension()==1))
406
    {
Francois Bertel's avatar
Francois Bertel committed
407
408
409
      // no parent=>no sibling=>no sibling which are not leaves=>easy
      // just create a voxel/pixel/line and clip it.

410
      vtkCell *cell=nullptr;
Francois Bertel's avatar
Francois Bertel committed
411
      vtkIdType cellId=cursor->GetLeafId(); // only one cell.
412

Francois Bertel's avatar
Francois Bertel committed
413
414
415
      vtkDoubleArray *cellScalars;
      cellScalars=vtkDoubleArray::New();// scalar at each corner point.
      cellScalars->Allocate(VTK_CELL_SIZE);
416

Francois Bertel's avatar
Francois Bertel committed
417
418
419
420
421
422
423
      vtkIdType numPts=0;
      double pt[3];
      vtkVoxel *v;
      vtkPixel *p;
      vtkLine *l;
      int coord;
      switch(this->Input->GetDimension())
424
      {
Francois Bertel's avatar
Francois Bertel committed
425
426
427
428
429
430
        case 3:
          v=vtkVoxel::New();
          cell=v;
          numPts=8;
          coord=0;
          while(coord<3)
431
          {
Francois Bertel's avatar
Francois Bertel committed
432
433
            pt[0]=bounds[coord*2];
            ++coord;
434
          }
Francois Bertel's avatar
Francois Bertel committed
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
          v->GetPoints()->SetPoint(0,pt);
          pt[0]=bounds[1];
          v->GetPoints()->SetPoint(1,pt);
          pt[0]=bounds[0];
          pt[1]=bounds[3];
          v->GetPoints()->SetPoint(2,pt);
          pt[0]=bounds[1];
          v->GetPoints()->SetPoint(3,pt);
          pt[0]=bounds[0];
          pt[1]=bounds[2];
          pt[2]=bounds[5];
          v->GetPoints()->SetPoint(4,pt);
          pt[0]=bounds[1];
          v->GetPoints()->SetPoint(5,pt);
          pt[0]=bounds[0];
          pt[1]=bounds[3];
          v->GetPoints()->SetPoint(6,pt);
          pt[0]=bounds[1];
          v->GetPoints()->SetPoint(7,pt);
          break;
        case 2:
          p=vtkPixel::New();
          cell=p;
          numPts=4;
          coord=0;
          while(coord<3)
461
          {
Francois Bertel's avatar
Francois Bertel committed
462
463
            pt[0]=bounds[coord*2];
            ++coord;
464
          }
Francois Bertel's avatar
Francois Bertel committed
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
          p->GetPoints()->SetPoint(0,pt);
          pt[0]=bounds[1];
          p->GetPoints()->SetPoint(1,pt);
          pt[0]=bounds[0];
          pt[1]=bounds[3];
          p->GetPoints()->SetPoint(2,pt);
          pt[0]=bounds[1];
          p->GetPoints()->SetPoint(3,pt);
          break;
        case 1:
          l=vtkLine::New();
          cell=l;
          numPts=2;
          pt[0]=bounds[0];
          pt[1]=bounds[2];
          pt[2]=bounds[4];
          l->GetPoints()->SetPoint(0,pt);
          pt[0]=bounds[1];
          l->GetPoints()->SetPoint(1,pt);
          break;
        default:
          assert("check: impossible" && 0);
          break;
488
      }
489

490
      vtkDataArray *clipScalars=nullptr;
491

Francois Bertel's avatar
Francois Bertel committed
492
      vtkPointData *inPD=this->Input->GetPointData();
493

494
      if(this->ClipFunction!=nullptr)
495
      {
Francois Bertel's avatar
Francois Bertel committed
496
497
498
499
500
501
        vtkDoubleArray *tmpScalars = vtkDoubleArray::New();
        tmpScalars->SetNumberOfTuples(numPts);
        tmpScalars->SetName("ClipDataSetScalars");
        inPD = vtkPointData::New();
        inPD->ShallowCopy(this->Input->GetPointData());//copies original
        if(this->GenerateClipScalars)
502
        {
Francois Bertel's avatar
Francois Bertel committed
503
          inPD->SetScalars(tmpScalars);
504
        }
Francois Bertel's avatar
Francois Bertel committed
505
        for ( int i=0; i < numPts; i++ )
506
        {
Francois Bertel's avatar
Francois Bertel committed
507
508
509
          double s = this->ClipFunction->FunctionValue(cell->GetPoints()->GetPoint(i));
          tmpScalars->SetTuple1(i,s);
        }
510
511
        clipScalars = tmpScalars;
      }
512

Francois Bertel's avatar
Francois Bertel committed
513
514
515
#if 0 // just to not break compilation
      outPD->InterpolateAllocate(inPD,estimatedSize,estimatedSize/2);
#endif
516
517
      int i;
      for (i=0; i < numPts; i++ )
518
      {
Francois Bertel's avatar
Francois Bertel committed
519
520
        double s = clipScalars->GetComponent(i, 0);
        cellScalars->InsertTuple(i, &s);
521
      }
522
523

      if ( this->ClipFunction )
524
      {
Francois Bertel's avatar
Francois Bertel committed
525
526
        clipScalars->UnRegister(this);
        inPD->UnRegister(this);
527
      }
528

Francois Bertel's avatar
Francois Bertel committed
529
530
531
      // perform clipping
      int num[2]; num[0]=num[1]=0;
      int numNew[2]; numNew[0]=numNew[1]=0;
532

Francois Bertel's avatar
Francois Bertel committed
533
534
535
536
537
      cell->Clip(this->Value,cellScalars,this->Locator,this->Conn[0],inPD,
                 this->OutPD[0],this->InCD,cellId,this->OutCD[0],
                 this->InsideOut);
      numNew[0] = this->Conn[0]->GetNumberOfCells() - num[0];
      num[0] = this->Conn[0]->GetNumberOfCells();
538

Francois Bertel's avatar
Francois Bertel committed
539
      if(this->GenerateClippedOutput)
540
      {
Francois Bertel's avatar
Francois Bertel committed
541
542
543
544
545
        cell->Clip(this->Value, cellScalars, this->Locator2, this->Conn[1],
                   inPD, this->OutPD[1], this->InCD, cellId, this->OutCD[1],
                   !this->InsideOut);
        numNew[1] = this->Conn[1]->GetNumberOfCells() - num[1];
        num[1] = this->Conn[1]->GetNumberOfCells();
546
      }
547

Francois Bertel's avatar
Francois Bertel committed
548
      cellScalars->UnRegister(this);
549

Francois Bertel's avatar
Francois Bertel committed
550
551
      int numOutputs;
      if(this->GenerateClippedOutput)
552
      {
Francois Bertel's avatar
Francois Bertel committed
553
        numOutputs=2;
554
      }
Francois Bertel's avatar
Francois Bertel committed
555
      else
556
      {
Francois Bertel's avatar
Francois Bertel committed
557
        numOutputs=1;
558
      }
559
560

      vtkIdType npts=0;
Francois Bertel's avatar
Francois Bertel committed
561
562
      vtkIdType *pts;
      int cellType;
563

564
      for (i=0; i<numOutputs; i++) //for both outputs
565
      {
566
        for (int j=0; j < numNew[i]; j++)
567
        {
Francois Bertel's avatar
Francois Bertel committed
568
569
570
          this->Locs[i]->InsertNextValue(
            this->Conn[i]->GetTraversalLocation());
          this->Conn[i]->GetNextCell(npts,pts);
571

Francois Bertel's avatar
Francois Bertel committed
572
573
          //For each new cell added, got to set the type of the cell
          switch ( cell->GetCellDimension() )
574
          {
Francois Bertel's avatar
Francois Bertel committed
575
576
577
            case 1: //lines are generated---------------------------------
              cellType = (npts > 2 ? VTK_POLY_LINE : VTK_LINE);
              break;
578

Francois Bertel's avatar
Francois Bertel committed
579
            case 2: //polygons are generated------------------------------
580
              cellType = (npts == 3 ? VTK_TRIANGLE :
Francois Bertel's avatar
Francois Bertel committed
581
582
                          (npts == 4 ? VTK_QUAD : VTK_POLYGON));
              break;
583

Francois Bertel's avatar
Francois Bertel committed
584
585
586
587
588
589
590
591
            case 3: //tetrahedra or wedges are generated------------------
              cellType = (npts == 4 ? VTK_TETRA : VTK_WEDGE);
              break;
            default:
              assert("check: impossible case" && 0);
              cellType=0; // useless, only for removing warning about
              // unitialized function.
              break;
592
          } //switch
593

Francois Bertel's avatar
Francois Bertel committed
594
          this->Types[i]->InsertNextValue(cellType);
595
596
        } //for each new cell
      } //for both outputs
Francois Bertel's avatar
Francois Bertel committed
597
598
      cell->UnRegister(this);

599
    }
Francois Bertel's avatar
Francois Bertel committed
600
    else
601
    {
Francois Bertel's avatar
Francois Bertel committed
602
603
604
605
606
607
608
609
610
611
612
      // some parent=>have sibling=>some sibling may have children
      // => those children may create points on some face of cursor
      // => difficult case
      //
      // Even worst, if the siblings don't have children, the
      // sibling of the parent may have children that create points
      // on some face.
      //
      // Even if there is no children, the neighbor cell tessellation
      // has to be compatible with the current cell tessellation.
      // In any case, we need the ordered triangulator.
613
614


Francois Bertel's avatar
Francois Bertel committed
615
616
617
618
619
620
621
      // Add the points of the current leaf
      // use the bounds

      vtkIdType ptId;
      int i;
      int j;
      int k;
622

Francois Bertel's avatar
Francois Bertel committed
623
624
625
      int pi;
      int pj;
      int pk;
626

Francois Bertel's avatar
Francois Bertel committed
627
628
629
      int x;
      int y;
      int z;
630

Francois Bertel's avatar
Francois Bertel committed
631
632
      // resolution in point along each axis.
      int resolution=(1<<(this->Input->GetNumberOfLevels()-1))+1;
633

Francois Bertel's avatar
Francois Bertel committed
634
635
      int deltaLevel=this->Input->GetNumberOfLevels()-1-level;
      assert("check: positive_deltaLevel" && deltaLevel>=0);
636

Francois Bertel's avatar
Francois Bertel committed
637
      double ratio=1.0/(resolution-1);
638

Francois Bertel's avatar
Francois Bertel committed
639
640
      double pt[3];
      double pcoords[3];
641

Francois Bertel's avatar
Francois Bertel committed
642
643
644
645
      int allOut=1; // bool
      int allIn=1; // bool
      int clipPoint; // bool
      double s;
646

Francois Bertel's avatar
Francois Bertel committed
647
648
      // index of the node
      if(this->Input->GetDimension()==3)
649
      {
Francois Bertel's avatar
Francois Bertel committed
650
651
        vtkIdType nbpts=this->Input->GetMaxNumberOfPointsOnBoundary(level);
        double pbounds[6]={0,1,0,1,0,1};
652

Francois Bertel's avatar
Francois Bertel committed
653
        this->Triangulator->InitTriangulation(pbounds,nbpts);
654

Francois Bertel's avatar
Francois Bertel committed
655
656
        this->Triangulator->PreSortedOff();
        this->Grabber->InitPointInsertion();
657

Francois Bertel's avatar
Francois Bertel committed
658
659
660
661
662
663
        i=cursor->GetIndex(0);
        j=cursor->GetIndex(1);
        k=cursor->GetIndex(2);
        z=0;
        pk=k;
        while(z<2)
664
        {
Francois Bertel's avatar
Francois Bertel committed
665
666
667
          pj=j;
          y=0;
          while(y<2)
668
          {
Francois Bertel's avatar
Francois Bertel committed
669
670
671
            pi=i;
            x=0;
            while(x<2)
672
            {
Francois Bertel's avatar
Francois Bertel committed
673
674
675
              pt[0]=bounds[x];
              pt[1]=bounds[2+y];
              pt[2]=bounds[4+z];
676

Francois Bertel's avatar
Francois Bertel committed
677
678
679
680
681
682
              assert("check: in_bounds" && pt[0]>=this->Input->GetBounds()[0] && pt[0]<=this->Input->GetBounds()[1] && pt[1]>=this->Input->GetBounds()[2] && pt[1]<=this->Input->GetBounds()[3] && pt[2]>=this->Input->GetBounds()[4] && pt[2]<=this->Input->GetBounds()[5]);
              // Get some parametric coords in [0,1]
              // [0,1] covers the whole dataset axis.
              pcoords[0]=(pi<<deltaLevel)*ratio;
              pcoords[1]=(pj<<deltaLevel)*ratio;
              pcoords[2]=(pk<<deltaLevel)*ratio;
683

Francois Bertel's avatar
Francois Bertel committed
684
685
686
              ptId=((pk<<deltaLevel)*resolution+(pj<<deltaLevel))*resolution
                +(pi<<deltaLevel);
              this->Triangulator->InsertPoint(ptId,pt,pcoords,0);
687
688


Francois Bertel's avatar
Francois Bertel committed
689
690
691
692
              // Test if the point is out or in the clipped part.
              // We have to put this code in the insertion loop of the
              // point because there is no method in vtkOrderedTriangulator
              // to access to inserted points.
693

Francois Bertel's avatar
Francois Bertel committed
694
695
              s=this->ClipFunction->FunctionValue(pt);
              if(this->InsideOut)
696
              {
Francois Bertel's avatar
Francois Bertel committed
697
                clipPoint=s<=this->Value; // keep point if true
698
              }
Francois Bertel's avatar
Francois Bertel committed
699
              else
700
              {
Francois Bertel's avatar
Francois Bertel committed
701
                clipPoint=s>=this->Value; // keep point if true
702
              }
Francois Bertel's avatar
Francois Bertel committed
703
              if(clipPoint)
704
              {
Francois Bertel's avatar
Francois Bertel committed
705
                allOut=0;
706
              }
Francois Bertel's avatar
Francois Bertel committed
707
              else
708
              {
Francois Bertel's avatar
Francois Bertel committed
709
                allIn=0;
710
              }
711

Francois Bertel's avatar
Francois Bertel committed
712
713
              ++x;
              ++pi;
714
            }
Francois Bertel's avatar
Francois Bertel committed
715
716
            ++y;
            ++pj;
717
          }
Francois Bertel's avatar
Francois Bertel committed
718
719
720
          ++z;
          ++pk;
        }
721
      }
Francois Bertel's avatar
Francois Bertel committed
722
      else
723
      {
Francois Bertel's avatar
Francois Bertel committed
724
725
726
727
        // this->Input->GetDimension()==2
        pt[2]=this->Input->GetOrigin()[2];
        y=0;
        while(y<2)
728
        {
Francois Bertel's avatar
Francois Bertel committed
729
730
          x=0;
          while(x<2)
731
          {
Francois Bertel's avatar
Francois Bertel committed
732
733
            pt[0]=bounds[x];
            pt[1]=bounds[2+y];
734

Francois Bertel's avatar
Francois Bertel committed
735
736
737
738
            // Test if the point is out or in the clipped part.
            // We have to put this code in the insertion loop of the
            // point because there is no method in vtkOrderedTriangulator
            // to access to inserted points.
739

Francois Bertel's avatar
Francois Bertel committed
740
741
            s=this->ClipFunction->FunctionValue(pt);
            if(this->InsideOut)
742
            {
Francois Bertel's avatar
Francois Bertel committed
743
              clipPoint=s<=this->Value; // keep point if true
744
            }
Francois Bertel's avatar
Francois Bertel committed
745
            else
746
            {
Francois Bertel's avatar
Francois Bertel committed
747
              clipPoint=s>=this->Value; // keep point if true
748
            }
Francois Bertel's avatar
Francois Bertel committed
749
            if(clipPoint)
750
            {
Francois Bertel's avatar
Francois Bertel committed
751
              allOut=0;
752
            }
Francois Bertel's avatar
Francois Bertel committed
753
            else
754
            {
Francois Bertel's avatar
Francois Bertel committed
755
756
              allIn=0;
            }
757
            ++x;
Francois Bertel's avatar
Francois Bertel committed
758
          }
759
          ++y;
Francois Bertel's avatar
Francois Bertel committed
760
        }
761
      }
762
763


Francois Bertel's avatar
Francois Bertel committed
764
765
766
767
768
769
770
771
772
      // see if we got a chance to either
      // 1. remove the leaf (!this->GenerateClippedOutput && allOut), no need
      // for triangulation, nor clipping, just skip the leaf.
      // 2. triangulate and passing the result without clipping each
      //    sub-tetra (allIn==1). Can work also if this->GenerateClippedOutput
      // is true. For one output, the sub-tetra will be passed (allIn), for the
      // other there will be nothing to pass or the clip. Or, if allIn is false
      // but allOut is true, there is nothing to do with the first output,
      // and passing everything to the second output.
773
774


Francois Bertel's avatar
Francois Bertel committed
775
      if(!this->GenerateClippedOutput && allOut)
776
      {
Kunda's avatar
Kunda committed
777
//        cout<<"this child is all out and we don't need to generate the other output"<<endl;
Francois Bertel's avatar
Francois Bertel committed
778
        return; // we've just save a lot of useless computation
779
      }
780

Francois Bertel's avatar
Francois Bertel committed
781
      int lastLevelLeaf=level>=(this->Input->GetNumberOfLevels()-1);
782

Francois Bertel's avatar
Francois Bertel committed
783
      if(this->Input->GetDimension()==3)
784
      {
Francois Bertel's avatar
Francois Bertel committed
785
        if(!lastLevelLeaf)
786
        {
Francois Bertel's avatar
Francois Bertel committed
787
788
789
790
791
          // Ok, now ask my parent if I have sibling with children on my
          // faces and even worst, if my parent has sibling with children
          // that have children on my face, or if the parent of my parent
          // has sibling with children that have children, that have children
          // on my face, until I reach the root...
792

Bill Lorensen's avatar
Bill Lorensen committed
793
          // list the 3 faces of the parent, the current node is laying on.
Francois Bertel's avatar
Francois Bertel committed
794
          int faces[3];
795

Francois Bertel's avatar
Francois Bertel committed
796
          int child=cursor->GetChildIndex();
797

Francois Bertel's avatar
Francois Bertel committed
798
799
800
          faces[0]=(child&1)==1; // false: -x, true: +x
          faces[1]=(child&2)==2; // false: -y, true: +y
          faces[2]=(child&4)==4; // false: -z, true: +z
801

Francois Bertel's avatar
Francois Bertel committed
802
803
          // sibling on faces that are not on a parent face
          int siblings[3];
804

Francois Bertel's avatar
Francois Bertel committed
805
806
807
          int inc=1;
          i=0;
          while(i<3)
808
          {
Francois Bertel's avatar
Francois Bertel committed
809
            if(faces[i])
810
            {
Francois Bertel's avatar
Francois Bertel committed
811
              siblings[i]=child-inc;
812
            }
Francois Bertel's avatar
Francois Bertel committed
813
            else
814
            {
Francois Bertel's avatar
Francois Bertel committed
815
              siblings[i]=child+inc;
816
            }
Francois Bertel's avatar
Francois Bertel committed
817
818
            ++i;
            inc<<=1;
819
          }
820

Francois Bertel's avatar
Francois Bertel committed
821
822
823
          this->Sibling->ToSameNode(cursor);
          this->Sibling->ToParent();
          // ask the 3 sibling, one on each face of the current node
824
          i=0;
Francois Bertel's avatar
Francois Bertel committed
825
826
          int faceOffset=0;
          while(i<3)
827
          {
Francois Bertel's avatar
Francois Bertel committed
828
829
830
            this->Sibling->ToChild(siblings[i]);
            assert("check: we are not visiting ourselves" && this->Sibling->GetChildIndex()!=child);
            if(!this->Sibling->CurrentIsLeaf())
831
            {
Francois Bertel's avatar
Francois Bertel committed
832
              assert("check: if the sibling is not a leaf we cannot be at the last level" && level<(this->Input->GetNumberOfLevels()-1));
833

Francois Bertel's avatar
Francois Bertel committed
834
835
836
              // get the points of this sibling on some given face
              int siblingFace=faceOffset;
              if(faces[i])
837
              {
Francois Bertel's avatar
Francois Bertel committed
838
839
                ++siblingFace;
              }
840
841
              this->Input->GetPointsOnFace(this->Sibling,siblingFace,level,this->Grabber);
            }
Francois Bertel's avatar
Francois Bertel committed
842
843
844
            this->Sibling->ToParent();
            ++i;
            faceOffset+=2;
845
          }
846

Francois Bertel's avatar
Francois Bertel committed
847
848
          // Get points on faces shared with the parent node.
          this->Input->GetPointsOnParentFaces(faces,level,cursor,this->Grabber);
849

Francois Bertel's avatar
Francois Bertel committed
850
851
852
853
854
          // Get the points from the edge-only neighbors.
          int childIndices[3];
          childIndices[2]=(child&4)>>2;
          childIndices[1]=(child&2)>>1;
          childIndices[0]=child&1;
855

Francois Bertel's avatar
Francois Bertel committed
856
857
858
859
860
861
          assert("check valid_range_c2" && childIndices[2]>=0 &&
                 childIndices[2]<=1);
          assert("check valid_range_c1" && childIndices[1]>=0 &&
                 childIndices[1]<=1);
          assert("check valid_range_c0" && childIndices[0]>=0 &&
                 childIndices[0]<=1);
862

Francois Bertel's avatar
Francois Bertel committed
863
864
865
866
          // First the edges aligned on X axis
          int axis=0;
          int a=2;
          int b=1;
867

Francois Bertel's avatar
Francois Bertel committed
868
869
          this->Sibling->ToSameNode(cursor);
          this->Sibling->ToParent();
870

Francois Bertel's avatar
Francois Bertel committed
871
          while(axis<3)
872
          {
Francois Bertel's avatar
Francois Bertel committed
873
874
            k=0;
            while(k<2)
875
            {
Francois Bertel's avatar
Francois Bertel committed
876
877
              j=0;
              while(j<2)
878
              {
Francois Bertel's avatar
Francois Bertel committed
879
                if(k!=childIndices[a]&&j!=childIndices[b])
880
                {
Francois Bertel's avatar
Francois Bertel committed
881
882
883
                  this->Sibling->ToChild((k<<a)+(j<<b)+
                                         (childIndices[axis]<<axis));
                  if(!this->Sibling->CurrentIsLeaf())
884
                  {
Francois Bertel's avatar
Francois Bertel committed
885
886
887
                    this->Input->GetPointsOnEdge(this->Sibling,level,axis,!k,
                                                 !j,this->Grabber);
                  }
888
889
                  this->Sibling->ToParent();
                }
Francois Bertel's avatar
Francois Bertel committed
890
                else
891
                {
Francois Bertel's avatar
Francois Bertel committed
892
893
894
                  this->Input->GetPointsOnParentEdge(cursor,level,axis,k,j,
                    this->Grabber);
                }
895
                ++j;
Francois Bertel's avatar
Francois Bertel committed
896
              }
897
898
              ++k;
            }
Francois Bertel's avatar
Francois Bertel committed
899
900
901
            ++axis;
            ++a;
            if(a>2)
902
            {
Francois Bertel's avatar
Francois Bertel committed
903
              a-=3;
904
            }
Francois Bertel's avatar
Francois Bertel committed
905
906
            ++b;
            if(b>2)
907
            {
Francois Bertel's avatar
Francois Bertel committed
908
909
              b-=3;
            }
910
911
912
          }
        } // if not leaf at last level
      }
Francois Bertel's avatar
Francois Bertel committed
913
      else
914
      {
Francois Bertel's avatar
Francois Bertel committed
915
916
        // this->Input->GetDimension()==2
        // counter- clockwise direction matters here.
917

Francois Bertel's avatar
Francois Bertel committed
918
919
        int edges[2];
        int child=cursor->GetChildIndex();
920

Francois Bertel's avatar
Francois Bertel committed
921
922
        this->Polygon->GetPointIds()->SetNumberOfIds(0);
        this->Polygon->GetPoints()->SetNumberOfPoints(0);
923

Francois Bertel's avatar
Francois Bertel committed
924
        if(!lastLevelLeaf)
925
        {
Francois Bertel's avatar
Francois Bertel committed
926
927
          this->Sibling->ToSameNode(cursor);
          this->Sibling->ToParent();
928

Bill Lorensen's avatar
Bill Lorensen committed
929
          // list the 2 edges of the parent, the current node is laying on.
Francois Bertel's avatar
Francois Bertel committed
930
931
          edges[0]=(child&1)==1; // false: -x, true: +x
          edges[1]=(child&