source: branches/1.2/packages/vizservers/nanovis/ReaderCommon.cpp @ 3589

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

Fix normalizeScalar

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