Ignore:
Timestamp:
Mar 6, 2008 4:01:20 PM (16 years ago)
Author:
gah
Message:

collect limits for axes

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/vizservers/nanovis/dxReader.cpp

    r911 r927  
    1 
     1 
    22/*
    33 * ----------------------------------------------------------------------
    … …  
    55 *
    66 *  dxReader.cpp
    7  *      This module contains openDX readers for 2D and 3D volumes.
     7 *      This module contains openDX readers for 2D and 3D volumes.
    88 *
    99 * ======================================================================
    … …  
    5151    int ngen = 0, sindex = 0;
    5252    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;
     53        data[ngen++] = scalar[sindex];
     54        data[ngen++] = g[sindex].x;
     55        data[ngen++] = g[sindex].y;
     56        data[ngen++] = g[sindex].z;
    5757    }
    5858    return data;
    … …  
    6666        for (int i = 0; i < count; ++i) {
    6767            fdata[i] = fdata[i] / v;
    68         }
     68        }
    6969    }
    7070}
    … …  
    7272static float*
    7373computeGradient(float* fdata, int width, int height, int depth,
    74                 float min, float max)
     74                float min, float max)
    7575{
    7676    float* gradients = (float *)malloc(width * height * depth * 3 *
    77                                        sizeof(float));
     77                                       sizeof(float));
    7878    float* tempGradients = (float *)malloc(width * height * depth * 3 *
    79                                            sizeof(float));
     79                                           sizeof(float));
    8080    int sizes[3] = { width, height, depth };
    8181    computeGradients(tempGradients, fdata, sizes, DATRAW_FLOAT);
    … …  
    9797    char line[128], type[128], *start;
    9898
    99     dx = dy = dz = 0.0;         // Suppress compiler warning.
     99    dx = dy = dz = 0.0;         // Suppress compiler warning.
    100100    while (!fin.eof()) {
    101101        fin.getline(line, sizeof(line) - 1);
    102         if (fin.fail()) {
    103             //return result.error("error in data stream");
    104             return;
    105         }
     102        if (fin.fail()) {
     103            //return result.error("error in data stream");
     104            return;
     105        }
    106106        for (start=&line[0]; *start == ' ' || *start == '\t'; start++)
    107107            ;  // skip leading blanks
    … …  
    115115                // found one of the delta lines
    116116                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                 }
     117                    dx = ddx;
     118                } else if (ddy != 0.0) {
     119                    dy = ddy;
     120                } else if (ddz != 0.0) {
     121                    dz = ddz;
     122                }
    123123            } else if (sscanf(start, "object %d class array type %s shape 3 rank 1 items %d data follows", &dummy, type, &npts) == 3) {
    124124                if (npts != nx*ny*nz) {
    … …  
    185185#ifndef NV40
    186186        // must be an even power of 2 for older cards
    187             nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
    188             ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
    189             nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
     187        nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
     188        ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
     189        nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
    190190#endif
    191191
    … …  
    211211                    vz = zfield.value(xval,yval,zval);
    212212
    213                     double vm = sqrt(vx*vx + vy*vy + vz*vz);
    214 
    215                     if (vm < vmin) { vmin = vm; }
    216                     if (vm > vmax) { vmax = vm; }
    217                     if (vm != 0.0f && vm < nzero_min)
    218                     {
     213                    double vm;
     214                    vm = sqrt(vx*vx + vy*vy + vz*vz);
     215                    if (vm < vmin) {
     216                        vmin = vm;
     217                    } else if (vm > vmax) {
     218                        vmax = vm;
     219                    }
     220                    if ((vm != 0.0f) && (vm < nzero_min)) {
    219221                        nzero_min = vm;
    220222                    }
    221 
    222223                    data[ngen++] = vx;
    223224                    data[ngen++] = vy;
    … …  
    233234            data[ngen] = (data[ngen]/(2.0*vmax) + 0.5);
    234235        }
    235         NanoVis::load_volume(index, nx, ny, nz, 3, data, vmin, vmax, nzero_min);
     236        Volume *volPtr;
     237        volPtr = NanoVis::load_volume(index, nx, ny, nz, 3, data, vmin, vmax,
     238                nzero_min);
     239        volPtr->set_limits(0, x0, x0 + (nx * ddx));
     240        volPtr->set_limits(1, y0, y0 + (ny * ddy));
     241        volPtr->set_limits(2, z0, z0 + (nz * ddz));
    236242        delete [] data;
    237243    } else {
    … …  
    254260
    255261    int isrect = 1;
    256     dx = dy = dz = 0.0;         // Suppress compiler warning.
     262    dx = dy = dz = 0.0;         // Suppress compiler warning.
    257263    do {
    258264        fin.getline(line,sizeof(line)-1);
    … …  
    265271                // found grid size
    266272                isrect = 1;
    267             }
    268             else if (sscanf(start, "object %d class array type float rank 1 shape 3 items %d data follows", &dummy, &nxy) == 2) {
     273            } else if (sscanf(start, "object %d class array type float rank 1 shape 3 items %d data follows", &dummy, &nxy) == 2) {
    269274                isrect = 0;
    270275                double xx, yy, zz;
    … …  
    280285                result.error(mesg);
    281286                return result;
     287
     288                char fpts[128];
     289                sprintf(fpts, "/tmp/tmppts%d", getpid());
     290                char fcells[128];
     291                sprintf(fcells, "/tmp/tmpcells%d", getpid());
     292
     293                std::ofstream ftmp(fpts);
     294                // save corners of bounding box first, to work around meshing
     295                // problems in voronoi utility
     296                ftmp << xymesh.rangeMin(Rappture::xaxis) << " "
     297                     << xymesh.rangeMin(Rappture::yaxis) << std::endl;
     298                ftmp << xymesh.rangeMax(Rappture::xaxis) << " "
     299                     << xymesh.rangeMin(Rappture::yaxis) << std::endl;
     300                ftmp << xymesh.rangeMax(Rappture::xaxis) << " "
     301                     << xymesh.rangeMax(Rappture::yaxis) << std::endl;
     302                ftmp << xymesh.rangeMin(Rappture::xaxis) << " "
     303                     << xymesh.rangeMax(Rappture::yaxis) << std::endl;
     304                for (int i=0; i < nxy; i++) {
     305                    ftmp << xymesh.atNode(i).x() << " " << xymesh.atNode(i).y() << std::endl;
     306               
     307                }
     308                ftmp.close();
     309
     310                char cmdstr[512];
     311                sprintf(cmdstr, "voronoi -t < %s > %s", fpts, fcells);
     312                if (system(cmdstr) == 0) {
     313                    int cx, cy, cz;
     314                    std::ifstream ftri(fcells);
     315                    while (!ftri.eof()) {
     316                        ftri.getline(line,sizeof(line)-1);
     317                        if (sscanf(line, "%d %d %d", &cx, &cy, &cz) == 3) {
     318                            if (cx >= 4 && cy >= 4 && cz >= 4) {
     319                                // skip first 4 boundary points
     320                                xymesh.addCell(cx-4, cy-4, cz-4);
     321                            }
     322                        }
     323                    }
     324                    ftri.close();
     325                } else {
     326                    return result.error("triangularization failed");
     327                }
     328
     329                sprintf(cmdstr, "rm -f %s %s", fpts, fcells);
     330                system(cmdstr);
     331            } else if (sscanf(start, "object %d class regulararray count %d", &dummy, &nz) == 2) {
     332                // found z-grid
     333            } else if (sscanf(start, "origin %lg %lg %lg", &x0, &y0, &z0) == 3) {
     334                // found origin
     335            } else if (sscanf(start, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {
     336                int count = 0;
     337                // 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                }
     354            } else if (sscanf(start, "object %d class array type %s rank 0 items %d data follows", &dummy, type, &npts) == 3) {
     355                if (isrect && (npts != nx*ny*nz)) {
     356                    char mesg[256];
     357                    sprintf(mesg,"inconsistent data: expected %d points but found %d points", nx*ny*nz, npts);
     358                    return result.error(mesg);
     359                } else if (!isrect && (npts != nxy*nz)) {
     360                    char mesg[256];
     361                    sprintf(mesg,"inconsistent data: expected %d points but found %d points", nxy*nz, npts);
     362                    return result.error(mesg);
     363                }
     364                break;
     365            } else if (sscanf(start, "object %d class array type %s rank 0 times %d data follows", &dummy, type, &npts) == 3) {
     366                if (npts != nx*ny*nz) {
     367                    char mesg[256];
     368                    sprintf(mesg,"inconsistent data: expected %d points but found %d points", nx*ny*nz, npts);
     369                    return result.error(mesg);
     370                }
     371                break;
     372            }
     373        }
     374    } while (!fin.eof());
     375
     376    // read data points
     377    if (!fin.eof()) {
     378        if (isrect) {
     379            double dval[6];
     380            int nread = 0;
     381            int ix = 0;
     382            int iy = 0;
     383            int iz = 0;
     384            float* data = new float[nx *  ny *  nz * 4];
     385            memset(data, 0, nx*ny*nz*4);
     386            double vmin = 1e21;
     387            double nzero_min = 1e21;
     388            double vmax = -1e21;
     389
     390
     391            while (!fin.eof() && nread < npts) {
     392                fin.getline(line,sizeof(line)-1);
     393                int n = sscanf(line, "%lg %lg %lg %lg %lg %lg", &dval[0], &dval[1], &dval[2], &dval[3], &dval[4], &dval[5]);
     394
     395                for (int p=0; p < n; p++) {
     396                    int nindex = (iz*nx*ny + iy*nx + ix) * 4;
     397                    data[nindex] = dval[p];
     398
     399                    if (dval[p] < vmin) {
     400                        vmin = dval[p];
     401                    } else if (dval[p] > vmax) {
     402                        vmax = dval[p];
     403                    }
     404                    if (dval[p] != 0.0f && dval[p] < nzero_min) {
     405                        nzero_min = dval[p];
     406                    }
     407
     408                    nread++;
     409                    if (++iz >= nz) {
     410                        iz = 0;
     411                        if (++iy >= ny) {
     412                            iy = 0;
     413                            ++ix;
     414                        }
     415                    }
     416                }
     417            }
     418
     419            // make sure that we read all of the expected points
     420            if (nread != nx*ny*nz) {
     421                char mesg[256];
     422                sprintf(mesg,"inconsistent data: expected %d points but found %d points", nx*ny*nz, nread);
     423                result.error(mesg);
     424                return result;
     425            }
     426
     427            double dv = vmax - vmin;
     428            int count = nx*ny*nz;
     429            int ngen = 0;
     430            double v;
     431            printf("test2\n");
     432            fflush(stdout);
     433            if (dv == 0.0) {
     434                dv = 1.0;
     435            }
     436            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            }
     443            // Compute the gradient of this data.  BE CAREFUL: center
     444            // calculation on each node to avoid skew in either direction.
     445            ngen = 0;
     446            for (int iz=0; iz < nz; iz++) {
     447                for (int iy=0; iy < ny; iy++) {
     448                    for (int ix=0; ix < nx; ix++) {
     449                        // gradient in x-direction
     450                        double valm1 = (ix == 0) ? 0.0 : data[ngen - 4];
     451                        double valp1 = (ix == nx-1) ? 0.0 : data[ngen + 4];
     452                        if (valm1 < 0 || valp1 < 0) {
     453                            data[ngen+1] = 0.0;
     454                        } else {
     455                            data[ngen+1] = valp1-valm1; // assume dx=1
     456                            //data[ngen+1] = ((valp1-valm1) + 1) *  0.5; // assume dx=1
     457                        }
     458
     459                        // gradient in y-direction
     460                        valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
     461                        valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
     462                        if (valm1 < 0 || valp1 < 0) {
     463                            data[ngen+2] = 0.0;
     464                        } else {
     465                            data[ngen+2] = valp1-valm1; // assume dx=1
     466                            //data[ngen+2] = ((valp1-valm1) + 1) *  0.5; // assume dy=1
     467                        }
     468
     469                        // gradient in z-direction
     470                        valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
     471                        valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
     472                        if (valm1 < 0 || valp1 < 0) {
     473                            data[ngen+3] = 0.0;
     474                        } else {
     475                            data[ngen+3] = valp1-valm1; // assume dx=1
     476                            //data[ngen+3] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1
     477                        }
     478
     479                        ngen += 4;
     480                    }
     481                }
     482            }
     483
     484            dx = nx;
     485            dy = ny;
     486            dz = nz;
     487            Volume *volPtr;
     488            volPtr = NanoVis::load_volume(index, nx, ny, nz, 4, data,
     489                                          vmin, vmax, nzero_min);
     490            volPtr->set_limits(NanoVis::X, x0, x0 + (nx * ddx));
     491            volPtr->set_limits(NanoVis::Y, y0, y0 + (ny * ddy));
     492            volPtr->set_limits(NanoVis::Z, z0, z0 + (nz * ddz));
     493            delete [] data;
     494
     495        } else {
     496            Rappture::Mesh1D zgrid(z0, z0+nz*dz, nz);
     497            Rappture::FieldPrism3D field(xymesh, zgrid);
     498
     499            double dval;
     500            int nread = 0;
     501            int ixy = 0;
     502            int iz = 0;
     503            while (!fin.eof() && nread < npts) {
     504                fin >> dval;
     505                if (fin.fail()) {
     506                    char mesg[256];
     507                    sprintf(mesg,"after %d of %d points: can't read number",
     508                            nread, npts);
     509                    return result.error(mesg);
     510                } else {
     511                    int nid = nxy*iz + ixy;
     512                    field.define(nid, dval);
     513
     514                    nread++;
     515                    if (++iz >= nz) {
     516                        iz = 0;
     517                        ixy++;
     518                    }
     519                }
     520            }
     521
     522            // make sure that we read all of the expected points
     523            if (nread != nxy*nz) {
     524                char mesg[256];
     525                sprintf(mesg,"inconsistent data: expected %d points but found %d points", nxy*nz, nread);
     526                return result.error(mesg);
     527            }
     528
     529            // figure out a good mesh spacing
     530            int nsample = 30;
     531            x0 = field.rangeMin(Rappture::xaxis);
     532            dx = field.rangeMax(Rappture::xaxis) - field.rangeMin(Rappture::xaxis);
     533            y0 = field.rangeMin(Rappture::yaxis);
     534            dy = field.rangeMax(Rappture::yaxis) - field.rangeMin(Rappture::yaxis);
     535            z0 = field.rangeMin(Rappture::zaxis);
     536            dz = field.rangeMax(Rappture::zaxis) - field.rangeMin(Rappture::zaxis);
     537            double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);
     538
     539            nx = (int)ceil(dx/dmin);
     540            ny = (int)ceil(dy/dmin);
     541            nz = (int)ceil(dz/dmin);
     542#ifndef NV40
     543            // must be an even power of 2 for older cards
     544            nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
     545            ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
     546            nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
     547#endif
     548            float *data = new float[4*nx*ny*nz];
     549
     550            double vmin = field.valueMin();
     551            double dv = field.valueMax() - field.valueMin();
     552            if (dv == 0.0) {
     553                dv = 1.0;
     554            }
     555            // generate the uniformly sampled data that we need for a volume
     556            int ngen = 0;
     557            double nzero_min = 0.0;
     558            for (iz=0; iz < nz; iz++) {
     559                double zval = z0 + iz*dmin;
     560                for (int iy=0; iy < ny; iy++) {
     561                    double yval = y0 + iy*dmin;
     562                    for (int ix=0; ix < nx; ix++) {
     563                        double xval = x0 + ix*dmin;
     564                        double v = field.value(xval,yval,zval);
     565
     566                        if (v != 0.0f && v < nzero_min)
     567                            {
     568                                nzero_min = v;
     569                            }
     570                        // scale all values [0-1], -1 => out of bounds
     571                        v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
     572                        data[ngen] = v;
     573
     574                        ngen += 4;
     575                    }
     576                }
     577            }
     578
     579            // Compute the gradient of this data.  BE CAREFUL: center
     580            // calculation on each node to avoid skew in either direction.
     581            ngen = 0;
     582            for (int iz=0; iz < nz; iz++) {
     583                for (int iy=0; iy < ny; iy++) {
     584                    for (int ix=0; ix < nx; ix++) {
     585                        // gradient in x-direction
     586                        double valm1 = (ix == 0) ? 0.0 : data[ngen-4];
     587                        double valp1 = (ix == nx-1) ? 0.0 : data[ngen+4];
     588                        if (valm1 < 0 || valp1 < 0) {
     589                            data[ngen+1] = 0.0;
     590                        } else {
     591                            //data[ngen+1] = valp1-valm1; // assume dx=1
     592                            data[ngen+1] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1
     593                        }
     594
     595                        // gradient in y-direction
     596                        valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
     597                        valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
     598                        if (valm1 < 0 || valp1 < 0) {
     599                            data[ngen+2] = 0.0;
     600                        } else {
     601                            //data[ngen+2] = valp1-valm1; // assume dy=1
     602                            data[ngen+2] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1
     603                        }
     604
     605                        // gradient in z-direction
     606                        valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
     607                        valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
     608                        if (valm1 < 0 || valp1 < 0) {
     609                            data[ngen+3] = 0.0;
     610                        } else {
     611                            //data[ngen+3] = valp1-valm1; // assume dz=1
     612                            data[ngen+3] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1
     613                        }
     614
     615                        ngen += 4;
     616                    }
     617                }
     618            }
     619
     620            Volume *volPtr;
     621            volPtr = NanoVis::load_volume(index, nx, ny, nz, 4, data,
     622                field.valueMin(), field.valueMax(), nzero_min);
     623            volPtr->set_limits(0, field.rangeMin(Rappture::xaxis),
     624                               field.rangeMax(Rappture::xaxis));
     625            volPtr->set_limits(1, field.rangeMin(Rappture::yaxis),
     626                               field.rangeMax(Rappture::yaxis));
     627            volPtr->set_limits(2, field.rangeMin(Rappture::zaxis),
     628                               field.rangeMax(Rappture::zaxis));
     629            delete [] data;
     630        }
     631    } else {
     632        return result.error("data not found in stream");
     633    }
     634
     635    //
     636    // Center this new volume on the origin.
     637    //
     638    float dx0 = -0.5;
     639    float dy0 = -0.5*dy/dx;
     640    float dz0 = -0.5*dz/dx;
     641    NanoVis::volume[index]->move(Vector3(dx0, dy0, dz0));
     642
     643    return result;
     644}
     645
     646Rappture::Outcome
     647load_volume_stream(int index, std::iostream& fin)
     648{
     649    Rappture::Outcome result;
     650
     651    Rappture::MeshTri2D xymesh;
     652    int dummy, nx, ny, nz, nxy, npts;
     653    double x0, y0, z0, dx, dy, dz, ddx, ddy, ddz;
     654    char line[128], type[128], *start;
     655
     656    int isrect = 1;
     657
     658    dx = dy = dz = 0.0;         // Suppress compiler warning.
     659    while (!fin.eof()) {
     660        fin.getline(line, sizeof(line) - 1);
     661        if (fin.fail()) {
     662            return result.error("error in data stream");
     663        }
     664        for (start=line; *start == ' ' || *start == '\t'; start++)
     665            ;  // skip leading blanks
     666
     667        if (*start != '#') {  // skip comment lines
     668            if (sscanf(start, "object %d class gridpositions counts %d %d %d", &dummy, &nx, &ny, &nz) == 4) {
     669                // found grid size
     670                isrect = 1;
     671            } else if (sscanf(start, "object %d class array type float rank 1 shape 3 items %d data follows", &dummy, &nxy) == 2) {
     672                isrect = 0;
     673
     674                double xx, yy, zz;
     675                for (int i=0; i < nxy; i++) {
     676                    fin.getline(line,sizeof(line)-1);
     677                    if (sscanf(line, "%lg %lg %lg", &xx, &yy, &zz) == 3) {
     678                        xymesh.addNode( Rappture::Node2D(xx,yy) );
     679                    }
     680                }
    282681
    283682                char fpts[128];
    … …  
    353752            }
    354753        }
    355     } while (!fin.eof());
     754    }
    356755
    357756    // read data points
    358757    if (!fin.eof()) {
    359758        if (isrect) {
     759            Rappture::Mesh1D xgrid(x0, x0+nx*dx, nx);
     760            Rappture::Mesh1D ygrid(y0, y0+ny*dy, ny);
     761            Rappture::Mesh1D zgrid(z0, z0+nz*dz, nz);
     762            Rappture::FieldRect3D field(xgrid, ygrid, zgrid);
     763
    360764            double dval[6];
    361765            int nread = 0;
    … …  
    363767            int iy = 0;
    364768            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 
    372769            while (!fin.eof() && nread < npts) {
    373770                fin.getline(line,sizeof(line)-1);
     771                if (fin.fail()) {
     772                    return result.error("error reading data points");
     773                }
    374774                int n = sscanf(line, "%lg %lg %lg %lg %lg %lg", &dval[0], &dval[1], &dval[2], &dval[3], &dval[4], &dval[5]);
    375775
    376776                for (int p=0; p < n; 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 
     777                    int nindex = iz*nx*ny + iy*nx + ix;
     778                    field.define(nindex, dval[p]);
    386779                    nread++;
    387780                    if (++iz >= nz) {
    … …  
    403796            }
    404797
    405             double dv = vmax - vmin;
    406             int count = nx*ny*nz;
     798            // figure out a good mesh spacing
     799            int nsample = 30;
     800            dx = field.rangeMax(Rappture::xaxis) - field.rangeMin(Rappture::xaxis);
     801            dy = field.rangeMax(Rappture::yaxis) - field.rangeMin(Rappture::yaxis);
     802            dz = field.rangeMax(Rappture::zaxis) - field.rangeMin(Rappture::zaxis);
     803            double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);
     804
     805            nx = (int)ceil(dx/dmin);
     806            ny = (int)ceil(dy/dmin);
     807            nz = (int)ceil(dz/dmin);
     808
     809#ifndef NV40
     810            // must be an even power of 2 for older cards
     811            nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
     812            ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
     813            nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
     814#endif
     815
     816            float *cdata = new float[nx*ny*nz];
    407817            int ngen = 0;
    408             double v;
    409             printf("test2\n");
    410                         fflush(stdout);
    411             if (dv == 0.0) { dv = 1.0; }
    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             }
     818            double nzero_min = 0.0;
     819            for (int iz=0; iz < nz; iz++) {
     820                double zval = z0 + iz*dmin;
     821                for (int iy=0; iy < ny; iy++) {
     822                    double yval = y0 + iy*dmin;
     823                    for (int ix=0; ix < nx; ix++) {
     824                        double xval = x0 + ix*dmin;
     825                        double v = field.value(xval,yval,zval);
     826
     827                        if (v != 0.0f && v < nzero_min) {
     828                            nzero_min = v;
     829                        }
     830
     831                        // scale all values [0-1], -1 => out of bounds
     832                        v = (isnan(v)) ? -1.0 : v;
     833
     834                        cdata[ngen] = v;
     835                        ++ngen;
     836                    }
     837                }
     838            }
     839
     840            float* data = computeGradient(cdata, nx, ny, nz, field.valueMin(),
     841                                          field.valueMax());
    420842
    421843            // Compute the gradient of this data.  BE CAREFUL: center
    422             // calculation on each node to avoid skew in either direction.
    423             ngen = 0;
    424             for (int iz=0; iz < nz; iz++) {
    425                 for (int iy=0; iy < ny; iy++) {
    426                     for (int ix=0; ix < nx; ix++) {
    427                         // gradient in x-direction
    428                         double valm1 = (ix == 0) ? 0.0 : data[ngen - 4];
    429                         double valp1 = (ix == nx-1) ? 0.0 : data[ngen + 4];
    430                         if (valm1 < 0 || valp1 < 0) {
    431                             data[ngen+1] = 0.0;
    432                         } else {
    433                             data[ngen+1] = valp1-valm1; // assume dx=1
    434                             //data[ngen+1] = ((valp1-valm1) + 1) *  0.5; // assume dx=1
    435                         }
    436 
    437                         // gradient in y-direction
    438                         valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
    439                         valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
    440                         if (valm1 < 0 || valp1 < 0) {
    441                             data[ngen+2] = 0.0;
    442                         } else {
    443                             data[ngen+2] = valp1-valm1; // assume dx=1
    444                             //data[ngen+2] = ((valp1-valm1) + 1) *  0.5; // assume dy=1
    445                         }
    446 
    447                         // gradient in z-direction
    448                         valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
    449                         valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
    450                         if (valm1 < 0 || valp1 < 0) {
    451                             data[ngen+3] = 0.0;
    452                         } else {
    453                             data[ngen+3] = valp1-valm1; // assume dx=1
    454                             //data[ngen+3] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1
    455                         }
    456 
    457                         ngen += 4;
    458                     }
    459                 }
    460             }
    461 
    462             dx = nx;
    463             dy = ny;
    464             dz = nz;
    465             NanoVis::load_volume(index, nx, ny, nz, 4, data,
    466                 vmin, vmax, nzero_min);
    467 
     844            /*
     845              float *data = new float[4*nx*ny*nz];
     846
     847              double vmin = field.valueMin();
     848              double dv = field.valueMax() - field.valueMin();
     849              if (dv == 0.0) { dv = 1.0; }
     850
     851              // generate the uniformly sampled data that we need for a volume
     852              int ngen = 0;
     853              double nzero_min = 0.0;
     854              for (int iz=0; iz < nz; iz++) {
     855              double zval = z0 + iz*dmin;
     856              for (int iy=0; iy < ny; iy++) {
     857              double yval = y0 + iy*dmin;
     858              for (int ix=0; ix < nx; ix++) {
     859              double xval = x0 + ix*dmin;
     860              double v = field.value(xval,yval,zval);
     861
     862              if (v != 0.0f && v < nzero_min)
     863              {
     864              nzero_min = v;
     865              }
     866
     867              // scale all values [0-1], -1 => out of bounds
     868              v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
     869
     870              data[ngen] = v;
     871              ngen += 4;
     872              }
     873              }
     874              }
     875              // Compute the gradient of this data.  BE CAREFUL: center
     876              // calculation on each node to avoid skew in either direction.
     877              ngen = 0;
     878              for (int iz=0; iz < nz; iz++) {
     879              for (int iy=0; iy < ny; iy++) {
     880              for (int ix=0; ix < nx; ix++) {
     881              // gradient in x-direction
     882              double valm1 = (ix == 0) ? 0.0 : data[ngen-4];
     883              double valp1 = (ix == nx-1) ? 0.0 : data[ngen+4];
     884              if (valm1 < 0 || valp1 < 0) {
     885              data[ngen+1] = 0.0;
     886              } else {
     887              data[ngen+1] = ((valp1-valm1) + 1) *  0.5; // assume dx=1
     888              }
     889
     890              // gradient in y-direction
     891              valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
     892              valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
     893              if (valm1 < 0 || valp1 < 0) {
     894              data[ngen+2] = 0.0;
     895              } else {
     896              //data[ngen+2] = valp1-valm1; // assume dy=1
     897              data[ngen+2] = ((valp1-valm1) + 1) *  0.5; // assume dx=1
     898              }
     899
     900              // gradient in z-direction
     901              valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
     902              valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
     903              if (valm1 < 0 || valp1 < 0) {
     904              data[ngen+3] = 0.0;
     905              } else {
     906              //data[ngen+3] = valp1-valm1; // assume dz=1
     907              data[ngen+3] = ((valp1-valm1) + 1) *  0.5; // assume dz=1
     908              }
     909
     910              ngen += 4;
     911              }
     912              }
     913              }
     914            */
     915
     916            Volume *volPtr;
     917            volPtr = NanoVis::load_volume(index, nx, ny, nz, 4, data,
     918                field.valueMin(), field.valueMax(), nzero_min);
     919            volPtr->set_limits(0, field.rangeMin(Rappture::xaxis),
     920                               field.rangeMax(Rappture::xaxis));
     921            volPtr->set_limits(1, field.rangeMin(Rappture::yaxis),
     922                               field.rangeMax(Rappture::yaxis));
     923            volPtr->set_limits(2, field.rangeMin(Rappture::zaxis),
     924                               field.rangeMax(Rappture::zaxis));
     925            // TBD..
     926            // POINTSET
     927            /*
     928              PointSet* pset = new PointSet();
     929              pset->initialize(volume[index], (float*) data);
     930              pset->setVisible(true);
     931              NanoVis::pointSet.push_back(pset);
     932              updateColor(pset);
     933              NanoVis::volume[index]->pointsetIndex = NanoVis::pointSet.size() - 1;
     934            */
     935 
    468936            delete [] data;
    469937
    … …  
    477945            int iz = 0;
    478946            while (!fin.eof() && nread < npts) {
    479                 fin >> dval;
     947                fin >> dval;
    480948                if (fin.fail()) {
    481                     char mesg[256];
    482                     sprintf(mesg,"after %d of %d points: can't read number",
    483                             nread, npts);
    484                     return result.error(mesg);
    485                 } else {
     949                    char mesg[256];
     950                    sprintf(mesg,"after %d of %d points: can't read number",
     951                            nread, npts);
     952                    return result.error(mesg);
     953                } else {
    486954                    int nid = nxy*iz + ixy;
    487955                    field.define(nid, dval);
    … …  
    517985#ifndef NV40
    518986            // must be an even power of 2 for older cards
    519                 nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
    520                 ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
    521                 nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
     987            nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
     988            ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
     989            nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
    522990#endif
    523991            float *data = new float[4*nx*ny*nz];
    … …  
    5391007
    5401008                        if (v != 0.0f && v < nzero_min)
    541                         {
    542                             nzero_min = v;
    543                         }
    544                         // scale all values [0-1], -1 => out of bounds
    545                         v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
    546                         data[ngen] = v;
    547 
    548                         ngen += 4;
    549                     }
    550                 }
    551             }
    552 
    553             // Compute the gradient of this data.  BE CAREFUL: center
    554             // calculation on each node to avoid skew in either direction.
    555             ngen = 0;
    556             for (int iz=0; iz < nz; iz++) {
    557                 for (int iy=0; iy < ny; iy++) {
    558                     for (int ix=0; ix < nx; ix++) {
    559                         // gradient in x-direction
    560                         double valm1 = (ix == 0) ? 0.0 : data[ngen-4];
    561                         double valp1 = (ix == nx-1) ? 0.0 : data[ngen+4];
    562                         if (valm1 < 0 || valp1 < 0) {
    563                             data[ngen+1] = 0.0;
    564                         } else {
    565                             //data[ngen+1] = valp1-valm1; // assume dx=1
    566                             data[ngen+1] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1
    567                         }
    568 
    569                         // gradient in y-direction
    570                         valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
    571                         valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
    572                         if (valm1 < 0 || valp1 < 0) {
    573                             data[ngen+2] = 0.0;
    574                         } else {
    575                             //data[ngen+2] = valp1-valm1; // assume dy=1
    576                             data[ngen+2] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1
    577                         }
    578 
    579                         // gradient in z-direction
    580                         valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
    581                         valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
    582                         if (valm1 < 0 || valp1 < 0) {
    583                             data[ngen+3] = 0.0;
    584                         } else {
    585                             //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 
    614 Rappture::Outcome
    615 load_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);
     1009                            {
     1010                                nzero_min = v;
    6831011                            }
    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
    875                             data[ngen+3] = ((valp1-valm1) + 1) *  0.5; // assume dz=1
    876                         }
    877 
    878                         ngen += 4;
    879                     }
    880                 }
    881             }
    882             */
    883 
    884             NanoVis::load_volume(index, nx, ny, nz, 4, data,
    885                 field.valueMin(), field.valueMax(), nzero_min);
    886 
    887             // TBD..
    888             // POINTSET
    889             /*
    890             PointSet* pset = new PointSet();
    891             pset->initialize(volume[index], (float*) data);
    892             pset->setVisible(true);
    893             NanoVis::pointSet.push_back(pset);
    894             updateColor(pset);
    895             NanoVis::volume[index]->pointsetIndex = NanoVis::pointSet.size() - 1;
    896             */
    897  
    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                         }
    9741012                        // scale all values [0-1], -1 => out of bounds
    9751013                        v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
    … …  
    10221060            }
    10231061
    1024                 NanoVis::load_volume(index, nx, ny, nz, 4, data,
    1025                 field.valueMin(), field.valueMax(), nzero_min);
     1062            Volume *volPtr;
     1063            volPtr = NanoVis::load_volume(index, nx, ny, nz, 4, data,
     1064                field.valueMin(), field.valueMax(), nzero_min);
     1065            volPtr->set_limits(0, field.rangeMin(Rappture::xaxis),
     1066                               field.rangeMax(Rappture::xaxis));
     1067            volPtr->set_limits(1, field.rangeMin(Rappture::yaxis),
     1068                               field.rangeMax(Rappture::yaxis));
     1069            volPtr->set_limits(2, field.rangeMin(Rappture::zaxis),
     1070                               field.rangeMax(Rappture::zaxis));
    10261071
    10271072            // TBD..
    10281073            // POINTSET
    10291074            /*
    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;
     1075              PointSet* pset = new PointSet();
     1076              pset->initialize(volume[index], (float*) data);
     1077              pset->setVisible(true);
     1078              NanoVis::pointSet.push_back(pset);
     1079              updateColor(pset);
     1080              NanoVis::volume[index]->pointsetIndex = NanoVis::pointSet.size() - 1;
    10361081            */
    10371082 
Note: See TracChangeset for help on using the changeset viewer.