vtkSimpleScalarTree.cxx 12.1 KB
Newer Older
1
2
3
4
5
/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkSimpleScalarTree.cxx

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

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

=========================================================================*/
#include "vtkSimpleScalarTree.h"
16

17
#include "vtkCell.h"
18
#include "vtkDataSet.h"
Ken Martin's avatar
Ken Martin committed
19
#include "vtkDoubleArray.h"
20
#include "vtkIdList.h"
21
#include "vtkObjectFactory.h"
22
#include "vtkPointData.h"
23

24
25
vtkStandardNewMacro(vtkSimpleScalarTree);

26
27
28
29
30
31
32
33
34
35
class vtkScalarNode {};

template <class TScalar>
class vtkScalarRange : public vtkScalarNode
{
public:
  TScalar min;
  TScalar max;
};

36
//-----------------------------------------------------------------------------
37
// Instantiate scalar tree with maximum level of 20 and branching
38
// factor of 3.
39
40
41
vtkSimpleScalarTree::vtkSimpleScalarTree()
{
  this->MaxLevel = 20;
42
  this->Level = 0;
43
44
45
  this->BranchingFactor = 3;
  this->Tree = NULL;
  this->TreeSize = 0;
46
47
48
49
50
51
52
53
54
55

  // Variables that support serial traversal
  this->NumCells = 0;
  this->TreeIndex = 0;
  this->ChildNumber = 0;
  this->CellId = 0;

  // For supporting parallel computing, list of possible candidates
  this->CandidateCells = NULL;
  this->NumCandidates = 0;
56
57
}

58
//-----------------------------------------------------------------------------
59
60
vtkSimpleScalarTree::~vtkSimpleScalarTree()
{
61
  delete [] this->Tree;
62
  if ( this->CandidateCells )
63
  {
64
65
    delete [] this->CandidateCells;
    this->CandidateCells = NULL;
66
  }
67
68
}

69
//-----------------------------------------------------------------------------
70
71
72
// Initialize locator. Frees memory and resets object as appropriate.
void vtkSimpleScalarTree::Initialize()
{
73
  delete [] this->Tree;
74
75
76
  this->Tree = NULL;
}

