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

Last change on this file since 3452 was 3452, checked in by ldelgass, 6 years ago

Remove XINETD define from nanovis. We only support server mode now, no glut
interaction loop with mouse/keyboard handlers. Fixes for trace logging to make
output closer to vtkvis: inlcude function name for trace messages, remove
newlines from format strings in macros since newlines get added by syslog.

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