[2798] | 1 | /* -*- mode: c++; c-basic-offset: 4; indent-tabs-mode: nil -*- */ |
---|
[3502] | 2 | /* |
---|
| 3 | * Copyright (C) 2004-2013 HUBzero Foundation, LLC |
---|
| 4 | * |
---|
| 5 | */ |
---|
[931] | 6 | |
---|
[4889] | 7 | #include <cstdlib> |
---|
| 8 | |
---|
[3492] | 9 | #include <vrmath/Vector3f.h> |
---|
| 10 | |
---|
[4889] | 11 | #include "ReaderCommon.h" |
---|
| 12 | #include "GradientFilter.h" |
---|
[931] | 13 | |
---|
[4889] | 14 | using namespace nv; |
---|
[3492] | 15 | using namespace vrmath; |
---|
| 16 | |
---|
[931] | 17 | float * |
---|
[4889] | 18 | nv::merge(float *scalar, float *gradient, int size) |
---|
[931] | 19 | { |
---|
[2844] | 20 | float *data = (float *)malloc(sizeof(float) * 4 * size); |
---|
[931] | 21 | |
---|
[3492] | 22 | Vector3f *g = (Vector3f *)gradient; |
---|
[931] | 23 | |
---|
| 24 | int ngen = 0, sindex = 0; |
---|
[2844] | 25 | for (sindex = 0; sindex < size; ++sindex) { |
---|
| 26 | data[ngen++] = scalar[sindex]; |
---|
[2877] | 27 | data[ngen++] = 1.0 - g[sindex].x; |
---|
| 28 | data[ngen++] = 1.0 - g[sindex].y; |
---|
| 29 | data[ngen++] = 1.0 - g[sindex].z; |
---|
[931] | 30 | } |
---|
| 31 | return data; |
---|
| 32 | } |
---|
| 33 | |
---|
[3575] | 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 | */ |
---|
[931] | 46 | void |
---|
[4889] | 47 | nv::normalizeScalar(float *data, int count, int stride, double vmin, double vmax) |
---|
[931] | 48 | { |
---|
[3574] | 49 | double dv = vmax - vmin; |
---|
[3580] | 50 | dv = (dv == 0.0) ? 1.0 : dv; |
---|
[3574] | 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; |
---|
[3575] | 55 | } else if (v < vmin) { |
---|
| 56 | data[i] = 0.0f; |
---|
| 57 | } else if (v > vmax) { |
---|
| 58 | data[i] = 1.0f; |
---|
| 59 | } else { |
---|
[3574] | 60 | data[i] = (float)((v - vmin)/ dv); |
---|
[2844] | 61 | } |
---|
[931] | 62 | } |
---|
| 63 | } |
---|
| 64 | |
---|
[2877] | 65 | /** |
---|
[4904] | 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 | */ |
---|
| 73 | void |
---|
| 74 | nv::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[i ] = 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 | /** |
---|
[2877] | 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. |
---|
[3574] | 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 |
---|
[2877] | 104 | */ |
---|
[2844] | 105 | float * |
---|
[4889] | 106 | nv::computeGradient(float *data, |
---|
| 107 | int nx, int ny, int nz, |
---|
| 108 | float dx, float dy, float dz, |
---|
| 109 | float min, float max) |
---|
[931] | 110 | { |
---|
[3574] | 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 }; |
---|
[2877] | 115 | float spacing[3] = { dx, dy, dz }; |
---|
[3574] | 116 | computeGradients(tempGradients, data, sizes, spacing, DATRAW_FLOAT); |
---|
[931] | 117 | filterGradients(tempGradients, sizes); |
---|
| 118 | quantizeGradients(tempGradients, gradients, sizes, DATRAW_FLOAT); |
---|
[2874] | 119 | free(tempGradients); |
---|
[3574] | 120 | normalizeScalar(data, npts, 1, min, max); |
---|
| 121 | float *dataOut = merge(data, gradients, npts); |
---|
[2874] | 122 | free(gradients); |
---|
[3574] | 123 | return dataOut; |
---|
[931] | 124 | } |
---|
| 125 | |
---|
[2877] | 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 | */ |
---|
[956] | 146 | void |
---|
[4889] | 147 | nv::computeSimpleGradient(float *data, |
---|
| 148 | int nx, int ny, int nz, |
---|
| 149 | float dx, float dy, float dz) |
---|
[956] | 150 | { |
---|
[2877] | 151 | bool clampToEdge = true; |
---|
| 152 | double borderVal = 0.0; |
---|
| 153 | |
---|
| 154 | #define BORDER ((clampToEdge ? data[ngen] : borderVal)) |
---|
| 155 | |
---|
[956] | 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; |
---|
[2844] | 159 | for (int iz = 0; iz < nz; iz++) { |
---|
| 160 | for (int iy = 0; iy < ny; iy++) { |
---|
| 161 | for (int ix = 0; ix < nx; ix++) { |
---|
[956] | 162 | // gradient in x-direction |
---|
[2877] | 163 | double valm1 = (ix == 0) ? BORDER : data[ngen - 4]; |
---|
| 164 | double valp1 = (ix == nx-1) ? BORDER : data[ngen + 4]; |
---|
[2844] | 165 | if (valm1 < 0.0 || valp1 < 0.0) { |
---|
[956] | 166 | data[ngen+1] = 0.0; |
---|
| 167 | } else { |
---|
[2877] | 168 | data[ngen+1] = -(valp1-valm1)/(2. * dx); |
---|
| 169 | } |
---|
[956] | 170 | |
---|
| 171 | // gradient in y-direction |
---|
[2877] | 172 | valm1 = (iy == 0) ? BORDER : data[ngen - 4*nx]; |
---|
| 173 | valp1 = (iy == ny-1) ? BORDER : data[ngen + 4*nx]; |
---|
[2844] | 174 | if (valm1 < 0.0 || valp1 < 0.0) { |
---|
[956] | 175 | data[ngen+2] = 0.0; |
---|
| 176 | } else { |
---|
[2877] | 177 | data[ngen+2] = -(valp1-valm1)/(2. * dy); |
---|
[956] | 178 | } |
---|
| 179 | |
---|
| 180 | // gradient in z-direction |
---|
[2877] | 181 | valm1 = (iz == 0) ? BORDER : data[ngen - 4*nx*ny]; |
---|
| 182 | valp1 = (iz == nz-1) ? BORDER : data[ngen + 4*nx*ny]; |
---|
[2844] | 183 | if (valm1 < 0.0 || valp1 < 0.0) { |
---|
[956] | 184 | data[ngen+3] = 0.0; |
---|
| 185 | } else { |
---|
[2877] | 186 | data[ngen+3] = -(valp1-valm1)/(2. * dz); |
---|
[956] | 187 | } |
---|
[2877] | 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 | } |
---|
[956] | 202 | ngen += 4; |
---|
| 203 | } |
---|
| 204 | } |
---|
| 205 | } |
---|
[2877] | 206 | |
---|
| 207 | #undef BORDER |
---|
[956] | 208 | } |
---|