source: trunk/vizservers/nanovis/GradientFilter.cpp @ 884

Last change on this file since 884 was 883, checked in by vrinside, 16 years ago

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