77
//-----------------------------------------------------------------------------
78
79
80
81
// Construct the scalar tree from the dataset provided. Checks build times
// and modified time from input and reconstructs the tree if necessaery.
void vtkSimpleScalarTree::BuildTree()
{
82
  vtkIdType cellId, i, j, numScalars;
83
84
85
86
  int level, offset, parentOffset, prod;
  vtkIdType numNodes, node, numLeafs, leaf, numParentLeafs;
  vtkCell *cell;
  vtkIdList *cellPts;
Ken Martin's avatar
Ken Martin committed
87
88
89
  vtkScalarRange<double> *tree, *parent;
  double *s;
  vtkDoubleArray *cellScalars;
90
91
92

  // Check input...see whether we have to rebuild
  //
93
94
  if ( !this->DataSet ||
       (this->NumCells = this->DataSet->GetNumberOfCells()) < 1 )
95
  {
96
97
    vtkErrorMacro( << "No data to build tree with");
    return;
98
  }
99

100
  if ( this->Tree != NULL && this->BuildTime > this->MTime
101
    && this->BuildTime > this->DataSet->GetMTime() )
102
  {
103
    return;
104
  }
105
106
107

  vtkDebugMacro( << "Building scalar tree..." );

108
109
  // If no scalars set then try and grab them from dataset
  if ( ! this->Scalars )
110
  {
111
    this->SetScalars(this->DataSet->GetPointData()->GetScalars());
112
  }
113
  if ( ! this->Scalars )
114
  {
115
116
    vtkErrorMacro( << "No scalar data to build trees with");
    return;
117
  }
118
119

  this->Initialize();
Ken Martin's avatar
Ken Martin committed
120
  cellScalars = vtkDoubleArray::New();
121
  cellScalars->Allocate(100);
122

123
124
  // Compute the number of levels in the tree
  //
Francois Bertel's avatar
Francois Bertel committed
125
  numLeafs = static_cast<int>(
126
    ceil(static_cast<double>(this->NumCells)/this->BranchingFactor));
127
  for (prod=1, numNodes=1, this->Level=0;
128
       prod < numLeafs && this->Level <= this->MaxLevel; this->Level++ )
129
  {
130
131
    prod *= this->BranchingFactor;
    numNodes += prod;
132
  }
133
134

  this->LeafOffset = offset = numNodes - prod;
Ken Martin's avatar
Ken Martin committed
135
  vtkScalarRange<double> *TTree;
136
  this->TreeSize = numNodes - (prod - numLeafs);
Ken Martin's avatar
Ken Martin committed
137
  this->Tree = TTree = new vtkScalarRange<double>[this->TreeSize];
138
  for ( i=0; i < this->TreeSize; i++ )
139
  {
Ken Martin's avatar
Ken Martin committed
140
141
    TTree[i].min = VTK_DOUBLE_MAX;
    TTree[i].max = -VTK_DOUBLE_MAX;
142
  }
143
144
145
146

  // Loop over all cells getting range of scalar data and place into leafs
  //
  for ( cellId=0, node=0; node < numLeafs; node++ )
147
  {
148
    tree = TTree + offset + node;
149
    for ( i=0; i < this->BranchingFactor && cellId < this->NumCells; i++, cellId++ )
150
    {
151
152
153
154
155
156
157
158
      cell = this->DataSet->GetCell(cellId);
      cellPts = cell->GetPointIds();
      numScalars = cellPts->GetNumberOfIds();
      cellScalars->SetNumberOfTuples(numScalars);
      this->Scalars->GetTuples(cellPts, cellScalars);
      s = cellScalars->GetPointer(0);

      for ( j=0; j < numScalars; j++ )
159
      {
160
        if ( s[j] < tree->min )
161
        {
162
          tree->min = s[j];
163
        }
164
        if ( s[j] > tree->max )
165
        {
166
167
168
169
          tree->max = s[j];
        }
      }
    }
170
  }
171
172
173
174

  // Now build top levels of tree in bottom-up fashion
  //
  for ( level=this->Level; level > 0; level-- )
175
  {
176
177
    parentOffset = offset - prod/this->BranchingFactor;
    prod /= this->BranchingFactor;
Francois Bertel's avatar
Francois Bertel committed
178
179
    numParentLeafs = static_cast<int>(
      ceil(static_cast<double>(numLeafs)/this->BranchingFactor));
180
181

    for ( leaf=0, node=0; node < numParentLeafs; node++ )
182
    {
183
184
      parent = TTree + parentOffset + node;
      for ( i=0; i < this->BranchingFactor && leaf < numLeafs; i++, leaf++ )
185
      {
186
187
        tree = TTree + offset + leaf;
        if ( tree->min < parent->min )
188
        {
189
          parent->min = tree->min;
190
        }
191
        if ( tree->max > parent->max )
192
        {
193
194
195
          parent->max = tree->max;
        }
      }
196
    }
197
198
199

    numLeafs = numParentLeafs;
    offset = parentOffset;
200
  }
201
202
203
204
205

  this->BuildTime.Modified();
  cellScalars->Delete();
}

206
//-----------------------------------------------------------------------------
207
208
// Begin to traverse the cells based on a scalar value. Returned cells
// will have scalar values that span the scalar value specified.
Ken Martin's avatar
Ken Martin committed
209
void vtkSimpleScalarTree::InitTraversal(double scalarValue)
210
{
Will Schroeder's avatar
Will Schroeder committed
211
  this->BuildTree();
212
  vtkScalarRange<double> *TTree =
Ken Martin's avatar
Ken Martin committed
213
    static_cast< vtkScalarRange<double> * > (this->Tree);
214
215
216
217
218
219
220

  this->ScalarValue = scalarValue;
  this->TreeIndex = this->TreeSize;

  // Check root of tree for overlap with scalar value
  //
  if ( TTree[0].min > scalarValue || TTree[0].max < scalarValue )
221
  {
222
    return;
223
  }
224

225
  else //find leaf that overlaps the isocontour value
226
  {
227
    this->FindStartLeaf(0,0); //recursive function call
228
  }
229
230
}

