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

Last change on this file since 1984 was 1984, checked in by gah, 13 years ago

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