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

Last change on this file since 3502 was 3502, checked in by ldelgass, 7 years ago

Add basic VTK structured points reader to nanovis, update copyright dates.

  • Property svn:eol-style set to native
File size: 16.7 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, di, vdi, 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    di = 0;
199    vdi = 0;
200    gp = gradients;
201    for (idz = 0; idz < sizes[2]; idz++) {
202        for (idy = 0; idy < sizes[1]; idy++) {
203            for (idx = 0; idx < sizes[0]; idx++) {
204#if SOBEL == 1
205                if (idx > 0 && idx < sizes[0] - 1 &&
206                    idy > 0 && idy < sizes[1] - 1 &&
207                    idz > 0 && idz < sizes[2] - 1) {
208
209                    for (dir = 0; dir < 3; dir++) {
210                        gp[dir] = 0.0;
211                        for (i = -1; i < 2; i++) {
212                            for (j = -1; j < 2; j++) {
213                                for (k = -1; k < 2; k++) {
214                                    gp[dir] += weights[dir][i + 1]
215                                                           [j + 1]
216                                                           [k + 1] *
217                                               getVoxel(idx + i,
218                                                      idy + j,
219                                                      idz + k,
220                                                      dataType);
221                                }
222                            }
223                        }
224
225                        gp[dir] /= 2.0 * g_sliceDists[dir];
226                    }
227                } else {
228                    /* X-direction */
229                    if (idx < 1) {
230                        gp[0] = (getVoxel(idx + 1, idy, idz, dataType) -
231                                 getVoxel(idx, idy, idz, dataType))/
232                                (g_sliceDists[0]);
233                    } else {
234                        gp[0] = (getVoxel(idx, idy, idz, dataType) -
235                                 getVoxel(idx - 1, idy, idz, dataType))/
236                                (g_sliceDists[0]);
237                    }
238
239                    /* Y-direction */
240                    if (idy < 1) {
241                        gp[1] = (getVoxel(idx, idy + 1, idz, dataType) -
242                                   getVoxel(idx, idy, idz, dataType))/
243                                   (g_sliceDists[1]);
244                    } else {
245                        gp[1] = (getVoxel(idx, idy, idz, dataType) -
246                                   getVoxel(idx, idy - 1, idz, dataType))/
247                                   (g_sliceDists[1]);
248                    }
249
250                    /* Z-direction */
251                    if (idz < 1) {
252                        gp[2] = (getVoxel(idx, idy, idz + 1, dataType) -
253                                 getVoxel(idx, idy, idz, dataType))/
254                                (g_sliceDists[2]);
255                    } else {
256                        gp[2] = (getVoxel(idx, idy, idz, dataType) -
257                                 getVoxel(idx, idy, idz - 1, dataType))/
258                                (g_sliceDists[2]);
259                    }
260                }
261#else
262                /* X-direction */
263                if (idx < 1) {
264                    gp[0] = (getVoxel(idx + 1, idy, idz, dataType) -
265                             getVoxel(idx, idy, idz, dataType))/
266                            (g_sliceDists[0]);
267                } else if (idx > g_numOfSlices[0] - 1) {
268                    gp[0] = (getVoxel(idx, idy, idz, dataType) -
269                             getVoxel(idx - 1, idy, idz, dataType))/
270                            (g_sliceDists[0]);
271                } else {
272                    gp[0] = (getVoxel(idx + 1, idy, idz, dataType) -
273                             getVoxel(idx - 1, idy, idz, dataType))/
274                            (2.0 * g_sliceDists[0]);
275                }
276
277                /* Y-direction */
278                if (idy < 1) {
279                    gp[1] = (getVoxel(idx, idy + 1, idz, dataType) -
280                             getVoxel(idx, idy, idz, dataType))/
281                            (g_sliceDists[1]);
282                } else if (idy > g_numOfSlices[1] - 1) {
283                    gp[1] = (getVoxel(idx, idy, idz, dataType) -
284                             getVoxel(idx, idy - 1, idz, dataType))/
285                            (g_sliceDists[1]);
286                } else {
287                    gp[1] = (getVoxel(idx, idy + 1, idz, dataType) -
288                             getVoxel(idx, idy - 1, idz, dataType))/
289                            (2.0 * g_sliceDists[1]);
290                }
291
292                /* Z-direction */
293                if (idz < 1) {
294                    gp[2] = (getVoxel(idx, idy, idz + 1, dataType) -
295                             getVoxel(idx, idy, idz, dataType))/
296                            (g_sliceDists[2]);
297                } else if (idz > g_numOfSlices[2] - 1) {
298                    gp[2] = (getVoxel(idx, idy, idz, dataType) -
299                             getVoxel(idx, idy, idz - 1, dataType))/
300                            (g_sliceDists[2]);
301                } else {
302                    gp[2] = (getVoxel(idx, idy, idz + 1, dataType) -
303                             getVoxel(idx, idy, idz - 1, dataType))/
304                            (2.0 * g_sliceDists[2]);
305                }
306#endif
307                gp += 3;
308            }
309        }
310    }
311}
312
313void filterGradients(float *gradients, int *sizes)
314{
315    int i, j, k, idz, idy, idx, gi, ogi, filterWidth, n, borderDist[3];
316    float sum, *filteredGradients, ****filter;
317
318    TRACE("filtering gradients ... may also take a while");
319
320    if (! (filteredGradients = (float *)malloc(sizes[0] * sizes[1] * sizes[2]
321                                               * 3 * sizeof(float)))) {
322        ERROR("not enough memory for filtered gradients");
323        exit(1);
324    }
325
326    /* Allocate storage for filter kernels */
327    if (! (filter = (float ****)malloc((GRAD_FILTER_SIZE/2 + 1) * 
328                                       sizeof(float ***)))) {
329        ERROR("failed to allocate gradient filter");
330        exit(1);
331    }
332
333    for (i = 0; i < GRAD_FILTER_SIZE/2 + 1; i++) {
334        if (! (filter[i] = (float ***)malloc((GRAD_FILTER_SIZE) * 
335                                             sizeof(float **)))) {
336            ERROR("failed to allocate gradient filter");
337            exit(1);
338        }
339    }
340    for (i = 0; i < GRAD_FILTER_SIZE/2 + 1; i++) {
341        for (j = 0; j < GRAD_FILTER_SIZE; j++) {
342            if (! (filter[i][j] = (float **)malloc((GRAD_FILTER_SIZE) * 
343                                                   sizeof(float *)))) {
344                ERROR("failed to allocate gradient filter");
345                exit(1);
346            }
347        }
348    }
349    for (i = 0; i < GRAD_FILTER_SIZE/2 + 1; i++) {
350        for (j = 0; j < GRAD_FILTER_SIZE; j++) {
351            for (k = 0; k < GRAD_FILTER_SIZE; k++) {
352                if (! (filter[i][j][k] = (float *)malloc((GRAD_FILTER_SIZE) *
353                                                         sizeof(float)))) {
354                    ERROR("failed to allocate gradient filter");
355                    exit(1);
356                }
357            }
358        }
359    }
360
361    filterWidth = GRAD_FILTER_SIZE/2;
362
363    /* Compute the filter kernels */
364    for (n = 0; n <= filterWidth; n++) {
365        sum = 0.0f;
366        for (k = -filterWidth; k <= filterWidth; k++) {
367            for (j = -filterWidth; j <= filterWidth; j++) {
368                for (i = -filterWidth; i <= filterWidth; i++) {
369                    sum += (filter[n][filterWidth + k]
370                                     [filterWidth + j]
371                                     [filterWidth + i] =
372                            exp(-(SQR(i) + SQR(j) + SQR(k)) / SIGMA2));
373                }
374            }
375        }
376        for (k = -filterWidth; k <= filterWidth; k++) {
377            for (j = -filterWidth; j <= filterWidth; j++) {
378                for (i = -filterWidth; i <= filterWidth; i++) {
379                    filter[n][filterWidth + k]
380                             [filterWidth + j]
381                             [filterWidth + i] /= sum;
382                }
383            }
384        }
385    }
386
387    gi = 0;
388    /* Filter the gradients */
389    for (idz = 0; idz < sizes[2]; idz++) {
390        for (idy = 0; idy < sizes[1]; idy++) {
391            for (idx = 0; idx < sizes[0]; idx++) {
392                borderDist[0] = MIN(idx, sizes[0] - idx - 1);
393                borderDist[1] = MIN(idy, sizes[1] - idy - 1);
394                borderDist[2] = MIN(idz, sizes[2] - idz - 1);
395
396                filterWidth = MIN(GRAD_FILTER_SIZE/2,
397                                  MIN(MIN(borderDist[0], 
398                                          borderDist[1]),
399                                          borderDist[2]));
400
401                for (n = 0; n < 3; n++) {
402                    filteredGradients[gi] = 0.0;
403                    for (k = -filterWidth; k <= filterWidth; k++) {
404                        for (j = -filterWidth; j <= filterWidth; j++) {
405                            for (i = -filterWidth; i <= filterWidth; i++) {
406                                ogi = (((idz + k) * sizes[1]  + (idy + j)) *
407                                      sizes[0] + (idx + i)) * 3 + 
408                                      n;
409                                filteredGradients[gi] += filter[filterWidth]
410                                                       [filterWidth + k]
411                                                       [filterWidth + j]
412                                                       [filterWidth + i] * 
413                                                       gradients[ogi];
414                            }
415                        }
416                    }
417                    gi++;
418                }
419            }
420        }
421    }
422
423    /* Replace the orignal gradients by the filtered gradients */
424    memcpy(gradients, filteredGradients,
425           sizes[0] * sizes[1] * sizes[2] * 3 * sizeof(float));
426
427    free(filteredGradients);
428
429    /* free storage of filter kernel(s) */
430    for (i = 0; i < GRAD_FILTER_SIZE/2 + 1; i++) {
431        for (j = 0; j < GRAD_FILTER_SIZE; j++) {
432            for (k = 0; k < GRAD_FILTER_SIZE; k++) {
433                free(filter[i][j][k]);
434            }
435        }
436    }
437    for (i = 0; i < GRAD_FILTER_SIZE/2 + 1; i++) {
438        for (j = 0; j < GRAD_FILTER_SIZE; j++) {
439            free(filter[i][j]);
440        }
441    }
442    for (i = 0; i < GRAD_FILTER_SIZE/2 + 1; i++) {
443        free(filter[i]);
444    }
445    free(filter);
446}
447
448void quantize8(float *grad, unsigned char *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 char)((grad[i] + 1.0)/2.0 * UCHAR_MAX);
465    }
466}
467
468void quantize16(float *grad, unsigned short *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] = (unsigned short)((grad[i] + 1.0)/2.0 * USHRT_MAX);
485    }
486}
487
488void quantizeFloat(float *grad, float *data)
489{
490    float len;
491    int i;
492
493    len = sqrt(SQR(grad[0]) + SQR(grad[1]) + SQR(grad[2]));
494
495    if (len < EPS) {
496        grad[0] = grad[1] = grad[2] = 0.0;
497    } else {
498        grad[0] /= len;
499        grad[1] /= len;
500        grad[2] /= len;
501    }
502
503    for (i = 0; i < 3; i++) {
504        data[i] = (float)((grad[i] + 1.0)/2.0);
505    }
506}
507
508void quantizeGradients(float *gradientsIn, void *gradientsOut,
509                       int *sizes, DataType dataType)
510{
511    int idx, idy, idz, di;
512
513    di = 0;
514    for (idz = 0; idz < sizes[2]; idz++) {
515        for (idy = 0; idy < sizes[1]; idy++) {
516            for (idx = 0; idx < sizes[0]; idx++) {
517                switch (dataType) {
518                case DATRAW_UCHAR:
519                    quantize8(&gradientsIn[di], 
520                              &((unsigned char*)gradientsOut)[di]);
521                    break;
522                case DATRAW_USHORT:
523                    quantize16(&gradientsIn[di], 
524                               &((unsigned short*)gradientsOut)[di]);
525                    break;
526                case DATRAW_FLOAT:
527                    quantizeFloat(&gradientsIn[di], 
528                               &((float*)gradientsOut)[di]);
529                    break;
530                default:
531                    ERROR("unsupported data type");
532                    break;
533                }
534                di += 3;
535            }
536        }
537    }
538}
Note: See TracBrowser for help on using the repository browser.