source: trunk/packages/vizservers/nanovis/VtkDataSetReader.cpp @ 3943

Last change on this file since 3943 was 3943, checked in by ldelgass, 6 years ago

Fix vector magnitude computation in nanovis VTK reader

File size: 5.2 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 <vtkDataArray.h>
15#include <vtkDataSetReader.h>
16#include <vtkCharArray.h>
17
18#include <RpOutcome.h>
19
20#include <vrmath/Vector3f.h>
21
22#include "ReaderCommon.h"
23#include "DataSetResample.h"
24
25#include "config.h"
26#include "nanovis.h"
27#include "VtkDataSetReader.h"
28#include "Volume.h"
29#include "Trace.h"
30
31using namespace nv;
32
33Volume *
34nv::load_vtk_volume_stream(Rappture::Outcome& result, const char *tag, const char *bytes, int nBytes)
35{
36    TRACE("Enter tag:%s", tag);
37
38    vtkSmartPointer<vtkImageData> resampledDataSet;
39    {
40        vtkSmartPointer<vtkDataSet> dataSet;
41        {
42            vtkSmartPointer<vtkDataSetReader> reader = vtkSmartPointer<vtkDataSetReader>::New();
43            vtkSmartPointer<vtkCharArray> dataSetString = vtkSmartPointer<vtkCharArray>::New();
44
45            dataSetString->SetArray(const_cast<char *>(bytes), nBytes, 1);
46            reader->SetInputArray(dataSetString);
47            reader->ReadFromInputStringOn();
48            //reader->ReadAllNormalsOn();
49            //reader->ReadAllTCoordsOn();
50            reader->ReadAllScalarsOn();
51            //reader->ReadAllColorScalarsOn();
52            reader->ReadAllVectorsOn();
53            //reader->ReadAllTensorsOn();
54            //reader->ReadAllFieldsOn();
55            reader->SetLookupTableName("");
56            reader->Update();
57
58            dataSet = reader->GetOutput();
59        }
60
61        if (dataSet == NULL)
62            return NULL;
63
64        int maxDim = 64;
65
66        resampledDataSet = vtkImageData::SafeDownCast(dataSet.GetPointer());
67        if (resampledDataSet != NULL) {
68            // Have a uniform grid, check if we need to resample
69#ifdef DOWNSAMPLE_DATA
70            if (resampledDataSet->GetDimensions()[0] > maxDim ||
71                resampledDataSet->GetDimensions()[1] > maxDim ||
72                resampledDataSet->GetDimensions()[2] > maxDim) {
73                resampledDataSet = resampleVTKDataSet(dataSet, maxDim);
74            }
75#endif
76        } else {
77            resampledDataSet = resampleVTKDataSet(dataSet, maxDim);
78        }
79    }
80
81    int nx, ny, nz, npts;
82    // origin
83    double x0, y0, z0;
84    // deltas (cell dimensions)
85    double dx, dy, dz;
86    // length of volume edges
87    double lx, ly, lz;
88
89    nx = resampledDataSet->GetDimensions()[0];
90    ny = resampledDataSet->GetDimensions()[1];
91    nz = resampledDataSet->GetDimensions()[2];
92    npts = nx * ny * nz;
93    resampledDataSet->GetSpacing(dx, dy, dz);
94    resampledDataSet->GetOrigin(x0, y0, z0);
95
96    lx = (nx - 1) * dx;
97    ly = (ny - 1) * dy;
98    lz = (nz - 1) * dz;
99
100    Volume *volume = NULL;
101    double vmin = DBL_MAX;
102    double nzero_min = DBL_MAX;
103    double vmax = -DBL_MAX;
104
105    float *data = new float[npts * 4];
106    memset(data, 0, npts * 4);
107
108    bool isVectorData = (resampledDataSet->GetPointData()->GetVectors() != NULL);
109
110    int ix = 0;
111    int iy = 0;
112    int iz = 0;
113    for (int p = 0; p < npts; p++) {
114        int nindex = p * 4;
115        double val;
116        if (isVectorData) {
117            double vec[3];
118            int loc[3];
119            loc[0] = ix; loc[1] = iy; loc[2] = iz;
120            vtkIdType idx = resampledDataSet->ComputePointId(loc);
121            resampledDataSet->GetPointData()->GetVectors()->GetTuple(idx, vec);
122            val = sqrt(vec[0]*vec[0] + vec[1]*vec[1] + vec[2]*vec[2]);
123            data[nindex] = (float)val;
124            data[nindex+1] = (float)vec[0];
125            data[nindex+2] = (float)vec[1];
126            data[nindex+3] = (float)vec[2];
127        } else {
128            val = resampledDataSet->GetScalarComponentAsDouble(ix, iy, iz, 0);
129            data[nindex] = (float)val;
130        }
131         if (val < vmin) {
132            vmin = val;
133        } else if (val > vmax) {
134            vmax = val;
135        }
136        if (val != 0.0 && val < nzero_min) {
137            nzero_min = val;
138        }
139
140        if (++ix >= nx) {
141            ix = 0;
142            if (++iy >= ny) {
143                iy = 0;
144                ++iz;
145            }
146        }
147    }
148
149    if (isVectorData) {
150        // Normalize magnitude [0,1] and vector components [0,1]
151        normalizeVector(data, npts, vmin, vmax);
152    } else {
153        // scale all values [0-1], -1 => out of bounds
154        normalizeScalar(data, npts, 4, vmin, vmax);
155        computeSimpleGradient(data, nx, ny, nz,
156                              dx, dy, dz);
157    }
158
159    TRACE("isVectorData: %s", isVectorData ? "yes" : "no");
160    TRACE("nx = %i ny = %i nz = %i", nx, ny, nz);
161    TRACE("x0 = %lg y0 = %lg z0 = %lg", x0, y0, z0);
162    TRACE("lx = %lg ly = %lg lz = %lg", lx, ly, lz);
163    TRACE("dx = %lg dy = %lg dz = %lg", dx, dy, dz);
164    TRACE("dataMin = %lg dataMax = %lg nzero_min = %lg",
165          vmin, vmax, nzero_min);
166
167    volume = NanoVis::loadVolume(tag, nx, ny, nz, 4, data,
168                                 vmin, vmax, nzero_min);
169    volume->xAxis.setRange(x0, x0 + lx);
170    volume->yAxis.setRange(y0, y0 + ly);
171    volume->zAxis.setRange(z0, z0 + lz);
172    volume->updatePending = true;
173
174    delete [] data;
175
176    TRACE("Leave");
177    return volume;
178}
Note: See TracBrowser for help on using the repository browser.