Changeset 2838 for trunk


Ignore:
Timestamp:
Mar 10, 2012, 7:25:22 PM (13 years ago)
Author:
ldelgass
Message:

Remove unused function

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/packages/vizservers/nanovis/dxReader.cpp

    r2831 r2838  
    833833    return volPtr;
    834834}
    835 
    836 
    837 Volume *
    838 load_volume_stream_insoo(Rappture::Outcome &result, const char *tag,
    839                          std::iostream& fin)
    840 {
    841     TRACE("load_volume_stream\n");
    842  
    843     Rappture::MeshTri2D xymesh;
    844     int dummy, nx, ny, nz, nxy, npts;
    845     double x0, y0, z0, dx, dy, dz, ddx, ddy, ddz;
    846     char line[128], type[128], *start;
    847 
    848     int isrect = 1;
    849 
    850     dx = dy = dz = 0.0;         // Suppress compiler warning.
    851     x0 = y0 = z0 = 0.0;         // May not have an origin line.
    852     nx = ny = nz = npts = nxy = 0;
    853     while (!fin.eof()) {
    854         fin.getline(line, sizeof(line) - 1);
    855         if (fin.fail()) {
    856             result.addError("line \"%s\"error in data stream");
    857             return NULL;
    858         }
    859         for (start=line; *start == ' ' || *start == '\t'; start++)
    860             ;  // skip leading blanks
    861 
    862         if (*start != '#') {  // skip comment lines
    863             if (sscanf(start, "object %d class gridpositions counts %d %d %d", &dummy, &nx, &ny, &nz) == 4) {
    864                 // found grid size
    865                 isrect = 1;
    866             } else if (sscanf(start, "object %d class array type float rank 1 shape 3 items %d data follows", &dummy, &nxy) == 2) {
    867                 isrect = 0;
    868 
    869                 double xx, yy, zz;
    870                 for (int i=0; i < nxy; i++) {
    871                     fin.getline(line,sizeof(line)-1);
    872                     if (sscanf(line, "%lg %lg %lg", &xx, &yy, &zz) == 3) {
    873                         xymesh.addNode( Rappture::Node2D(xx,yy) );
    874                     }
    875                 }
    876 
    877                 char fpts[128];
    878                 sprintf(fpts, "/tmp/tmppts%d", getpid());
    879                 char fcells[128];
    880                 sprintf(fcells, "/tmp/tmpcells%d", getpid());
    881 
    882                 std::ofstream ftmp(fpts);
    883                 // save corners of bounding box first, to work around meshing
    884                 // problems in voronoi utility
    885                 ftmp << xymesh.rangeMin(Rappture::xaxis) << " "
    886                      << xymesh.rangeMin(Rappture::yaxis) << std::endl;
    887                 ftmp << xymesh.rangeMax(Rappture::xaxis) << " "
    888                      << xymesh.rangeMin(Rappture::yaxis) << std::endl;
    889                 ftmp << xymesh.rangeMax(Rappture::xaxis) << " "
    890                      << xymesh.rangeMax(Rappture::yaxis) << std::endl;
    891                 ftmp << xymesh.rangeMin(Rappture::xaxis) << " "
    892                      << xymesh.rangeMax(Rappture::yaxis) << std::endl;
    893                 for (int i=0; i < nxy; i++) {
    894                     ftmp << xymesh.atNode(i).x() << " " << xymesh.atNode(i).y() << std::endl;
    895 
    896                 }
    897                 ftmp.close();
    898 
    899                 char cmdstr[512];
    900                 sprintf(cmdstr, "voronoi -t < %s > %s", fpts, fcells);
    901                 if (system(cmdstr) == 0) {
    902                     int cx, cy, cz;
    903                     std::ifstream ftri(fcells);
    904                     while (!ftri.eof()) {
    905                         ftri.getline(line,sizeof(line)-1);
    906                         if (sscanf(line, "%d %d %d", &cx, &cy, &cz) == 3) {
    907                             if (cx >= 4 && cy >= 4 && cz >= 4) {
    908                                 // skip first 4 boundary points
    909                                 xymesh.addCell(cx-4, cy-4, cz-4);
    910                             }
    911                         }
    912                     }
    913                     ftri.close();
    914                 } else {
    915                     result.error("triangularization failed");
    916                     return NULL;
    917                 }
    918                 unlink(fpts), unlink(fcells);
    919             } else if (sscanf(start, "object %d class regulararray count %d", &dummy, &nz) == 2) {
    920                 // found z-grid
    921             } else if (sscanf(start, "origin %lg %lg %lg", &x0, &y0, &z0) == 3) {
    922                 // found origin
    923             } else if (sscanf(start, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {
    924                 // found one of the delta lines
    925                 if (ddx != 0.0) { dx = ddx; }
    926                 else if (ddy != 0.0) { dy = ddy; }
    927                 else if (ddz != 0.0) { dz = ddz; }
    928             } else if (sscanf(start, "object %d class array type %s rank 0 items %d data follows", &dummy, type, &npts) == 3) {
    929                 if (isrect && (npts != nx*ny*nz)) {
    930                     result.addError("inconsistent data: expected %d points"
    931                                     " but found %d points", nx*ny*nz, npts);
    932                     return NULL;
    933                 } else if (!isrect && (npts != nxy*nz)) {
    934                     result.addError("inconsistent data: expected %d points"
    935                                     " but found %d points", nx*ny*nz, npts);
    936                     return NULL;
    937                 }
    938                 break;
    939             } else if (sscanf(start, "object %d class array type %s rank 0 times %d data follows", &dummy, type, &npts) == 3) {
    940                 if (npts != nx*ny*nz) {
    941                     result.addError("inconsistent data: expected %d points"
    942                                     " but found %d points", nx*ny*nz, npts);
    943                     return NULL;
    944                 }
    945                 break;
    946             }
    947         }
    948     }
    949 
    950     // read data points
    951     if (fin.eof()) {
    952         result.error("data not found in stream");
    953         return NULL;
    954     }
    955 
    956     Volume* volPtr = 0;
    957     if (isrect) {
    958         Rappture::Mesh1D xgrid(x0, x0+nx*dx, nx);
    959         Rappture::Mesh1D ygrid(y0, y0+ny*dy, ny);
    960         Rappture::Mesh1D zgrid(z0, z0+nz*dz, nz);
    961         Rappture::FieldRect3D field(xgrid, ygrid, zgrid);
    962        
    963         double dval[6];
    964         int nread = 0;
    965         int ix = 0;
    966         int iy = 0;
    967         int iz = 0;
    968         while (!fin.eof() && nread < npts) {
    969             fin.getline(line,sizeof(line)-1);
    970             if (fin.fail()) {
    971                 result.error("error reading data points");
    972                 return NULL;
    973             }
    974             int n = sscanf(line, "%lg %lg %lg %lg %lg %lg", &dval[0], &dval[1], &dval[2], &dval[3], &dval[4], &dval[5]);
    975            
    976             for (int p=0; p < n; p++) {
    977                 int nindex = iz*nx*ny + iy*nx + ix;
    978                 field.define(nindex, dval[p]);
    979                 TRACE("nindex = %i\tdval[%i] = %lg\n", nindex, p, dval[p]);
    980                 nread++;
    981                 if (++iz >= nz) {
    982                     iz = 0;
    983                     if (++iy >= ny) {
    984                         iy = 0;
    985                         ++ix;
    986                     }
    987                 }
    988             }
    989         }
    990        
    991         // make sure that we read all of the expected points
    992         if (nread != nx*ny*nz) {
    993             result.addError("inconsistent data: expected %d points"
    994                             " but found %d points", nx*ny*nz, npts);
    995             return NULL;
    996         }
    997        
    998         // figure out a good mesh spacing
    999         int nsample = 30;
    1000         dx = field.rangeMax(Rappture::xaxis) - field.rangeMin(Rappture::xaxis);
    1001         dy = field.rangeMax(Rappture::yaxis) - field.rangeMin(Rappture::yaxis);
    1002         dz = field.rangeMax(Rappture::zaxis) - field.rangeMin(Rappture::zaxis);
    1003         double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);
    1004        
    1005         nx = (int)ceil(dx/dmin);
    1006         ny = (int)ceil(dy/dmin);
    1007         nz = (int)ceil(dz/dmin);
    1008        
    1009 #ifndef NV40
    1010         // must be an even power of 2 for older cards
    1011         nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
    1012         ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
    1013         nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
    1014 #endif
    1015        
    1016         float *data = new float[4*nx*ny*nz];
    1017        
    1018         double vmin = field.valueMin();
    1019         double dv = field.valueMax() - field.valueMin();
    1020         if (dv == 0.0) { dv = 1.0; }
    1021        
    1022         int ngen = 0;
    1023         double nzero_min = 0.0;
    1024         for (iz=0; iz < nz; iz++) {
    1025             double zval = z0 + iz*dmin;
    1026             for (int iy=0; iy < ny; iy++) {
    1027                 double yval = y0 + iy*dmin;
    1028                 for (int ix=0; ix < nx; ix++) {
    1029                     double xval = x0 + ix*dmin;
    1030                     double v = field.value(xval,yval,zval);
    1031                    
    1032                     if (v != 0.0f && v < nzero_min) {
    1033                         nzero_min = v;
    1034                     }
    1035                     // scale all values [0-1], -1 => out of bounds
    1036                     v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
    1037                     data[ngen] = v;
    1038                    
    1039                     ngen += 4;
    1040                 }
    1041             }
    1042         }
    1043        
    1044         // Compute the gradient of this data.  BE CAREFUL: center
    1045         // calculation on each node to avoid skew in either direction.
    1046         ngen = 0;
    1047         for (int iz=0; iz < nz; iz++) {
    1048             for (int iy=0; iy < ny; iy++) {
    1049                 for (int ix=0; ix < nx; ix++) {
    1050                     // gradient in x-direction
    1051                     double valm1 = (ix == 0) ? 0.0 : data[ngen-1];
    1052                     double valp1 = (ix == nx-1) ? 0.0 : data[ngen+1];
    1053                     if (valm1 < 0 || valp1 < 0) {
    1054                         data[ngen+1] = 0.0;
    1055                     } else {
    1056                         data[ngen+1] = valp1-valm1; // assume dx=1
    1057                         //data[ngen+1] = ((valp1-valm1) + 1) *  0.5; // assume dx=1 (ISO)
    1058                     }
    1059                    
    1060                     // gradient in y-direction
    1061                     valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
    1062                     valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
    1063                     if (valm1 < 0 || valp1 < 0) {
    1064                         data[ngen+2] = 0.0;
    1065                     } else {
    1066                         data[ngen+2] = valp1-valm1; // assume dy=1
    1067                         //data[ngen+2] = ((valp1-valm1) + 1) *  0.5; // assume dy=1 (ISO)
    1068                     }
    1069                    
    1070                     // gradient in z-direction
    1071                     valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
    1072                     valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
    1073                     if (valm1 < 0 || valp1 < 0) {
    1074                         data[ngen+3] = 0.0;
    1075                     } else {
    1076                         data[ngen+3] = valp1-valm1; // assume dz=1
    1077                         //data[ngen+3] = ((valp1-valm1) + 1) *  0.5; // assume dz=1 (ISO)
    1078                     }
    1079                    
    1080                     ngen += 4;
    1081                 }
    1082             }
    1083         }
    1084        
    1085         /*
    1086           float *cdata = new float[nx*ny*nz];
    1087           int ngen = 0;
    1088           double nzero_min = 0.0;
    1089           for (int iz=0; iz < nz; iz++) {
    1090           double zval = z0 + iz*dmin;
    1091           for (int iy=0; iy < ny; iy++) {
    1092           double yval = y0 + iy*dmin;
    1093           for (int ix=0; ix < nx; ix++) {
    1094           double xval = x0 + ix*dmin;
    1095           double v = field.value(xval,yval,zval);
    1096          
    1097           if (v != 0.0f && v < nzero_min) {
    1098           nzero_min = v;
    1099           }
    1100          
    1101           // scale all values [0-1], -1 => out of bounds
    1102           v = (isnan(v)) ? -1.0 : v;
    1103          
    1104           cdata[ngen] = v;
    1105           ++ngen;
    1106           }
    1107                 }
    1108             }
    1109            
    1110             float* data = computeGradient(cdata, nx, ny, nz, field.valueMin(),
    1111                                           field.valueMax());
    1112                                          
    1113             for (int i=0; i<nx*ny*nz; i++) {
    1114                 TRACE("enddata[%i] = %lg\n",i,data[i]);
    1115             }
    1116 
    1117             TRACE("End Data Stats volDataID = %i\n",volDataID);
    1118             TRACE("nx = %i ny = %i nz = %i\n",nx,ny,nz);
    1119             TRACE("dx = %lg dy = %lg dz = %lg\n",dx,dy,dz);
    1120             TRACE("dataMin = %lg\tdataMax = %lg\tnzero_min = %lg\n",
    1121             field.valueMin(),field.valueMax(),nzero_min);
    1122             */
    1123 
    1124         volPtr = NanoVis::load_volume(tag, nx, ny, nz, 4, data,
    1125                 field.valueMin(), field.valueMax(), nzero_min);
    1126         volPtr->xAxis.SetRange(field.rangeMin(Rappture::xaxis),
    1127                                field.rangeMax(Rappture::xaxis));
    1128         volPtr->yAxis.SetRange(field.rangeMin(Rappture::yaxis),
    1129                                field.rangeMax(Rappture::yaxis));
    1130         volPtr->zAxis.SetRange(field.rangeMin(Rappture::zaxis),
    1131                                field.rangeMax(Rappture::zaxis));
    1132         volPtr->wAxis.SetRange(field.valueMin(), field.valueMax());
    1133         volPtr->update_pending = true;
    1134         // TBD..
    1135         // POINTSET
    1136         /*
    1137           PointSet* pset = new PointSet();
    1138           pset->initialize(volume[index], (float*) data);
    1139           pset->setVisible(true);
    1140           NanoVis::pointSet.push_back(pset);
    1141           updateColor(pset);
    1142           NanoVis::volume[index]->pointsetIndex = NanoVis::pointSet.size() - 1;
    1143         */
    1144        
    1145         delete [] data;
    1146        
    1147     } else {
    1148         Rappture::Mesh1D zgrid(z0, z0+nz*dz, nz);
    1149         Rappture::FieldPrism3D field(xymesh, zgrid);
    1150        
    1151         double dval;
    1152         int nread = 0;
    1153         int ixy = 0;
    1154         int iz = 0;
    1155         while (!fin.eof() && nread < npts) {
    1156             fin >> dval;
    1157             if (fin.fail()) {
    1158                 char mesg[256];
    1159                 sprintf(mesg,"after %d of %d points: can't read number",
    1160                         nread, npts);
    1161                 result.error(mesg);
    1162                 return NULL;
    1163             } else {
    1164                 int nid = nxy*iz + ixy;
    1165                 field.define(nid, dval);
    1166                
    1167                 nread++;
    1168                 if (++iz >= nz) {
    1169                     iz = 0;
    1170                     ixy++;
    1171                 }
    1172             }
    1173         }
    1174        
    1175         // make sure that we read all of the expected points
    1176         if (nread != nxy*nz) {
    1177             result.addError("inconsistent data: expected %d points"
    1178                             " but found %d points", nx*ny*nz, npts);
    1179             return NULL;
    1180         }
    1181        
    1182         // figure out a good mesh spacing
    1183         int nsample = 30;
    1184         x0 = field.rangeMin(Rappture::xaxis);
    1185         dx = field.rangeMax(Rappture::xaxis) - field.rangeMin(Rappture::xaxis);
    1186         y0 = field.rangeMin(Rappture::yaxis);
    1187         dy = field.rangeMax(Rappture::yaxis) - field.rangeMin(Rappture::yaxis);
    1188         z0 = field.rangeMin(Rappture::zaxis);
    1189         dz = field.rangeMax(Rappture::zaxis) - field.rangeMin(Rappture::zaxis);
    1190         double dmin = pow((dx*dy*dz)/(nsample*nsample*nsample), 0.333);
    1191        
    1192         nx = (int)ceil(dx/dmin);
    1193         ny = (int)ceil(dy/dmin);
    1194         nz = (int)ceil(dz/dmin);
    1195 #ifndef NV40
    1196         // must be an even power of 2 for older cards
    1197         nx = (int)pow(2.0, ceil(log10((double)nx)/log10(2.0)));
    1198         ny = (int)pow(2.0, ceil(log10((double)ny)/log10(2.0)));
    1199         nz = (int)pow(2.0, ceil(log10((double)nz)/log10(2.0)));
    1200 #endif
    1201         float *data = new float[4*nx*ny*nz];
    1202        
    1203         double vmin = field.valueMin();
    1204         double dv = field.valueMax() - field.valueMin();
    1205         if (dv == 0.0) { dv = 1.0; }
    1206        
    1207         // generate the uniformly sampled data that we need for a volume
    1208         int ngen = 0;
    1209         double nzero_min = 0.0;
    1210         for (iz=0; iz < nz; iz++) {
    1211             double zval = z0 + iz*dmin;
    1212             for (int iy=0; iy < ny; iy++) {
    1213                 double yval = y0 + iy*dmin;
    1214                 for (int ix=0; ix < nx; ix++) {
    1215                     double xval = x0 + ix*dmin;
    1216                     double v = field.value(xval,yval,zval);
    1217                    
    1218                     if (v != 0.0f && v < nzero_min) {
    1219                         nzero_min = v;
    1220                     }
    1221                     // scale all values [0-1], -1 => out of bounds
    1222                     v = (isnan(v)) ? -1.0 : (v - vmin)/dv;
    1223                     data[ngen] = v;
    1224                    
    1225                     ngen += 4;
    1226                 }
    1227             }
    1228         }
    1229        
    1230         // Compute the gradient of this data.  BE CAREFUL: center
    1231         // calculation on each node to avoid skew in either direction.
    1232         ngen = 0;
    1233         for (int iz=0; iz < nz; iz++) {
    1234             for (int iy=0; iy < ny; iy++) {
    1235                 for (int ix=0; ix < nx; ix++) {
    1236                     // gradient in x-direction
    1237                     double valm1 = (ix == 0) ? 0.0 : data[ngen-1];
    1238                     double valp1 = (ix == nx-1) ? 0.0 : data[ngen+1];
    1239                     if (valm1 < 0 || valp1 < 0) {
    1240                         data[ngen+1] = 0.0;
    1241                     } else {
    1242                         data[ngen+1] = valp1-valm1; // assume dx=1
    1243                         //data[ngen+1] = ((valp1-valm1) + 1) *  0.5; // assume dx=1 (ISO)
    1244                     }
    1245                    
    1246                     // gradient in y-direction
    1247                     valm1 = (iy == 0) ? 0.0 : data[ngen-4*nx];
    1248                     valp1 = (iy == ny-1) ? 0.0 : data[ngen+4*nx];
    1249                     if (valm1 < 0 || valp1 < 0) {
    1250                         data[ngen+2] = 0.0;
    1251                     } else {
    1252                         data[ngen+2] = valp1-valm1; // assume dy=1
    1253                         //data[ngen+2] = ((valp1-valm1) + 1) *  0.5; // assume dy=1 (ISO)
    1254                     }
    1255                    
    1256                     // gradient in z-direction
    1257                     valm1 = (iz == 0) ? 0.0 : data[ngen-4*nx*ny];
    1258                     valp1 = (iz == nz-1) ? 0.0 : data[ngen+4*nx*ny];
    1259                     if (valm1 < 0 || valp1 < 0) {
    1260                         data[ngen+3] = 0.0;
    1261                     } else {
    1262                         data[ngen+3] = valp1-valm1; // assume dz=1
    1263                         //data[ngen+3] = ((valp1-valm1) + 1) *  0.5; // assume dz=1 (ISO)
    1264                     }
    1265                    
    1266                     ngen += 4;
    1267                 }
    1268             }
    1269         }
    1270        
    1271         volPtr = NanoVis::load_volume(tag, nx, ny, nz, 4, data,
    1272                 field.valueMin(), field.valueMax(), nzero_min);
    1273         volPtr->xAxis.SetRange(field.rangeMin(Rappture::xaxis),
    1274                                field.rangeMax(Rappture::xaxis));
    1275         volPtr->yAxis.SetRange(field.rangeMin(Rappture::yaxis),
    1276                                field.rangeMax(Rappture::yaxis));
    1277         volPtr->zAxis.SetRange(field.rangeMin(Rappture::zaxis),
    1278                                field.rangeMax(Rappture::zaxis));
    1279         volPtr->wAxis.SetRange(field.valueMin(), field.valueMax());
    1280         volPtr->update_pending = true;
    1281         // TBD..
    1282         // POINTSET
    1283         /*
    1284           PointSet* pset = new PointSet();
    1285           pset->initialize(volume[index], (float*) data);
    1286           pset->setVisible(true);
    1287           NanoVis::pointSet.push_back(pset);
    1288           updateColor(pset);
    1289           NanoVis::volume[index]->pointsetIndex = NanoVis::pointSet.size() - 1;
    1290         */
    1291        
    1292        
    1293         delete [] data;
    1294     }
    1295 
    1296     //
    1297     // Center this new volume on the origin.
    1298     //
    1299     float dx0 = -0.5;
    1300     float dy0 = -0.5*dy/dx;
    1301     float dz0 = -0.5*dz/dx;
    1302     if (volPtr) {
    1303         volPtr->location(Vector3(dx0, dy0, dz0));
    1304     }
    1305     return volPtr;
    1306 }
Note: See TracChangeset for help on using the changeset viewer.