Changeset 931 for trunk/vizservers


Ignore:
Timestamp:
Mar 7, 2008, 4:52:41 PM (17 years ago)
Author:
dkearney
Message:

adding new Rappture::DX class which interacts with libdx to produce positions and data as read from the dx file. moved some common functions needed by open dx and nanovis dx to the dxReaderCommon header file. this version allows the user to visualize some dx files from abinit.

Location:
trunk/vizservers/nanovis
Files:
4 added
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/vizservers/nanovis/CmdProc.cpp

    r913 r931  
    1313 *
    1414 * Results:
    15  *      If the string matches unambiguously the index of the specification in
    16  *      the array is returned.  If the string does not match, even as an
    17  *      abbreviation, any operation, -1 is returned.  If the string matches,
    18  *      but ambiguously -2 is returned.
     15 *  If the string matches unambiguously the index of the specification in
     16 *  the array is returned.  If the string does not match, even as an
     17 *  abbreviation, any operation, -1 is returned.  If the string matches,
     18 *  but ambiguously -2 is returned.
    1919 *
    2020 *-------------------------------------------------------------------------------
    … …  
    2424    Rappture::CmdSpec *specs,
    2525    int nSpecs,
    26     char *string)               /* Name of minor operation to search for */
     26    char *string)       /* Name of minor operation to search for */
    2727{
    2828    char c;
    … …  
    3535    length = strlen(string);
    3636    while (low <= high) {
    37         Rappture::CmdSpec *specPtr;
    38         int compare;
    39         int median;
    40        
    41         median = (low + high) >> 1;
    42         specPtr = specs + median;
    43 
    44         /* Test the first character */
    45         compare = c - specPtr->name[0];
    46         if (compare == 0) {
    47             /* Now test the entire string */
    48             compare = strncmp(string, specPtr->name, length);
    49             if (compare == 0) {
    50                 if ((int)length < specPtr->minChars) {
    51                     return -2;  /* Ambiguous operation name */
    52                 }
    53             }
    54         }
    55         if (compare < 0) {
    56             high = median - 1;
    57         } else if (compare > 0) {
    58             low = median + 1;
    59         } else {
    60             return median;      /* Op found. */
    61         }
    62     }
    63     return -1;                  /* Can't find operation */
     37    Rappture::CmdSpec *specPtr;
     38    int compare;
     39    int median;
     40   
     41    median = (low + high) >> 1;
     42    specPtr = specs + median;
     43
     44    /* Test the first character */
     45    compare = c - specPtr->name[0];
     46    if (compare == 0) {
     47        /* Now test the entire string */
     48        compare = strncmp(string, specPtr->name, length);
     49        if (compare == 0) {
     50        if ((int)length < specPtr->minChars) {
     51            return -2;  /* Ambiguous operation name */
     52        }
     53        }
     54    }
     55    if (compare < 0) {
     56        high = median - 1;
     57    } else if (compare > 0) {
     58        low = median + 1;
     59    } else {
     60        return median;  /* Op found. */
     61    }
     62    }
     63    return -1;          /* Can't find operation */
    6464}
    6565
    … …  
    7575 *
    7676 * Results:
    77  *      If the string matches unambiguously the index of the specification in
    78  *      the array is returned.  If the string does not match, even as an
    79  *      abbreviation, any operation, -1 is returned.  If the string matches,
    80  *      but ambiguously -2 is returned.
     77 *  If the string matches unambiguously the index of the specification in
     78 *  the array is returned.  If the string does not match, even as an
     79 *  abbreviation, any operation, -1 is returned.  If the string matches,
     80 *  but ambiguously -2 is returned.
    8181 *
    8282 *-------------------------------------------------------------------------------
    … …  
    8686    Rappture::CmdSpec *specs,
    8787    int nSpecs,
    88     char *string)               /* Name of minor operation to search for */
     88    char *string)       /* Name of minor operation to search for */
    8989{
    9090    Rappture::CmdSpec *specPtr;
    … …  
    9999    last = -1;
    100100    for (specPtr = specs, i = 0; i < nSpecs; i++, specPtr++) {
    101         if ((c == specPtr->name[0]) &&
    102             (strncmp(string, specPtr->name, length) == 0)) {
    103             last = i;
    104             nMatches++;
    105             if ((int)length == specPtr->minChars) {
    106                 break;
    107             }
    108         }
     101    if ((c == specPtr->name[0]) &&
     102        (strncmp(string, specPtr->name, length) == 0)) {
     103        last = i;
     104        nMatches++;
     105        if ((int)length == specPtr->minChars) {
     106        break;
     107        }
     108    }
    109109    }
    110110    if (nMatches > 1) {
    111         return -2;              /* Ambiguous operation name */
     111    return -2;      /* Ambiguous operation name */
    112112    }
    113113    if (nMatches == 0) {
    114         return -1;              /* Can't find operation */
     114    return -1;      /* Can't find operation */
    115115    }
    116     return last;                /* Op found. */
     116    return last;        /* Op found. */
    117117}
    118118
    … …  
    134134Tcl_ObjCmdProc *
    135135Rappture::GetOpFromObj(
    136     Tcl_Interp *interp,         /* Interpreter to report errors to */
    137     int nSpecs,                 /* Number of specifications in array */
    138     CmdSpec *specs,             /* Op specification array */
    139     int operPos,                /* Position of operation in argument list. */
    140     int objc,                   /* Number of arguments in the argument vector.
    141                                 * This includes any prefixed arguments */
    142     Tcl_Obj *CONST *objv,       /* Argument vector */
     136    Tcl_Interp *interp,     /* Interpreter to report errors to */
     137    int nSpecs,         /* Number of specifications in array */
     138    CmdSpec *specs,     /* Op specification array */
     139    int operPos,        /* Position of operation in argument list. */
     140    int objc,           /* Number of arguments in the argument vector.
     141                * This includes any prefixed arguments */
     142    Tcl_Obj *CONST *objv,   /* Argument vector */
    143143    int flags)
    144144{
    … …  
    147147    int n;
    148148
    149     if (objc <= operPos) {      /* No operation argument */
    150         Tcl_AppendResult(interp, "wrong # args: ", (char *)NULL);
     149    if (objc <= operPos) {  /* No operation argument */
     150    Tcl_AppendResult(interp, "wrong # args: ", (char *)NULL);
    151151      usage:
    152         Tcl_AppendResult(interp, "should be one of...", (char *)NULL);
    153         for (n = 0; n < nSpecs; n++) {
    154             int i;
    155 
    156             Tcl_AppendResult(interp, "\n  ", (char *)NULL);
    157             for (i = 0; i < operPos; i++) {
    158                 Tcl_AppendResult(interp, Tcl_GetString(objv[i]), " ",
    159                         (char *)NULL);
    160             }
    161             specPtr = specs + n;
    162             Tcl_AppendResult(interp, specPtr->name, " ", specPtr->usage,
    163                 (char *)NULL);
    164         }
    165         return NULL;
     152    Tcl_AppendResult(interp, "should be one of...", (char *)NULL);
     153    for (n = 0; n < nSpecs; n++) {
     154        int i;
     155
     156        Tcl_AppendResult(interp, "\n  ", (char *)NULL);
     157        for (i = 0; i < operPos; i++) {
     158        Tcl_AppendResult(interp, Tcl_GetString(objv[i]), " ",
     159            (char *)NULL);
     160        }
     161        specPtr = specs + n;
     162        Tcl_AppendResult(interp, specPtr->name, " ", specPtr->usage,
     163        (char *)NULL);
     164    }
     165    return NULL;
    166166    }
    167167    string = Tcl_GetString(objv[operPos]);
    168168    if (flags & CMDSPEC_LINEAR_SEARCH) {
    169         n = LinearOpSearch(specs, nSpecs, string);
     169    n = LinearOpSearch(specs, nSpecs, string);
    170170    } else {
    171         n = BinaryOpSearch(specs, nSpecs, string);
     171    n = BinaryOpSearch(specs, nSpecs, string);
    172172    }
    173173    if (n == -2) {
    174         char c;
    175         size_t length;
    176 
    177         Tcl_AppendResult(interp, "ambiguous", (char *)NULL);
    178         if (operPos > 2) {
    179             Tcl_AppendResult(interp, " ", Tcl_GetString(objv[operPos - 1]),
    180                 (char *)NULL);
    181         }
    182         Tcl_AppendResult(interp, " operation \"", string, "\" matches: ",
    183             (char *)NULL);
    184 
    185         c = string[0];
    186         length = strlen(string);
    187         for (n = 0; n < nSpecs; n++) {
    188             specPtr = specs + n;
    189             if ((c == specPtr->name[0]) &&
    190                 (strncmp(string, specPtr->name, length) == 0)) {
    191                 Tcl_AppendResult(interp, " ", specPtr->name, (char *)NULL);
    192             }
    193         }
    194         return NULL;
    195 
    196     } else if (n == -1) {       /* Can't find operation, display help */
    197         Tcl_AppendResult(interp, "bad", (char *)NULL);
    198         if (operPos > 2) {
    199             Tcl_AppendResult(interp, " ", Tcl_GetString(objv[operPos - 1]),
    200                 (char *)NULL);
    201         }
    202         Tcl_AppendResult(interp, " operation \"", string, "\": ", (char *)NULL);
    203         goto usage;
     174    char c;
     175    size_t length;
     176
     177    Tcl_AppendResult(interp, "ambiguous", (char *)NULL);
     178    if (operPos > 2) {
     179        Tcl_AppendResult(interp, " ", Tcl_GetString(objv[operPos - 1]),
     180        (char *)NULL);
     181    }
     182    Tcl_AppendResult(interp, " operation \"", string, "\" matches: ",
     183        (char *)NULL);
     184
     185    c = string[0];
     186    length = strlen(string);
     187    for (n = 0; n < nSpecs; n++) {
     188        specPtr = specs + n;
     189        if ((c == specPtr->name[0]) &&
     190        (strncmp(string, specPtr->name, length) == 0)) {
     191        Tcl_AppendResult(interp, " ", specPtr->name, (char *)NULL);
     192        }
     193    }
     194    return NULL;
     195
     196    } else if (n == -1) {   /* Can't find operation, display help */
     197    Tcl_AppendResult(interp, "bad", (char *)NULL);
     198    if (operPos > 2) {
     199        Tcl_AppendResult(interp, " ", Tcl_GetString(objv[operPos - 1]),
     200        (char *)NULL);
     201    }
     202    Tcl_AppendResult(interp, " operation \"", string, "\": ", (char *)NULL);
     203    goto usage;
    204204    }
    205205    specPtr = specs + n;
    206206    if ((objc < specPtr->minArgs) ||
    207         ((specPtr->maxArgs > 0) && (objc > specPtr->maxArgs))) {
    208         int i;
    209 
    210         Tcl_AppendResult(interp, "wrong # args: should be \"", (char *)NULL);
    211         for (i = 0; i < operPos; i++) {
    212             Tcl_AppendResult(interp, Tcl_GetString(objv[i]), " ",
    213                 (char *)NULL);
    214         }
    215         Tcl_AppendResult(interp, specPtr->name, " ", specPtr->usage, "\"",
    216             (char *)NULL);
    217         return NULL;
     207    ((specPtr->maxArgs > 0) && (objc > specPtr->maxArgs))) {
     208    int i;
     209
     210    Tcl_AppendResult(interp, "wrong # args: should be \"", (char *)NULL);
     211    for (i = 0; i < operPos; i++) {
     212        Tcl_AppendResult(interp, Tcl_GetString(objv[i]), " ",
     213        (char *)NULL);
     214    }
     215    Tcl_AppendResult(interp, specPtr->name, " ", specPtr->usage, "\"",
     216        (char *)NULL);
     217    return NULL;
    218218    }
    219219    return specPtr->proc;

  • trunk/vizservers/nanovis/Command.cpp

    r930 r931  
    13061306        printf("Loading DX using OpenDX library...\n");
    13071307        fflush(stdout);
    1308         //err = load_volume_stream_odx(n, buf.bytes()+5, buf.size()-5);
     1308        err = load_volume_stream_odx(n, buf.bytes()+5, buf.size()-5);
    13091309        //err = load_volume_stream2(n, fdata);
    13101310        if (err) {

  • trunk/vizservers/nanovis/GradientFilter.cpp

    r887 r931  
    1212#endif
    1313
    14 #define GRADIENTS_EXT           ".grd"
     14#define GRADIENTS_EXT       ".grd"
    1515int g_numOfSlices[3] = { 256, 256, 256 };
    1616void* g_volData = 0;
    1717float g_sliceDists[3];
    1818
    19 #define SOBEL                           1
     19#define SOBEL               1
    2020#define GRAD_FILTER_SIZE    5
    21 #define SIGMA2                          5.0
     21#define SIGMA2              5.0
    2222#define MAX(a,b) ((a) > (b) ? (a) : (b))
    2323#define MIN(a,b) ((a) < (b) ? (a) : (b))
    … …  
    2626static char *getFloatGradientsFilename(void)
    2727{
    28         char floatExt[] = "_float";
     28    char floatExt[] = "_float";
    2929    char *filename;
    3030
    31         /*
     31    /*
    3232    if (! (filename = (char *)malloc(strlen(g.basename) +
    33                                                                         strlen(floatExt) +
     33                                    strlen(floatExt) +
    3434                                     strlen(GRADIENTS_EXT) + 1))) {
    35                                                                                 */
    36         if (! (filename = (char *)malloc(strlen("base") +
    37                                                                         strlen(floatExt) +
     35                                        */
     36    if (! (filename = (char *)malloc(strlen("base") +
     37                                    strlen(floatExt) +
    3838                                     strlen(GRADIENTS_EXT) + 1))) {
    3939        fprintf(stderr, "not enough memory for filename\n");
    … …  
    4242
    4343    //strcpy(filename, g.basename);
    44         strcpy(filename, "base");
    45        
    46         strcat(filename, floatExt);
     44    strcpy(filename, "base");
     45
     46    strcat(filename, floatExt);
    4747    strcat(filename, GRADIENTS_EXT);
    4848
    … …  
    6262    }
    6363
    64     if (fwrite(gradients, 3 * sizes[0] * sizes[1] * sizes[2] * sizeof(float), 
    65                            1, fp) != 1) {
     64    if (fwrite(gradients, 3 * sizes[0] * sizes[1] * sizes[2] * sizeof(float),
     65               1, fp) != 1) {
    6666        fprintf(stderr, "writing float gradients failed\n");
    6767        exit(1);
    … …  
    7474int getNextPowerOfTwo(int n)
    7575{
    76         int i;
    77 
    78         i = 1;
    79         while (i < n) {
    80                 i *= 2;
    81         }
    82 
    83         return i;
     76    int i;
     77
     78    i = 1;
     79    while (i < n) {
     80        i *= 2;
     81    }
     82
     83    return i;
    8484}
    8585
    8686static unsigned char getVoxel8(int x, int y, int z)
    8787{
    88         return ((unsigned char*)g_volData)[z * g_numOfSlices[0] * g_numOfSlices[1] +
    89                                                                            y * g_numOfSlices[0] +
    90                                                                            x];
     88    return ((unsigned char*)g_volData)[z * g_numOfSlices[0] * g_numOfSlices[1] +
     89                                       y * g_numOfSlices[0] +
     90                                       x];
    9191}
    9292
    9393static unsigned short getVoxel16(int x, int y, int z)
    9494{
    95         return ((unsigned short*)g_volData)[z * g_numOfSlices[0] * g_numOfSlices[1] +
    96                                                                                 y * g_numOfSlices[0] +
    97                                                                                 x];
     95    return ((unsigned short*)g_volData)[z * g_numOfSlices[0] * g_numOfSlices[1] +
     96                                        y * g_numOfSlices[0] +
     97                                        x];
    9898}
    9999
    100100static unsigned char getVoxelFloat(int x, int y, int z)
    101101{
    102         return ((float*)g_volData)[z * g_numOfSlices[0] * g_numOfSlices[1] +
    103                                                                            y * g_numOfSlices[0] +
    104                                                                            x];
     102    return ((float*)g_volData)[z * g_numOfSlices[0] * g_numOfSlices[1] +
     103                                       y * g_numOfSlices[0] +
     104                                       x];
    105105}
    106106
    107107static float getVoxel(int x, int y, int z, DataType dataType)
    108108{
    109         switch (dataType) {
    110                 case DATRAW_UCHAR:
    111                         return (float)getVoxel8(x, y, z);
    112                         break;
    113                 case DATRAW_USHORT:
    114                         return (float)getVoxel16(x, y, z);
    115                         break;
    116                 case DATRAW_FLOAT :
    117                         return (float)getVoxelFloat(x, y, z);
    118                         break;
    119                 default:
    120                         fprintf(stderr, "Unsupported data type\n");
    121                         exit(1);
    122                         break;
    123         }
    124         return 0.0;
     109    switch (dataType) {
     110        case DATRAW_UCHAR:
     111            return (float)getVoxel8(x, y, z);
     112            break;
     113        case DATRAW_USHORT:
     114            return (float)getVoxel16(x, y, z);
     115            break;
     116        case DATRAW_FLOAT :
     117            return (float)getVoxelFloat(x, y, z);
     118            break;
     119        default:
     120            fprintf(stderr, "Unsupported data type\n");
     121            exit(1);
     122            break;
     123    }
     124    return 0.0;
    125125}
    126126
    … …  
    128128void computeGradients(float *gradients, void* volData, int *sizes, DataType dataType)
    129129{
    130         ::g_volData = volData;
    131         g_numOfSlices[0] = sizes[0];
    132         g_numOfSlices[1] = sizes[1];
    133         g_numOfSlices[2] = sizes[2];
    134 
    135         g_sliceDists[0] = 1.0f / sizes[0];
    136         g_sliceDists[1] = 1.0f / sizes[1];
    137         g_sliceDists[2] = 1.0f / sizes[2];
     130    ::g_volData = volData;
     131    g_numOfSlices[0] = sizes[0];
     132    g_numOfSlices[1] = sizes[1];
     133    g_numOfSlices[2] = sizes[2];
     134
     135    g_sliceDists[0] = 1.0f / sizes[0];
     136    g_sliceDists[1] = 1.0f / sizes[1];
     137    g_sliceDists[2] = 1.0f / sizes[2];
    138138
    139139    int i, j, k, dir, di, vdi, idz, idy, idx;
    … …  
    170170    };
    171171
    172         fprintf(stderr, "computing gradients ... may take a while\n");
     172    fprintf(stderr, "computing gradients ... may take a while\n");
    173173
    174174    di = 0;
    175175    vdi = 0;
    176         gp = gradients;
     176    gp = gradients;
    177177    for (idz = 0; idz < sizes[2]; idz++) {
    178178        for (idy = 0; idy < sizes[1]; idy++) {
    … …  
    194194                                                      idy + j,
    195195                                                      idz + k,
    196                                                                                                           dataType);
     196                                                      dataType);
    197197                                }
    198198                            }
    199199                        }
    200200
    201                                                 gp[dir] /= 2.0 * g_sliceDists[dir];
     201                        gp[dir] /= 2.0 * g_sliceDists[dir];
    202202                    }
    203203                } else {
    … …  
    293293    float sum, *filteredGradients, ****filter;
    294294
    295         fprintf(stderr, "filtering gradients ... may also take a while\n");
     295    fprintf(stderr, "filtering gradients ... may also take a while\n");
    296296
    297297    if (! (filteredGradients = (float *)malloc(sizes[0] * sizes[1] * sizes[2]
    298                                                                                            * 3 * sizeof(float)))) {
     298                                               * 3 * sizeof(float)))) {
    299299        fprintf(stderr, "not enough memory for filtered gradients\n");
    300300        exit(1);
    301301    }
    302302
    303         /* Allocate storage for filter kernels */
    304         if (! (filter = (float ****)malloc((GRAD_FILTER_SIZE/2 + 1) *
    305                                                                            sizeof(float ***)))) {
    306                 fprintf(stderr, "failed to allocate gradient filter\n");
    307                 exit(1);
    308         }
    309        
    310         for (i = 0; i < GRAD_FILTER_SIZE/2 + 1; i++) {
    311                 if (! (filter[i] = (float ***)malloc((GRAD_FILTER_SIZE) *
    312                                                                                         sizeof(float **)))) {
    313                         fprintf(stderr, "failed to allocate gradient filter\n");
    314                         exit(1);
    315                 }
    316         }
    317         for (i = 0; i < GRAD_FILTER_SIZE/2 + 1; i++) {
    318                 for (j = 0; j < GRAD_FILTER_SIZE; j++) {
    319                         if (! (filter[i][j] = (float **)malloc((GRAD_FILTER_SIZE) *
    320                                                                                                    sizeof(float *)))) {
    321                                 fprintf(stderr, "failed to allocate gradient filter\n");
    322                                 exit(1);
    323                         }
    324                 }
    325         }
    326         for (i = 0; i < GRAD_FILTER_SIZE/2 + 1; i++) {
    327                 for (j = 0; j < GRAD_FILTER_SIZE; j++) {
    328                         for (k = 0; k < GRAD_FILTER_SIZE; k++) {
    329                                 if (! (filter[i][j][k] = (float *)malloc((GRAD_FILTER_SIZE) *
    330                                                                                                                 sizeof(float)))) {
    331                                         fprintf(stderr, "failed to allocate gradient filter\n");
    332                                         exit(1);
    333                                 }
    334                         }
    335                 }
    336         }
     303    /* Allocate storage for filter kernels */
     304    if (! (filter = (float ****)malloc((GRAD_FILTER_SIZE/2 + 1) *
     305                                       sizeof(float ***)))) {
     306        fprintf(stderr, "failed to allocate gradient filter\n");
     307        exit(1);
     308    }
     309
     310    for (i = 0; i < GRAD_FILTER_SIZE/2 + 1; i++) {
     311        if (! (filter[i] = (float ***)malloc((GRAD_FILTER_SIZE) *
     312                                            sizeof(float **)))) {
     313            fprintf(stderr, "failed to allocate gradient filter\n");
     314            exit(1);
     315        }
     316    }
     317    for (i = 0; i < GRAD_FILTER_SIZE/2 + 1; i++) {
     318        for (j = 0; j < GRAD_FILTER_SIZE; j++) {
     319            if (! (filter[i][j] = (float **)malloc((GRAD_FILTER_SIZE) *
     320                                                   sizeof(float *)))) {
     321                fprintf(stderr, "failed to allocate gradient filter\n");
     322                exit(1);
     323            }
     324        }
     325    }
     326    for (i = 0; i < GRAD_FILTER_SIZE/2 + 1; i++) {
     327        for (j = 0; j < GRAD_FILTER_SIZE; j++) {
     328            for (k = 0; k < GRAD_FILTER_SIZE; k++) {
     329                if (! (filter[i][j][k] = (float *)malloc((GRAD_FILTER_SIZE) *
     330                                                        sizeof(float)))) {
     331                    fprintf(stderr, "failed to allocate gradient filter\n");
     332                    exit(1);
     333                }
     334            }
     335        }
     336    }
    337337
    338338    filterWidth = GRAD_FILTER_SIZE/2;
    339        
    340         /* Compute the filter kernels */
     339
     340    /* Compute the filter kernels */
    341341    for (n = 0; n <= filterWidth; n++) {
    342342        sum = 0.0f;
    … …  
    345345                for (i = -filterWidth; i <= filterWidth; i++) {
    346346                    sum += (filter[n][filterWidth + k]
    347                                                                         [filterWidth + j]
    348                                                                         [filterWidth + i] =
    349                                                         exp(-(SQR(i) + SQR(j) + SQR(k)) / SIGMA2));
     347                                    [filterWidth + j]
     348                                    [filterWidth + i] =
     349                            exp(-(SQR(i) + SQR(j) + SQR(k)) / SIGMA2));
    350350                }
    351351            }
    … …  
    355355                for (i = -filterWidth; i <= filterWidth; i++) {
    356356                    filter[n][filterWidth + k]
    357                                                         [filterWidth + j]
    358                                                         [filterWidth + i] /= sum;
     357                            [filterWidth + j]
     358                            [filterWidth + i] /= sum;
    359359                }
    360360            }
    … …  
    373373                filterWidth = MIN(GRAD_FILTER_SIZE/2,
    374374                                  MIN(MIN(borderDist[0],
    375                                                                                   borderDist[1]),
    376                                                                                   borderDist[2]));
     375                                          borderDist[1]),
     376                                          borderDist[2]));
    377377
    378378                for (n = 0; n < 3; n++) {
    … …  
    382382                            for (i = -filterWidth; i <= filterWidth; i++) {
    383383                                ogi = (((idz + k) * sizes[1]  + (idy + j)) *
    384                                                                           sizes[0] + (idx + i)) * 3 +
    385                                                                           n;
     384                                      sizes[0] + (idx + i)) * 3 +
     385                                      n;
    386386                                filteredGradients[gi] += filter[filterWidth]
    387                                                                                                            [filterWidth + k]
    388                                                                                                            [filterWidth + j]
    389                                                                                                            [filterWidth + i] *
    390                                                                                                            gradients[ogi];
     387                                                       [filterWidth + k]
     388                                                       [filterWidth + j]
     389                                                       [filterWidth + i] *
     390                                                       gradients[ogi];
    391391                            }
    392392                        }
    … …  
    398398    }
    399399
    400         /* Replace the orignal gradients by the filtered gradients */
    401         memcpy(gradients, filteredGradients,
    402                    sizes[0] * sizes[1] * sizes[2] * 3 * sizeof(float));
     400    /* Replace the orignal gradients by the filtered gradients */
     401    memcpy(gradients, filteredGradients,
     402           sizes[0] * sizes[1] * sizes[2] * 3 * sizeof(float));
    403403
    404404    free(filteredGradients);
    405405
    406         /* free storage of filter kernel(s) */
    407         for (i = 0; i < GRAD_FILTER_SIZE/2 + 1; i++) {
    408                 for (j = 0; j < GRAD_FILTER_SIZE; j++) {
    409                         for (k = 0; k < GRAD_FILTER_SIZE; k++) {
    410                                 free(filter[i][j][k]);
    411                         }
    412                 }
    413         }
    414         for (i = 0; i < GRAD_FILTER_SIZE/2 + 1; i++) {
    415                 for (j = 0; j < GRAD_FILTER_SIZE; j++) {
    416                         free(filter[i][j]);
    417                 }
    418         }
    419         for (i = 0; i < GRAD_FILTER_SIZE/2 + 1; i++) {
    420                 free(filter[i]);
    421         }
    422         free(filter);
     406    /* free storage of filter kernel(s) */
     407    for (i = 0; i < GRAD_FILTER_SIZE/2 + 1; i++) {
     408        for (j = 0; j < GRAD_FILTER_SIZE; j++) {
     409            for (k = 0; k < GRAD_FILTER_SIZE; k++) {
     410                free(filter[i][j][k]);
     411            }
     412        }
     413    }
     414    for (i = 0; i < GRAD_FILTER_SIZE/2 + 1; i++) {
     415        for (j = 0; j < GRAD_FILTER_SIZE; j++) {
     416            free(filter[i][j]);
     417        }
     418    }
     419    for (i = 0; i < GRAD_FILTER_SIZE/2 + 1; i++) {
     420        free(filter[i]);
     421    }
     422    free(filter);
    423423}
    424424
    … …  
    426426void quantize8(float *grad, unsigned char *data)
    427427{
    428         float len;
    429         int i;
    430 
    431         len = sqrt(SQR(grad[0]) + SQR(grad[1]) + SQR(grad[2]));
    432 
    433         if (len < EPS) {
    434                 grad[0] = grad[1] = grad[2] = 0.0;
    435         } else {
    436                 grad[0] /= len;
    437                 grad[1] /= len;
    438                 grad[2] /= len;
    439         }
    440 
    441         for (i = 0; i < 3; i++) {
    442                 data[i] = (unsigned char)((grad[i] + 1.0)/2.0 * UCHAR_MAX);
    443         }
     428    float len;
     429    int i;
     430
     431    len = sqrt(SQR(grad[0]) + SQR(grad[1]) + SQR(grad[2]));
     432
     433    if (len < EPS) {
     434        grad[0] = grad[1] = grad[2] = 0.0;
     435    } else {
     436        grad[0] /= len;
     437        grad[1] /= len;
     438        grad[2] /= len;
     439    }
     440
     441    for (i = 0; i < 3; i++) {
     442        data[i] = (unsigned char)((grad[i] + 1.0)/2.0 * UCHAR_MAX);
     443    }
    444444}
    445445
    446446void quantize16(float *grad, unsigned short *data)
    447447{
    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 short)((grad[i] + 1.0)/2.0 * USHRT_MAX);
    463         }
     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 short)((grad[i] + 1.0)/2.0 * USHRT_MAX);
     463    }
    464464}
    465465
    466466void quantizeFloat(float *grad, float *data)
    467467{
    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] = (float)((grad[i] + 1.0)/2.0);
    483         }
     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] = (float)((grad[i] + 1.0)/2.0);
     483    }
    484484}
    485485
    486486void quantizeGradients(float *gradientsIn, void *gradientsOut,
    487                                            int *sizes, DataType dataType)
    488 {
    489 
    490         int idx, idy, idz, di;
    491        
     487                       int *sizes, DataType dataType)
     488{
     489
     490    int idx, idy, idz, di;
     491
    492492    di = 0;
    493493    for (idz = 0; idz < sizes[2]; idz++) {
    494494        for (idy = 0; idy < sizes[1]; idy++) {
    495495            for (idx = 0; idx < sizes[0]; idx++) {
    496                                 switch (dataType) {
    497                                 case DATRAW_UCHAR:
    498                                         quantize8(&gradientsIn[di],
    499                                                           &((unsigned char*)gradientsOut)[di]);
    500                                         break;
    501                                 case DATRAW_USHORT:
    502                                         quantize16(&gradientsIn[di],
    503                                                            &((unsigned short*)gradientsOut)[di]);
    504                                         break;
    505                                 case DATRAW_FLOAT:
    506                                         quantizeFloat(&gradientsIn[di],
    507                                                            &((float*)gradientsOut)[di]);
    508                                         break;
    509                                 default:
    510                                         fprintf(stderr, "unsupported data type\n");
    511                                         break;
    512                                 }
     496                switch (dataType) {
     497                case DATRAW_UCHAR:
     498                    quantize8(&gradientsIn[di],
     499                              &((unsigned char*)gradientsOut)[di]);
     500                    break;
     501                case DATRAW_USHORT:
     502                    quantize16(&gradientsIn[di],
     503                               &((unsigned short*)gradientsOut)[di]);
     504                    break;
     505                case DATRAW_FLOAT:
     506                    quantizeFloat(&gradientsIn[di],
     507                               &((float*)gradientsOut)[di]);
     508                    break;
     509                default:
     510                    fprintf(stderr, "unsupported data type\n");
     511                    break;
     512                }
    513513                di += 3;
    514514            }
    … …  
    521521{
    522522    char *filename;
    523         int size;
     523    int size;
    524524    FILE *fp;
    525525
    … …  
    530530    }
    531531
    532         size = 3 * sizes[0] * sizes[1] * sizes[2]
    533                    * getDataTypeSize(dataType);
    534        
     532    size = 3 * sizes[0] * sizes[1] * sizes[2]
     533           * getDataTypeSize(dataType);
     534
    535535    if (fwrite(gradients, size, 1, fp) != 1) {
    536536        fprintf(stderr, "writing gradients failed\n");

  • trunk/vizservers/nanovis/Makefile.in

    r915 r931  
    114114        RenderVertexArray.o \
    115115        Renderable.o \
     116        RpDX.o \
    116117        ScreenSnapper.o \
    117118        Socket.o \
    … …  
    132133        dxReader.o \
    133134        dxReader2.o \
    134         nanovis.o
     135        dxReaderCommon.o \
     136        nanovis.o
    135137
    136138
    … …  
    179181dxReader2.o: $(srcdir)/dxReader2.cpp
    180182        $(CC) -c $(CC_SWITCHES) $(DX_INC_SPEC) $?
     183RpDX.o: $(srcdir)/RpDX.cpp
     184        $(CC) -c $(CC_SWITCHES) $(DX_INC_SPEC) $?
     185
    181186ColorGradient.o: transfer-function/ColorGradient.cpp
    182187        $(CC) $(CC_SWITCHES) -o $@ -c $<

  • trunk/vizservers/nanovis/dxReader.cpp

    r929 r931  
    1919 */
    2020
     21// common dx functions
     22#include "dxReaderCommon.h"
     23
    2124#include <stdio.h>
    2225#include <math.h>
    … …  
    3841#include "ZincBlendeVolume.h"
    3942#include "NvZincBlendeReconstructor.h"
    40 #include "GradientFilter.h"
    4143
    4244#define  _LOCAL_ZINC_TEST_ 0
    4345
    44 static float *
    45 merge(float* scalar, float* gradient, int size)
    46 {
    47     float* data = (float*) malloc(sizeof(float) * 4 * size);
    48 
    49     Vector3* g = (Vector3*) gradient;
    50    
    51     int ngen = 0, sindex = 0;
    52     for (sindex = 0; sindex <size; ++sindex) {
    53         data[ngen++] = scalar[sindex];
    54         data[ngen++] = g[sindex].x;
    55         data[ngen++] = g[sindex].y;
    56         data[ngen++] = g[sindex].z;
    57     }
    58     return data;
    59 }
    60 
    61 static void
    62 normalizeScalar(float* fdata, int count, float min, float max)
    63 {
    64     float v = max - min;
    65     if (v != 0.0f) {
    66         for (int i = 0; i < count; ++i) {
    67             fdata[i] = fdata[i] / v;
    68         }
    69     }
    70 }
    71 
    72 static float*
    73 computeGradient(float* fdata, int width, int height, int depth,
    74                 float min, float max)
    75 {
    76     float* gradients = (float *)malloc(width * height * depth * 3 *
    77                                        sizeof(float));
    78     float* tempGradients = (float *)malloc(width * height * depth * 3 *
    79                                            sizeof(float));
    80     int sizes[3] = { width, height, depth };
    81     computeGradients(tempGradients, fdata, sizes, DATRAW_FLOAT);
    82     filterGradients(tempGradients, sizes);
    83     quantizeGradients(tempGradients, gradients, sizes, DATRAW_FLOAT);
    84     normalizeScalar(fdata, width * height * depth, min, max);
    85     float* data = merge(fdata, gradients, width * height * depth);
    86     return data;
    87 }
    88 
    89 /*
     46/*
    9047 * Load a 3D vector field from a dx-format file
    9148 */
    9249void
    93 load_vector_stream(int index, std::istream& fin) 
     50load_vector_stream(int index, std::istream& fin)
    9451{
    9552    int dummy, nx, ny, nz, npts;
    … …  
    11471            } else if (sscanf(start, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {
    11572                // found one of the delta lines
    116                 if (ddx != 0.0) { 
    117                     dx = ddx; 
    118                 } else if (ddy != 0.0) { 
    119                     dy = ddy; 
    120                 } else if (ddz != 0.0) { 
    121                     dz = ddz; 
     73                if (ddx != 0.0) {
     74                    dx = ddx;
     75                } else if (ddy != 0.0) {
     76                    dy = ddy;
     77                } else if (ddz != 0.0) {
     78                    dz = ddz;
    12279                }
    12380            } else if (sscanf(start, "object %d class array type %s shape 3 rank 1 items %d data follows", &dummy, type, &npts) == 3) {
    … …  
    13592            }
    13693        }
    137     } 
     94    }
    13895
    13996    // read data points
    … …  
    234191            data[ngen] = (data[ngen]/(2.0*vmax) + 0.5);
    235192        }
    236         Volume *volPtr;
    237         volPtr = NanoVis::load_volume(index, nx, ny, nz, 3, data, vmin, vmax, 
    238                 nzero_min);
    239         volPtr->SetRange(AxisRange::X, x0, x0 + (nx * ddx));
    240         volPtr->SetRange(AxisRange::Y, y0, y0 + (ny * ddy));
    241         volPtr->SetRange(AxisRange::Z, z0, z0 + (nz * ddz));
     193        Volume *volPtr;
     194        volPtr = NanoVis::load_volume(index, nx, ny, nz, 3, data, vmin, vmax,
     195                    nzero_min);
     196        volPtr->SetRange(AxisRange::X, x0, x0 + (nx * ddx));
     197        volPtr->SetRange(AxisRange::Y, y0, y0 + (ny * ddy));
     198        volPtr->SetRange(AxisRange::Z, z0, z0 + (nz * ddz));
    242199        delete [] data;
    243200    } else {
    … …  
    250207 */
    251208Rappture::Outcome
    252 load_volume_stream2(int index, std::iostream& fin) 
     209load_volume_stream2(int index, std::iostream& fin)
    253210{
    254211    Rappture::Outcome result;
    … …  
    304261                for (int i=0; i < nxy; i++) {
    305262                    ftmp << xymesh.atNode(i).x() << " " << xymesh.atNode(i).y() << std::endl;
    306                
    307263                }
    308264                ftmp.close();
    … …  
    334290                // found origin
    335291            } else if (sscanf(start, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {
    336                 int count = 0;
     292        int count = 0;
    337293                // found one of the delta lines
    338                 if (ddx != 0.0) { 
    339                     dx = ddx;
    340                     count++;
    341                 }
    342                 if (ddy != 0.0) {
    343                     dy = ddy;
    344                     count++;
    345                 }
    346                 if (ddz != 0.0) {
    347                     dz = ddz;
    348                     count++;
    349                 }
    350                 if (count > 1) {
    351                     return result.error(
    352                         "don't know how to handle multiple non-zero delta values");
    353                 }
     294                if (ddx != 0.0) {
     295            dx = ddx;
     296            count++;
     297        }
     298        if (ddy != 0.0) {
     299            dy = ddy;
     300            count++;
     301        }
     302        if (ddz != 0.0) {
     303            dz = ddz;
     304            count++;
     305        }
     306        if (count > 1) {
     307            return result.error(
     308            "don't know how to handle multiple non-zero delta values");
     309        }
    354310            } else if (sscanf(start, "object %d class array type %s rank 0 items %d data follows", &dummy, type, &npts) == 3) {
    355311                if (isrect && (npts != nx*ny*nz)) {
    … …  
    398354
    399355                    if (dval[p] < vmin) {
    400                         vmin = dval[p];
    401                     } else if (dval[p] > vmax) {
    402                         vmax = dval[p];
    403                     }
     356            vmin = dval[p];
     357            } else if (dval[p] > vmax) {
     358            vmax = dval[p];
     359            }
    404360                    if (dval[p] != 0.0f && dval[p] < nzero_min) {
    405361                        nzero_min = dval[p];
    … …  
    431387            printf("test2\n");
    432388            fflush(stdout);
    433             if (dv == 0.0) { 
    434                 dv = 1.0;
    435             }
     389            if (dv == 0.0) {
     390        dv = 1.0;
     391        }
    436392            for (int i = 0; i < count; ++i) {
    437                 v = data[ngen];
    438                 // scale all values [0-1], -1 => out of bounds
    439                 v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
    440                 data[ngen] = v;
    441                 ngen += 4;
    442             }
     393        v = data[ngen];
     394        // scale all values [0-1], -1 => out of bounds
     395        v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
     396        data[ngen] = v;
     397        ngen += 4;
     398        }
    443399            // Compute the gradient of this data.  BE CAREFUL: center
    444400            // calculation on each node to avoid skew in either direction.
    … …  
    485441            dy = ny;
    486442            dz = nz;
    487             Volume *volPtr;
     443            Volume *volPtr;
    488444            volPtr = NanoVis::load_volume(index, nx, ny, nz, 4, data,
    489                                           vmin, vmax, nzero_min);
    490             volPtr->SetRange(AxisRange::X, x0, x0 + (nx * ddx));
    491             volPtr->SetRange(AxisRange::Y, y0, y0 + (ny * ddy));
    492             volPtr->SetRange(AxisRange::Z, z0, z0 + (nz * ddz));
     445                                          vmin, vmax, nzero_min);
     446            volPtr->SetRange(AxisRange::X, x0, x0 + (nx * ddx));
     447            volPtr->SetRange(AxisRange::Y, y0, y0 + (ny * ddy));
     448            volPtr->SetRange(AxisRange::Z, z0, z0 + (nz * ddz));
    493449            delete [] data;
    494450
    … …  
    550506            double vmin = field.valueMin();
    551507            double dv = field.valueMax() - field.valueMin();
    552             if (dv == 0.0) { 
    553                 dv = 1.0;
    554             }
     508            if (dv == 0.0) {
     509                dv = 1.0;
     510            }
    555511            // generate the uniformly sampled data that we need for a volume
    556512            int ngen = 0;
    … …  
    618574            }
    619575
    620             Volume *volPtr;
     576            Volume *volPtr;
    621577            volPtr = NanoVis::load_volume(index, nx, ny, nz, 4, data,
    622                 field.valueMin(), field.valueMax(), nzero_min);
    623             volPtr->SetRange(AxisRange::X, field.rangeMin(Rappture::xaxis),
    624                                field.rangeMax(Rappture::xaxis));
    625             volPtr->SetRange(AxisRange::Y, field.rangeMin(Rappture::yaxis),
    626                                field.rangeMax(Rappture::yaxis));
    627             volPtr->SetRange(AxisRange::Z, field.rangeMin(Rappture::zaxis),
    628                                field.rangeMax(Rappture::zaxis));
     578            field.valueMin(), field.valueMax(), nzero_min);
     579            volPtr->SetRange(AxisRange::X, field.rangeMin(Rappture::xaxis),
     580                   field.rangeMax(Rappture::xaxis));
     581            volPtr->SetRange(AxisRange::Y, field.rangeMin(Rappture::yaxis),
     582                   field.rangeMax(Rappture::yaxis));
     583            volPtr->SetRange(AxisRange::Z, field.rangeMin(Rappture::zaxis),
     584                   field.rangeMax(Rappture::zaxis));
    629585            delete [] data;
    630586        }
    … …  
    645601
    646602Rappture::Outcome
    647 load_volume_stream(int index, std::iostream& fin) 
     603load_volume_stream(int index, std::iostream& fin)
    648604{
    649605    Rappture::Outcome result;
    … …  
    698654                for (int i=0; i < nxy; i++) {
    699655                    ftmp << xymesh.atNode(i).x() << " " << xymesh.atNode(i).y() << std::endl;
    700                
     656
    701657                }
    702658                ftmp.close();
    … …  
    752708            }
    753709        }
    754     } 
     710    }
    755711
    756712    // read data points
    … …  
    777733                    int nindex = iz*nx*ny + iy*nx + ix;
    778734                    field.define(nindex, dval[p]);
     735                    fprintf(stdout,"nindex = %i\tdval[%i] = %lg\n",nindex,p,dval[p]);
     736                    fflush(stdout);
    779737                    nread++;
    780738                    if (++iz >= nz) {
    … …  
    838796            }
    839797
    840             float* data = computeGradient(cdata, nx, ny, nz, field.valueMin(), 
     798            float* data = computeGradient(cdata, nx, ny, nz, field.valueMin(),
    841799                                          field.valueMax());
    842800
     801            for (int i=0; i<nx*ny*nz; i++) {
     802                fprintf(stdout,"enddata[%i] = %lg\n",i,data[i]);
     803                fflush(stdout);
     804            }
    843805            // Compute the gradient of this data.  BE CAREFUL: center
    844806            /*
    … …  
    914876            */
    915877
    916             Volume *volPtr;
     878            Volume *volPtr;
    917879            volPtr = NanoVis::load_volume(index, nx, ny, nz, 4, data,
    918                 field.valueMin(), field.valueMax(), nzero_min);
    919             volPtr->SetRange(AxisRange::X, field.rangeMin(Rappture::xaxis),
    920                                field.rangeMax(Rappture::xaxis));
    921             volPtr->SetRange(AxisRange::Y, field.rangeMin(Rappture::yaxis),
    922                                field.rangeMax(Rappture::yaxis));
    923             volPtr->SetRange(AxisRange::Z, field.rangeMin(Rappture::zaxis),
    924                                field.rangeMax(Rappture::zaxis));
     880            field.valueMin(), field.valueMax(), nzero_min);
     881            volPtr->SetRange(AxisRange::X, field.rangeMin(Rappture::xaxis),
     882                       field.rangeMax(Rappture::xaxis));
     883            volPtr->SetRange(AxisRange::Y, field.rangeMin(Rappture::yaxis),
     884                       field.rangeMax(Rappture::yaxis));
     885            volPtr->SetRange(AxisRange::Z, field.rangeMin(Rappture::zaxis),
     886                       field.rangeMax(Rappture::zaxis));
    925887            // TBD..
    926888            // POINTSET
    … …  
    933895              NanoVis::volume[index]->pointsetIndex = NanoVis::pointSet.size() - 1;
    934896            */
    935  
     897
    936898            delete [] data;
    937899
    … …  
    10601022            }
    10611023
    1062             Volume *volPtr;
     1024            Volume *volPtr;
    10631025            volPtr = NanoVis::load_volume(index, nx, ny, nz, 4, data,
    1064                 field.valueMin(), field.valueMax(), nzero_min);
    1065             volPtr->SetRange(AxisRange::X, field.rangeMin(Rappture::xaxis),
    1066                                field.rangeMax(Rappture::xaxis));
    1067             volPtr->SetRange(AxisRange::Y, field.rangeMin(Rappture::yaxis),
    1068                                field.rangeMax(Rappture::yaxis));
    1069             volPtr->SetRange(AxisRange::Z, field.rangeMin(Rappture::zaxis),
    1070                                field.rangeMax(Rappture::zaxis));
     1026            field.valueMin(), field.valueMax(), nzero_min);
     1027            volPtr->SetRange(AxisRange::X, field.rangeMin(Rappture::xaxis),
     1028                       field.rangeMax(Rappture::xaxis));
     1029            volPtr->SetRange(AxisRange::Y, field.rangeMin(Rappture::yaxis),
     1030                       field.rangeMax(Rappture::yaxis));
     1031            volPtr->SetRange(AxisRange::Z, field.rangeMin(Rappture::zaxis),
     1032                       field.rangeMax(Rappture::zaxis));
    10711033
    10721034            // TBD..
    … …  
    10801042              NanoVis::volume[index]->pointsetIndex = NanoVis::pointSet.size() - 1;
    10811043            */
    1082  
     1044
    10831045
    10841046            delete [] data;
    … …  
    10971059    return result;
    10981060}
     1061

  • trunk/vizservers/nanovis/dxReader2.cpp

    r886 r931  
    11
    2 //opendx (libdx4-dev) headers
    3 #include <dx/dx.h>
     2#include "dxReaderCommon.h"
    43
    54#include <stdio.h>
    … …  
    1716#include "RpFieldPrism3D.h"
    1817
     18
     19// FIXME: move newCoolInterplation to the Rappture::DX object
     20// this definition is kept only so we can test newCoolInterpolation
    1921typedef struct {
    20     int n;
    21     int rank;
    22     int shape[3];
    23     Category category;
    24     Type type;
    25 } arrayInfo;
    26 
    27 int initArrayInfo(arrayInfo* a)
    28 {
    29     a->n = 0;
    30     a->rank = 0;
    31     a->shape[0] = 0;
    32     a->shape[1] = 0;
    33     a->shape[2] = 0;
    34     return 0;
    35 }
    36 
    37 int getInterpPos2(  float* numPts, float *origin,
     22    int nx;
     23    int ny;
     24    int nz;
     25    double dataMin;
     26    double dataMax;
     27    double nzero_min;
     28    float* data;
     29} dataInfo;
     30
     31// This is legacy code that can be removed during code cleanup
     32int getInterpPos(  float* numPts, float *origin,
    3833                    float* max, int rank, float* arr)
    3934{
    … …  
    4641    // dn holds the delta you want for interpolation
    4742    for (i = 0; i < rank; i++) {
    48         dn[i] = (max[i] - origin[i]) / (numPts[i] - 1);
     43        // dn[i] = (max[i] - origin[i]) / (numPts[i] - 1);
     44        dn[i] = (max[i] - origin[i]) / (numPts[i]);
    4945    }
    5046
    … …  
    6258    }
    6359
    64     /*
    65     for (z = origin[2]-(numPts[2]/2.0); z < origin[2]+(numPts[2]/2.0); z++) {
    66         for (y = origin[1]-(numPts[1]/2.0); y < origin[1]+(numPts[1]/2.0); y++) {
    67             for (x = origin[0]-(numPts[0]/2.0); x < origin[0]+(numPts[0]/2.0); x++) {
    68                 arr[i++] = x*dn[0];
    69                 arr[i++] = y*dn[1];
    70                 arr[i++] = z*dn[2];
    71                 fprintf(stderr, "(%f,%f,%f)\n",arr[i-3],arr[i-2],arr[i-1]);
    72             }
    73         }
    74     }
    75     */
    76 
    77     return 0;
    78 }
    79 
    80 
    81 int getDataStats(   int n, int rank, int shape,
    82                     float* data, float* delta, float* max)
    83 {
    84     int lcv = 0;
    85     int r = rank*shape;
    86 
    87     // initialize the max array and delta array
    88     // max holds the maximum value found for each index
    89     // delta holds the difference between each entry's value
    90     for (lcv = 0; lcv < r; lcv++) {
    91         max[lcv] = data[lcv];
    92         delta[lcv] = data[lcv];
    93     }
    94 
    95     for (lcv=lcv; lcv < n*r; lcv++) {
    96         if (data[lcv] > max[lcv%r]) {
    97             max[lcv%r] = data[lcv];
    98         }
    99         if (delta[lcv%r] == data[lcv%r]) {
    100             if (data[lcv] != data[lcv-r]) {
    101                 delta[lcv%r] = data[lcv] - data[lcv-r];
    102             }
    103         }
    104     }
    105 
    106     return lcv;
     60    return 0;
     61}
     62
     63// This is legacy code that can be removed during code cleanup
     64int getInterpData(  int rank, float* numPts, float* max, float* dxpos, Object* dxobj,
     65                    double* dataMin, double* dataMax, float** result)
     66{
     67    int pts = int(numPts[0]*numPts[1]*numPts[2]);
     68    int interppts = pts;
     69    int arrSize = interppts*rank;
     70    float interppos[arrSize];
     71    Interpolator interpolator;
     72
     73//     commented this out to see what happens when we use the dxpos array
     74//     which holds the positions dx generated for interpolation.
     75//
     76    // generate the positions we want to interpolate on
     77//    getInterpPos(numPts,dxpos,max,rank,interppos);
     78//    fprintf(stdout, "after getInterpPos\n");
     79//    fflush(stdout);
     80
     81    // build the interpolator and interpolate
     82    interpolator = DXNewInterpolator(*dxobj,INTERP_INIT_IMMEDIATE,-1.0);
     83    fprintf(stdout, "after DXNewInterpolator=%x\n", (unsigned int)interpolator);
     84    fflush(stdout);
     85//    DXInterpolate(interpolator,&interppts,interppos,*result);
     86    DXInterpolate(interpolator,&interppts,dxpos,*result);
     87    fprintf(stdout,"interppts = %i\n",interppts);
     88    fflush(stdout);
     89
     90    // print debug info
     91    for (int lcv = 0, pt = 0; lcv < pts; lcv+=3,pt+=9) {
     92        fprintf(stdout,
     93            "(%f,%f,%f)|->% 8e\n(%f,%f,%f)|->% 8e\n(%f,%f,%f)|->% 8e\n",
     94//            interppos[pt],interppos[pt+1],interppos[pt+2], (*result)[lcv],
     95//            interppos[pt+3],interppos[pt+4],interppos[pt+5],(*result)[lcv+1],
     96//            interppos[pt+6],interppos[pt+7],interppos[pt+8],(*result)[lcv+2]);
     97            dxpos[pt],dxpos[pt+1],dxpos[pt+2], (*result)[lcv],
     98            dxpos[pt+3],dxpos[pt+4],dxpos[pt+5],(*result)[lcv+1],
     99            dxpos[pt+6],dxpos[pt+7],dxpos[pt+8],(*result)[lcv+2]);
     100        fflush(stdout);
     101    }
     102
     103    for (int lcv = 0; lcv < pts; lcv=lcv+3) {
     104        if ((*result)[lcv] < *dataMin) {
     105            *dataMin = (*result)[lcv];
     106        }
     107        if ((*result)[lcv] > *dataMax) {
     108            *dataMax = (*result)[lcv];
     109        }
     110    }
     111
     112    return 0;
    107113}
    108114
    … …  
    118124 * xyz - array containing the values of the x, y, and z axis
    119125 */
    120 int idx2xyz(int *axisLen, int idx, int *xyz) {
     126int idx2xyz(const int *axisLen, int idx, int *xyz) {
    121127    int mult = 1;
    122128    int axis = 0;
    … …  
    141147 * idx - the index that represents xyz in a linear array
    142148 */
    143 int xyz2idx(int *axisLen, int *xyz, int *idx) {
     149int xyz2idx(const int *axisLen, int *xyz, int *idx) {
    144150    int mult = 1;
    145151    int axis = 0;
    … …  
    150156        *idx = *idx + xyz[axis]*mult;
    151157        mult = mult*axisLen[axis];
     158    }
     159
     160    return 0;
     161}
     162
     163Rappture::FieldRect3D* fillRectMesh(const int* numPts, const float* delta, const float* dxpos, const float* data)
     164{
     165    Rappture::Mesh1D xgrid(dxpos[0], dxpos[0]+numPts[0]*delta[0], (int)(ceil(numPts[0])));
     166    Rappture::Mesh1D ygrid(dxpos[1], dxpos[1]+numPts[1]*delta[1], (int)(ceil(numPts[1])));
     167    Rappture::Mesh1D zgrid(dxpos[2], dxpos[2]+numPts[2]*delta[2], (int)(ceil(numPts[2])));
     168    Rappture::FieldRect3D* field = new Rappture::FieldRect3D(xgrid, ygrid, zgrid);
     169
     170    int i = 0;
     171    int idx = 0;
     172    int xyz[] = {0,0,0};
     173
     174    //FIXME: This can be switched to 3 for loops to avoid the if checks
     175    for (i=0; i<numPts[0]*numPts[1]*numPts[2]; i++,xyz[2]++) {
     176        if (xyz[2] >= numPts[2]) {
     177            xyz[2] = 0;
     178            xyz[1]++;
     179            if (xyz[1] >= numPts[1]) {
     180                xyz[1] = 0;
     181                xyz[0]++;
     182                if (xyz[0] >= numPts[0]) {
     183                    xyz[0] = 0;
     184                }
     185            }
     186        }
     187        xyz2idx(numPts,xyz,&idx);
     188        field->define(idx, data[i]);
     189        fprintf(stdout,"xyz = {%i,%i,%i}\tidx = %i\tdata[%i] = % 8e\n",xyz[0],xyz[1],xyz[2],idx,i,data[i]);
     190        fflush(stdout);
     191    }
     192
     193    return field;
     194}
     195
     196float* oldCrustyInterpolation(int* nx, int* ny, int* nz, double* dx, double* dy, double* dz, const float* dxpos, Rappture::FieldRect3D* field, double* nzero_min)
     197{
     198    // figure out a good mesh spacing
     199    int nsample = 30;
     200    *dx = field->rangeMax(Rappture::xaxis) - field->rangeMin(Rappture::xaxis);
     201    *dy = field->rangeMax(Rappture::yaxis) - field->rangeMin(Rappture::yaxis);
     202    *dz = field->rangeMax(Rappture::zaxis) - field->rangeMin(Rappture::zaxis);
     203    double dmin = pow((*dx**dy**dz)/(nsample*nsample*nsample), 0.333);
     204
     205    *nx = (int)ceil(*dx/dmin);
     206    *ny = (int)ceil(*dy/dmin);
     207    *nz = (int)ceil(*dz/dmin);
     208
     209#ifndef NV40
     210    // must be an even power of 2 for older cards
     211    *nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
     212    *ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
     213    *nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
     214#endif
     215
     216    float *cdata = new float[*nx**ny**nz];
     217    int ngen = 0;
     218    *nzero_min = 0.0;
     219    for (int iz=0; iz < *nz; iz++) {
     220        double zval = dxpos[2] + iz*dmin;
     221        for (int iy=0; iy < *ny; iy++) {
     222            double yval = dxpos[1] + iy*dmin;
     223            for (int ix=0; ix < *nx; ix++)
     224            {
     225                double xval = dxpos[0] + ix*dmin;
     226                double v = field->value(xval,yval,zval);
     227
     228                if (v != 0.0f && v < *nzero_min)
     229                {
     230                    *nzero_min = v;
     231                }
     232
     233                // scale all values [0-1], -1 => out of bounds
     234                v = (isnan(v)) ? -1.0 : v;
     235
     236                cdata[ngen] = v;
     237                ++ngen;
     238            }
     239        }
     240    }
     241
     242    float* data = computeGradient(cdata, *nx, *ny, *nz, field->valueMin(), field->valueMax());
     243
     244    delete[] cdata;
     245    return data;
     246}
     247
     248// FIXME: move newCoolInterplation to the Rappture::DX object
     249int
     250newCoolInterpolation(int pts, dataInfo* inData, dataInfo* outData, float* max, float* dxpos, float* delta, float* numPts, double* dx, double* dy, double* dz)
     251{
     252    // wouldnt dataMin have the non-zero min of the data?
     253    //double nzero_min = 0.0;
     254    outData->nzero_min = inData->dataMin;
     255
     256    // need to rewrite this starting here!!! dsk
     257    // figure out a good mesh spacing
     258    // dxpos[0,1,2] are the origins for x,y,z axis
     259    *dx = max[0] - dxpos[0];
     260    *dy = max[1] - dxpos[1];
     261    *dz = max[2] - dxpos[2];
     262    double dmin = pow((*dx**dy**dz)/(pts), (double)(1.0/3.0));
     263
     264    fprintf(stdout, "dx = %f    dy = %f    dz = %f\n",*dx,*dy,*dz);
     265    fprintf(stdout, "dmin = %f\n",dmin);
     266    fprintf(stdout, "max = [%f,%f,%f]\n",max[0],max[1],max[2]);
     267    fprintf(stdout, "delta = [%f,%f,%f]\n",delta[0],delta[1],delta[2]);
     268    fprintf(stdout, "numPts = [%f,%f,%f]\n",numPts[0],numPts[1],numPts[2]);
     269    fflush(stdout);
     270
     271
     272    outData->nx = (int)ceil(numPts[0]);
     273    outData->ny = (int)ceil(numPts[1]);
     274    outData->nz = (int)ceil(numPts[2]);
     275
     276    fprintf(stdout, "nx = %i    ny = %i    nz = %i\n",outData->nx,outData->ny,outData->nz);
     277
     278    outData->data = new float[outData->nx*outData->ny*outData->nz];
     279    for (int i=0; i<(outData->nx*outData->ny*outData->nz); i+=1) {
     280        // scale all values [0,1], -1 => out of bounds
     281        outData->data[i] = (isnan(inData->data[i])) ? -1.0 :
     282            (inData->data[i] - inData->dataMin)/(inData->dataMax-inData->dataMin);
     283        fprintf(stdout, "outData[%i] = % 8e\tinData[%i] = % 8e\n",i,outData->data[i],i,inData->data[i]);
    152284    }
    153285
    … …  
    161293    Rappture::Outcome outcome;
    162294
    163     int nx, ny, nz;
    164     double dx, dy, dz;
    165295    char dxfilename[128];
    166 
    167     Object dxobj;
    168     Array dxarr;
    169     arrayInfo dataInfo;
    170     arrayInfo posInfo;
    171296
    172297    if (nBytes == 0) {
    … …  
    178303    sprintf(dxfilename, "/tmp/dx%d.dx", getpid());
    179304
    180 #ifdef notdef
    181     std::ofstream ftmp(dxfilename);
    182     fin.seekg(0,std::ios::end);
    183     int finLen = fin.tellg();
    184     fin.seekg(0,std::ios::beg);
    185     char *finBuf = new char[finLen+1];
    186     finBuf[finLen] = '\0';
    187     fprintf(stderr, "length = %d\n",finLen);
    188     fin.read(finBuf,finLen);
    189     ftmp << finBuf;
    190     ftmp.close();
    191 #else
    192305    FILE *f;
    193    
     306
    194307    f = fopen(dxfilename, "w");
    195308    fwrite(buf, sizeof(char), nBytes, f);
    196309    fclose(f);
    197 #endif
    198 
    199     fprintf(stdout, "Calling DXImportDX(%s)\n", dxfilename);
     310
     311    Rappture::DX dxObj(dxfilename);
     312
     313    if (remove(dxfilename) != 0) {
     314        fprintf(stderr, "Error deleting dx file: %s\n", dxfilename);
     315        fflush(stderr);
     316    }
     317
     318    // store our data in a rappture structure to be interpolated
     319    Rappture::FieldRect3D* field = fillRectMesh(dxObj.axisLen(),
     320                                                dxObj.delta(),
     321                                                dxObj.positions(),
     322                                                dxObj.data());
     323
     324    int nx = 0;
     325    int ny = 0;
     326    int nz = 0;
     327    double dx = 0.0;
     328    double dy = 0.0;
     329    double dz = 0.0;
     330    double nzero_min = 0.0;
     331    float* data = oldCrustyInterpolation(&nx,&ny,&nz,&dx,&dy,&dz,
     332                                        dxObj.positions(),field,&nzero_min);
     333
     334    double dataMin = field->valueMin();
     335    double dataMax = field->valueMax();
     336
     337
     338    for (int i=0; i<nx*ny*nz; i++) {
     339        fprintf(stdout,"enddata[%i] = %lg\n",i,data[i]);
     340        fflush(stdout);
     341    }
     342
     343    fprintf(stdout,"End Data Stats index = %i\n",index);
     344    fprintf(stdout,"nx = %i ny = %i nz = %i\n",nx,ny,nz);
     345    fprintf(stdout,"dx = %lg dy = %lg dz = %lg\n",dx,dy,dz);
     346    fprintf(stdout,"dataMin = %lg\tdataMax = %lg\tnzero_min = %lg\n", dataMin,dataMax,nzero_min);
    200347    fflush(stdout);
    201     dxobj = DXImportDX(dxfilename,NULL,NULL,NULL,NULL);
    202     fprintf(stdout, "dxobj=%x\n", dxobj);
    203     fflush(stdout);
    204     initArrayInfo(&dataInfo);
    205     dxarr = (Array)DXGetComponentValue((Field) dxobj, (char *)"data");
    206 
    207     DXGetArrayInfo( dxarr, &dataInfo.n, &dataInfo.type,
    208                     &dataInfo.category, &dataInfo.rank, dataInfo.shape);
    209     fprintf(stdout, "after DXGetArrayInfo\n");
    210     fflush(stdout);
    211 
    212     initArrayInfo(&posInfo);
    213     dxarr = (Array) DXGetComponentValue((Field) dxobj, (char *)"positions");
    214     DXGetArrayInfo( dxarr, &posInfo.n, &posInfo.type,
    215                     &posInfo.category, &posInfo.rank, posInfo.shape);
    216     fprintf(stdout, "after DXGetArrayInfo\n");
    217     fflush(stdout);
    218     float* dxpos = (float*) DXGetArrayData(dxarr);
    219 
    220     float delta[] = {0.0,0.0,0.0};
    221     float max[] = {0.0,0.0,0.0};
    222     getDataStats(posInfo.n,posInfo.rank,posInfo.shape[0],dxpos,delta,max);
    223     fprintf(stdout, "after getDataStats\n");
    224 
    225     int rank = 3;
    226     int nsample = 30;
    227     float numPts[] = { nsample, nsample, nsample};
    228     int pts = int(numPts[0]*numPts[1]*numPts[2]);
    229     int interppts = pts;
    230     int arrSize = interppts*rank;
    231     float interppos[arrSize];
    232     float result[interppts];
    233     Interpolator interpolator;
    234     double dataMin = 0.0;
    235     double dataMax = 0.0;
    236     double nzero_min = 0.0;
    237    
    238     getInterpPos2(numPts,dxpos,max,rank,interppos);
    239     fprintf(stdout, "after getInterpPos2\n");
    240     fflush(stdout);
    241     interpolator = DXNewInterpolator(dxobj,INTERP_INIT_IMMEDIATE,-1.0);
    242     fprintf(stdout, "after DXNewInterpolator=%x\n", interpolator);
    243     fflush(stdout);
    244     DXInterpolate(interpolator,&interppts,interppos,result);
    245     fprintf(stdout,"interppts = %i\n",interppts);
    246     fflush(stdout);
    247     for (int lcv = 0, pt = 0; lcv < pts; lcv+=3,pt+=9) {
    248         fprintf(stdout,
    249                 "(%f,%f,%f)|->% 8e\n(%f,%f,%f)|->% 8e\n(%f,%f,%f)|->% 8e\n",
    250                 interppos[pt],interppos[pt+1],interppos[pt+2], result[lcv],
    251                 interppos[pt+3],interppos[pt+4],interppos[pt+5],result[lcv+1],
    252                 interppos[pt+6],interppos[pt+7],interppos[pt+8],result[lcv+2]);
    253     fflush(stdout);
    254     }
    255    
    256     for (int lcv = 0; lcv < pts; lcv=lcv+3) {
    257         if (result[lcv] < dataMin) {
    258             dataMin = result[lcv];
    259         }
    260         if (result[lcv] > dataMax) {
    261             dataMax = result[lcv];
    262         }
    263     }
    264    
    265     // i dont understand what nzero_min is for
    266     // i assume it means non-zero min, but its the non-zero min of what?
    267     // wouldnt dataMin have the non-zero min of the data?
    268     nzero_min = dataMin;
    269    
    270     // need to rewrite this starting here!!! dsk
    271     // figure out a good mesh spacing
    272     // dxpos[0,1,2] are the origins for x,y,z axis
    273     dx = max[0] - dxpos[0];
    274     dy = max[1] - dxpos[1];
    275     dz = max[2] - dxpos[2];
    276     double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);
    277    
    278     fprintf(stdout, "dx = %f    dy = %f    dz = %f\n",dx,dy,dz);
    279     fprintf(stdout, "dmin = %f\n",dmin);
    280     fflush(stdout);
    281     nx = (int)ceil(dx/dmin);
    282     ny = (int)ceil(dy/dmin);
    283     nz = (int)ceil(dz/dmin);
    284    
    285     fprintf(stdout, "nx = %i    ny = %i    nz = %i\n",nx,ny,nz);
    286 #ifndef NV40
    287     // must be an even power of 2 for older cards
    288     nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
    289     ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
    290     nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
    291 #endif
    292    
    293     fprintf(stdout, "nx = %i    ny = %i    nz = %i\n",nx,ny,nz);
    294     float *data = new float[nx*ny*nz];
    295     for (int i=0; i<(nx*ny*nz); i+=1) {
    296         // scale all values [0,1], -1 => out of bounds
    297         data[i] = (isnan(result[i])) ? -1.0 :
    298             (result[i] - dataMin)/(dataMax-dataMin);
    299         fprintf(stdout, "data[%i] = % 8e\n",i,data[i]);
    300     }
    301     NanoVis::load_volume(index, nx, ny, nz, 1, data, dataMin, dataMax,
    302                          nzero_min);
    303 
     348    NanoVis::load_volume(index, nx, ny, nz, 4, data, dataMin, dataMax,
     349                             nzero_min);
     350
     351    // free the result array
     352    // delete[] data;
     353    // delete[] result;
     354    // delete field;
     355
     356    //
     357    // Center this new volume on the origin.
     358    //
    304359    float dx0 = -0.5;
    305360    float dy0 = -0.5*dy/dx;
    306361    float dz0 = -0.5*dz/dx;
    307362    NanoVis::volume[index]->move(Vector3(dx0, dy0, dz0));
     363
    308364    return outcome;
    309365}
    310 
Note: See TracChangeset for help on using the changeset viewer.