source: trunk/packages/vizservers/nanovis/dxReaderCommon.cpp @ 3492

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

Fix camera reset for nanovis. Includes refactoring of vector/matrix classes
in nanovis to consolidate into vrmath library. Also add preliminary canonical
view control to clients for testing.

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