- Timestamp:
- Feb 26, 2008, 7:27:47 PM (17 years ago)
- Location:
- trunk/vizservers/nanovis
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/vizservers/nanovis/Command.cpp
r906 r911 36 36 */ 37 37 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 39 43 40 44 #include "Command.h" … … 70 74 #include <NvLIC.h> 71 75 72 // FOR testing new functions 73 //#define _LOCAL_ZINC_TEST_ 74 #ifdef _LOCAL_ZINC_TEXT_ 76 #if _LOCAL_ZINC_TEST_ 75 77 #include "Test.h" 76 78 #endif 77 #include "Test.h" 79 78 80 // EXTERN DECLARATIONS 79 81 // in Nv.cpp … … 83 85 // in nanovis.cpp 84 86 extern vector<PointSet*> g_pointSet; 85 86 extern float live_rot_x; //object rotation angles87 extern float live_rot_y;88 extern float live_rot_z;89 extern float live_obj_x; //object translation location from the origin90 extern float live_obj_y;91 extern float live_obj_z;92 87 93 88 extern PlaneRenderer* plane_render; … … 139 134 static Tcl_CmdProc GridCmd; 140 135 static Tcl_CmdProc LegendCmd; 136 #if PLANE_CMDS 141 137 static Tcl_CmdProc PlaneEnableCmd; 142 138 static Tcl_CmdProc PlaneLinkCmd; 143 139 static Tcl_CmdProc PlaneNewCmd; 140 #endif 144 141 static Tcl_CmdProc ScreenCmd; 145 142 static Tcl_CmdProc ScreenShotCmd; … … 163 160 vector<unsigned int>* vectorPtr); 164 161 static int GetAxis(Tcl_Interp *interp, const char *string, int *valPtr); 165 static int GetColor(Tcl_Interp *interp, const char *string, float *rgbPtr); 162 static int GetColor(Tcl_Interp *interp, int argc, const char *argv[], 163 float *rgbPtr); 166 164 static int GetDataStream(Tcl_Interp *interp, Rappture::Buffer &buf, 167 165 int nBytes); … … 670 668 } 671 669 672 #if defISO_TEST670 #if ISO_TEST 673 671 /* 674 672 for (i=0; i < nslots; i++) { … … 850 848 header[5] = '\0'; 851 849 852 #if def_LOCAL_ZINC_TEST_850 #if _LOCAL_ZINC_TEST_ 853 851 //FILE* fp = fopen("/home/nanohub/vrinside/nv/data/HOON/QDWL_100_100_50_strain_8000i.nd_zatom_12_1", "rb"); 854 852 FILE* fp; … … 875 873 //vol = NvZincBlendeReconstructor::getInstance()->loadFromStream(fdata); 876 874 877 #if def_LOCAL_ZINC_TEST_875 #if _LOCAL_ZINC_TEST_ 878 876 vol = NvZincBlendeReconstructor::getInstance()->loadFromMemory(b); 879 877 #else … … 901 899 NanoVis::volume[n] = vol; 902 900 } 903 #ifdef __TEST_CODE__ 904 /* 901 #if __TEST_CODE__ 905 902 } else if (strcmp(header, "<FET>") == 0) { 906 903 printf("FET loading...\n"); … … 914 911 return TCL_ERROR; 915 912 } 916 */917 913 #endif /*__TEST_CODE__*/ 918 914 } else if (strcmp(header, "<ODX>") == 0) { … … 934 930 std::stringstream fdata; 935 931 fdata.write(buf.bytes(),buf.size()); 936 #ifndef ISO_TEST 932 #if ISO_TEST 933 err = load_volume_stream2(n, fdata); 934 #else 937 935 err = load_volume_stream(n, fdata); 938 #else939 err = load_volume_stream2(n, fdata);940 936 #endif 941 937 if (err) { … … 1033 1029 } 1034 1030 } else if ((c == 'c') && (strcmp(argv[2],"color") == 0)) { 1035 if (argc < 3) {1031 if (argc < 6) { 1036 1032 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); 1038 1034 return TCL_ERROR; 1039 1035 } 1040 1036 float rgb[3]; 1041 if (GetColor(interp, arg v[3], rgb) != TCL_OK) {1037 if (GetColor(interp, argc - 3, argv + 3, rgb) != TCL_OK) { 1042 1038 return TCL_ERROR; 1043 1039 } 1044 1040 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) { 1046 1042 return TCL_ERROR; 1047 1043 } … … 2093 2089 */ 2094 2090 static 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); 2091 GetColor(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); 2106 2096 return TCL_ERROR; 2107 2097 } … … 2109 2099 (GetFloat(interp, argv[1], rgbPtr + 1) != TCL_OK) || 2110 2100 (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 2119 2108 static int 2120 2109 PlaneNewCmd(ClientData cdata, Tcl_Interp *interp, int argc, … … 2208 2197 return TCL_OK; 2209 2198 } 2210 2199 #endif /*PLANE_CMDS*/ 2211 2200 2212 2201 void … … 2281 2270 (ClientData)NULL, (Tcl_CmdDeleteProc*)NULL); 2282 2271 2283 #if def__TEST_CODE__2272 #if __TEST_CODE__ 2284 2273 Tcl_CreateCommand(interp, "test", TestCmd, 2285 2274 (ClientData)NULL, (Tcl_CmdDeleteProc*)NULL); -
trunk/vizservers/nanovis/Trace.cpp
r848 r911 1 1 #include <Trace.h> 2 2 #include <stdio.h> 3 #include <stdarg.h> 3 4 4 void Trace(c har* format, ...)5 void Trace(const char* format, ...) 5 6 { 6 7 char buff[1024]; 8 va_list lst; 7 9 8 va_list lst;9 10 va_start(lst, format); 10 vs printf(buff, format, lst);11 11 vsnprintf(buff, 1023, format, lst); 12 buff[1023] = '\0'; 12 13 printf(buff); 13 14 fflush(stdout); -
trunk/vizservers/nanovis/Trace.h
r848 r911 2 2 #define __TRACE_H__ 3 3 4 #include <stdarg.h> 5 #include <stdio.h> 6 7 extern void Trace(char* format, ...); 4 extern void Trace(const char* format, ...); 8 5 9 6 #endif -
trunk/vizservers/nanovis/Vector3.h
r900 r911 61 61 v.y = y + scalar; 62 62 v.z = z + scalar; 63 return v; 63 64 } 64 65 … … 69 70 v.y = y - scalar; 70 71 v.z = z - scalar; 72 return v; 71 73 } 72 74 #endif -
trunk/vizservers/nanovis/dxReader.cpp
r883 r911 40 40 #include "GradientFilter.h" 41 41 42 //#define _LOCAL_ZINC_TEST_ 43 float* merge(float* scalar, float* gradient, int size) 42 #define _LOCAL_ZINC_TEST_ 0 43 44 static float * 45 merge(float* scalar, float* gradient, int size) 44 46 { 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; 58 59 } 59 60 60 void normalizeScalar(float* fdata, int count, float min, float max) 61 static void 62 normalizeScalar(float* fdata, int count, float min, float max) 61 63 { 62 64 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) { 66 67 fdata[i] = fdata[i] / v; 68 } 67 69 } 68 70 } 69 71 70 float* computeGradient(float* fdata, int width, int height, int depth, float min, float max) 72 static float* 73 computeGradient(float* fdata, int width, int height, int depth, 74 float min, float max) 71 75 { 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; 82 87 } 83 88 … … 88 93 load_vector_stream(int index, std::istream& fin) 89 94 { 90 int dummy, nx, ny, nz, n xy, npts;95 int dummy, nx, ny, nz, npts; 91 96 double x0, y0, z0, dx, dy, dz, ddx, ddy, ddz; 92 97 char line[128], type[128], *start; 93 98 99 dx = dy = dz = 0.0; // Suppress compiler warning. 94 100 while (!fin.eof()) { 95 101 fin.getline(line, sizeof(line) - 1); … … 104 110 if (sscanf(start, "object %d class gridpositions counts %d %d %d", &dummy, &nx, &ny, &nz) == 4) { 105 111 // 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) { 108 113 // 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) { 111 115 // 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) { 117 124 if (npts != nx*ny*nz) { 118 125 std::cerr << "inconsistent data: expected " << nx*ny*nz << " points but found " << npts << " points" << std::endl; … … 120 127 } 121 128 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) { 124 130 if (npts != nx*ny*nz) { 125 131 std::cerr << "inconsistent data: expected " << nx*ny*nz << " points but found " << npts << " points" << std::endl; … … 224 230 225 231 // scale should be accounted. 226 for (ngen=0; ngen < npts; ngen++) 227 { 232 for (ngen=0; ngen < npts; ngen++) { 228 233 data[ngen] = (data[ngen]/(2.0*vmax) + 0.5); 229 234 } 230 231 235 NanoVis::load_volume(index, nx, ny, nz, 3, data, vmin, vmax, nzero_min); 232 233 236 delete [] data; 234 } 235 else { 237 } else { 236 238 std::cerr << "WARNING: data not found in stream" << std::endl; 237 239 } … … 252 254 253 255 int isrect = 1; 254 255 float* voldata = 0; 256 dx = dy = dz = 0.0; // Suppress compiler warning. 256 257 do { 257 258 fin.getline(line,sizeof(line)-1); … … 279 280 result.error(mesg); 280 281 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 meshing289 // problems in voronoi utility290 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 points314 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-grid328 }329 else if (sscanf(start, "origin %lg %lg %lg", &x0, &y0, &z0) == 3) {330 // found origin331 }332 else if (sscanf(start, "delta %lg %lg %lg", &ddx, &ddy, &ddz) == 3) {333 // found one of the delta lines334 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 points363 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 points404 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 bounds422 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: center428 // 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-direction434 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=1440 //data[ngen+1] = ((valp1-valm1) + 1) * 0.5; // assume dx=1441 }442 443 // gradient in y-direction444 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=1450 //data[ngen+2] = ((valp1-valm1) + 1) * 0.5; // assume dy=1451 }452 453 // gradient in z-direction454 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=1460 //data[ngen+3] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1461 }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 points504 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 spacing511 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 NV40524 // must be an even power of 2 for older cards525 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 #endif529 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 volume536 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 bounds551 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: center560 // 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-direction566 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=1572 data[ngen+1] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1573 }574 575 // gradient in y-direction576 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=1582 data[ngen+2] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1583 }584 585 // gradient in z-direction586 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=1592 data[ngen+3] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1593 }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::Outcome621 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 blanks639 640 if (*start != '#') { // skip comment lines641 if (sscanf(start, "object %d class gridpositions counts %d %d %d", &dummy, &nx, &ny, &nz) == 4) {642 // found grid size643 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 }654 282 655 283 char fpts[128]; … … 725 353 } 726 354 } 727 } 355 } while (!fin.eof()); 728 356 729 357 // read data points 730 358 if (!fin.eof()) { 731 359 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 737 360 double dval[6]; 738 361 int nread = 0; … … 740 363 int iy = 0; 741 364 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 742 372 while (!fin.eof() && nread < npts) { 743 373 fin.getline(line,sizeof(line)-1); 744 if (fin.fail()) {745 return result.error("error reading data points");746 }747 374 int n = sscanf(line, "%lg %lg %lg %lg %lg %lg", &dval[0], &dval[1], &dval[2], &dval[3], &dval[4], &dval[5]); 748 375 749 376 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 752 386 nread++; 753 387 if (++iz >= nz) { … … 769 403 } 770 404 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; 790 407 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); 823 411 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 849 421 // Compute the gradient of this data. BE CAREFUL: center 850 422 // calculation on each node to avoid skew in either direction. … … 854 426 for (int ix=0; ix < nx; ix++) { 855 427 // 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]; 858 430 if (valm1 < 0 || valp1 < 0) { 859 431 data[ngen+1] = 0.0; 860 432 } 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 862 435 } 863 436 … … 868 441 data[ngen+2] = 0.0; 869 442 } else { 870 //data[ngen+2] = valp1-valm1; // assume dy=1871 data[ngen+2] = ((valp1-valm1) + 1) * 0.5; // assume dx=1443 data[ngen+2] = valp1-valm1; // assume dx=1 444 //data[ngen+2] = ((valp1-valm1) + 1) * 0.5; // assume dy=1 872 445 } 873 446 … … 878 451 data[ngen+3] = 0.0; 879 452 } else { 880 //data[ngen+3] = valp1-valm1; // assume dz=1881 data[ngen+3] = ((valp1-valm1) + 1) *0.5; // assume dz=1453 data[ngen+3] = valp1-valm1; // assume dx=1 454 //data[ngen+3] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1 882 455 } 883 456 … … 886 459 } 887 460 } 888 */ 889 461 462 dx = nx; 463 dy = ny; 464 dz = nz; 890 465 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 904 468 delete [] data; 905 469 … … 994 558 for (int ix=0; ix < nx; ix++) { 995 559 // 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]; 998 562 if (valm1 < 0 || valp1 < 0) { 999 563 data[ngen+1] = 0.0; 1000 564 } else { 1001 565 //data[ngen+1] = valp1-valm1; // assume dx=1 1002 data[ngen+1] = ((valp1-valm1) + 1 ) * 0.5; // assume dx=1566 data[ngen+1] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1 1003 567 } 1004 568 … … 1010 574 } else { 1011 575 //data[ngen+2] = valp1-valm1; // assume dy=1 1012 data[ngen+2] = ((valp1-valm1) + 1 ) * 0.5; // assume dy=1576 data[ngen+2] = ((valp1-valm1) + 1.0) * 0.5; // assume dz=1 1013 577 } 1014 578 … … 1020 584 } else { 1021 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); 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 1022 875 data[ngen+3] = ((valp1-valm1) + 1) * 0.5; // assume dz=1 1023 876 } … … 1027 880 } 1028 881 } 1029 1030 NanoVis::load_volume(index, nx, ny, nz, 4, data, 882 */ 883 884 NanoVis::load_volume(index, nx, ny, nz, 4, data, 1031 885 field.valueMin(), field.valueMax(), nzero_min); 1032 886 … … 1042 896 */ 1043 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 } 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 1044 1038 1045 1039 delete [] data;
Note: See TracChangeset
for help on using the changeset viewer.