source: trunk/packages/vizservers/nanovis/ReaderCommon.cpp @ 3611

Last change on this file since 3611 was 3611, checked in by ldelgass, 7 years ago

Use nv namespace for classes in nanovis rather than prefixing class names with
Nv (still need to convert shader classes).

  • Property svn:eol-style set to native
File size: 6.3 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 Compute Sobel filtered gradients for a 3D volume
67 *
68 * This technique is fairly expensive in terms of memory and
69 * running time due to the filter extent.
70 *
71 * \param data Data array with X the fastest running, stride of
72 * 1 float
73 * \param nx number of values in X direction
74 * \param ny number of values in Y direction
75 * \param nz number of values in Z direction
76 * \param dx Size of voxels in X direction
77 * \param dy Size of voxels in Y direction
78 * \param dz Size of voxels in Z direction
79 * \param min Minimum value in data
80 * \param max Maximum value in data
81 * \return Returns a float array with stride of 4 floats
82 * containing the normalized scalar and the x,y,z components of
83 * the (normalized) gradient vector
84 */
85float *
86nv::computeGradient(float *data,
87                    int nx, int ny, int nz,
88                    float dx, float dy, float dz,
89                    float min, float max)
90{
91    int npts = nx * ny * nz;
92    float *gradients = (float *)malloc(npts * 3 * sizeof(float));
93    float *tempGradients = (float *)malloc(npts * 3 * sizeof(float));
94    int sizes[3] = { nx, ny, nz };
95    float spacing[3] = { dx, dy, dz };
96    computeGradients(tempGradients, data, sizes, spacing, DATRAW_FLOAT);
97    filterGradients(tempGradients, sizes);
98    quantizeGradients(tempGradients, gradients, sizes, DATRAW_FLOAT);
99    free(tempGradients);
100    normalizeScalar(data, npts, 1, min, max);
101    float *dataOut = merge(data, gradients, npts);
102    free(gradients);
103    return dataOut;
104}
105
106/**
107 * \brief Compute gradients for a 3D volume with cubic cells
108 *
109 * The gradients are estimated using the central difference
110 * method.  This function assumes the data are normalized
111 * to [0,1] with missing data/NaNs represented by a negative
112 * value.
113 *
114 * \param data Data array with X the fastest running. There
115 * should be 4 floats allocated for each node, with the
116 * first float containing the scalar value.  The subsequent
117 * 3 floats will be filled with the x,y,z components of the
118 * gradient vector
119 * \param nx The number of nodes in the X direction
120 * \param ny The number of nodes in the Y direction
121 * \param nz The number of nodes in the Z direction
122 * \param dx The spacing (cell length) in the X direction
123 * \param dy The spacing (cell length) in the Y direction
124 * \param dz The spacing (cell length) in the Z direction
125 */
126void
127nv::computeSimpleGradient(float *data,
128                          int nx, int ny, int nz,
129                          float dx, float dy, float dz)
130{
131    bool clampToEdge = true;
132    double borderVal = 0.0;
133
134#define BORDER ((clampToEdge ? data[ngen] : borderVal))
135
136    // Compute the gradient of this data.  BE CAREFUL: center
137    // calculation on each node to avoid skew in either direction.
138    int ngen = 0;
139    for (int iz = 0; iz < nz; iz++) {
140        for (int iy = 0; iy < ny; iy++) {
141            for (int ix = 0; ix < nx; ix++) {
142                // gradient in x-direction
143                double valm1 = (ix == 0) ? BORDER : data[ngen - 4];
144                double valp1 = (ix == nx-1) ? BORDER : data[ngen + 4];
145                if (valm1 < 0.0 || valp1 < 0.0) {
146                    data[ngen+1] = 0.0;
147                } else {
148                    data[ngen+1] = -(valp1-valm1)/(2. * dx);
149                 }
150
151                // gradient in y-direction
152                valm1 = (iy == 0) ? BORDER : data[ngen - 4*nx];
153                valp1 = (iy == ny-1) ? BORDER : data[ngen + 4*nx];
154                if (valm1 < 0.0 || valp1 < 0.0) {
155                    data[ngen+2] = 0.0;
156                } else {
157                    data[ngen+2] = -(valp1-valm1)/(2. * dy);
158                }
159
160                // gradient in z-direction
161                valm1 = (iz == 0) ? BORDER : data[ngen - 4*nx*ny];
162                valp1 = (iz == nz-1) ? BORDER : data[ngen + 4*nx*ny];
163                if (valm1 < 0.0 || valp1 < 0.0) {
164                    data[ngen+3] = 0.0;
165                } else {
166                    data[ngen+3] = -(valp1-valm1)/(2. * dz);
167                }
168                // Normalize and scale/bias to [0,1] range
169                // The volume shader will expand to [-1,1]
170                double len = sqrt(data[ngen+1]*data[ngen+1] +
171                                  data[ngen+2]*data[ngen+2] +
172                                  data[ngen+3]*data[ngen+3]);
173                if (len < 1.0e-6) {
174                    data[ngen+1] = 0.0;
175                    data[ngen+2] = 0.0;
176                    data[ngen+3] = 0.0;
177                } else {
178                    data[ngen+1] = (data[ngen+1]/len + 1.0) * 0.5;
179                    data[ngen+2] = (data[ngen+2]/len + 1.0) * 0.5;
180                    data[ngen+3] = (data[ngen+3]/len + 1.0) * 0.5;
181                }
182                ngen += 4;
183            }
184        }
185    }
186
187#undef BORDER
188}
Note: See TracBrowser for help on using the repository browser.