source: nanovis/trunk/GradientFilter.cpp @ 5841

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

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