1 | /* -*- mode: c++; c-basic-offset: 4; indent-tabs-mode: nil -*- */ |
---|
2 | #include "dxReaderCommon.h" |
---|
3 | #include "GradientFilter.h" |
---|
4 | |
---|
5 | #include "Vector3.h" |
---|
6 | #include "stdlib.h" |
---|
7 | |
---|
8 | float * |
---|
9 | merge(float *scalar, float *gradient, int size) |
---|
10 | { |
---|
11 | float *data = (float *)malloc(sizeof(float) * 4 * size); |
---|
12 | |
---|
13 | Vector3 *g = (Vector3 *)gradient; |
---|
14 | |
---|
15 | int ngen = 0, sindex = 0; |
---|
16 | for (sindex = 0; sindex < size; ++sindex) { |
---|
17 | data[ngen++] = scalar[sindex]; |
---|
18 | data[ngen++] = g[sindex].x; |
---|
19 | data[ngen++] = g[sindex].y; |
---|
20 | data[ngen++] = g[sindex].z; |
---|
21 | } |
---|
22 | return data; |
---|
23 | } |
---|
24 | |
---|
25 | void |
---|
26 | normalizeScalar(float *fdata, int count, float min, float max) |
---|
27 | { |
---|
28 | float v = max - min; |
---|
29 | if (v != 0.0f) { |
---|
30 | for (int i = 0; i < count; ++i) { |
---|
31 | fdata[i] = fdata[i] / v; |
---|
32 | } |
---|
33 | } |
---|
34 | } |
---|
35 | |
---|
36 | float * |
---|
37 | computeGradient(float *fdata, int width, int height, int depth, |
---|
38 | float min, float max) |
---|
39 | { |
---|
40 | float *gradients = (float *)malloc(width * height * depth * 3 * |
---|
41 | sizeof(float)); |
---|
42 | float *tempGradients = (float *)malloc(width * height * depth * 3 * |
---|
43 | sizeof(float)); |
---|
44 | int sizes[3] = { width, height, depth }; |
---|
45 | computeGradients(tempGradients, fdata, sizes, DATRAW_FLOAT); |
---|
46 | filterGradients(tempGradients, sizes); |
---|
47 | quantizeGradients(tempGradients, gradients, sizes, DATRAW_FLOAT); |
---|
48 | normalizeScalar(fdata, width * height * depth, min, max); |
---|
49 | float *data = merge(fdata, gradients, width * height * depth); |
---|
50 | return data; |
---|
51 | } |
---|
52 | |
---|
53 | void |
---|
54 | computeSimpleGradient(float *data, int nx, int ny, int nz) |
---|
55 | { |
---|
56 | // Compute the gradient of this data. BE CAREFUL: center |
---|
57 | // calculation on each node to avoid skew in either direction. |
---|
58 | int ngen = 0; |
---|
59 | for (int iz = 0; iz < nz; iz++) { |
---|
60 | for (int iy = 0; iy < ny; iy++) { |
---|
61 | for (int ix = 0; ix < nx; ix++) { |
---|
62 | // gradient in x-direction |
---|
63 | double valm1 = (ix == 0) ? 0.0 : data[ngen - 4]; |
---|
64 | double valp1 = (ix == nx-1) ? 0.0 : data[ngen + 4]; |
---|
65 | if (valm1 < 0.0 || valp1 < 0.0) { |
---|
66 | data[ngen+1] = 0.0; |
---|
67 | } else { |
---|
68 | data[ngen+1] = valp1-valm1; // assume dx=1 |
---|
69 | //data[ngen+1] = ((valp1-valm1) + 1.0) * 0.5; // assume dx=1 (ISO) |
---|
70 | } |
---|
71 | |
---|
72 | // gradient in y-direction |
---|
73 | valm1 = (iy == 0) ? 0.0 : data[ngen - 4*nx]; |
---|
74 | valp1 = (iy == ny-1) ? 0.0 : data[ngen + 4*nx]; |
---|
75 | if (valm1 < 0.0 || valp1 < 0.0) { |
---|
76 | data[ngen+2] = 0.0; |
---|
77 | } else { |
---|
78 | data[ngen+2] = valp1-valm1; // assume dy=1 |
---|
79 | //data[ngen+2] = ((valp1-valm1) + .01) * 0.5; // assume dy=1 (ISO) |
---|
80 | } |
---|
81 | |
---|
82 | // gradient in z-direction |
---|
83 | valm1 = (iz == 0) ? 0.0 : data[ngen - 4*nx*ny]; |
---|
84 | valp1 = (iz == nz-1) ? 0.0 : data[ngen + 4*nx*ny]; |
---|
85 | if (valm1 < 0.0 || valp1 < 0.0) { |
---|
86 | data[ngen+3] = 0.0; |
---|
87 | } else { |
---|
88 | data[ngen+3] = valp1-valm1; // assume dz=1 |
---|
89 | //data[ngen+3] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1 (ISO) |
---|
90 | } |
---|
91 | |
---|
92 | ngen += 4; |
---|
93 | } |
---|
94 | } |
---|
95 | } |
---|
96 | } |
---|