231
//-----------------------------------------------------------------------------
232
233
234
int vtkSimpleScalarTree::FindStartLeaf(vtkIdType index, int level)
{
  if ( level < this->Level )
235
  {
236
237
    int i;
    vtkIdType childIndex=this->BranchingFactor*index+1;
238

239
    level++;
240
    for ( i=0; i < this->BranchingFactor; i++ )
241
    {
242
243
      index = childIndex + i;
      if ( index >= this->TreeSize )
244
      {
245
246
        this->TreeIndex = this->TreeSize;
        return 0;
247
      }
248
      else if ( this->FindStartLeaf(childIndex+i, level) )
249
      {
250
251
        return 1;
      }
252
    }
253
254

    return 0;
255
  }
256
257

  else //recursion terminated
258
  {
259
    vtkScalarRange<double> *tree = static_cast<
Ken Martin's avatar
Ken Martin committed
260
      vtkScalarRange<double>*>(this->Tree) + index;
261
262

    if ( tree->min > this->ScalarValue || tree->max < this->ScalarValue )
263
    {
264
      return 0;
265
    }
266
    else
267
    {
268
269
270
271
272
      this->ChildNumber = 0;
      this->TreeIndex = index;
      this->CellId = (index - this->LeafOffset) * this->BranchingFactor;
      return 1;
    }
273
  }
274
275
}

276
//-----------------------------------------------------------------------------
277
278
279
280
281
282
283
284
285
286
287
int vtkSimpleScalarTree::FindNextLeaf(vtkIdType childIndex, int childLevel)
{
  vtkIdType myIndex=(childIndex-1)/this->BranchingFactor;
  int myLevel=childLevel-1;
  vtkIdType firstChildIndex, childNum, index;

  //Find which child invoked this method
  firstChildIndex = myIndex*this->BranchingFactor + 1;
  childNum = childIndex - firstChildIndex;

  for ( childNum++; childNum < this->BranchingFactor; childNum++ )
288
  {
289
    index = firstChildIndex + childNum;
290
    if ( index >= this->TreeSize )
291
    {
292
293
      this->TreeIndex = this->TreeSize;
      return 0;
294
    }
295
    else if ( this->FindStartLeaf(index, childLevel) )
296
    {
297
298
      return 1;
    }
299
  }
300
301
302

  //If here, didn't find anything yet
  if ( myLevel <= 0 ) //at root, can't go any higher in tree
303
  {
304
305
    this->TreeIndex = this->TreeSize;
    return 0;
306
  }
307
  else
308
  {
309
    return this->FindNextLeaf(myIndex,myLevel);
310
  }
311
312
313
}


314
//-----------------------------------------------------------------------------
315
316
317
318
// Return the next cell that may contain scalar value specified to
// initialize traversal. The value NULL is returned if the list is
// exhausted. Make sure that InitTraversal() has been invoked first or
// you'll get erratic behavior.
319
vtkCell *vtkSimpleScalarTree::GetNextCell(vtkIdType& cellId,
Ken Martin's avatar
Ken Martin committed
320
321
                                          vtkIdList* &cellPts,
                                          vtkDataArray *cellScalars)
