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

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

Add basic VTK structured points reader to nanovis, update copyright dates.

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