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

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

Use normalizeScalar common function where possible. Also, don't set NaNs? to
-1 in unnormalized data since -1 could be a valid data value. Need to make
sure gradient filtering properly handles NaNs?.

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