source: nanovis/trunk/VtkDataSetReader.cpp @ 6295

Last change on this file since 6295 was 5717, checked in by ldelgass, 9 years ago

Downsample volume if dims exceed hardware limits
File size: 9.7 KB
Line 
1/* -*- mode: c++; c-basic-offset: 4; indent-tabs-mode: nil -*- */
2/*
3 * Copyright (C) 2004-2013  HUBzero Foundation, LLC
4 *
5 */
6
7#include <iostream>
8#include <float.h>
9
10#include <vtkSmartPointer.h>
11#include <vtkDataSet.h>
12#include <vtkImageData.h>
13#include <vtkPointData.h>
14#include <vtkCellData.h>
15#include <vtkDataArray.h>
16#include <vtkDataSetReader.h>
17#include <vtkCharArray.h>
18
19#include <vrmath/Vector3f.h>
20
21#include "ReaderCommon.h"
22#include "DataSetResample.h"
23
24#include "config.h"
25#include "nanovis.h"
26#include "VtkDataSetReader.h"
27#include "Volume.h"
28#include "Trace.h"
29
30using namespace nv;
31
32static void print(vtkDataSet *ds)
33{
34    if (ds == NULL)
35        return;
36
37    TRACE("DataSet class: %s", ds->GetClassName());
38
39    TRACE("DataSet memory: %g MiB", ds->GetActualMemorySize()/1024.);
40
41    double bounds[6];
42    ds->GetBounds(bounds);
43
44    // Topology
45    TRACE("DataSet bounds: %g %g %g %g %g %g",
46          bounds[0], bounds[1],
47          bounds[2], bounds[3],
48          bounds[4], bounds[5]);
49    TRACE("Points: %d Cells: %d", ds->GetNumberOfPoints(), ds->GetNumberOfCells());
50
51    double dataRange[2];
52    if (ds->GetPointData() != NULL) {
53        TRACE("PointData arrays: %d", ds->GetPointData()->GetNumberOfArrays());
54        for (int i = 0; i < ds->GetPointData()->GetNumberOfArrays(); i++) {
55            if (ds->GetPointData()->GetArray(i) != NULL) {
56                ds->GetPointData()->GetArray(i)->GetRange(dataRange, -1);
57                TRACE("PointData[%d]: '%s' comp:%d, (%g,%g)", i,
58                      ds->GetPointData()->GetArrayName(i),
59                      ds->GetPointData()->GetAbstractArray(i)->GetNumberOfComponents(),
60                      dataRange[0], dataRange[1]);
61            } else {
62                TRACE("PointData[%d]: '%s' comp:%d", i,
63                      ds->GetPointData()->GetArrayName(i),
64                      ds->GetPointData()->GetAbstractArray(i)->GetNumberOfComponents());
65            }
66        }
67        if (ds->GetPointData()->GetScalars() != NULL) {
68            TRACE("Active point scalars: %s", ds->GetPointData()->GetScalars()->GetName());
69        }
70        if (ds->GetPointData()->GetVectors() != NULL) {
71            TRACE("Active point vectors: %s", ds->GetPointData()->GetVectors()->GetName());
72        }
73    }
74    if (ds->GetCellData() != NULL) {
75        TRACE("CellData arrays: %d", ds->GetCellData()->GetNumberOfArrays());
76        for (int i = 0; i < ds->GetCellData()->GetNumberOfArrays(); i++) {
77            if (ds->GetCellData()->GetArray(i) != NULL) {
78                ds->GetCellData()->GetArray(i)->GetRange(dataRange, -1);
79                TRACE("CellData[%d]: '%s' comp:%d, (%g,%g)", i,
80                      ds->GetCellData()->GetArrayName(i),
81                      ds->GetCellData()->GetAbstractArray(i)->GetNumberOfComponents(),
82                      dataRange[0], dataRange[1]);
83            } else {
84                TRACE("CellData[%d]: '%s' comp:%d", i,
85                      ds->GetCellData()->GetArrayName(i),
86                      ds->GetCellData()->GetAbstractArray(i)->GetNumberOfComponents());
87            }
88        }
89        if (ds->GetCellData()->GetScalars() != NULL) {
90            TRACE("Active cell scalars: %s", ds->GetCellData()->GetScalars()->GetName());
91        }
92        if (ds->GetCellData()->GetVectors() != NULL) {
93            TRACE("Active cell vectors: %s", ds->GetCellData()->GetVectors()->GetName());
94        }
95    }
96    if (ds->GetFieldData() != NULL) {
97        TRACE("FieldData arrays: %d", ds->GetFieldData()->GetNumberOfArrays());
98        for (int i = 0; i < ds->GetFieldData()->GetNumberOfArrays(); i++) {
99            if (ds->GetFieldData()->GetArray(i) != NULL) {
100                ds->GetFieldData()->GetArray(i)->GetRange(dataRange, -1);
101                TRACE("FieldData[%d]: '%s' comp:%d, tuples:%d (%g,%g)", i,
102                      ds->GetFieldData()->GetArrayName(i),
103                      ds->GetFieldData()->GetAbstractArray(i)->GetNumberOfComponents(),
104                      ds->GetFieldData()->GetAbstractArray(i)->GetNumberOfTuples(),
105                      dataRange[0], dataRange[1]);
106            } else {
107                TRACE("FieldData[%d]: '%s' comp:%d, tuples:%d", i,
108                      ds->GetFieldData()->GetArrayName(i),
109                      ds->GetFieldData()->GetAbstractArray(i)->GetNumberOfComponents(),
110                      ds->GetFieldData()->GetAbstractArray(i)->GetNumberOfTuples());
111            }
112        }
113    }
114}
115
116Volume *
117nv::load_vtk_volume_stream(const char *tag, const char *bytes, int nBytes)
118{
119    TRACE("Enter tag:%s", tag);
120
121    vtkSmartPointer<vtkImageData> resampledDataSet;
122    {
123        vtkSmartPointer<vtkDataSet> dataSet;
124        {
125            vtkSmartPointer<vtkDataSetReader> reader = vtkSmartPointer<vtkDataSetReader>::New();
126            vtkSmartPointer<vtkCharArray> dataSetString = vtkSmartPointer<vtkCharArray>::New();
127
128            dataSetString->SetArray(const_cast<char *>(bytes), nBytes, 1);
129            reader->SetInputArray(dataSetString);
130            reader->ReadFromInputStringOn();
131            //reader->ReadAllNormalsOn();
132            //reader->ReadAllTCoordsOn();
133            reader->ReadAllScalarsOn();
134            //reader->ReadAllColorScalarsOn();
135            reader->ReadAllVectorsOn();
136            //reader->ReadAllTensorsOn();
137            reader->ReadAllFieldsOn();
138            reader->SetLookupTableName("");
139            reader->Update();
140
141            dataSet = reader->GetOutput();
142        }
143
144        if (dataSet == NULL)
145            return NULL;
146
147        print(dataSet);
148
149        resampledDataSet = vtkImageData::SafeDownCast(dataSet.GetPointer());
150        if (resampledDataSet != NULL) {
151            // Have a uniform grid, check if we need to resample
152#ifdef DOWNSAMPLE_DATA
153            int maxDim = 64;
154#else
155            // This is the hardware limit
156            int maxDim = NanoVis::max3dTextureSize;
157#endif
158            if (resampledDataSet->GetDimensions()[0] > maxDim ||
159                resampledDataSet->GetDimensions()[1] > maxDim ||
160                resampledDataSet->GetDimensions()[2] > maxDim) {
161                resampledDataSet = resampleVTKDataSet(dataSet, maxDim);
162            }
163        } else {
164            int maxDim = 64;
165            resampledDataSet = resampleVTKDataSet(dataSet, maxDim);
166        }
167    }
168
169    int nx, ny, nz, npts;
170    // origin
171    double x0, y0, z0;
172    // deltas (cell dimensions)
173    double dx, dy, dz;
174    // length of volume edges
175    double lx, ly, lz;
176
177    nx = resampledDataSet->GetDimensions()[0];
178    ny = resampledDataSet->GetDimensions()[1];
179    nz = resampledDataSet->GetDimensions()[2];
180    npts = nx * ny * nz;
181    resampledDataSet->GetSpacing(dx, dy, dz);
182    resampledDataSet->GetOrigin(x0, y0, z0);
183
184    lx = (nx - 1) * dx;
185    ly = (ny - 1) * dy;
186    lz = (nz - 1) * dz;
187
188    Volume *volume = NULL;
189    double vmin = DBL_MAX;
190    double nzero_min = DBL_MAX;
191    double vmax = -DBL_MAX;
192
193    float *data = new float[npts * 4];
194    memset(data, 0, npts * 4);
195
196    bool isVectorData = (resampledDataSet->GetPointData()->GetVectors() != NULL);
197
198    int ix = 0;
199    int iy = 0;
200    int iz = 0;
201    vtkDataArray *mask = resampledDataSet->GetPointData()->GetArray("vtkValidPointMask");
202    int numValidPoints = 0;
203    for (int p = 0; p < npts; p++) {
204        int nindex = p * 4;
205        double val;
206        int loc[3];
207        loc[0] = ix; loc[1] = iy; loc[2] = iz;
208        vtkIdType idx = resampledDataSet->ComputePointId(loc);
209        bool valid = (mask == NULL) ? true : (mask->GetComponent(idx, 0) != 0.0);
210        if (isVectorData) {
211            double vec[3];
212            resampledDataSet->GetPointData()->GetVectors()->GetTuple(idx, vec);
213            val = sqrt(vec[0]*vec[0] + vec[1]*vec[1] + vec[2]*vec[2]);
214            data[nindex] = (float)val;
215            data[nindex+1] = (float)vec[0];
216            data[nindex+2] = (float)vec[1];
217            data[nindex+3] = (float)vec[2];
218        } else {
219            //val = resampledDataSet->GetScalarComponentAsDouble(ix, iy, iz, 0);
220            val = resampledDataSet->GetPointData()->GetScalars()->GetComponent(idx, 0);
221            data[nindex] = valid ? (float)val : -FLT_MAX;
222        }
223        if (valid) {
224            numValidPoints++;
225            if (val < vmin) {
226                vmin = val;
227            } else if (val > vmax) {
228                vmax = val;
229            }
230            if (val != 0.0 && val < nzero_min) {
231                nzero_min = val;
232            }
233        }
234
235        if (++ix >= nx) {
236            ix = 0;
237            if (++iy >= ny) {
238                iy = 0;
239                ++iz;
240            }
241        }
242    }
243
244    if (numValidPoints < 1) {
245        ERROR("No valid points in resample volume");
246        vmin = 0.0;
247        vmax = 1.0;
248    } else {
249        TRACE("Valid points: %d of %d", numValidPoints, resampledDataSet->GetNumberOfPoints());
250    }
251
252    if (isVectorData) {
253        // Normalize magnitude [0,1] and vector components [0,1]
254        normalizeVector(data, npts, vmin, vmax);
255    } else {
256        // Normalize all values [0,1], -1 => out of bounds
257        normalizeScalar(data, npts, 4, vmin, vmax);
258        computeSimpleGradient(data, nx, ny, nz,
259                              dx, dy, dz);
260    }
261
262    TRACE("isVectorData: %s", isVectorData ? "yes" : "no");
263    TRACE("nx = %i ny = %i nz = %i", nx, ny, nz);
264    TRACE("x0 = %g y0 = %g z0 = %g", x0, y0, z0);
265    TRACE("lx = %g ly = %g lz = %g", lx, ly, lz);
266    TRACE("dx = %g dy = %g dz = %g", dx, dy, dz);
267    TRACE("dataMin = %g dataMax = %g nzero_min = %g",
268          vmin, vmax, nzero_min);
269
270    volume = NanoVis::loadVolume(tag, nx, ny, nz, 4, data,
271                                 vmin, vmax, nzero_min);
272    volume->xAxis.setRange(x0, x0 + lx);
273    volume->yAxis.setRange(y0, y0 + ly);
274    volume->zAxis.setRange(z0, z0 + lz);
275    volume->updatePending = true;
276
277    delete [] data;
278
279    TRACE("Leave");
280    return volume;
281}
Note: See TracBrowser for help on using the repository browser.