[2798] | 1 | /* -*- mode: c++; c-basic-offset: 4; indent-tabs-mode: nil -*- */ |
---|
[3502] | 2 | /* |
---|
| 3 | * Copyright (C) 2004-2013 HUBzero Foundation, LLC |
---|
| 4 | * |
---|
| 5 | */ |
---|
| 6 | #include "ReaderCommon.h" |
---|
[931] | 7 | #include "GradientFilter.h" |
---|
| 8 | |
---|
[3492] | 9 | #include <vrmath/Vector3f.h> |
---|
| 10 | |
---|
[931] | 11 | #include "stdlib.h" |
---|
| 12 | |
---|
[3492] | 13 | using namespace vrmath; |
---|
| 14 | |
---|
[931] | 15 | float * |
---|
[2844] | 16 | merge(float *scalar, float *gradient, int size) |
---|
[931] | 17 | { |
---|
[2844] | 18 | float *data = (float *)malloc(sizeof(float) * 4 * size); |
---|
[931] | 19 | |
---|
[3492] | 20 | Vector3f *g = (Vector3f *)gradient; |
---|
[931] | 21 | |
---|
| 22 | int ngen = 0, sindex = 0; |
---|
[2844] | 23 | for (sindex = 0; sindex < size; ++sindex) { |
---|
| 24 | data[ngen++] = scalar[sindex]; |
---|
[2877] | 25 | data[ngen++] = 1.0 - g[sindex].x; |
---|
| 26 | data[ngen++] = 1.0 - g[sindex].y; |
---|
| 27 | data[ngen++] = 1.0 - g[sindex].z; |
---|
[931] | 28 | } |
---|
| 29 | return data; |
---|
| 30 | } |
---|
| 31 | |
---|
[3575] | 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 | */ |
---|
[931] | 44 | void |
---|
[3574] | 45 | normalizeScalar(float *data, int count, int stride, double vmin, double vmax) |
---|
[931] | 46 | { |
---|
[3574] | 47 | double dv = vmax - vmin; |
---|
[3580] | 48 | dv = (dv == 0.0) ? 1.0 : dv; |
---|
[3574] | 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; |
---|
[3575] | 53 | } else if (v < vmin) { |
---|
| 54 | data[i] = 0.0f; |
---|
| 55 | } else if (v > vmax) { |
---|
| 56 | data[i] = 1.0f; |
---|
| 57 | } else { |
---|
[3574] | 58 | data[i] = (float)((v - vmin)/ dv); |
---|
[2844] | 59 | } |
---|
[931] | 60 | } |
---|
| 61 | } |
---|
| 62 | |
---|
[2877] | 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. |
---|
[3574] | 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 |
---|
[2877] | 82 | */ |
---|
[2844] | 83 | float * |
---|
[3574] | 84 | computeGradient(float *data, |
---|
| 85 | int nx, int ny, int nz, |
---|
[2877] | 86 | float dx, float dy, float dz, |
---|
[2844] | 87 | float min, float max) |
---|
[931] | 88 | { |
---|
[3574] | 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 }; |
---|
[2877] | 93 | float spacing[3] = { dx, dy, dz }; |
---|
[3574] | 94 | computeGradients(tempGradients, data, sizes, spacing, DATRAW_FLOAT); |
---|
[931] | 95 | filterGradients(tempGradients, sizes); |
---|
| 96 | quantizeGradients(tempGradients, gradients, sizes, DATRAW_FLOAT); |
---|
[2874] | 97 | free(tempGradients); |
---|
[3574] | 98 | normalizeScalar(data, npts, 1, min, max); |
---|
| 99 | float *dataOut = merge(data, gradients, npts); |
---|
[2874] | 100 | free(gradients); |
---|
[3574] | 101 | return dataOut; |
---|
[931] | 102 | } |
---|
| 103 | |
---|
[2877] | 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 | */ |
---|
[956] | 124 | void |
---|
[3574] | 125 | computeSimpleGradient(float *data, |
---|
| 126 | int nx, int ny, int nz, |
---|
[2877] | 127 | float dx, float dy, float dz) |
---|
[956] | 128 | { |
---|
[2877] | 129 | bool clampToEdge = true; |
---|
| 130 | double borderVal = 0.0; |
---|
| 131 | |
---|
| 132 | #define BORDER ((clampToEdge ? data[ngen] : borderVal)) |
---|
| 133 | |
---|
[956] | 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; |
---|
[2844] | 137 | for (int iz = 0; iz < nz; iz++) { |
---|
| 138 | for (int iy = 0; iy < ny; iy++) { |
---|
| 139 | for (int ix = 0; ix < nx; ix++) { |
---|
[956] | 140 | // gradient in x-direction |
---|
[2877] | 141 | double valm1 = (ix == 0) ? BORDER : data[ngen - 4]; |
---|
| 142 | double valp1 = (ix == nx-1) ? BORDER : data[ngen + 4]; |
---|
[2844] | 143 | if (valm1 < 0.0 || valp1 < 0.0) { |
---|
[956] | 144 | data[ngen+1] = 0.0; |
---|
| 145 | } else { |
---|
[2877] | 146 | data[ngen+1] = -(valp1-valm1)/(2. * dx); |
---|
| 147 | } |
---|
[956] | 148 | |
---|
| 149 | // gradient in y-direction |
---|
[2877] | 150 | valm1 = (iy == 0) ? BORDER : data[ngen - 4*nx]; |
---|
| 151 | valp1 = (iy == ny-1) ? BORDER : data[ngen + 4*nx]; |
---|
[2844] | 152 | if (valm1 < 0.0 || valp1 < 0.0) { |
---|
[956] | 153 | data[ngen+2] = 0.0; |
---|
| 154 | } else { |
---|
[2877] | 155 | data[ngen+2] = -(valp1-valm1)/(2. * dy); |
---|
[956] | 156 | } |
---|
| 157 | |
---|
| 158 | // gradient in z-direction |
---|
[2877] | 159 | valm1 = (iz == 0) ? BORDER : data[ngen - 4*nx*ny]; |
---|
| 160 | valp1 = (iz == nz-1) ? BORDER : data[ngen + 4*nx*ny]; |
---|
[2844] | 161 | if (valm1 < 0.0 || valp1 < 0.0) { |
---|
[956] | 162 | data[ngen+3] = 0.0; |
---|
| 163 | } else { |
---|
[2877] | 164 | data[ngen+3] = -(valp1-valm1)/(2. * dz); |
---|
[956] | 165 | } |
---|
[2877] | 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 | } |
---|
[956] | 180 | ngen += 4; |
---|
| 181 | } |
---|
| 182 | } |
---|
| 183 | } |
---|
[2877] | 184 | |
---|
| 185 | #undef BORDER |
---|
[956] | 186 | } |
---|