322
{
Ken Martin's avatar
Ken Martin committed
323
  double s, min=VTK_DOUBLE_MAX, max=(-VTK_DOUBLE_MAX);
324
325
  vtkIdType i, numScalars;
  vtkCell *cell;
326
  vtkIdType numCells = this->NumCells;
327
328

  while ( this->TreeIndex < this->TreeSize )
329
  {
330
    for ( ; this->ChildNumber<this->BranchingFactor && this->CellId<numCells;
331
    this->ChildNumber++, this->CellId++ )
332
    {
333
334
335
      cell = this->DataSet->GetCell(this->CellId);
      cellPts = cell->GetPointIds();
      numScalars = cellPts->GetNumberOfIds();
Ken Martin's avatar
Ken Martin committed
336
337
      cellScalars->SetNumberOfTuples(numScalars);
      this->Scalars->GetTuples(cellPts, cellScalars);
338
      for (i=0; i < numScalars; i++)
339
      {
Ken Martin's avatar
Ken Martin committed
340
341
        s = cellScalars->GetTuple1(i);
        if ( s < min )
342
        {
Ken Martin's avatar
Ken Martin committed
343
          min = s;
344
        }
Ken Martin's avatar
Ken Martin committed
345
        if ( s > max )
346
        {
Ken Martin's avatar
Ken Martin committed
347
          max = s;
348
        }
349
      }
350
      if ( this->ScalarValue >= min && this->ScalarValue <= max )
351
      {
352
353
354
355
        cellId = this->CellId;
        this->ChildNumber++; //prepare for next time
        this->CellId++;
        return cell;
356
357
      }
    } //for each cell in this leaf
358
359
360

    // If here, must have not found anything in this leaf
    this->FindNextLeaf(this->TreeIndex, this->Level);
361
  } //while not all leafs visited
362
363
364
365

  return NULL;
}

366
367
368
369
370
371
372
373
//-----------------------------------------------------------------------------
// Return the number of cell batches.
vtkIdType vtkSimpleScalarTree::GetNumberOfCellBatches()
{
  // Basically we do a traversal of the tree and identify potential candidates.
  // It is essential that InitTraversal() has been called first.
  this->NumCandidates = 0;
  if ( this->CandidateCells )
374
  {
375
376
    delete [] this->CandidateCells;
    this->CandidateCells = NULL;
377
  }
378
  if ( this->NumCells < 1 )
379
  {
380
    return 0;
381
  }
382
383
384
385
  this->CandidateCells = new vtkIdType [this->NumCells];

  // Now begin traversing tree. InitTraversal() sets the first CellId.
  while ( this->TreeIndex < this->TreeSize )
386
  {
387
388
    for ( ; this->ChildNumber < this->BranchingFactor && this->CellId < this->NumCells;
          this->ChildNumber++, this->CellId++ )
389
    {
390
      this->CandidateCells[this->NumCandidates++] = this->CellId;
391
    } //for each cell in this leaf
392
393
394

    // If here, must have not found anything in this leaf
    this->FindNextLeaf(this->TreeIndex, this->Level);
395
  } //while not all leafs visited
396
397
398

  // Watch for boundary conditions
  if ( this->NumCandidates < 1 )
399
  {
400
    return 0;
401
  }
402
  else
403
  {
404
    return ( ((this->NumCandidates-1)/this->BranchingFactor) + 1);
405
  }
406
407
408
409
410
411
412
413
414
415
416
417
}

//-----------------------------------------------------------------------------
// Return the next batch of cells to process. Here we just use the branching
// factor (i.e., group size) as the batch size.
const vtkIdType* vtkSimpleScalarTree::
GetCellBatch(vtkIdType batchNum, vtkIdType& numCells)
{
  // Make sure that everything is hunky dory
  vtkIdType pos = batchNum * this->BranchingFactor;
  if ( this->NumCells < 1 || ! this->CandidateCells ||
       pos > this->NumCandidates )
418
  {
419
420
    numCells = 0;
    return NULL;
421
  }
422
423

  if ( (this->NumCandidates - pos) >= this->BranchingFactor )
424
  {
425
    numCells = this->BranchingFactor;
426
  }
427
  else
428
  {
429
    numCells = this->NumCandidates % this->BranchingFactor;
430
  }
431
432
433
434
  return this->CandidateCells + pos;
}

//-----------------------------------------------------------------------------
435
436
void vtkSimpleScalarTree::PrintSelf(ostream& os, vtkIndent indent)
{
437
  this->Superclass::PrintSelf(os,indent);
438

439
440
441
  os << indent << "Level: " << this->GetLevel() << "\n" ;
  os << indent << "Max Level: " << this->GetMaxLevel() << "\n" ;
  os << indent << "Branching Factor: " << this->GetBranchingFactor() << "\n" ;
442
}