Changeset 911 for trunk


Ignore:
Timestamp:
Feb 26, 2008 7:27:47 PM (16 years ago)
Author:
gah
Message:

turned on ISO_TEST

Location:
trunk/vizservers/nanovis
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/vizservers/nanovis/Command.cpp

    r906 r911  
    3636 */
    3737
    38 //#define ISO_TEST
     38#define ISO_TEST                0
     39#define PLANE_CMDS              0
     40#define __TEST_CODE__           0
     41// FOR testing new functions
     42#define _LOCAL_ZINC_TEST_       0
    3943
    4044#include "Command.h"
    … …  
    7074#include <NvLIC.h>
    7175
    72 // FOR testing new functions
    73 //#define  _LOCAL_ZINC_TEST_
    74 #ifdef _LOCAL_ZINC_TEXT_
     76#if _LOCAL_ZINC_TEST_
    7577#include "Test.h"
    7678#endif
    77 #include "Test.h"
     79
    7880// EXTERN DECLARATIONS
    7981// in Nv.cpp
    … …  
    8385// in nanovis.cpp
    8486extern vector<PointSet*> g_pointSet;
    85 
    86 extern float live_rot_x;                //object rotation angles
    87 extern float live_rot_y;
    88 extern float live_rot_z;
    89 extern float live_obj_x;        //object translation location from the origin
    90 extern float live_obj_y;
    91 extern float live_obj_z;
    9287
    9388extern PlaneRenderer* plane_render;
    … …  
    139134static Tcl_CmdProc GridCmd;
    140135static Tcl_CmdProc LegendCmd;
     136#if PLANE_CMDS
    141137static Tcl_CmdProc PlaneEnableCmd;
    142138static Tcl_CmdProc PlaneLinkCmd;
    143139static Tcl_CmdProc PlaneNewCmd;
     140#endif
    144141static Tcl_CmdProc ScreenCmd;
    145142static Tcl_CmdProc ScreenShotCmd;
    … …  
    163160        vector<unsigned int>* vectorPtr);
    164161static int GetAxis(Tcl_Interp *interp, const char *string, int *valPtr);
    165 static int GetColor(Tcl_Interp *interp, const char *string, float *rgbPtr);
     162static int GetColor(Tcl_Interp *interp, int argc, const char *argv[],
     163        float *rgbPtr);
    166164static int GetDataStream(Tcl_Interp *interp, Rappture::Buffer &buf,
    167165        int nBytes);
    … …  
    670668        }
    671669
    672 #ifdef ISO_TEST
     670#if ISO_TEST
    673671/*
    674672        for (i=0; i < nslots; i++) {
    … …  
    850848            header[5] = '\0';
    851849
    852 #ifdef _LOCAL_ZINC_TEST_
     850#if _LOCAL_ZINC_TEST_
    853851            //FILE* fp = fopen("/home/nanohub/vrinside/nv/data/HOON/QDWL_100_100_50_strain_8000i.nd_zatom_12_1", "rb");
    854852            FILE* fp;
    … …  
    875873                //vol = NvZincBlendeReconstructor::getInstance()->loadFromStream(fdata);
    876874               
    877 #ifdef _LOCAL_ZINC_TEST_
     875#if _LOCAL_ZINC_TEST_
    878876                vol = NvZincBlendeReconstructor::getInstance()->loadFromMemory(b);
    879877#else
    … …  
    901899                    NanoVis::volume[n] = vol;
    902900                }
    903 #ifdef __TEST_CODE__
    904 /*
     901#if __TEST_CODE__
    905902            } else if (strcmp(header, "<FET>") == 0) {
    906903                printf("FET loading...\n");
    … …  
    914911                    return TCL_ERROR;
    915912                }
    916 */
    917913#endif  /*__TEST_CODE__*/
    918914            } else if (strcmp(header, "<ODX>") == 0) {
    … …  
    934930                std::stringstream fdata;
    935931                fdata.write(buf.bytes(),buf.size());
    936 #ifndef ISO_TEST
     932#if ISO_TEST
     933                err = load_volume_stream2(n, fdata);
     934#else
    937935                err = load_volume_stream(n, fdata);
    938 #else
    939                 err = load_volume_stream2(n, fdata);
    940936#endif
    941937                if (err) {
    … …  
    10331029            }
    10341030        } else if ((c == 'c') && (strcmp(argv[2],"color") == 0)) {
    1035             if (argc < 3) {
     1031            if (argc < 6) {
    10361032                Tcl_AppendResult(interp, "wrong # args: should be \"", argv[0],
    1037                     " outline color {R G B} ?volume ...? \"", (char*)NULL);
     1033                    " outline color R G B ?volume ...? \"", (char*)NULL);
    10381034                return TCL_ERROR;
    10391035            }
    10401036            float rgb[3];
    1041             if (GetColor(interp, argv[3], rgb) != TCL_OK) {
     1037            if (GetColor(interp, argc - 3, argv + 3, rgb) != TCL_OK) {
    10421038                return TCL_ERROR;
    10431039            }
    10441040            vector<Volume *> ivol;
    1045             if (GetVolumes(interp, argc - 4, argv + 4, &ivol) != TCL_OK) {
     1041            if (GetVolumes(interp, argc - 6, argv + 6, &ivol) != TCL_OK) {
    10461042                return TCL_ERROR;
    10471043            }
    … …  
    20932089 */
    20942090static int
    2095 GetColor(Tcl_Interp *interp, const char *string, float *rgbPtr)
    2096 {
    2097     int argc;
    2098     const char **argv;
    2099 
    2100     if (Tcl_SplitList(interp, string, &argc, &argv) != TCL_OK) {
    2101         return TCL_ERROR;
    2102     }
    2103     if (argc != 3) {
    2104         Tcl_AppendResult(interp, "bad color \"", string,
    2105             "\": should list of R G B values 0.0 - 1.0", (char*)NULL);
     2091GetColor(Tcl_Interp *interp, int argc, const char *argv[], float *rgbPtr)
     2092{
     2093    if (argc < 3) {
     2094        Tcl_AppendResult(interp, "missing color values\": ",
     2095                "should list of R G B values 0.0 - 1.0", (char*)NULL);
    21062096        return TCL_ERROR;
    21072097    }
    … …  
    21092099        (GetFloat(interp, argv[1], rgbPtr + 1) != TCL_OK) ||
    21102100        (GetFloat(interp, argv[2], rgbPtr + 2) != TCL_OK)) {
    2111         Tcl_Free((char*)argv);
    2112         return TCL_ERROR;
    2113     }
    2114     Tcl_Free((char*)argv);
    2115     return TCL_OK;
    2116 }
    2117 
    2118 
     2101        return TCL_ERROR;
     2102    }
     2103    return TCL_OK;
     2104}
     2105
     2106
     2107#if PLANE_CMDS
    21192108static int
    21202109PlaneNewCmd(ClientData cdata, Tcl_Interp *interp, int argc,
    … …  
    22082197    return TCL_OK;
    22092198}
    2210 
     2199#endif  /*PLANE_CMDS*/
    22112200
    22122201void
    … …  
    22812270        (ClientData)NULL, (Tcl_CmdDeleteProc*)NULL);
    22822271
    2283 #ifdef __TEST_CODE__
     2272#if __TEST_CODE__
    22842273    Tcl_CreateCommand(interp, "test", TestCmd,
    22852274        (ClientData)NULL, (Tcl_CmdDeleteProc*)NULL);

  • trunk/vizservers/nanovis/Trace.cpp

    r848 r911  
    11#include <Trace.h>
    22#include <stdio.h>
     3#include <stdarg.h>
    34
    4 void Trace(char* format, ...)
     5void Trace(const char* format, ...)
    56{
    67    char buff[1024];
     8    va_list lst;
    79
    8     va_list lst;
    910    va_start(lst, format);
    10     vsprintf(buff, format, lst);
    11 
     11    vsnprintf(buff, 1023, format, lst);
     12    buff[1023] = '\0';
    1213    printf(buff);
    1314    fflush(stdout);

  • trunk/vizservers/nanovis/Trace.h

    r848 r911  
    22#define __TRACE_H__
    33
    4 #include <stdarg.h>
    5 #include <stdio.h>
    6 
    7 extern void Trace(char* format, ...);
     4extern void Trace(const char* format, ...);
    85
    96#endif

  • trunk/vizservers/nanovis/Vector3.h

    r900 r911  
    6161    v.y = y + scalar;
    6262    v.z = z + scalar;
     63    return v;
    6364}
    6465
    … …  
    6970    v.y = y - scalar;
    7071    v.z = z - scalar;
     72    return v;
    7173}
    7274#endif

  • trunk/vizservers/nanovis/dxReader.cpp

    r883 r911  
    4040#include "GradientFilter.h"
    4141
    42 //#define  _LOCAL_ZINC_TEST_
    43 float* merge(float* scalar, float* gradient, int size)
     42#define  _LOCAL_ZINC_TEST_ 0
     43
     44static float *
     45merge(float* scalar, float* gradient, int size)
    4446{
    45         float* data = (float*) malloc(sizeof(float) * 4 * size);
    46 
    47         Vector3* g = (Vector3*) gradient;
    48 
    49         int ngen = 0, sindex = 0;
    50         for (sindex = 0; sindex <size; ++sindex)
    51         {
    52                 data[ngen++] = scalar[sindex];
    53                 data[ngen++] = g[sindex].x;
    54                 data[ngen++] = g[sindex].y;
    55                 data[ngen++] = g[sindex].z;
    56         }
    57         return data;
     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;
    5859}
    5960
    60 void normalizeScalar(float* fdata, int count, float min, float max)
     61static void
     62normalizeScalar(float* fdata, int count, float min, float max)
    6163{
    6264    float v = max - min;
    63     if (v != 0.0f)
    64     {
    65         for (int i = 0; i < count; ++i)
     65    if (v != 0.0f) {
     66        for (int i = 0; i < count; ++i) {
    6667            fdata[i] = fdata[i] / v;
     68        }
    6769    }
    6870}
    6971
    70 float* computeGradient(float* fdata, int width, int height, int depth, float min, float max)
     72static float*
     73computeGradient(float* fdata, int width, int height, int depth,
     74                float min, float max)
    7175{
    72                 float* gradients = (float *)malloc(width * height * depth * 3 * sizeof(float));
    73                 float* tempGradients = (float *)malloc(width * height * depth * 3 * sizeof(float));
    74                 int sizes[3] = { width, height, depth };
    75                 computeGradients(tempGradients, fdata, sizes, DATRAW_FLOAT);
    76                 filterGradients(tempGradients, sizes);
    77                 quantizeGradients(tempGradients, gradients, sizes, DATRAW_FLOAT);
    78                 normalizeScalar(fdata, width * height * depth, min, max);
    79                 float* data = merge(fdata, gradients, width * height * depth);
    80 
    81         return data;
     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;
    8287}
    8388
    … …  
    8893load_vector_stream(int index, std::istream& fin)
    8994{
    90     int dummy, nx, ny, nz, nxy, npts;
     95    int dummy, nx, ny, nz, npts;
    9196    double x0, y0, z0, dx, dy, dz, ddx, ddy, ddz;
    9297    char line[128], type[128], *start;
    9398
     99    dx = dy = dz = 0.0;         // Suppress compiler warning.
    94100    while (!fin.eof()) {
    95101        fin.getline(line, sizeof(line) - 1);
    … …  
    104110            if (sscanf(start, "object %d class gridpositions counts %d %d %d", &dummy, &nx, &ny, &nz) == 4) {
    105111                // found grid size
    106             }
    107             else if (sscanf(start, "origin %lg %lg %lg", &x0, &y0, &z0) == 3) {
     112            } else if (sscanf(start, "origin %lg %lg %lg", &x0, &y0, &z0) == 3) {
    108113                // found origin
    109             }
    110             else if (sscanf(start, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {
     114            } else if (sscanf(start, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {
    111115                // found one of the delta lines
    112                 if (ddx != 0.0) { dx = ddx; }
    113                 else if (ddy != 0.0) { dy = ddy; }
    114                 else if (ddz != 0.0) { dz = ddz; }
    115             }
    116             else if (sscanf(start, "object %d class array type %s shape 3 rank 1 items %d data follows", &dummy, type, &npts) == 3) {
     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;
     122                }
     123            } else if (sscanf(start, "object %d class array type %s shape 3 rank 1 items %d data follows", &dummy, type, &npts) == 3) {
    117124                if (npts != nx*ny*nz) {
    118125                    std::cerr << "inconsistent data: expected " << nx*ny*nz << " points but found " << npts << " points" << std::endl;
    … …  
    120127                }
    121128                break;
    122             }
    123             else if (sscanf(start, "object %d class array type %s rank 0 times %d data follows", &dummy, type, &npts) == 3) {
     129            } else if (sscanf(start, "object %d class array type %s rank 0 times %d data follows", &dummy, type, &npts) == 3) {
    124130                if (npts != nx*ny*nz) {
    125131                    std::cerr << "inconsistent data: expected " << nx*ny*nz << " points but found " << npts << " points" << std::endl;
    … …  
    224230
    225231        // scale should be accounted.
    226         for (ngen=0; ngen < npts; ngen++)
    227         {
     232        for (ngen=0; ngen < npts; ngen++) {
    228233            data[ngen] = (data[ngen]/(2.0*vmax) + 0.5);
    229234        }
    230 
    231235        NanoVis::load_volume(index, nx, ny, nz, 3, data, vmin, vmax, nzero_min);
    232 
    233236        delete [] data;
    234     }
    235     else {
     237    } else {
    236238        std::cerr << "WARNING: data not found in stream" << std::endl;
    237239    }
    … …  
    252254
    253255    int isrect = 1;
    254 
    255     float* voldata = 0;
     256    dx = dy = dz = 0.0;         // Suppress compiler warning.
    256257    do {
    257258        fin.getline(line,sizeof(line)-1);
    … …  
    279280                result.error(mesg);
    280281                return result;
    281 
    282                 char fpts[128];
    283                 sprintf(fpts, "/tmp/tmppts%d", getpid());
    284                 char fcells[128];
    285                 sprintf(fcells, "/tmp/tmpcells%d", getpid());
    286 
    287                 std::ofstream ftmp(fpts);
    288                 // save corners of bounding box first, to work around meshing
    289                 // problems in voronoi utility
    290                 ftmp << xymesh.rangeMin(Rappture::xaxis) << " "
    291                      << xymesh.rangeMin(Rappture::yaxis) << std::endl;
    292                 ftmp << xymesh.rangeMax(Rappture::xaxis) << " "
    293                      << xymesh.rangeMin(Rappture::yaxis) << std::endl;
    294                 ftmp << xymesh.rangeMax(Rappture::xaxis) << " "
    295                      << xymesh.rangeMax(Rappture::yaxis) << std::endl;
    296                 ftmp << xymesh.rangeMin(Rappture::xaxis) << " "
    297                      << xymesh.rangeMax(Rappture::yaxis) << std::endl;
    298                 for (int i=0; i < nxy; i++) {
    299                     ftmp << xymesh.atNode(i).x() << " " << xymesh.atNode(i).y() << std::endl;
    300                
    301                 }
    302                 ftmp.close();
    303 
    304                 char cmdstr[512];
    305                 sprintf(cmdstr, "voronoi -t < %s > %s", fpts, fcells);
    306                 if (system(cmdstr) == 0) {
    307                     int cx, cy, cz;
    308                     std::ifstream ftri(fcells);
    309                     while (!ftri.eof()) {
    310                         ftri.getline(line,sizeof(line)-1);
    311                         if (sscanf(line, "%d %d %d", &cx, &cy, &cz) == 3) {
    312                             if (cx >= 4 && cy >= 4 && cz >= 4) {
    313                                 // skip first 4 boundary points
    314                                 xymesh.addCell(cx-4, cy-4, cz-4);
    315                             }
    316                         }
    317                     }
    318                     ftri.close();
    319                 } else {
    320                     return result.error("triangularization failed");
    321                 }
    322 
    323                 sprintf(cmdstr, "rm -f %s %s", fpts, fcells);
    324                 system(cmdstr);
    325             }
    326             else if (sscanf(start, "object %d class regulararray count %d", &dummy, &nz) == 2) {
    327                 // found z-grid
    328             }
    329             else if (sscanf(start, "origin %lg %lg %lg", &x0, &y0, &z0) == 3) {
    330                 // found origin
    331             }
    332             else if (sscanf(start, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {
    333                 // found one of the delta lines
    334                 if (ddx != 0.0) { dx = ddx; }
    335                 else if (ddy != 0.0) { dy = ddy; }
    336                 else if (ddz != 0.0) { dz = ddz; }
    337             }
    338             else if (sscanf(start, "object %d class array type %s rank 0 items %d data follows", &dummy, type, &npts) == 3) {
    339                 if (isrect && (npts != nx*ny*nz)) {
    340                     char mesg[256];
    341                     sprintf(mesg,"inconsistent data: expected %d points but found %d points", nx*ny*nz, npts);
    342                     return result.error(mesg);
    343                 }
    344                 else if (!isrect && (npts != nxy*nz)) {
    345                     char mesg[256];
    346                     sprintf(mesg,"inconsistent data: expected %d points but found %d points", nxy*nz, npts);
    347                     return result.error(mesg);
    348                 }
    349                 break;
    350             }
    351             else if (sscanf(start, "object %d class array type %s rank 0 times %d data follows", &dummy, type, &npts) == 3) {
    352                 if (npts != nx*ny*nz) {
    353                     char mesg[256];
    354                     sprintf(mesg,"inconsistent data: expected %d points but found %d points", nx*ny*nz, npts);
    355                     return result.error(mesg);
    356                 }
    357                 break;
    358             }
    359         }
    360     } while (!fin.eof());
    361 
    362     // read data points
    363     if (!fin.eof()) {
    364         if (isrect) {
    365             double dval[6];
    366             int nread = 0;
    367             int ix = 0;
    368             int iy = 0;
    369             int iz = 0;
    370             float* data = new float[nx *  ny *  nz * 4];
    371             memset(data, 0, nx*ny*nz*4);
    372             double vmin = 1e21;
    373             double nzero_min = 1e21;
    374             double vmax = -1e21;
    375 
    376 
    377             while (!fin.eof() && nread < npts) {
    378                 fin.getline(line,sizeof(line)-1);
    379                 int n = sscanf(line, "%lg %lg %lg %lg %lg %lg", &dval[0], &dval[1], &dval[2], &dval[3], &dval[4], &dval[5]);
    380 
    381                 for (int p=0; p < n; p++) {
    382                     int nindex = (iz*nx*ny + iy*nx + ix) * 4;
    383                     data[nindex] = dval[p];
    384 
    385                     if (dval[p] < vmin) vmin = dval[p];
    386                     if (dval[p] > vmax) vmax = dval[p];
    387                     if (dval[p] != 0.0f && dval[p] < nzero_min)
    388                     {
    389                          nzero_min = dval[p];
    390                     }
    391 
    392                     nread++;
    393                     if (++iz >= nz) {
    394                         iz = 0;
    395                         if (++iy >= ny) {
    396                             iy = 0;
    397                             ++ix;
    398                         }
    399                     }
    400                 }
    401             }
    402 
    403             // make sure that we read all of the expected points
    404             if (nread != nx*ny*nz) {
    405                 char mesg[256];
    406                 sprintf(mesg,"inconsistent data: expected %d points but found %d points", nx*ny*nz, nread);
    407                 result.error(mesg);
    408                 return result;
    409             }
    410 
    411             double dv = vmax - vmin;
    412             int count = nx*ny*nz;
    413             int ngen = 0;
    414             double v;
    415             printf("test2\n");
    416                         fflush(stdout);
    417             if (dv == 0.0) { dv = 1.0; }
    418             for (int i = 0; i < count; ++i)
    419             {
    420                 v = data[ngen];
    421                 // scale all values [0-1], -1 => out of bounds
    422                 v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
    423                 data[ngen] = v;
    424                 ngen += 4;
    425             }
    426 
    427             // Compute the gradient of this data.  BE CAREFUL: center
    428             // calculation on each node to avoid skew in either direction.
    429             ngen = 0;
    430             for (int iz=0; iz < nz; iz++) {
    431                 for (int iy=0; iy < ny; iy++) {
    432                     for (int ix=0; ix < nx; ix++) {
    433                         // gradient in x-direction
    434                         double valm1 = (ix == 0) ? 0.0 : data[ngen - 4];
    435                         double valp1 = (ix == nx-1) ? 0.0 : data[ngen + 4];
    436                         if (valm1 < 0 || valp1 < 0) {
    437                             data[ngen+1] = 0.0;
    438                         } else {
    439                             data[ngen+1] = valp1-valm1; // assume dx=1
    440                             //data[ngen+1] = ((valp1-valm1) + 1) *  0.5; // assume dx=1
    441                         }
    442 
    443                         // gradient in y-direction
    444                         valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
    445                         valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
    446                         if (valm1 < 0 || valp1 < 0) {
    447                             data[ngen+2] = 0.0;
    448                         } else {
    449                             data[ngen+2] = valp1-valm1; // assume dx=1
    450                             //data[ngen+2] = ((valp1-valm1) + 1) *  0.5; // assume dy=1
    451                         }
    452 
    453                         // gradient in z-direction
    454                         valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
    455                         valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
    456                         if (valm1 < 0 || valp1 < 0) {
    457                             data[ngen+3] = 0.0;
    458                         } else {
    459                             data[ngen+3] = valp1-valm1; // assume dx=1
    460                             //data[ngen+3] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1
    461                         }
    462 
    463                         ngen += 4;
    464                     }
    465                 }
    466             }
    467 
    468             dx = nx;
    469             dy = ny;
    470             dz = nz;
    471             NanoVis::load_volume(index, nx, ny, nz, 4, data,
    472                 vmin, vmax, nzero_min);
    473 
    474             delete [] data;
    475 
    476         } else {
    477             Rappture::Mesh1D zgrid(z0, z0+nz*dz, nz);
    478             Rappture::FieldPrism3D field(xymesh, zgrid);
    479 
    480             double dval;
    481             int nread = 0;
    482             int ixy = 0;
    483             int iz = 0;
    484             while (!fin.eof() && nread < npts) {
    485                 fin >> dval;
    486                 if (fin.fail()) {
    487                     char mesg[256];
    488                     sprintf(mesg,"after %d of %d points: can't read number",
    489                             nread, npts);
    490                     return result.error(mesg);
    491                 } else {
    492                     int nid = nxy*iz + ixy;
    493                     field.define(nid, dval);
    494 
    495                     nread++;
    496                     if (++iz >= nz) {
    497                         iz = 0;
    498                         ixy++;
    499                     }
    500                 }
    501             }
    502 
    503             // make sure that we read all of the expected points
    504             if (nread != nxy*nz) {
    505                 char mesg[256];
    506                 sprintf(mesg,"inconsistent data: expected %d points but found %d points", nxy*nz, nread);
    507                 return result.error(mesg);
    508             }
    509 
    510             // figure out a good mesh spacing
    511             int nsample = 30;
    512             x0 = field.rangeMin(Rappture::xaxis);
    513             dx = field.rangeMax(Rappture::xaxis) - field.rangeMin(Rappture::xaxis);
    514             y0 = field.rangeMin(Rappture::yaxis);
    515             dy = field.rangeMax(Rappture::yaxis) - field.rangeMin(Rappture::yaxis);
    516             z0 = field.rangeMin(Rappture::zaxis);
    517             dz = field.rangeMax(Rappture::zaxis) - field.rangeMin(Rappture::zaxis);
    518             double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);
    519 
    520             nx = (int)ceil(dx/dmin);
    521             ny = (int)ceil(dy/dmin);
    522             nz = (int)ceil(dz/dmin);
    523 #ifndef NV40
    524             // must be an even power of 2 for older cards
    525                 nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
    526                 ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
    527                 nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
    528 #endif
    529             float *data = new float[4*nx*ny*nz];
    530 
    531             double vmin = field.valueMin();
    532             double dv = field.valueMax() - field.valueMin();
    533             if (dv == 0.0) { dv = 1.0; }
    534 
    535             // generate the uniformly sampled data that we need for a volume
    536             int ngen = 0;
    537             double nzero_min = 0.0;
    538             for (iz=0; iz < nz; iz++) {
    539                 double zval = z0 + iz*dmin;
    540                 for (int iy=0; iy < ny; iy++) {
    541                     double yval = y0 + iy*dmin;
    542                     for (int ix=0; ix < nx; ix++) {
    543                         double xval = x0 + ix*dmin;
    544                         double v = field.value(xval,yval,zval);
    545 
    546                         if (v != 0.0f && v < nzero_min)
    547                         {
    548                             nzero_min = v;
    549                         }
    550                         // scale all values [0-1], -1 => out of bounds
    551                         v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
    552                         data[ngen] = v;
    553 
    554                         ngen += 4;
    555                     }
    556                 }
    557             }
    558 
    559             // Compute the gradient of this data.  BE CAREFUL: center
    560             // calculation on each node to avoid skew in either direction.
    561             ngen = 0;
    562             for (int iz=0; iz < nz; iz++) {
    563                 for (int iy=0; iy < ny; iy++) {
    564                     for (int ix=0; ix < nx; ix++) {
    565                         // gradient in x-direction
    566                         double valm1 = (ix == 0) ? 0.0 : data[ngen-4];
    567                         double valp1 = (ix == nx-1) ? 0.0 : data[ngen+4];
    568                         if (valm1 < 0 || valp1 < 0) {
    569                             data[ngen+1] = 0.0;
    570                         } else {
    571                             //data[ngen+1] = valp1-valm1; // assume dx=1
    572                             data[ngen+1] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1
    573                         }
    574 
    575                         // gradient in y-direction
    576                         valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
    577                         valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
    578                         if (valm1 < 0 || valp1 < 0) {
    579                             data[ngen+2] = 0.0;
    580                         } else {
    581                             //data[ngen+2] = valp1-valm1; // assume dy=1
    582                             data[ngen+2] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1
    583                         }
    584 
    585                         // gradient in z-direction
    586                         valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
    587                         valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
    588                         if (valm1 < 0 || valp1 < 0) {
    589                             data[ngen+3] = 0.0;
    590                         } else {
    591                             //data[ngen+3] = valp1-valm1; // assume dz=1
    592                             data[ngen+3] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1
    593                         }
    594 
    595                         ngen += 4;
    596                     }
    597                 }
    598             }
    599 
    600             NanoVis::load_volume(index, nx, ny, nz, 4, data,
    601                 field.valueMin(), field.valueMax(), nzero_min);
    602 
    603             delete [] data;
    604         }
    605     } else {
    606         return result.error("data not found in stream");
    607     }
    608 
    609     //
    610     // Center this new volume on the origin.
    611     //
    612     float dx0 = -0.5;
    613     float dy0 = -0.5*dy/dx;
    614     float dz0 = -0.5*dz/dx;
    615     NanoVis::volume[index]->move(Vector3(dx0, dy0, dz0));
    616 
    617     return result;
    618 }
    619 
    620 Rappture::Outcome
    621 load_volume_stream(int index, std::iostream& fin)
    622 {
    623     Rappture::Outcome result;
    624 
    625     Rappture::MeshTri2D xymesh;
    626     int dummy, nx, ny, nz, nxy, npts;
    627     double x0, y0, z0, dx, dy, dz, ddx, ddy, ddz;
    628     char line[128], type[128], *start;
    629 
    630     int isrect = 1;
    631 
    632     while (!fin.eof()) {
    633         fin.getline(line, sizeof(line) - 1);
    634         if (fin.fail()) {
    635             return result.error("error in data stream");
    636         }
    637         for (start=line; *start == ' ' || *start == '\t'; start++)
    638             ;  // skip leading blanks
    639 
    640         if (*start != '#') {  // skip comment lines
    641             if (sscanf(start, "object %d class gridpositions counts %d %d %d", &dummy, &nx, &ny, &nz) == 4) {
    642                 // found grid size
    643                 isrect = 1;
    644             } else if (sscanf(start, "object %d class array type float rank 1 shape 3 items %d data follows", &dummy, &nxy) == 2) {
    645                 isrect = 0;
    646 
    647                 double xx, yy, zz;
    648                 for (int i=0; i < nxy; i++) {
    649                     fin.getline(line,sizeof(line)-1);
    650                     if (sscanf(line, "%lg %lg %lg", &xx, &yy, &zz) == 3) {
    651                         xymesh.addNode( Rappture::Node2D(xx,yy) );
    652                     }
    653                 }
    654282
    655283                char fpts[128];
    … …  
    725353            }
    726354        }
    727     }
     355    } while (!fin.eof());
    728356
    729357    // read data points
    730358    if (!fin.eof()) {
    731359        if (isrect) {
    732             Rappture::Mesh1D xgrid(x0, x0+nx*dx, nx);
    733             Rappture::Mesh1D ygrid(y0, y0+ny*dy, ny);
    734             Rappture::Mesh1D zgrid(z0, z0+nz*dz, nz);
    735             Rappture::FieldRect3D field(xgrid, ygrid, zgrid);
    736 
    737360            double dval[6];
    738361            int nread = 0;
    … …  
    740363            int iy = 0;
    741364            int iz = 0;
     365            float* data = new float[nx *  ny *  nz * 4];
     366            memset(data, 0, nx*ny*nz*4);
     367            double vmin = 1e21;
     368            double nzero_min = 1e21;
     369            double vmax = -1e21;
     370
     371
    742372            while (!fin.eof() && nread < npts) {
    743373                fin.getline(line,sizeof(line)-1);
    744                 if (fin.fail()) {
    745                     return result.error("error reading data points");
    746                 }
    747374                int n = sscanf(line, "%lg %lg %lg %lg %lg %lg", &dval[0], &dval[1], &dval[2], &dval[3], &dval[4], &dval[5]);
    748375
    749376                for (int p=0; p < n; p++) {
    750                     int nindex = iz*nx*ny + iy*nx + ix;
    751                     field.define(nindex, dval[p]);
     377                    int nindex = (iz*nx*ny + iy*nx + ix) * 4;
     378                    data[nindex] = dval[p];
     379
     380                    if (dval[p] < vmin) vmin = dval[p];
     381                    if (dval[p] > vmax) vmax = dval[p];
     382                    if (dval[p] != 0.0f && dval[p] < nzero_min) {
     383                         nzero_min = dval[p];
     384                    }
     385
    752386                    nread++;
    753387                    if (++iz >= nz) {
    … …  
    769403            }
    770404
    771             // figure out a good mesh spacing
    772             int nsample = 30;
    773             dx = field.rangeMax(Rappture::xaxis) - field.rangeMin(Rappture::xaxis);
    774             dy = field.rangeMax(Rappture::yaxis) - field.rangeMin(Rappture::yaxis);
    775             dz = field.rangeMax(Rappture::zaxis) - field.rangeMin(Rappture::zaxis);
    776             double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);
    777 
    778             nx = (int)ceil(dx/dmin);
    779             ny = (int)ceil(dy/dmin);
    780             nz = (int)ceil(dz/dmin);
    781 
    782 #ifndef NV40
    783             // must be an even power of 2 for older cards
    784                 nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
    785                 ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
    786                 nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
    787 #endif
    788 
    789             float *cdata = new float[nx*ny*nz];
     405            double dv = vmax - vmin;
     406            int count = nx*ny*nz;
    790407            int ngen = 0;
    791             double nzero_min = 0.0;
    792             for (int iz=0; iz < nz; iz++) {
    793                 double zval = z0 + iz*dmin;
    794                 for (int iy=0; iy < ny; iy++) {
    795                     double yval = y0 + iy*dmin;
    796                     for (int ix=0; ix < nx; ix++)
    797                     {
    798                         double xval = x0 + ix*dmin;
    799                         double v = field.value(xval,yval,zval);
    800 
    801                         if (v != 0.0f && v < nzero_min)
    802                         {
    803                             nzero_min = v;
    804                         }
    805 
    806                         // scale all values [0-1], -1 => out of bounds
    807                         v = (isnan(v)) ? -1.0 : v;
    808 
    809                         cdata[ngen] = v;
    810                         ++ngen;
    811                     }
    812                 }
    813             }
    814 
    815             float* data = computeGradient(cdata, nx, ny, nz, field.valueMin(), field.valueMax());
    816 
    817             // Compute the gradient of this data.  BE CAREFUL: center
    818             /*
    819             float *data = new float[4*nx*ny*nz];
    820 
    821             double vmin = field.valueMin();
    822             double dv = field.valueMax() - field.valueMin();
     408            double v;
     409            printf("test2\n");
     410                        fflush(stdout);
    823411            if (dv == 0.0) { dv = 1.0; }
    824 
    825             // generate the uniformly sampled data that we need for a volume
    826             int ngen = 0;
    827             double nzero_min = 0.0;
    828             for (int iz=0; iz < nz; iz++) {
    829                 double zval = z0 + iz*dmin;
    830                 for (int iy=0; iy < ny; iy++) {
    831                     double yval = y0 + iy*dmin;
    832                     for (int ix=0; ix < nx; ix++) {
    833                         double xval = x0 + ix*dmin;
    834                         double v = field.value(xval,yval,zval);
    835 
    836                         if (v != 0.0f && v < nzero_min)
    837                         {
    838                             nzero_min = v;
    839                         }
    840 
    841                         // scale all values [0-1], -1 => out of bounds
    842                         v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
    843 
    844                         data[ngen] = v;
    845                         ngen += 4;
    846                     }
    847                 }
    848             }
     412            for (int i = 0; i < count; ++i)
     413            {
     414                v = data[ngen];
     415                // scale all values [0-1], -1 => out of bounds
     416                v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
     417                data[ngen] = v;
     418                ngen += 4;
     419            }
     420
    849421            // Compute the gradient of this data.  BE CAREFUL: center
    850422            // calculation on each node to avoid skew in either direction.
    … …  
    854426                    for (int ix=0; ix < nx; ix++) {
    855427                        // gradient in x-direction
    856                         double valm1 = (ix == 0) ? 0.0 : data[ngen-4];
    857                         double valp1 = (ix == nx-1) ? 0.0 : data[ngen+4];
     428                        double valm1 = (ix == 0) ? 0.0 : data[ngen - 4];
     429                        double valp1 = (ix == nx-1) ? 0.0 : data[ngen + 4];
    858430                        if (valm1 < 0 || valp1 < 0) {
    859431                            data[ngen+1] = 0.0;
    860432                        } else {
    861                             data[ngen+1] = ((valp1-valm1) + 1) *  0.5; // assume dx=1
     433                            data[ngen+1] = valp1-valm1; // assume dx=1
     434                            //data[ngen+1] = ((valp1-valm1) + 1) *  0.5; // assume dx=1
    862435                        }
    863436
    … …  
    868441                            data[ngen+2] = 0.0;
    869442                        } else {
    870                             //data[ngen+2] = valp1-valm1; // assume dy=1
    871                             data[ngen+2] = ((valp1-valm1) + 1) *  0.5; // assume dx=1
     443                            data[ngen+2] = valp1-valm1; // assume dx=1
     444                            //data[ngen+2] = ((valp1-valm1) + 1) *  0.5; // assume dy=1
    872445                        }
    873446
    … …  
    878451                            data[ngen+3] = 0.0;
    879452                        } else {
    880                             //data[ngen+3] = valp1-valm1; // assume dz=1
    881                             data[ngen+3] = ((valp1-valm1) + 1) * 0.5; // assume dz=1
     453                            data[ngen+3] = valp1-valm1; // assume dx=1
     454                            //data[ngen+3] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1
    882455                        }
    883456
    … …  
    886459                }
    887460            }
    888             */
    889 
     461
     462            dx = nx;
     463            dy = ny;
     464            dz = nz;
    890465            NanoVis::load_volume(index, nx, ny, nz, 4, data,
    891                 field.valueMin(), field.valueMax(), nzero_min);
    892 
    893             // TBD..
    894             // POINTSET
    895             /*
    896             PointSet* pset = new PointSet();
    897             pset->initialize(volume[index], (float*) data);
    898             pset->setVisible(true);
    899             NanoVis::pointSet.push_back(pset);
    900             updateColor(pset);
    901             NanoVis::volume[index]->pointsetIndex = NanoVis::pointSet.size() - 1;
    902             */
    903  
     466                vmin, vmax, nzero_min);
     467
    904468            delete [] data;
    905469
    … …  
    994558                    for (int ix=0; ix < nx; ix++) {
    995559                        // gradient in x-direction
    996                         double valm1 = (ix == 0) ? 0.0 : data[ngen-1];
    997                         double valp1 = (ix == nx-1) ? 0.0 : data[ngen+1];
     560                        double valm1 = (ix == 0) ? 0.0 : data[ngen-4];
     561                        double valp1 = (ix == nx-1) ? 0.0 : data[ngen+4];
    998562                        if (valm1 < 0 || valp1 < 0) {
    999563                            data[ngen+1] = 0.0;
    1000564                        } else {
    1001565                            //data[ngen+1] = valp1-valm1; // assume dx=1
    1002                             data[ngen+1] = ((valp1-valm1) + 1) *  0.5; // assume dx=1
     566                            data[ngen+1] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1
    1003567                        }
    1004568
    … …  
    1010574                        } else {
    1011575                            //data[ngen+2] = valp1-valm1; // assume dy=1
    1012                             data[ngen+2] = ((valp1-valm1) + 1) *  0.5; // assume dy=1
     576                            data[ngen+2] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1
    1013577                        }
    1014578
    … …  
    1020584                        } else {
    1021585                            //data[ngen+3] = valp1-valm1; // assume dz=1
     586                            data[ngen+3] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1
     587                        }
     588
     589                        ngen += 4;
     590                    }
     591                }
     592            }
     593
     594            NanoVis::load_volume(index, nx, ny, nz, 4, data,
     595                field.valueMin(), field.valueMax(), nzero_min);
     596
     597            delete [] data;
     598        }
     599    } else {
     600        return result.error("data not found in stream");
     601    }
     602
     603    //
     604    // Center this new volume on the origin.
     605    //
     606    float dx0 = -0.5;
     607    float dy0 = -0.5*dy/dx;
     608    float dz0 = -0.5*dz/dx;
     609    NanoVis::volume[index]->move(Vector3(dx0, dy0, dz0));
     610
     611    return result;
     612}
     613
     614Rappture::Outcome
     615load_volume_stream(int index, std::iostream& fin)
     616{
     617    Rappture::Outcome result;
     618
     619    Rappture::MeshTri2D xymesh;
     620    int dummy, nx, ny, nz, nxy, npts;
     621    double x0, y0, z0, dx, dy, dz, ddx, ddy, ddz;
     622    char line[128], type[128], *start;
     623
     624    int isrect = 1;
     625
     626    dx = dy = dz = 0.0;         // Suppress compiler warning.
     627    while (!fin.eof()) {
     628        fin.getline(line, sizeof(line) - 1);
     629        if (fin.fail()) {
     630            return result.error("error in data stream");
     631        }
     632        for (start=line; *start == ' ' || *start == '\t'; start++)
     633            ;  // skip leading blanks
     634
     635        if (*start != '#') {  // skip comment lines
     636            if (sscanf(start, "object %d class gridpositions counts %d %d %d", &dummy, &nx, &ny, &nz) == 4) {
     637                // found grid size
     638                isrect = 1;
     639            } else if (sscanf(start, "object %d class array type float rank 1 shape 3 items %d data follows", &dummy, &nxy) == 2) {
     640                isrect = 0;
     641
     642                double xx, yy, zz;
     643                for (int i=0; i < nxy; i++) {
     644                    fin.getline(line,sizeof(line)-1);
     645                    if (sscanf(line, "%lg %lg %lg", &xx, &yy, &zz) == 3) {
     646                        xymesh.addNode( Rappture::Node2D(xx,yy) );
     647                    }
     648                }
     649
     650                char fpts[128];
     651                sprintf(fpts, "/tmp/tmppts%d", getpid());
     652                char fcells[128];
     653                sprintf(fcells, "/tmp/tmpcells%d", getpid());
     654
     655                std::ofstream ftmp(fpts);
     656                // save corners of bounding box first, to work around meshing
     657                // problems in voronoi utility
     658                ftmp << xymesh.rangeMin(Rappture::xaxis) << " "
     659                     << xymesh.rangeMin(Rappture::yaxis) << std::endl;
     660                ftmp << xymesh.rangeMax(Rappture::xaxis) << " "
     661                     << xymesh.rangeMin(Rappture::yaxis) << std::endl;
     662                ftmp << xymesh.rangeMax(Rappture::xaxis) << " "
     663                     << xymesh.rangeMax(Rappture::yaxis) << std::endl;
     664                ftmp << xymesh.rangeMin(Rappture::xaxis) << " "
     665                     << xymesh.rangeMax(Rappture::yaxis) << std::endl;
     666                for (int i=0; i < nxy; i++) {
     667                    ftmp << xymesh.atNode(i).x() << " " << xymesh.atNode(i).y() << std::endl;
     668               
     669                }
     670                ftmp.close();
     671
     672                char cmdstr[512];
     673                sprintf(cmdstr, "voronoi -t < %s > %s", fpts, fcells);
     674                if (system(cmdstr) == 0) {
     675                    int cx, cy, cz;
     676                    std::ifstream ftri(fcells);
     677                    while (!ftri.eof()) {
     678                        ftri.getline(line,sizeof(line)-1);
     679                        if (sscanf(line, "%d %d %d", &cx, &cy, &cz) == 3) {
     680                            if (cx >= 4 && cy >= 4 && cz >= 4) {
     681                                // skip first 4 boundary points
     682                                xymesh.addCell(cx-4, cy-4, cz-4);
     683                            }
     684                        }
     685                    }
     686                    ftri.close();
     687                } else {
     688                    return result.error("triangularization failed");
     689                }
     690
     691                sprintf(cmdstr, "rm -f %s %s", fpts, fcells);
     692                system(cmdstr);
     693            } else if (sscanf(start, "object %d class regulararray count %d", &dummy, &nz) == 2) {
     694                // found z-grid
     695            } else if (sscanf(start, "origin %lg %lg %lg", &x0, &y0, &z0) == 3) {
     696                // found origin
     697            } else if (sscanf(start, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {
     698                // found one of the delta lines
     699                if (ddx != 0.0) { dx = ddx; }
     700                else if (ddy != 0.0) { dy = ddy; }
     701                else if (ddz != 0.0) { dz = ddz; }
     702            } else if (sscanf(start, "object %d class array type %s rank 0 items %d data follows", &dummy, type, &npts) == 3) {
     703                if (isrect && (npts != nx*ny*nz)) {
     704                    char mesg[256];
     705                    sprintf(mesg,"inconsistent data: expected %d points but found %d points", nx*ny*nz, npts);
     706                    return result.error(mesg);
     707                } else if (!isrect && (npts != nxy*nz)) {
     708                    char mesg[256];
     709                    sprintf(mesg,"inconsistent data: expected %d points but found %d points", nxy*nz, npts);
     710                    return result.error(mesg);
     711                }
     712                break;
     713            } else if (sscanf(start, "object %d class array type %s rank 0 times %d data follows", &dummy, type, &npts) == 3) {
     714                if (npts != nx*ny*nz) {
     715                    char mesg[256];
     716                    sprintf(mesg,"inconsistent data: expected %d points but found %d points", nx*ny*nz, npts);
     717                    return result.error(mesg);
     718                }
     719                break;
     720            }
     721        }
     722    }
     723
     724    // read data points
     725    if (!fin.eof()) {
     726        if (isrect) {
     727            Rappture::Mesh1D xgrid(x0, x0+nx*dx, nx);
     728            Rappture::Mesh1D ygrid(y0, y0+ny*dy, ny);
     729            Rappture::Mesh1D zgrid(z0, z0+nz*dz, nz);
     730            Rappture::FieldRect3D field(xgrid, ygrid, zgrid);
     731
     732            double dval[6];
     733            int nread = 0;
     734            int ix = 0;
     735            int iy = 0;
     736            int iz = 0;
     737            while (!fin.eof() && nread < npts) {
     738                fin.getline(line,sizeof(line)-1);
     739                if (fin.fail()) {
     740                    return result.error("error reading data points");
     741                }
     742                int n = sscanf(line, "%lg %lg %lg %lg %lg %lg", &dval[0], &dval[1], &dval[2], &dval[3], &dval[4], &dval[5]);
     743
     744                for (int p=0; p < n; p++) {
     745                    int nindex = iz*nx*ny + iy*nx + ix;
     746                    field.define(nindex, dval[p]);
     747                    nread++;
     748                    if (++iz >= nz) {
     749                        iz = 0;
     750                        if (++iy >= ny) {
     751                            iy = 0;
     752                            ++ix;
     753                        }
     754                    }
     755                }
     756            }
     757
     758            // make sure that we read all of the expected points
     759            if (nread != nx*ny*nz) {
     760                char mesg[256];
     761                sprintf(mesg,"inconsistent data: expected %d points but found %d points", nx*ny*nz, nread);
     762                result.error(mesg);
     763                return result;
     764            }
     765
     766            // figure out a good mesh spacing
     767            int nsample = 30;
     768            dx = field.rangeMax(Rappture::xaxis) - field.rangeMin(Rappture::xaxis);
     769            dy = field.rangeMax(Rappture::yaxis) - field.rangeMin(Rappture::yaxis);
     770            dz = field.rangeMax(Rappture::zaxis) - field.rangeMin(Rappture::zaxis);
     771            double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);
     772
     773            nx = (int)ceil(dx/dmin);
     774            ny = (int)ceil(dy/dmin);
     775            nz = (int)ceil(dz/dmin);
     776
     777#ifndef NV40
     778            // must be an even power of 2 for older cards
     779                nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
     780                ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
     781                nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
     782#endif
     783
     784            float *cdata = new float[nx*ny*nz];
     785            int ngen = 0;
     786            double nzero_min = 0.0;
     787            for (int iz=0; iz < nz; iz++) {
     788                double zval = z0 + iz*dmin;
     789                for (int iy=0; iy < ny; iy++) {
     790                    double yval = y0 + iy*dmin;
     791                    for (int ix=0; ix < nx; ix++) {
     792                        double xval = x0 + ix*dmin;
     793                        double v = field.value(xval,yval,zval);
     794
     795                        if (v != 0.0f && v < nzero_min) {
     796                            nzero_min = v;
     797                        }
     798
     799                        // scale all values [0-1], -1 => out of bounds
     800                        v = (isnan(v)) ? -1.0 : v;
     801
     802                        cdata[ngen] = v;
     803                        ++ngen;
     804                    }
     805                }
     806            }
     807
     808            float* data = computeGradient(cdata, nx, ny, nz, field.valueMin(),
     809                                          field.valueMax());
     810
     811            // Compute the gradient of this data.  BE CAREFUL: center
     812            /*
     813            float *data = new float[4*nx*ny*nz];
     814
     815            double vmin = field.valueMin();
     816            double dv = field.valueMax() - field.valueMin();
     817            if (dv == 0.0) { dv = 1.0; }
     818
     819            // generate the uniformly sampled data that we need for a volume
     820            int ngen = 0;
     821            double nzero_min = 0.0;
     822            for (int iz=0; iz < nz; iz++) {
     823                double zval = z0 + iz*dmin;
     824                for (int iy=0; iy < ny; iy++) {
     825                    double yval = y0 + iy*dmin;
     826                    for (int ix=0; ix < nx; ix++) {
     827                        double xval = x0 + ix*dmin;
     828                        double v = field.value(xval,yval,zval);
     829
     830                        if (v != 0.0f && v < nzero_min)
     831                        {
     832                            nzero_min = v;
     833                        }
     834
     835                        // scale all values [0-1], -1 => out of bounds
     836                        v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
     837
     838                        data[ngen] = v;
     839                        ngen += 4;
     840                    }
     841                }
     842            }
     843            // Compute the gradient of this data.  BE CAREFUL: center
     844            // calculation on each node to avoid skew in either direction.
     845            ngen = 0;
     846            for (int iz=0; iz < nz; iz++) {
     847                for (int iy=0; iy < ny; iy++) {
     848                    for (int ix=0; ix < nx; ix++) {
     849                        // gradient in x-direction
     850                        double valm1 = (ix == 0) ? 0.0 : data[ngen-4];
     851                        double valp1 = (ix == nx-1) ? 0.0 : data[ngen+4];
     852                        if (valm1 < 0 || valp1 < 0) {
     853                            data[ngen+1] = 0.0;
     854                        } else {
     855                            data[ngen+1] = ((valp1-valm1) + 1) *  0.5; // assume dx=1
     856                        }
     857
     858                        // gradient in y-direction
     859                        valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
     860                        valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
     861                        if (valm1 < 0 || valp1 < 0) {
     862                            data[ngen+2] = 0.0;
     863                        } else {
     864                            //data[ngen+2] = valp1-valm1; // assume dy=1
     865                            data[ngen+2] = ((valp1-valm1) + 1) *  0.5; // assume dx=1
     866                        }
     867
     868                        // gradient in z-direction
     869                        valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
     870                        valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
     871                        if (valm1 < 0 || valp1 < 0) {
     872                            data[ngen+3] = 0.0;
     873                        } else {
     874                            //data[ngen+3] = valp1-valm1; // assume dz=1
    1022875                            data[ngen+3] = ((valp1-valm1) + 1) *  0.5; // assume dz=1
    1023876                        }
    … …  
    1027880                }
    1028881            }
    1029 
    1030                 NanoVis::load_volume(index, nx, ny, nz, 4, data,
     882            */
     883
     884            NanoVis::load_volume(index, nx, ny, nz, 4, data,
    1031885                field.valueMin(), field.valueMax(), nzero_min);
    1032886
    … …  
    1042896            */
    1043897 
     898            delete [] data;
     899
     900        } else {
     901            Rappture::Mesh1D zgrid(z0, z0+nz*dz, nz);
     902            Rappture::FieldPrism3D field(xymesh, zgrid);
     903
     904            double dval;
     905            int nread = 0;
     906            int ixy = 0;
     907            int iz = 0;
     908            while (!fin.eof() && nread < npts) {
     909                fin >> dval;
     910                if (fin.fail()) {
     911                    char mesg[256];
     912                    sprintf(mesg,"after %d of %d points: can't read number",
     913                            nread, npts);
     914                    return result.error(mesg);
     915                } else {
     916                    int nid = nxy*iz + ixy;
     917                    field.define(nid, dval);
     918
     919                    nread++;
     920                    if (++iz >= nz) {
     921                        iz = 0;
     922                        ixy++;
     923                    }
     924                }
     925            }
     926
     927            // make sure that we read all of the expected points
     928            if (nread != nxy*nz) {
     929                char mesg[256];
     930                sprintf(mesg,"inconsistent data: expected %d points but found %d points", nxy*nz, nread);
     931                return result.error(mesg);
     932            }
     933
     934            // figure out a good mesh spacing
     935            int nsample = 30;
     936            x0 = field.rangeMin(Rappture::xaxis);
     937            dx = field.rangeMax(Rappture::xaxis) - field.rangeMin(Rappture::xaxis);
     938            y0 = field.rangeMin(Rappture::yaxis);
     939            dy = field.rangeMax(Rappture::yaxis) - field.rangeMin(Rappture::yaxis);
     940            z0 = field.rangeMin(Rappture::zaxis);
     941            dz = field.rangeMax(Rappture::zaxis) - field.rangeMin(Rappture::zaxis);
     942            double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);
     943
     944            nx = (int)ceil(dx/dmin);
     945            ny = (int)ceil(dy/dmin);
     946            nz = (int)ceil(dz/dmin);
     947#ifndef NV40
     948            // must be an even power of 2 for older cards
     949                nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
     950                ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
     951                nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
     952#endif
     953            float *data = new float[4*nx*ny*nz];
     954
     955            double vmin = field.valueMin();
     956            double dv = field.valueMax() - field.valueMin();
     957            if (dv == 0.0) { dv = 1.0; }
     958
     959            // generate the uniformly sampled data that we need for a volume
     960            int ngen = 0;
     961            double nzero_min = 0.0;
     962            for (iz=0; iz < nz; iz++) {
     963                double zval = z0 + iz*dmin;
     964                for (int iy=0; iy < ny; iy++) {
     965                    double yval = y0 + iy*dmin;
     966                    for (int ix=0; ix < nx; ix++) {
     967                        double xval = x0 + ix*dmin;
     968                        double v = field.value(xval,yval,zval);
     969
     970                        if (v != 0.0f && v < nzero_min)
     971                        {
     972                            nzero_min = v;
     973                        }
     974                        // scale all values [0-1], -1 => out of bounds
     975                        v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
     976                        data[ngen] = v;
     977
     978                        ngen += 4;
     979                    }
     980                }
     981            }
     982
     983            // Compute the gradient of this data.  BE CAREFUL: center
     984            // calculation on each node to avoid skew in either direction.
     985            ngen = 0;
     986            for (int iz=0; iz < nz; iz++) {
     987                for (int iy=0; iy < ny; iy++) {
     988                    for (int ix=0; ix < nx; ix++) {
     989                        // gradient in x-direction
     990                        double valm1 = (ix == 0) ? 0.0 : data[ngen-1];
     991                        double valp1 = (ix == nx-1) ? 0.0 : data[ngen+1];
     992                        if (valm1 < 0 || valp1 < 0) {
     993                            data[ngen+1] = 0.0;
     994                        } else {
     995                            //data[ngen+1] = valp1-valm1; // assume dx=1
     996                            data[ngen+1] = ((valp1-valm1) + 1) *  0.5; // assume dx=1
     997                        }
     998
     999                        // gradient in y-direction
     1000                        valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
     1001                        valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
     1002                        if (valm1 < 0 || valp1 < 0) {
     1003                            data[ngen+2] = 0.0;
     1004                        } else {
     1005                            //data[ngen+2] = valp1-valm1; // assume dy=1
     1006                            data[ngen+2] = ((valp1-valm1) + 1) *  0.5; // assume dy=1
     1007                        }
     1008
     1009                        // gradient in z-direction
     1010                        valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
     1011                        valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
     1012                        if (valm1 < 0 || valp1 < 0) {
     1013                            data[ngen+3] = 0.0;
     1014                        } else {
     1015                            //data[ngen+3] = valp1-valm1; // assume dz=1
     1016                            data[ngen+3] = ((valp1-valm1) + 1) *  0.5; // assume dz=1
     1017                        }
     1018
     1019                        ngen += 4;
     1020                    }
     1021                }
     1022            }
     1023
     1024                NanoVis::load_volume(index, nx, ny, nz, 4, data,
     1025                field.valueMin(), field.valueMax(), nzero_min);
     1026
     1027            // TBD..
     1028            // POINTSET
     1029            /*
     1030            PointSet* pset = new PointSet();
     1031            pset->initialize(volume[index], (float*) data);
     1032            pset->setVisible(true);
     1033            NanoVis::pointSet.push_back(pset);
     1034            updateColor(pset);
     1035            NanoVis::volume[index]->pointsetIndex = NanoVis::pointSet.size() - 1;
     1036            */
     1037 
    10441038
    10451039            delete [] data;
Note: See TracChangeset for help on using the changeset viewer.