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

Last change on this file since 1358 was 1053, checked in by gah, 16 years ago

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