source: nanovis/branches/1.1/GradientFilter.cpp @ 4804

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