source: nanovis/trunk/ReaderCommon.cpp @ 4923

Last change on this file since 4923 was 4897, checked in by ldelgass, 5 years ago

add header deps

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