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

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

Cleanups, no functional changes
  • Property svn:eol-style set to native
File size: 16.9 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
16#define GRADIENTS_EXT       ".grd"
17static int g_numOfSlices[3] = { 256, 256, 256 };
18static void *g_volData = 0;
19static float g_sliceDists[3];
20
21#define SOBEL               1
22#define GRAD_FILTER_SIZE    5
23#define SIGMA2              5.0
24#define MAX(a,b) ((a) > (b) ? (a) : (b))
25#define MIN(a,b) ((a) < (b) ? (a) : (b))
26#define EPS 1e-5f
27
28#ifdef notused
29static char *
30getFloatGradientsFilename(void)
31{
32    char floatExt[] = "_float";
33    char *filename;
34
35    /*
36    if (! (filename = (char *)malloc(strlen(g.basename) +
37                                     strlen(floatExt) +
38                                     strlen(GRADIENTS_EXT) + 1))) {
39                                         */
40    if (! (filename = (char *)malloc(strlen("base") +
41                                     strlen(floatExt) +
42                                     strlen(GRADIENTS_EXT) + 1))) {
43        ERROR("not enough memory for filename\n");
44        exit(1);
45    }
46
47    //strcpy(filename, g.basename);
48    strcpy(filename, "base");
49
50    strcat(filename, floatExt);
51    strcat(filename, GRADIENTS_EXT);
52
53    return filename;
54}
55
56static void
57saveFloatGradients(float *gradients, int *sizes)
58{
59    char *filename;
60    FILE *fp;
61
62    filename = getFloatGradientsFilename();
63    if (! (fp = fopen(filename, "wb"))) {
64        perror("cannot open gradients file for writing");
65        exit(1);
66    }
67    if (fwrite(gradients, 3 * sizes[0] * sizes[1] * sizes[2] * sizeof(float),
68               1, fp) != 1) {
69        ERROR("writing float gradients failed\n");
70        exit(1);
71    }
72    fclose(fp);
73}
74#endif
75
76int
77getNextPowerOfTwo(int n)
78{
79    int i;
80
81    i = 1;
82    while (i < n) {
83        i *= 2;
84    }
85
86    return i;
87}
88
89static unsigned char getVoxel8(int x, int y, int z)
90{
91    return ((unsigned char*)g_volData)[z * g_numOfSlices[0] * g_numOfSlices[1] +
92                                       y * g_numOfSlices[0] +
93                                       x];
94}
95
96static unsigned short getVoxel16(int x, int y, int z)
97{
98    return ((unsigned short*)g_volData)[z * g_numOfSlices[0] * g_numOfSlices[1] +
99                                        y * g_numOfSlices[0] +
100                                        x];
101}
102
103static float getVoxelFloat(int x, int y, int z)
104{
105    return ((float*)g_volData)[z * g_numOfSlices[0] * g_numOfSlices[1] +
106                                       y * g_numOfSlices[0] +
107                                       x];
108}
109
110static float getVoxel(int x, int y, int z, DataType dataType)
111{
112    switch (dataType) {
113        case DATRAW_UCHAR:
114            return (float)getVoxel8(x, y, z);
115            break;
116        case DATRAW_USHORT:
117            return (float)getVoxel16(x, y, z);
118            break;
119        case DATRAW_FLOAT :
120            return (float)getVoxelFloat(x, y, z);
121            break;
122        default:
123            ERROR("Unsupported data type\n");
124            exit(1);
125            break;
126    }
127    return 0.0;
128}
129
130
131void computeGradients(float *gradients, void* volData, int *sizes, DataType dataType)
132{
133    ::g_volData = volData;
134    g_numOfSlices[0] = sizes[0];
135    g_numOfSlices[1] = sizes[1];
136    g_numOfSlices[2] = sizes[2];
137
138    g_sliceDists[0] = 1.0f / sizes[0];
139    g_sliceDists[1] = 1.0f / sizes[1];
140    g_sliceDists[2] = 1.0f / sizes[2];
141
142    int i, j, k, dir, di, vdi, idz, idy, idx;
143    float *gp;
144
145    static int weights[][3][3][3] = {
146        {{{-1, -3, -1},
147          {-3, -6, -3},
148          {-1, -3, -1}},
149         {{ 0,  0,  0},
150          { 0,  0,  0},
151          { 0,  0,  0}},
152         {{ 1,  3,  1},
153          { 3,  6,  3},
154          { 1,  3,  1}}},
155        {{{-1, -3, -1},
156          { 0,  0,  0},
157          { 1,  3,  1}},
158         {{-3, -6, -3},
159          { 0,  0,  0},
160          { 3,  6,  3}},
161         {{-1, -3, -1},
162          { 0,  0,  0},
163          { 1,  3,  1}}},
164        {{{-1,  0,  1},
165          {-3,  0,  3},
166          {-1,  0,  1}},
167         {{-3,  0,  3},
168          {-6,  0,  6},
169          {-3,  0,  3}},
170         {{-1,  0,  1},
171          {-3,  0,  3},
172          {-1,  0,  1}}}
173    };
174
175    TRACE("computing gradients ... may take a while\n");
176
177    di = 0;
178    vdi = 0;
179    gp = gradients;
180    for (idz = 0; idz < sizes[2]; idz++) {
181        for (idy = 0; idy < sizes[1]; idy++) {
182            for (idx = 0; idx < sizes[0]; idx++) {
183#if SOBEL == 1
184                if (idx > 0 && idx < sizes[0] - 1 &&
185                    idy > 0 && idy < sizes[1] - 1 &&
186                    idz > 0 && idz < sizes[2] - 1) {
187
188                    for (dir = 0; dir < 3; dir++) {
189                        gp[dir] = 0.0;
190                        for (i = -1; i < 2; i++) {
191                            for (j = -1; j < 2; j++) {
192                                for (k = -1; k < 2; k++) {
193                                    gp[dir] += weights[dir][i + 1]
194                                                           [j + 1]
195                                                           [k + 1] *
196                                               getVoxel(idx + i,
197                                                      idy + j,
198                                                      idz + k,
199                                                      dataType);
200                                }
201                            }
202                        }
203
204                        gp[dir] /= 2.0 * g_sliceDists[dir];
205                    }
206                } else {
207                    /* X-direction */
208                    if (idx < 1) {
209                        gp[0] = (getVoxel(idx + 1, idy, idz, dataType) -
210                                 getVoxel(idx, idy, idz, dataType))/
211                                (g_sliceDists[0]);
212                    } else {
213                        gp[0] = (getVoxel(idx, idy, idz, dataType) -
214                                 getVoxel(idx - 1, idy, idz, dataType))/
215                                (g_sliceDists[0]);
216                    }
217
218                    /* Y-direction */
219                    if (idy < 1) {
220                        gp[1] = (getVoxel(idx, idy + 1, idz, dataType) -
221                                   getVoxel(idx, idy, idz, dataType))/
222                                   (g_sliceDists[1]);
223                    } else {
224                        gp[1] = (getVoxel(idx, idy, idz, dataType) -
225                                   getVoxel(idx, idy - 1, idz, dataType))/
226                                   (g_sliceDists[1]);
227                    }
228
229                    /* Z-direction */
230                    if (idz < 1) {
231                        gp[2] = (getVoxel(idx, idy, idz + 1, dataType) -
232                                 getVoxel(idx, idy, idz, dataType))/
233                                (g_sliceDists[2]);
234                    } else {
235                        gp[2] = (getVoxel(idx, idy, idz, dataType) -
236                                 getVoxel(idx, idy, idz - 1, dataType))/
237                                (g_sliceDists[2]);
238                    }
239                }
240#else
241                /* X-direction */
242                if (idx < 1) {
243                    gp[0] = (getVoxel(idx + 1, idy, idz, dataType) -
244                             getVoxel(idx, idy, idz, dataType))/
245                            (g_sliceDists[0]);
246                } else if (idx > g_numOfSlices[0] - 1) {
247                    gp[0] = (getVoxel(idx, idy, idz, dataType) -
248                             getVoxel(idx - 1, idy, idz, dataType))/
249                            (g_sliceDists[0]);
250                } else {
251                    gp[0] = (getVoxel(idx + 1, idy, idz, dataType) -
252                             getVoxel(idx - 1, idy, idz, dataType))/
253                            (2.0 * g_sliceDists[0]);
254                }
255
256                /* Y-direction */
257                if (idy < 1) {
258                    gp[1] = (getVoxel(idx, idy + 1, idz, dataType) -
259                             getVoxel(idx, idy, idz, dataType))/
260                            (g_sliceDists[1]);
261                } else if (idy > g_numOfSlices[1] - 1) {
262                    gp[1] = (getVoxel(idx, idy, idz, dataType) -
263                             getVoxel(idx, idy - 1, idz, dataType))/
264                            (g_sliceDists[1]);
265                } else {
266                    gp[1] = (getVoxel(idx, idy + 1, idz, dataType) -
267                             getVoxel(idx, idy - 1, idz, dataType))/
268                            (2.0 * g_sliceDists[1]);
269                }
270
271                /* Z-direction */
272                if (idz < 1) {
273                    gp[2] = (getVoxel(idx, idy, idz + 1, dataType) -
274                             getVoxel(idx, idy, idz, dataType))/
275                            (g_sliceDists[2]);
276                } else if (idz > g_numOfSlices[2] - 1) {
277                    gp[2] = (getVoxel(idx, idy, idz, dataType) -
278                             getVoxel(idx, idy, idz - 1, dataType))/
279                            (g_sliceDists[2]);
280                } else {
281                    gp[2] = (getVoxel(idx, idy, idz + 1, dataType) -
282                             getVoxel(idx, idy, idz - 1, dataType))/
283                            (2.0 * g_sliceDists[2]);
284                }
285#endif
286                gp += 3;
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.