vtkVisItScalarTree.C 10.7 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.
//
85 86 87
//   Brad Whitlock, Thu Jul 23 16:01:46 PDT 2015
//   Support for non-standard memory layout.
//
hrchilds's avatar
hrchilds committed
88
// ***************************************************************************
hrchilds's avatar
hrchilds committed
89 90 91 92

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

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

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

112
    if ( this->tree != NULL && BuildTime > MTime && BuildTime > DataSet->GetMTime() )
hrchilds's avatar
hrchilds committed
113 114
        return;
    
115
    vtkDataArray *scalars = this->DataSet->GetPointData()->GetScalars();
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

194 195 196 197 198 199 200 201 202
    int accessMethod = 0;
    if(scalars->HasStandardMemoryLayout())
    {
        if(scalars->GetDataType() == VTK_FLOAT)
            accessMethod = 1;
        else if(scalars->GetDataType() == VTK_DOUBLE)
            accessMethod = 2;
    }

hrchilds's avatar
hrchilds committed
203
    // Now find the min/max for each bucket.
204
    for (vtkIdType cellId = 0 ; cellId < this->nCells ; cellId++)
hrchilds's avatar
hrchilds committed
205
    {
206
        vtkIdType*pts;
207
        vtkIdType npts;
208
        vtkIdType arr8[8];
hrchilds's avatar
hrchilds committed
209 210 211
        
        // Get the points
        if (meshType==VTK_RECTILINEAR_GRID || meshType==VTK_STRUCTURED_GRID)
hrchilds's avatar
hrchilds committed
212
        {
hrchilds's avatar
hrchilds committed
213
            if (structured2D)
hrchilds's avatar
hrchilds committed
214
            {
215 216
                vtkIdType cellI = cellId % cell_dims[0];
                vtkIdType cellJ = cellId / strideY;
hrchilds's avatar
hrchilds committed
217

218
                for (vtkIdType j = 0 ; j < 4 ; ++j)
hrchilds's avatar
hrchilds committed
219
                {
hrchilds's avatar
hrchilds committed
220 221 222 223 224 225
                    arr8[j] = (cellI + X_val[j]) +
                              (cellJ + Y_val[j]) * ptstrideY; 
                } 

                pts = arr8;
                npts = 4;
hrchilds's avatar
hrchilds committed
226
            }
hrchilds's avatar
hrchilds committed
227
            else
hrchilds's avatar
hrchilds committed
228
            {
229 230 231 232 233
                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
234 235 236 237 238 239 240
                {
                    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
241
            }
hrchilds's avatar
hrchilds committed
242 243 244 245 246 247 248
        } 
        else
        { 
            DataSet->GetCellPoints(cellId, lst);
            npts = lst->GetNumberOfIds();
            pts = lst->GetPointer(0);
        }
hrchilds's avatar
hrchilds committed
249
            
250
        ScalarRange &leaf = this->tree[this->leafOffset + (cellId / this->bucketSize)];
251
        if(accessMethod == 1)
hrchilds's avatar
hrchilds committed
252
        {
253
            const float *fscalars = static_cast<const float *>(scalars->GetVoidPointer(0));
254 255 256 257 258 259 260 261 262
            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;
            }
        }
263
        else if(accessMethod == 2)
264
        {
265
            const double *dscalars = static_cast<const double *>(scalars->GetVoidPointer(0));
266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284
            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
285 286 287 288 289 290
        }
    }
    lst->Delete();


    //
hrchilds's avatar
hrchilds committed
291
    // Now propagate up the changes to the branches of the tree.
hrchilds's avatar
hrchilds committed
292
    //
293
    for (vtkIdType lev = this->levels - 1; lev >= 0; --lev)
hrchilds's avatar
hrchilds committed
294 295 296 297
    {
        //
        // The first index of each level is ((b ^ lev) - 1) / (b - 1)
        //
298
        vtkIdType offset;
299
        offset = (vtkIdType)((pow(double(this->BranchingFactor), double(lev)) - 1) /
300
                 (this->BranchingFactor - 1));
hrchilds's avatar
hrchilds committed
301

302
        vtkIdType len = (vtkIdType)pow(double(this->BranchingFactor), double(lev));
hrchilds's avatar
hrchilds committed
303 304


305
        vtkIdType cRow;
306
        cRow = (vtkIdType)((pow(double(this->BranchingFactor), double(lev + 1)) - 1)
307
                         / (this->BranchingFactor - 1));
hrchilds's avatar
hrchilds committed
308
        
309
        for (vtkIdType i = 0; i < len; ++i)
hrchilds's avatar
hrchilds committed
310
        {
311 312
            ScalarRange *parent = &(this->tree[offset + i]);
            for (vtkIdType j = 0 ; j < this->BranchingFactor ; j++)
hrchilds's avatar
hrchilds committed
313
            {
314
                const ScalarRange *child = &(this->tree[cRow + i * this->BranchingFactor + j]);
hrchilds's avatar
hrchilds committed
315 316 317 318 319 320 321 322 323 324 325 326
                if (child->min < parent->min)
                    parent->min = child->min;
                if (child->max > parent->max)
                    parent->max = child->max;
            }
        }
    }

    BuildTime.Modified();
}

void
327
vtkVisItScalarTree::GetCellList(double scalarValue, std::vector<vtkIdType> &v)
hrchilds's avatar
hrchilds committed
328 329 330
{
    BuildTree();
    
331
    if (!this->tree)
hrchilds's avatar
hrchilds committed
332 333
        return;
    
334 335
    this->searchValue = scalarValue;
    if (this->searchValue >= this->tree->min && this->searchValue <= this->tree->max)
hrchilds's avatar
hrchilds committed
336 337 338 339
        RecursiveSearch(v, 0, 0); 
}

void
340
vtkVisItScalarTree::RecursiveSearch(std::vector<vtkIdType> &v, vtkIdType index, vtkIdType lev)
hrchilds's avatar
hrchilds committed
341 342
{
    // Reached a leaf node
343
    if (lev == this->levels)
hrchilds's avatar
hrchilds committed
344
    {
345 346
        vtkIdType cId = (index - this->leafOffset) * this->bucketSize;
        for (vtkIdType i = 0; i < this->bucketSize && cId < this->nCells; ++i)
hrchilds's avatar
hrchilds committed
347 348 349 350
            v.push_back(cId++); 
        return;
    }

351 352
    vtkIdType childIndex = index * this->BranchingFactor + 1;
    for (vtkIdType i = 0; i < this->BranchingFactor; ++i)
hrchilds's avatar
hrchilds committed
353
    {
354 355
        vtkIdType ix = childIndex + i;
        if (this->searchValue >= this->tree[ix].min && this->searchValue <= this->tree[ix].max)
hrchilds's avatar
hrchilds committed
356 357 358
            RecursiveSearch(v, ix, lev + 1);
    }
}