source: trunk/packages/vizservers/nanovis/GradientFilter.cpp @ 2884

Last change on this file since 2884 was 2877, checked in by ldelgass, 12 years ago

Some minor refactoring, also add some more fine grained config.h defines
(e.g. replace NV40 define with feature defines). Add tests for some required
OpenGL extensions (should always check for extensions or base version before
calling entry points from the extension). Also, clamp diffuse and specular
values on input and warn when they are out of range.
  • Property svn:eol-style set to native
File size: 16.6 KB
Line 
1/* -*- mode: c++; c-basic-offset: 4; indent-tabs-mode: nil -*- */
2#include <stdlib.h>
3#include <float.h>
4#include <math.h>
5#include <limits.h>
6#include <string.h>
7#include <stdio.h>
8
9#include "Trace.h"
10#include "GradientFilter.h"
11
12#ifndef SQR
13#define SQR(a) ((a) * (a))
14#endif
15
16static int g_numOfSlices[3] = { 256, 256, 256 };
17static void *g_volData = NULL;
18static float g_sliceDists[3];
19
20#define SOBEL               1
21#define GRAD_FILTER_SIZE    5
22#define SIGMA2              5.0
23#define MAX(a,b) ((a) > (b) ? (a) : (b))
24#define MIN(a,b) ((a) < (b) ? (a) : (b))
25#define EPS 1e-6f
26
27#ifdef notused
28static char *
29getFloatGradientsFilename()
30{
31    char floatExt[] = "_float";
32    char *filename;
33    char extension[] = ".grd";
34
35    if (! (filename = (char *)malloc(strlen("base") +
36                                     strlen(floatExt) +
37                                     strlen(extension) + 1))) {
38        ERROR("not enough memory for filename\n");
39        exit(1);
40    }
41
42    strcpy(filename, "base");
43    strcat(filename, floatExt);
44    strcat(filename, extension);
45
46    return filename;
47}
48
49static void
50saveFloatGradients(float *gradients, int *sizes)
51{
52    char *filename;
53    FILE *fp;
54
55    filename = getFloatGradientsFilename();
56    if (! (fp = fopen(filename, "wb"))) {
57        perror("cannot open gradients file for writing");
58        exit(1);
59    }
60    if (fwrite(gradients, 3 * sizes[0] * sizes[1] * sizes[2] * sizeof(float),
61               1, fp) != 1) {
62        ERROR("writing float gradients failed\n");
63        exit(1);
64    }
65    fclose(fp);
66}
67
68static void
69saveGradients(void *gradients, int *sizes, DataType dataType)
70{
71    char *filename;
72    int size;
73    FILE *fp;
74
75    filename = getGradientsFilename();
76    if (! (fp = fopen(filename, "wb"))) {
77        perror("cannot open gradients file for writing");
78        exit(1);
79    }
80
81    size = 3 * sizes[0] * sizes[1] * sizes[2]
82           * getDataTypeSize(dataType);
83
84    if (fwrite(gradients, size, 1, fp) != 1) {
85        ERROR("writing gradients failed\n");
86        exit(1);
87    }
88
89    fclose(fp);
90}
91#endif
92
93int
94getNextPowerOfTwo(int n)
95{
96    int i;
97
98    i = 1;
99    while (i < n) {
100        i *= 2;
101    }
102
103    return i;
104}
105
106static unsigned char getVoxel8(int x, int y, int z)
107{
108    return ((unsigned char*)g_volData)[z * g_numOfSlices[0] * g_numOfSlices[1] +
109                                       y * g_numOfSlices[0] +
110                                       x];
111}
112
113static unsigned short getVoxel16(int x, int y, int z)
114{
115    return ((unsigned short*)g_volData)[z * g_numOfSlices[0] * g_numOfSlices[1] +
116                                        y * g_numOfSlices[0] +
117                                        x];
118}
119
120static float getVoxelFloat(int x, int y, int z)
121{
122    return ((float*)g_volData)[z * g_numOfSlices[0] * g_numOfSlices[1] +
123                               y * g_numOfSlices[0] +
124                               x];
125}
126
127static float getVoxel(int x, int y, int z, DataType dataType)
128{
129    switch (dataType) {
130        case DATRAW_UCHAR:
131            return (float)getVoxel8(x, y, z);
132            break;
133        case DATRAW_USHORT:
134            return (float)getVoxel16(x, y, z);
135            break;
136        case DATRAW_FLOAT :
137            return (float)getVoxelFloat(x, y, z);
138            break;
139        default:
140            ERROR("Unsupported data type\n");
141            exit(1);
142            break;
143    }
144    return 0.0;
145}
146
147void computeGradients(float *gradients, void *volData, int *sizes,
148                      float *spacing, DataType dataType)
149{
150    ::g_volData = volData;
151    g_numOfSlices[0] = sizes[0];
152    g_numOfSlices[1] = sizes[1];
153    g_numOfSlices[2] = sizes[2];
154
155    g_sliceDists[0] = spacing[0];
156    g_sliceDists[1] = spacing[1];
157    g_sliceDists[2] = spacing[2];
158
159    int i, j, k, dir, di, vdi, idz, idy, idx;
160    float *gp;
161
162    static int weights[][3][3][3] = {
163        {{{-1, -3, -1},
164          {-3, -6, -3},
165          {-1, -3, -1}},
166         {{ 0,  0,  0},
167          { 0,  0,  0},
168          { 0,  0,  0}},
169         {{ 1,  3,  1},
170          { 3,  6,  3},
171          { 1,  3,  1}}},
172        {{{-1, -3, -1},
173          { 0,  0,  0},
174          { 1,  3,  1}},
175         {{-3, -6, -3},
176          { 0,  0,  0},
177          { 3,  6,  3}},
178         {{-1, -3, -1},
179          { 0,  0,  0},
180          { 1,  3,  1}}},
181        {{{-1,  0,  1},
182          {-3,  0,  3},
183          {-1,  0,  1}},
184         {{-3,  0,  3},
185          {-6,  0,  6},
186          {-3,  0,  3}},
187         {{-1,  0,  1},
188          {-3,  0,  3},
189          {-1,  0,  1}}}
190    };
191
192    TRACE("computing gradients ... may take a while\n");
193
194    di = 0;
195    vdi = 0;
196    gp = gradients;
197    for (idz = 0; idz < sizes[2]; idz++) {
198        for (idy = 0; idy < sizes[1]; idy++) {
199            for (idx = 0; idx < sizes[0]; idx++) {
200#if SOBEL == 1
201                if (idx > 0 && idx < sizes[0] - 1 &&
202                    idy > 0 && idy < sizes[1] - 1 &&
203                    idz > 0 && idz < sizes[2] - 1) {
204
205                    for (dir = 0; dir < 3; dir++) {
206                        gp[dir] = 0.0;
207                        for (i = -1; i < 2; i++) {
208                            for (j = -1; j < 2; j++) {
209                                for (k = -1; k < 2; k++) {
210                                    gp[dir] += weights[dir][i + 1]
211                                                           [j + 1]
212                                                           [k + 1] *
213                                               getVoxel(idx + i,
214                                                      idy + j,
215                                                      idz + k,
216                                                      dataType);
217                                }
218                            }
219                        }
220
221                        gp[dir] /= 2.0 * g_sliceDists[dir];
222                    }
223                } else {
224                    /* X-direction */
225                    if (idx < 1) {
226                        gp[0] = (getVoxel(idx + 1, idy, idz, dataType) -
227                                 getVoxel(idx, idy, idz, dataType))/
228                                (g_sliceDists[0]);
229                    } else {
230                        gp[0] = (getVoxel(idx, idy, idz, dataType) -
231                                 getVoxel(idx - 1, idy, idz, dataType))/
232                                (g_sliceDists[0]);
233                    }
234
235                    /* Y-direction */
236                    if (idy < 1) {
237                        gp[1] = (getVoxel(idx, idy + 1, idz, dataType) -
238                                   getVoxel(idx, idy, idz, dataType))/
239                                   (g_sliceDists[1]);
240                    } else {
241                        gp[1] = (getVoxel(idx, idy, idz, dataType) -
242                                   getVoxel(idx, idy - 1, idz, dataType))/
243                                   (g_sliceDists[1]);
244                    }
245
246                    /* Z-direction */
247                    if (idz < 1) {
248                        gp[2] = (getVoxel(idx, idy, idz + 1, dataType) -
249                                 getVoxel(idx, idy, idz, dataType))/
250                                (g_sliceDists[2]);
251                    } else {
252                        gp[2] = (getVoxel(idx, idy, idz, dataType) -
253                                 getVoxel(idx, idy, idz - 1, dataType))/
254                                (g_sliceDists[2]);
255                    }
256                }
257#else
258                /* X-direction */
259                if (idx < 1) {
260                    gp[0] = (getVoxel(idx + 1, idy, idz, dataType) -
261                             getVoxel(idx, idy, idz, dataType))/
262                            (g_sliceDists[0]);
263                } else if (idx > g_numOfSlices[0] - 1) {
264                    gp[0] = (getVoxel(idx, idy, idz, dataType) -
265                             getVoxel(idx - 1, idy, idz, dataType))/
266                            (g_sliceDists[0]);
267                } else {
268                    gp[0] = (getVoxel(idx + 1, idy, idz, dataType) -
269                             getVoxel(idx - 1, idy, idz, dataType))/
270                            (2.0 * g_sliceDists[0]);
271                }
272
273                /* Y-direction */
274                if (idy < 1) {
275                    gp[1] = (getVoxel(idx, idy + 1, idz, dataType) -
276                             getVoxel(idx, idy, idz, dataType))/
277                            (g_sliceDists[1]);
278                } else if (idy > g_numOfSlices[1] - 1) {
279                    gp[1] = (getVoxel(idx, idy, idz, dataType) -
280                             getVoxel(idx, idy - 1, idz, dataType))/
281                            (g_sliceDists[1]);
282                } else {
283                    gp[1] = (getVoxel(idx, idy + 1, idz, dataType) -
284                             getVoxel(idx, idy - 1, idz, dataType))/
285                            (2.0 * g_sliceDists[1]);
286                }
287
288                /* Z-direction */
289                if (idz < 1) {
290                    gp[2] = (getVoxel(idx, idy, idz + 1, dataType) -
291                             getVoxel(idx, idy, idz, dataType))/
292                            (g_sliceDists[2]);
293                } else if (idz > g_numOfSlices[2] - 1) {
294                    gp[2] = (getVoxel(idx, idy, idz, dataType) -
295                             getVoxel(idx, idy, idz - 1, dataType))/
296                            (g_sliceDists[2]);
297                } else {
298                    gp[2] = (getVoxel(idx, idy, idz + 1, dataType) -
299                             getVoxel(idx, idy, idz - 1, dataType))/
300                            (2.0 * g_sliceDists[2]);
301                }
302#endif
303                gp += 3;
304            }
305        }
306    }
307}
308
309void filterGradients(float *gradients, int *sizes)
310{
311    int i, j, k, idz, idy, idx, gi, ogi, filterWidth, n, borderDist[3];
312    float sum, *filteredGradients, ****filter;
313
314    TRACE("filtering gradients ... may also take a while\n");
315
316    if (! (filteredGradients = (float *)malloc(sizes[0] * sizes[1] * sizes[2]
317                                               * 3 * sizeof(float)))) {
318        ERROR("not enough memory for filtered gradients\n");
319        exit(1);
320    }
321
322    /* Allocate storage for filter kernels */
323    if (! (filter = (float ****)malloc((GRAD_FILTER_SIZE/2 + 1) *
324                                       sizeof(float ***)))) {
325        ERROR("failed to allocate gradient filter\n");
326        exit(1);
327    }
328
329    for (i = 0; i < GRAD_FILTER_SIZE/2 + 1; i++) {
330        if (! (filter[i] = (float ***)malloc((GRAD_FILTER_SIZE) *
331                                             sizeof(float **)))) {
332            ERROR("failed to allocate gradient filter\n");
333            exit(1);
334        }
335    }
336    for (i = 0; i < GRAD_FILTER_SIZE/2 + 1; i++) {
337        for (j = 0; j < GRAD_FILTER_SIZE; j++) {
338            if (! (filter[i][j] = (float **)malloc((GRAD_FILTER_SIZE) *
339                                                   sizeof(float *)))) {
340                ERROR("failed to allocate gradient filter\n");
341                exit(1);
342            }
343        }
344    }
345    for (i = 0; i < GRAD_FILTER_SIZE/2 + 1; i++) {
346        for (j = 0; j < GRAD_FILTER_SIZE; j++) {
347            for (k = 0; k < GRAD_FILTER_SIZE; k++) {
348                if (! (filter[i][j][k] = (float *)malloc((GRAD_FILTER_SIZE) *
349                                                         sizeof(float)))) {
350                    ERROR("failed to allocate gradient filter\n");
351                    exit(1);
352                }
353            }
354        }
355    }
356
357    filterWidth = GRAD_FILTER_SIZE/2;
358
359    /* Compute the filter kernels */
360    for (n = 0; n <= filterWidth; n++) {
361        sum = 0.0f;
362        for (k = -filterWidth; k <= filterWidth; k++) {
363            for (j = -filterWidth; j <= filterWidth; j++) {
364                for (i = -filterWidth; i <= filterWidth; i++) {
365                    sum += (filter[n][filterWidth + k]
366                                     [filterWidth + j]
367                                     [filterWidth + i] =
368                            exp(-(SQR(i) + SQR(j) + SQR(k)) / SIGMA2));
369                }
370            }
371        }
372        for (k = -filterWidth; k <= filterWidth; k++) {
373            for (j = -filterWidth; j <= filterWidth; j++) {
374                for (i = -filterWidth; i <= filterWidth; i++) {
375                    filter[n][filterWidth + k]
376                             [filterWidth + j]
377                             [filterWidth + i] /= sum;
378                }
379            }
380        }
381    }
382
383    gi = 0;
384    /* Filter the gradients */
385    for (idz = 0; idz < sizes[2]; idz++) {
386        for (idy = 0; idy < sizes[1]; idy++) {
387            for (idx = 0; idx < sizes[0]; idx++) {
388                borderDist[0] = MIN(idx, sizes[0] - idx - 1);
389                borderDist[1] = MIN(idy, sizes[1] - idy - 1);
390                borderDist[2] = MIN(idz, sizes[2] - idz - 1);
391
392                filterWidth = MIN(GRAD_FILTER_SIZE/2,
393                                  MIN(MIN(borderDist[0],
394                                          borderDist[1]),
395                                          borderDist[2]));
396
397                for (n = 0; n < 3; n++) {
398                    filteredGradients[gi] = 0.0;
399                    for (k = -filterWidth; k <= filterWidth; k++) {
400                        for (j = -filterWidth; j <= filterWidth; j++) {
401                            for (i = -filterWidth; i <= filterWidth; i++) {
402                                ogi = (((idz + k) * sizes[1]  + (idy + j)) *
403                                      sizes[0] + (idx + i)) * 3 +
404                                      n;
405                                filteredGradients[gi] += filter[filterWidth]
406                                                       [filterWidth + k]
407                                                       [filterWidth + j]
408                                                       [filterWidth + i] *
409                                                       gradients[ogi];
410                            }
411                        }
412                    }
413                    gi++;
414                }
415            }
416        }
417    }
418
419    /* Replace the orignal gradients by the filtered gradients */
420    memcpy(gradients, filteredGradients,
421           sizes[0] * sizes[1] * sizes[2] * 3 * sizeof(float));
422
423    free(filteredGradients);
424
425    /* free storage of filter kernel(s) */
426    for (i = 0; i < GRAD_FILTER_SIZE/2 + 1; i++) {
427        for (j = 0; j < GRAD_FILTER_SIZE; j++) {
428            for (k = 0; k < GRAD_FILTER_SIZE; k++) {
429                free(filter[i][j][k]);
430            }
431        }
432    }
433    for (i = 0; i < GRAD_FILTER_SIZE/2 + 1; i++) {
434        for (j = 0; j < GRAD_FILTER_SIZE; j++) {
435            free(filter[i][j]);
436        }
437    }
438    for (i = 0; i < GRAD_FILTER_SIZE/2 + 1; i++) {
439        free(filter[i]);
440    }
441    free(filter);
442}
443
444void quantize8(float *grad, unsigned char *data)
445{
446    float len;
447    int i;
448
449    len = sqrt(SQR(grad[0]) + SQR(grad[1]) + SQR(grad[2]));
450
451    if (len < EPS) {
452        grad[0] = grad[1] = grad[2] = 0.0;
453    } else {
454        grad[0] /= len;
455        grad[1] /= len;
456        grad[2] /= len;
457    }
458
459    for (i = 0; i < 3; i++) {
460        data[i] = (unsigned char)((grad[i] + 1.0)/2.0 * UCHAR_MAX);
461    }
462}
463
464void quantize16(float *grad, unsigned short *data)
465{
466    float len;
467    int i;
468
469    len = sqrt(SQR(grad[0]) + SQR(grad[1]) + SQR(grad[2]));
470
471    if (len < EPS) {
472        grad[0] = grad[1] = grad[2] = 0.0;
473    } else {
474        grad[0] /= len;
475        grad[1] /= len;
476        grad[2] /= len;
477    }
478
479    for (i = 0; i < 3; i++) {
480        data[i] = (unsigned short)((grad[i] + 1.0)/2.0 * USHRT_MAX);
481    }
482}
483
484void quantizeFloat(float *grad, float *data)
485{
486    float len;
487    int i;
488
489    len = sqrt(SQR(grad[0]) + SQR(grad[1]) + SQR(grad[2]));
490
491    if (len < EPS) {
492        grad[0] = grad[1] = grad[2] = 0.0;
493    } else {
494        grad[0] /= len;
495        grad[1] /= len;
496        grad[2] /= len;
497    }
498
499    for (i = 0; i < 3; i++) {
500        data[i] = (float)((grad[i] + 1.0)/2.0);
501    }
502}
503
504void quantizeGradients(float *gradientsIn, void *gradientsOut,
505                       int *sizes, DataType dataType)
506{
507    int idx, idy, idz, di;
508
509    di = 0;
510    for (idz = 0; idz < sizes[2]; idz++) {
511        for (idy = 0; idy < sizes[1]; idy++) {
512            for (idx = 0; idx < sizes[0]; idx++) {
513                switch (dataType) {
514                case DATRAW_UCHAR:
515                    quantize8(&gradientsIn[di],
516                              &((unsigned char*)gradientsOut)[di]);
517                    break;
518                case DATRAW_USHORT:
519                    quantize16(&gradientsIn[di],
520                               &((unsigned short*)gradientsOut)[di]);
521                    break;
522                case DATRAW_FLOAT:
523                    quantizeFloat(&gradientsIn[di],
524                               &((float*)gradientsOut)[di]);
525                    break;
526                default:
527                    ERROR("unsupported data type\n");
528                    break;
529                }
530                di += 3;
531            }
532        }
533    }
534}
Note: See TracBrowser for help on using the repository browser.