vtkVisItScalarTree.C 10.3 KB
Newer Older
hrchilds's avatar
hrchilds 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
/*=========================================================================

  Program:   Visualization Toolkit
  Module:    $RCSfile: vtkVisItScalarTree.cxx,v $
  Language:  C++
  Date:      $Date: 2002/02/22 21:16:54 $
  Version:   $Revision: 1.66 $

  Copyright (c) 1993-2002 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 "vtkVisItScalarTree.h"

#include <vector>

#include <math.h>

#include <vtkDataSet.h>
hrchilds's avatar
hrchilds committed
25 26 27
#include <vtkIdList.h>
#include <vtkObjectFactory.h>
#include <vtkPointData.h>
hrchilds's avatar
hrchilds committed
28
#include <vtkRectilinearGrid.h>
hrchilds's avatar
hrchilds committed
29
#include <vtkStructuredGrid.h>
hrchilds's avatar
hrchilds committed
30 31 32 33 34 35 36

using std::vector;

vtkStandardNewMacro(vtkVisItScalarTree);

vtkVisItScalarTree::vtkVisItScalarTree()
{
37 38 39 40
    this->DataSet = NULL;
    this->MaxLevel = 20;
    this->BranchingFactor = 3;
    this->tree = NULL;
hrchilds's avatar
hrchilds committed
41 42
}

hrchilds's avatar
hrchilds committed
43 44 45 46 47 48 49 50
// ****************************************************************************
//  Modifications:
//
//    Hank Childs, Mon Nov  1 09:55:59 PST 2004
//    Fixed memory leak.
//
// ****************************************************************************

hrchilds's avatar
hrchilds committed
51 52
vtkVisItScalarTree::~vtkVisItScalarTree()
{
53 54 55 56
    if (this->tree)
        delete[] this->tree;
    if (this->DataSet != NULL)
        this->DataSet->UnRegister(this);
hrchilds's avatar
hrchilds committed
57 58 59 60 61
}

void
vtkVisItScalarTree::Initialize()
{
62 63 64
    if (this->tree)
        delete[] this->tree;
    this->tree = NULL;
hrchilds's avatar
hrchilds committed
65 66
}

hrchilds's avatar
hrchilds committed
67 68 69 70
// ***************************************************************************
// Modifications:
//   Akira Haddox, Wed Jul 30 09:27:22 PDT 2003
//   Fixed a bug in optimization for 2D structured / rectlinear grids.
hrchilds's avatar
hrchilds committed
71 72 73 74
//
//   Eric Brugger, Tue Jul 27 07:59:13 PDT 2004
//   Add several casts to fix compile errors.
//
hrchilds's avatar
hrchilds committed
75 76 77
//   Hank Childs, Mon Aug  9 08:39:04 PDT 2004
//   Fix UMR.
//
hrchilds's avatar
hrchilds committed
78 79 80 81
//   Hank Childs, Wed Mar  9 07:28:55 PST 2005
//   Fix additional memory issues.  Tree created in slightly different manner
//   now.
//
82 83 84
//   Brad Whitlock, Tue Mar 13 13:55:36 PDT 2012
//   Added double implementation. Changed int to vtkIdType.
//
hrchilds's avatar
hrchilds committed
85
// ***************************************************************************
hrchilds's avatar
hrchilds committed
86 87 88 89

void
vtkVisItScalarTree::BuildTree()
{
90
    if (!this->DataSet)
hrchilds's avatar
hrchilds committed
91 92 93 94 95
    {
        vtkErrorMacro( << "No data to build tree with");
        return;
    }

96 97
    this->nCells = this->DataSet->GetNumberOfCells();
    if (this->nCells < 1)
hrchilds's avatar
hrchilds committed
98 99 100 101 102
    {
        vtkErrorMacro( << "No data to build tree with");
        return;
    }

103
    if (!this->DataSet->GetPointData()->GetScalars())
hrchilds's avatar
hrchilds committed
104 105 106 107 108
    {
        vtkErrorMacro( << "No scalar data to build tree with");
        return;
    }

109
    if ( this->tree != NULL && BuildTime > MTime && BuildTime > DataSet->GetMTime() )
hrchilds's avatar
hrchilds committed
110 111
        return;
    
112 113 114 115
    vtkDataArray *scalars = this->DataSet->GetPointData()->GetScalars();
    float *fscalars = (float *)scalars->GetVoidPointer(0);
    double *dscalars = (double *)scalars->GetVoidPointer(0);
    int stype = scalars->GetDataType();
hrchilds's avatar
hrchilds committed
116 117 118 119 120 121 122

    //
    // Information for speeding up rectilinear and structured grids.
    // Based on code in vtkVisItContourFilter.
    //
    int meshType = DataSet->GetDataObjectType();
    
hrchilds's avatar
hrchilds committed
123
    int pt_dims[3] = { 1, 1, 1 };
hrchilds's avatar
hrchilds committed
124 125 126 127
    if (meshType == VTK_RECTILINEAR_GRID)
        ((vtkRectilinearGrid *)DataSet)->GetDimensions(pt_dims);
    else if (meshType == VTK_STRUCTURED_GRID)
        ((vtkStructuredGrid *)DataSet)->GetDimensions(pt_dims);
hrchilds's avatar
hrchilds committed
128 129 130 131 132 133 134 135 136 137 138 139

    //
    // Just to ensure we don't crash, 2D structured datasets
    // not in the XY plane we'll not optimize for.
    // 
    if ((meshType == VTK_RECTILINEAR_GRID ||
         meshType == VTK_STRUCTURED_GRID ) &&
        (pt_dims[0] == 1 || pt_dims[1] == 1))
    {
        meshType = VTK_POLY_DATA;
    }

hrchilds's avatar
hrchilds committed
140 141 142
    bool structured2D = false;
    if (meshType == VTK_RECTILINEAR_GRID || meshType == VTK_STRUCTURED_GRID)
        structured2D = (pt_dims[2] == 1);
hrchilds's avatar
hrchilds committed
143

hrchilds's avatar
hrchilds committed
144
    int cell_dims[3] = { 0, 0, 0 };
hrchilds's avatar
hrchilds committed
145 146 147
    cell_dims[0] = pt_dims[0] - 1;
    cell_dims[1] = pt_dims[1] - 1;
    cell_dims[2] = pt_dims[2] - 1;
148 149 150 151 152 153 154 155 156 157
    vtkIdType strideY = cell_dims[0];
    vtkIdType strideZ = cell_dims[0]*cell_dims[1];
    vtkIdType ptstrideY = pt_dims[0];
    vtkIdType ptstrideZ = pt_dims[0]*pt_dims[1];
    const vtkIdType X_val[8] = { 0, 1, 1, 0, 0, 1, 1, 0 };
    const vtkIdType Y_val[8] = { 0, 0, 1, 1, 0, 0, 1, 1 };
    const vtkIdType Z_val[8] = { 0, 0, 0, 0, 1, 1, 1, 1 };

    if (this->tree)
        delete[] this->tree;
hrchilds's avatar
hrchilds committed
158
    
159 160 161
    vtkIdType numLeafs = (vtkIdType)ceil(double (nCells) / BranchingFactor);
    vtkIdType count = 1;
    for (this->levels = 0; count < numLeafs ; ++this->levels)
hrchilds's avatar
hrchilds committed
162
    {
163
        count *= this->BranchingFactor;
hrchilds's avatar
hrchilds committed
164 165 166
    }
   
    // If levels is too high, adjust appropriately.
167
    if (this->levels > this->MaxLevel)
hrchilds's avatar
hrchilds committed
168
    {
169 170
        this->levels = this->MaxLevel;
        this->bucketSize = (vtkIdType)ceil(double (nCells) /
171
                     pow((double)BranchingFactor, double(levels)));
hrchilds's avatar
hrchilds committed
172 173
    }
    else
174
        this->bucketSize = BranchingFactor;
hrchilds's avatar
hrchilds committed
175 176
    
    // Size of a full tree is a geometric series: (b^(n+1)-1) / (b-1)
177
    this->treeSize = (vtkIdType)((pow(double(BranchingFactor), double(levels + 1)) - 1)
hrchilds's avatar
hrchilds committed
178 179
                        / (BranchingFactor - 1));

180
    this->tree = new ScalarRange[this->treeSize];
hrchilds's avatar
hrchilds committed
181 182 183 184

    //
    // First, fill up the leaves of the tree.
    //
185
    this->leafOffset = this->treeSize - (vtkIdType)pow(double(this->BranchingFactor), double(this->levels));
hrchilds's avatar
hrchilds committed
186
    vtkIdList *lst = vtkIdList::New();
187
    for (vtkIdType t = 0 ; t < this->treeSize ; t++)
hrchilds's avatar
hrchilds committed
188
    {
189
        ScalarRange &leaf = this->tree[t];
bonnell's avatar
bonnell committed
190 191
        leaf.min = VTK_FLOAT_MAX;
        leaf.max = VTK_FLOAT_MIN;
hrchilds's avatar
hrchilds committed
192
    }
hrchilds's avatar
hrchilds committed
193

hrchilds's avatar
hrchilds committed
194
    // Now find the min/max for each bucket.
195
    for (vtkIdType cellId = 0 ; cellId < this->nCells ; cellId++)
hrchilds's avatar
hrchilds committed
196
    {
197
        vtkIdType*pts;
198
        vtkIdType npts;
199
        vtkIdType arr8[8];
hrchilds's avatar
hrchilds committed
200 201 202
        
        // Get the points
        if (meshType==VTK_RECTILINEAR_GRID || meshType==VTK_STRUCTURED_GRID)
hrchilds's avatar
hrchilds committed
203
        {
hrchilds's avatar
hrchilds committed
204
            if (structured2D)
hrchilds's avatar
hrchilds committed
205
            {
206 207
                vtkIdType cellI = cellId % cell_dims[0];
                vtkIdType cellJ = cellId / strideY;
hrchilds's avatar
hrchilds committed
208

209
                for (vtkIdType j = 0 ; j < 4 ; ++j)
hrchilds's avatar
hrchilds committed
210
                {
hrchilds's avatar
hrchilds committed
211 212 213 214 215 216
                    arr8[j] = (cellI + X_val[j]) +
                              (cellJ + Y_val[j]) * ptstrideY; 
                } 

                pts = arr8;
                npts = 4;
hrchilds's avatar
hrchilds committed
217
            }
hrchilds's avatar
hrchilds committed
218
            else
hrchilds's avatar
hrchilds committed
219
            {
220 221 222 223 224
                vtkIdType cellI = cellId % cell_dims[0];
                vtkIdType cellJ = (cellId/strideY) % cell_dims[1];
                vtkIdType cellK = (cellId/strideZ);

                for (vtkIdType j = 0 ; j < 8 ; ++j)
hrchilds's avatar
hrchilds committed
225 226 227 228 229 230 231
                {
                    arr8[j] = (cellI + X_val[j]) +
                              (cellJ + Y_val[j])*ptstrideY 
                                + (cellK + Z_val[j])*ptstrideZ;
                }
                pts = arr8;
                npts = 8;
hrchilds's avatar
hrchilds committed
232
            }
hrchilds's avatar
hrchilds committed
233 234 235 236 237 238 239
        } 
        else
        { 
            DataSet->GetCellPoints(cellId, lst);
            npts = lst->GetNumberOfIds();
            pts = lst->GetPointer(0);
        }
hrchilds's avatar
hrchilds committed
240
            
241 242
        ScalarRange &leaf = this->tree[this->leafOffset + (cellId / this->bucketSize)];
        if(stype == VTK_FLOAT)
hrchilds's avatar
hrchilds committed
243
        {
244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273
            for (vtkIdType p = 0; p < npts; ++p)
            {
                double val = fscalars[pts[p]];
                if (val < leaf.min)
                    leaf.min = val;
                if (val > leaf.max)
                    leaf.max = val;
            }
        }
        else if(stype == VTK_DOUBLE)
        {
            for (vtkIdType p = 0; p < npts; ++p)
            {
                double val = dscalars[pts[p]];
                if (val < leaf.min)
                    leaf.min = val;
                if (val > leaf.max)
                    leaf.max = val;
            }
        }
        else
        {
            for (vtkIdType p = 0; p < npts; ++p)
            {
                double val = scalars->GetTuple1(pts[p]);
                if (val < leaf.min)
                    leaf.min = val;
                if (val > leaf.max)
                    leaf.max = val;
            }
hrchilds's avatar
hrchilds committed
274 275 276 277 278 279
        }
    }
    lst->Delete();


    //
hrchilds's avatar
hrchilds committed
280
    // Now propagate up the changes to the branches of the tree.
hrchilds's avatar
hrchilds committed
281
    //
282
    for (vtkIdType lev = this->levels - 1; lev >= 0; --lev)
hrchilds's avatar
hrchilds committed
283 284 285 286
    {
        //
        // The first index of each level is ((b ^ lev) - 1) / (b - 1)
        //
287
        vtkIdType offset;
288
        offset = (vtkIdType)((pow(double(this->BranchingFactor), double(lev)) - 1) /
289
                 (this->BranchingFactor - 1));
hrchilds's avatar
hrchilds committed
290

291
        vtkIdType len = (vtkIdType)pow(double(this->BranchingFactor), double(lev));
hrchilds's avatar
hrchilds committed
292 293


294
        vtkIdType cRow;
295
        cRow = (vtkIdType)((pow(double(this->BranchingFactor), double(lev + 1)) - 1)
296
                         / (this->BranchingFactor - 1));
hrchilds's avatar
hrchilds committed
297
        
298
        for (vtkIdType i = 0; i < len; ++i)
hrchilds's avatar
hrchilds committed
299
        {
300 301
            ScalarRange *parent = &(this->tree[offset + i]);
            for (vtkIdType j = 0 ; j < this->BranchingFactor ; j++)
hrchilds's avatar
hrchilds committed
302
            {
303
                const ScalarRange *child = &(this->tree[cRow + i * this->BranchingFactor + j]);
hrchilds's avatar
hrchilds committed
304 305 306 307 308 309 310 311 312 313 314 315
                if (child->min < parent->min)
                    parent->min = child->min;
                if (child->max > parent->max)
                    parent->max = child->max;
            }
        }
    }

    BuildTime.Modified();
}

void
316
vtkVisItScalarTree::GetCellList(double scalarValue, std::vector<vtkIdType> &v)
hrchilds's avatar
hrchilds committed
317 318 319
{
    BuildTree();
    
320
    if (!this->tree)
hrchilds's avatar
hrchilds committed
321 322
        return;
    
323 324
    this->searchValue = scalarValue;
    if (this->searchValue >= this->tree->min && this->searchValue <= this->tree->max)
hrchilds's avatar
hrchilds committed
325 326 327 328
        RecursiveSearch(v, 0, 0); 
}

void
329
vtkVisItScalarTree::RecursiveSearch(std::vector<vtkIdType> &v, vtkIdType index, vtkIdType lev)
hrchilds's avatar
hrchilds committed
330 331
{
    // Reached a leaf node
332
    if (lev == this->levels)
hrchilds's avatar
hrchilds committed
333
    {
334 335
        vtkIdType cId = (index - this->leafOffset) * this->bucketSize;
        for (vtkIdType i = 0; i < this->bucketSize && cId < this->nCells; ++i)
hrchilds's avatar
hrchilds committed
336 337 338 339
            v.push_back(cId++); 
        return;
    }

340 341
    vtkIdType childIndex = index * this->BranchingFactor + 1;
    for (vtkIdType i = 0; i < this->BranchingFactor; ++i)
hrchilds's avatar
hrchilds committed
342
    {
343 344
        vtkIdType ix = childIndex + i;
        if (this->searchValue >= this->tree[ix].min && this->searchValue <= this->tree[ix].max)
hrchilds's avatar
hrchilds committed
345 346 347
            RecursiveSearch(v, ix, lev + 1);
    }
}