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

Last change on this file since 4385 was 4385, checked in by ldelgass, 10 years ago

Fix out-of-bounds points when resampling: don't consider out-of-boudnds points
in field min/max, set to normalized min.

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