source: nanovis/trunk/ReaderCommon.cpp @ 4895

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

First pass at loading VTK vector data for flows in nanovis

  • Property svn:eol-style set to native
File size: 7.0 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 <cstdlib>
8
9#include <vrmath/Vector3f.h>
10
11#include "ReaderCommon.h"
12#include "GradientFilter.h"
13
14using namespace nv;
15using namespace vrmath;
16
17float *
18nv::merge(float *scalar, float *gradient, int size)
19{
20    float *data = (float *)malloc(sizeof(float) * 4 * size);
21
22    Vector3f *g = (Vector3f *)gradient;
23
24    int ngen = 0, sindex = 0;
25    for (sindex = 0; sindex < size; ++sindex) {
26        data[ngen++] = scalar[sindex];
27        data[ngen++] = 1.0 - g[sindex].x;
28        data[ngen++] = 1.0 - g[sindex].y;
29        data[ngen++] = 1.0 - g[sindex].z;
30    }
31    return data;
32}
33
34/**
35 * \brief Normalize data to [0,1] based on vmin,vmax range
36 *
37 * Data outside of given range is clamped, and NaNs are set to
38 * -1 in the output
39 *
40 * \param data Float array of unnormalized data, will be normalized on return
41 * \param count Number of elts in array
42 * \param stride Stride between values in data array
43 * \param vmin Minimum value in data array
44 * \param vmax Maximum value in data array
45 */
46void
47nv::normalizeScalar(float *data, int count, int stride, double vmin, double vmax)
48{
49    double dv = vmax - vmin;
50    dv = (dv == 0.0) ? 1.0 : dv;
51    for (int pt = 0, i = 0; pt < count; ++pt, i += stride) {
52        double v = data[i];
53        if (isnan(v)) {
54            data[i] = -1.0f;
55        } else if (v < vmin) {
56            data[i] = 0.0f;
57        } else if (v > vmax) {
58            data[i] = 1.0f;
59        } else {
60            data[i] = (float)((v - vmin)/ dv);
61        }
62    }
63}
64
65/**
66 * \brief Normalize data to [0,1] based on vmin,vmax range
67 *
68 * Data outside of given range is clamped, and NaNs are set to
69 * -1 in the output
70 *
71 * \param data Float array of unnormalized data, will be normalized on return
72 * \param count Number of elts in array
73 * \param stride Stride between values in data array
74 * \param vmin Minimum value in data array
75 * \param vmax Maximum value in data array
76 */
77void
78nv::normalizeVector(float *data, int count, double vmin, double vmax)
79{
80    for (int p = 0; p < count; p++) {
81        int i = p * 4;
82        data[] = data[i]/vmax;
83        data[i+1] = data[i+1]/(2.0 * vmax) + 0.5;
84        data[i+2] = data[i+2]/(2.0 * vmax) + 0.5;
85        data[i+3] = data[i+3]/(2.0 * vmax) + 0.5;
86    }
87}
88
89/**
90 * \brief Compute Sobel filtered gradients for a 3D volume
91 *
92 * This technique is fairly expensive in terms of memory and
93 * running time due to the filter extent.
94 *
95 * \param data Data array with X the fastest running, stride of
96 * 1 float
97 * \param nx number of values in X direction
98 * \param ny number of values in Y direction
99 * \param nz number of values in Z direction
100 * \param dx Size of voxels in X direction
101 * \param dy Size of voxels in Y direction
102 * \param dz Size of voxels in Z direction
103 * \param min Minimum value in data
104 * \param max Maximum value in data
105 * \return Returns a float array with stride of 4 floats
106 * containing the normalized scalar and the x,y,z components of
107 * the (normalized) gradient vector
108 */
109float *
110nv::computeGradient(float *data,
111                    int nx, int ny, int nz,
112                    float dx, float dy, float dz,
113                    float min, float max)
114{
115    int npts = nx * ny * nz;
116    float *gradients = (float *)malloc(npts * 3 * sizeof(float));
117    float *tempGradients = (float *)malloc(npts * 3 * sizeof(float));
118    int sizes[3] = { nx, ny, nz };
119    float spacing[3] = { dx, dy, dz };
120    computeGradients(tempGradients, data, sizes, spacing, DATRAW_FLOAT);
121    filterGradients(tempGradients, sizes);
122    quantizeGradients(tempGradients, gradients, sizes, DATRAW_FLOAT);
123    free(tempGradients);
124    normalizeScalar(data, npts, 1, min, max);
125    float *dataOut = merge(data, gradients, npts);
126    free(gradients);
127    return dataOut;
128}
129
130/**
131 * \brief Compute gradients for a 3D volume with cubic cells
132 *
133 * The gradients are estimated using the central difference
134 * method.  This function assumes the data are normalized
135 * to [0,1] with missing data/NaNs represented by a negative
136 * value.
137 *
138 * \param data Data array with X the fastest running. There
139 * should be 4 floats allocated for each node, with the
140 * first float containing the scalar value.  The subsequent
141 * 3 floats will be filled with the x,y,z components of the
142 * gradient vector
143 * \param nx The number of nodes in the X direction
144 * \param ny The number of nodes in the Y direction
145 * \param nz The number of nodes in the Z direction
146 * \param dx The spacing (cell length) in the X direction
147 * \param dy The spacing (cell length) in the Y direction
148 * \param dz The spacing (cell length) in the Z direction
149 */
150void
151nv::computeSimpleGradient(float *data,
152                          int nx, int ny, int nz,
153                          float dx, float dy, float dz)
154{
155    bool clampToEdge = true;
156    double borderVal = 0.0;
157
158#define BORDER ((clampToEdge ? data[ngen] : borderVal))
159
160    // Compute the gradient of this data.  BE CAREFUL: center
161    // calculation on each node to avoid skew in either direction.
162    int ngen = 0;
163    for (int iz = 0; iz < nz; iz++) {
164        for (int iy = 0; iy < ny; iy++) {
165            for (int ix = 0; ix < nx; ix++) {
166                // gradient in x-direction
167                double valm1 = (ix == 0) ? BORDER : data[ngen - 4];
168                double valp1 = (ix == nx-1) ? BORDER : data[ngen + 4];
169                if (valm1 < 0.0 || valp1 < 0.0) {
170                    data[ngen+1] = 0.0;
171                } else {
172                    data[ngen+1] = -(valp1-valm1)/(2. * dx);
173                 }
174
175                // gradient in y-direction
176                valm1 = (iy == 0) ? BORDER : data[ngen - 4*nx];
177                valp1 = (iy == ny-1) ? BORDER : data[ngen + 4*nx];
178                if (valm1 < 0.0 || valp1 < 0.0) {
179                    data[ngen+2] = 0.0;
180                } else {
181                    data[ngen+2] = -(valp1-valm1)/(2. * dy);
182                }
183
184                // gradient in z-direction
185                valm1 = (iz == 0) ? BORDER : data[ngen - 4*nx*ny];
186                valp1 = (iz == nz-1) ? BORDER : data[ngen + 4*nx*ny];
187                if (valm1 < 0.0 || valp1 < 0.0) {
188                    data[ngen+3] = 0.0;
189                } else {
190                    data[ngen+3] = -(valp1-valm1)/(2. * dz);
191                }
192                // Normalize and scale/bias to [0,1] range
193                // The volume shader will expand to [-1,1]
194                double len = sqrt(data[ngen+1]*data[ngen+1] +
195                                  data[ngen+2]*data[ngen+2] +
196                                  data[ngen+3]*data[ngen+3]);
197                if (len < 1.0e-6) {
198                    data[ngen+1] = 0.0;
199                    data[ngen+2] = 0.0;
200                    data[ngen+3] = 0.0;
201                } else {
202                    data[ngen+1] = (data[ngen+1]/len + 1.0) * 0.5;
203                    data[ngen+2] = (data[ngen+2]/len + 1.0) * 0.5;
204                    data[ngen+3] = (data[ngen+3]/len + 1.0) * 0.5;
205                }
206                ngen += 4;
207            }
208        }
209    }
210
211#undef BORDER
212}
Note: See TracBrowser for help on using the repository browser.