00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #include <stdio.h>
00024 #include <stdlib.h>
00025 #include <string.h>
00026 #include <math.h>
00027
00028 #ifndef windows_platform
00029 #include <unistd.h>
00030 #endif
00031
00032 #ifdef windows_platform
00033 #define fmaxf max
00034 #define fminf min
00035 #endif
00036
00037 #include "clH5Handler.h"
00038 #include "h5File.h"
00039 #include "h5Dataset.h"
00040 #include "h5Group.h"
00041
00042 #define JPEG
00043
00044 #ifdef JPEG
00045 #include "jpeglib.h"
00046 #include "jerror.h"
00047 #endif
00048
00049 #define NODEBUG
00050 #define NAME_LENGTH 120
00051 #define LINE_LENGTH 120
00052 #define VAR_LENGTH 5
00053 #define NCOLOR 256
00054
00055 #define DNAME_REFINE_LEVEL "refine level"
00056 #define DNAME_NODE_TYPE "node type"
00057 #define DNAME_BOUNDING_BOX "bounding box"
00058
00059 #define PRINT_ERROR(_msg) { \
00060 fprintf(stderr, "%s\n", _msg);\
00061 ret_val = -1; \
00062 goto done; \
00063 }
00064
00065 typedef struct { float r, g, b; } color;
00066
00067
00068
00069
00070
00071
00072
00073
00074 static void
00075 help(char *name)
00076 {
00077 (void) printf("\n\nNAME:\n\tislice -- Extract a slice from a FLASH output file\n\n");
00078
00079 (void) printf("DESCRIPTION:\n\tislice -- Extract a slice, perpendicular to one of the coordinate axes,\n");
00080 (void) printf( "\tfrom a FLASH output file stored at iRODS server. Selection of a slice is\n");
00081 (void) printf( "\tbased on the tool, extract_slice_from_chkpnt, developed by Paul Ricker\n");
00082 (void) printf( "\t(UIUC/NCSA).\n\n");
00083
00084 (void) fprintf (stdout, "SYNOPSIS:");
00085 (void) fprintf (stdout, "\n\t%s -h[elp], OR", name);
00086 (void) fprintf (stdout, "\n\t%s [options] infile", name);
00087
00088 (void) fprintf (stdout, "\n\n\t-h[elp]:");
00089 (void) fprintf (stdout, "\n\t\tPrint this summary of usage, and exit.");
00090
00091 (void) fprintf (stdout, "\n\n\tinfile:");
00092 (void) fprintf (stdout, "\n\t\tName of the FLASH file to read");
00093
00094 (void) fprintf (stdout, "\n\nOPTIONS:");
00095 (void) fprintf (stdout, "\n\t-o outfile:");
00096 (void) fprintf (stdout, "\n\t\tName of the output file to create");
00097
00098 (void) fprintf (stdout, "\n\n\t-p position:");
00099 (void) fprintf (stdout, "\n\t\tOrientation of the slice (0 = xy, 1 = xz, 2 = yz)");
00100
00101 (void) fprintf (stdout, "\n\n\t-m mesh_variable:");
00102 (void) fprintf (stdout, "\n\t\tName of the FLASH mesh variable");
00103
00104 (void) fprintf (stdout, "\n\n\t-c coordinate:");
00105 (void) fprintf (stdout, "\n\t\tValue of the coordinate perpendicular to the slice");
00106
00107 (void) fprintf (stdout, "\n\n\t-j jpegfile:");
00108 (void) fprintf (stdout, "\n\t\tName of the output jpeg file");
00109
00110 (void) fprintf (stdout, "\n\n\t-t color_table:");
00111 (void) fprintf (stdout, "\n\t\tName of the file of the color table. The table must be a 256x3 2D array of RGB values\n\n");
00112 }
00113
00114
00115
00116
00117
00118
00119
00120
00121 static void
00122 usage(char *name)
00123 {
00124 (void) fprintf(stderr, "\nUsage:\t%s -h[elp], OR\n", name);
00125 (void) fprintf(stderr, "\t%s [options] infile\n", name);
00126 }
00127
00128
00129
00130
00131
00132
00133
00134
00135 static char *
00136 trim(char *str, int len)
00137 {
00138 int idx=0, n=len-1;
00139 char* str_idx;
00140
00141 str_idx = str;
00142 idx = 0;
00143 while (((int) *str_idx > 32) && (idx < n)) {
00144 idx++;
00145 str_idx++;
00146 }
00147
00148 *(str+idx) = '\0';
00149
00150 return str;
00151 }
00152
00153
00154
00155
00156
00157
00158
00159
00160 static rcComm_t *
00161 rodsConnect()
00162 {
00163 int status = 0;
00164 rodsEnv env;
00165 rcComm_t *conn=NULL;
00166 rErrMsg_t errMsg;
00167
00168 status = getRodsEnv (&env);
00169
00170 if (status < 0) {
00171 rodsLogError (LOG_ERROR, status, "main: getRodsEnv error. ");
00172 return NULL;
00173 }
00174
00175 conn = rcConnect (env.rodsHost, env.rodsPort, env.rodsUserName, env.rodsZone, 1, &errMsg);
00176
00177 if (conn == NULL) {
00178 rodsLogError (LOG_ERROR, errMsg.status, "rcConnect failure %s", errMsg.msg);
00179 return NULL;
00180 }
00181
00182 status = clientLogin(conn);
00183 if (status != 0) {
00184 rcDisconnect(conn);
00185 return NULL;
00186 }
00187
00188 return conn;
00189 }
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203 static void
00204 get_defualt_outfile_name(const char *infile, char *outfile, int pos,
00205 const char *var, float coor, int dims[], const char *ext)
00206 {
00207 char infile_cp[NAME_LENGTH];
00208
00209 if (outfile==NULL || strlen(outfile) > 0)
00210 return;
00211
00212 strcpy(infile_cp, infile);
00213
00214 const char *tmpfile = strrchr(infile, '/');
00215 if (tmpfile == NULL)
00216 tmpfile = infile_cp;
00217 else
00218 tmpfile++;
00219
00220 switch (pos) {
00221 case 1:
00222 sprintf(outfile, "%s_%s_xz_%f_%dx%d.%s", tmpfile, var, coor, dims[0], dims[1], ext);
00223 break;
00224 case 2:
00225 sprintf(outfile, "%s_%s_yz_%f_%dx%d.%s", tmpfile, var, coor, dims[0], dims[1], ext);
00226 break;
00227 default:
00228 sprintf(outfile, "%s_%s_xy_%f_%dx%d.%s", tmpfile, var, coor, dims[0], dims[1], ext);
00229 break;
00230 }
00231
00232 trim(outfile, NAME_LENGTH);
00233 }
00234
00235
00236 #ifdef JPEG
00237
00238
00239
00240 static void
00241 write_jpeg_file(const char *fname, JSAMPLE *buf, int height,
00242 int width, int ncomp, int quality)
00243 {
00244
00245 FILE *jpeg_file;
00246
00247
00248 struct jpeg_compress_struct cinfo;
00249
00250
00251 struct jpeg_error_mgr jerr;
00252
00253
00254 JSAMPROW row_pointer[1];
00255
00256
00257 int row_stride;
00258
00259
00260 cinfo.err = jpeg_std_error(&jerr);
00261 jpeg_create_compress(&cinfo);
00262
00263
00264 if ((jpeg_file = fopen(fname, "wb")) == NULL) {
00265 fprintf(stderr, "can't open %s\n", fname);
00266 exit(1);
00267 }
00268 jpeg_stdio_dest(&cinfo, jpeg_file);
00269
00270
00271 cinfo.image_width = width;
00272 cinfo.image_height = height;
00273 cinfo.input_components = ncomp;
00274 cinfo.in_color_space = JCS_RGB;
00275 jpeg_set_defaults(&cinfo);
00276 jpeg_set_quality(&cinfo, quality, TRUE);
00277
00278
00279 jpeg_start_compress(&cinfo, TRUE);
00280 row_stride = width * ncomp;
00281
00282 while (cinfo.next_scanline < cinfo.image_height) {
00283 row_pointer[0] = & buf[cinfo.next_scanline * row_stride];
00284 (void) jpeg_write_scanlines(&cinfo, row_pointer, 1);
00285 }
00286 jpeg_finish_compress(&cinfo);
00287
00288 fclose(jpeg_file);
00289 jpeg_destroy_compress(&cinfo);
00290 }
00291 #endif
00292
00293
00294
00295
00296
00297 static float *
00298 get_slice(const char *infile, int pos, const char *var,
00299 float coor, int *sizes)
00300 {
00301
00302 int c, b, select, is_float=0;
00303 int i1, i2, j1, j2, i, j, k, ii, jj;
00304 int tot_blocks, lrefmin, lrefmax;
00305 int Nx, Ny, N1, N2, Nz, N1b, N2b, N1r, N2r, nxb, nyb, nzb;
00306 float dx1, dx2, xm1, xm2;
00307 double xmin, xmax, ymin, ymax, zmin, zmax;
00308 float dx, dy, dz;
00309 int rebin_factor;
00310 int *lrefine, *nodetype;
00311 float *data_f, *slice=NULL;
00312 double *data_d, *bnd_box, *bnd_box_cp=NULL;
00313
00314 size_t *dims, *start, *stride, *count;
00315
00316
00317 rcComm_t *conn = NULL;
00318
00319
00320 H5File *h5file=NULL;
00321 H5Dataset *h5dset=NULL;
00322
00323 int ret_val = 0;
00324
00325 conn = rodsConnect();
00326 if (conn == NULL)
00327 PRINT_ERROR("Failed to make connection to iRODS server.");
00328
00329
00330 h5file = (H5File*)malloc(sizeof(H5File));
00331 H5File_ctor(h5file);
00332 h5file->filename = (char*)malloc(strlen(infile)+1);
00333 strcpy(h5file->filename, infile);
00334
00335
00336 h5file->opID = H5FILE_OP_OPEN;
00337 ret_val = h5ObjRequest(conn, h5file, H5OBJECT_FILE);
00338 if (ret_val < 0 || h5file->fid < 0)
00339 PRINT_ERROR ("H5FILE_OP_OPEN failed.");
00340
00341
00342 for (i=0; i<h5file->root->ndatasets; i++) {
00343 h5dset = (H5Dataset *) &(h5file->root->datasets[i]);
00344 if (strcmp(h5dset->fullpath+1, DNAME_REFINE_LEVEL) == 0)
00345 break;
00346 else
00347 h5dset = NULL;
00348 }
00349 if (h5dset == NULL)
00350 PRINT_ERROR("Failed to open dataset: REFINE_LEVEL.");
00351
00352 h5dset->opID = H5DATASET_OP_READ;
00353 ret_val = h5ObjRequest(conn, h5dset, H5OBJECT_DATASET);
00354 if (ret_val < 0)
00355 PRINT_ERROR("H5DATASET_OP_READ failed.");
00356
00357 tot_blocks = h5dset->space.dims[0];
00358 lrefine = (int *)h5dset->value;
00359
00360 #ifdef DEBUG
00361 printf("refine level:\n");
00362 for (i = 0; i < tot_blocks; i++)
00363 printf("%d, ", lrefine[i]);
00364 puts("\n");
00365 #endif
00366
00367
00368 for (i=0; i<h5file->root->ndatasets; i++) {
00369 h5dset = (H5Dataset *) &(h5file->root->datasets[i]);
00370 if (strcmp(h5dset->fullpath+1, DNAME_NODE_TYPE) == 0)
00371 break;
00372 else
00373 h5dset = NULL;
00374 }
00375 if (h5dset == NULL)
00376 PRINT_ERROR("Failed to open dataset: NODE_TYPE.");
00377
00378 h5dset->opID = H5DATASET_OP_READ;
00379 ret_val = h5ObjRequest(conn, h5dset, H5OBJECT_DATASET);
00380 if (ret_val < 0)
00381 PRINT_ERROR("H5DATASET_OP_READ failed.");
00382 nodetype = (int *)h5dset->value;
00383
00384 #ifdef DEBUG
00385 printf("node type:\n");
00386 for (i = 0; i < tot_blocks; i++)
00387 printf("%d, ", nodetype[i]);
00388 puts("\n");
00389 #endif
00390
00391
00392 for (i=0; i<h5file->root->ndatasets; i++) {
00393 h5dset = (H5Dataset *) &(h5file->root->datasets[i]);
00394 if (strcmp(h5dset->fullpath+1, DNAME_BOUNDING_BOX) == 0)
00395 break;
00396 else
00397 h5dset = NULL;
00398 }
00399 if (h5dset == NULL)
00400 PRINT_ERROR("Failed to open dataset: BOUNDING_BOX.");
00401
00402 h5dset->opID = H5DATASET_OP_READ;
00403 ret_val = h5ObjRequest(conn, h5dset, H5OBJECT_DATASET);
00404 if (ret_val < 0)
00405 PRINT_ERROR("H5DATASET_OP_READ failed.");
00406
00407 if (h5dset->type.size == 4) {
00408 is_float = 1;
00409 int npoints = h5dset->space.npoints;
00410 bnd_box = (double *)malloc(npoints*sizeof(double));
00411 float *bnd_box_f = (float *)h5dset->value;
00412
00413 for (i=0; i<npoints; i++)
00414 bnd_box[i] = (double)bnd_box_f[i];
00415 bnd_box_cp = bnd_box;
00416 } else
00417 bnd_box = (double *)h5dset->value;
00418
00419 xmin = 1.E99;
00420 xmax = -1.E99;
00421 ymin = 1.E99;
00422 ymax = -1.E99;
00423 zmin = 1.E99;
00424 zmax = -1.E99;
00425
00426 for (i = 0; i < 6*tot_blocks; i += 6) {
00427 if (bnd_box[i ] < xmin) { xmin = bnd_box[i ]; }
00428 if (bnd_box[i+1] > xmax) { xmax = bnd_box[i+1]; }
00429 if (bnd_box[i+2] < ymin) { ymin = bnd_box[i+2]; }
00430 if (bnd_box[i+3] > ymax) { ymax = bnd_box[i+3]; }
00431 if (bnd_box[i+4] < zmin) { zmin = bnd_box[i+4]; }
00432 if (bnd_box[i+5] > zmax) { zmax = bnd_box[i+5]; }
00433 }
00434
00435 #ifdef DEBUG
00436 printf("bounding box:\n");
00437 printf(" x = %e ... %e\n", xmin, xmax);
00438 printf(" y = %e ... %e\n", ymin, ymax);
00439 printf(" z = %e ... %e\n", zmin, zmax);
00440 #endif
00441
00442 lrefmin = 100000;
00443 lrefmax = 0;
00444
00445 for (i = 0; i < tot_blocks; i++) {
00446 if (lrefine[i] < lrefmin) { lrefmin = lrefine[i]; }
00447 if (lrefine[i] > lrefmax) { lrefmax = lrefine[i]; }
00448 }
00449
00450 #ifdef DEBUG
00451 printf("refinement levels: %d ... %d\n", lrefmin, lrefmax);
00452 #endif
00453
00454 for (i=0; i<h5file->root->ndatasets; i++) {
00455 h5dset = (H5Dataset *) &(h5file->root->datasets[i]);
00456 if (strncmp(h5dset->fullpath+1, var, strlen(var)) == 0)
00457 break;
00458 else
00459 h5dset = NULL;
00460 }
00461
00462 printf("\ndataset: %s\n", var);
00463
00464 if (h5dset == NULL)
00465 PRINT_ERROR("Failed to open dataset.");
00466
00467 dims = h5dset->space.dims;
00468 start = h5dset->space.start;
00469 stride = h5dset->space.stride;
00470 count = h5dset->space.count;
00471
00472 nxb = (int)dims[3];
00473 nyb = (int)dims[2];
00474 nzb = (int)dims[1];
00475
00476 start[0] = 0;
00477 start[1] = 0;
00478 start[2] = 0;
00479 start[3] = 0;
00480 stride[0] = 1;
00481 stride[1] = 1;
00482 stride[2] = 1;
00483 stride[3] = 1;
00484 count[0] = 1;
00485 count[1] = dims[1];
00486 count[2] = dims[2];
00487 count[3] = dims[3];
00488
00489 Nx = (int)dims[3];
00490 Ny = (int)dims[2];
00491 Nz = (int)dims[1];
00492 for (i = 0; i < lrefmax-1; i++) {
00493 Nx *= 2;
00494 Ny *= 2;
00495 Nz *= 2;
00496 }
00497
00498 if (pos == 0) {
00499 slice = (float *)malloc(Nx*Ny*sizeof(float));
00500 N1 = Nx;
00501 N2 = Ny;
00502 N1b = (int)dims[3];
00503 N2b = (int)dims[2];
00504 dx1 = (float) ((xmax - xmin) / ((float)N1));
00505 xm1 = (float)xmin;
00506 dx2 = (float) ((ymax - ymin) / ((float)N2));
00507 xm2 = (float)ymin;
00508 } else if (pos == 1) {
00509 slice = (float *)malloc(Nx*Nz*sizeof(float));
00510 N1 = Nx;
00511 N2 = Nz;
00512 N1b = (int)dims[3];
00513 N2b = (int)dims[1];
00514 dx1 = (float) ((xmax - xmin) / ((float)N1));
00515 xm1 = (float)xmin;
00516 dx2 = (float) ((zmax - zmin) / ((float)N2));
00517 xm2 = (float)zmin;
00518 } else {
00519 slice = (float *)malloc(Ny*Nz*sizeof(float));
00520 N1 = Ny;
00521 N2 = Nz;
00522 N1b = (int)dims[2];
00523 N2b = (int)dims[1];
00524 dx1 = (float) ((ymax - ymin) / ((float)N1));
00525 xm1 = (float)ymin;
00526 dx2 = (float) ((zmax - zmin) / ((float)N2));
00527 xm2 = (float)zmin;
00528 }
00529
00530 #ifdef DEBUG
00531 printf("slice grid size: %d x %d\n", N1, N2);
00532 printf("min zone spacing: %e x %e\n", dx1, dx2);
00533 #endif
00534
00535 c = 0;
00536 for (b = 0; b < tot_blocks; b++) {
00537 c = 6*b;
00538
00539 if (nodetype[b] != 1)
00540 continue;
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550 dx = (float)(bnd_box[c+1] - bnd_box[c ]) / ((float)dims[3]);
00551 dy = (float)(bnd_box[c+3] - bnd_box[c+2]) / ((float)dims[2]);
00552 dz = (float)(bnd_box[c+5] - bnd_box[c+4]) / ((float)dims[1]);
00553
00554 select = 0;
00555 switch (pos) {
00556 case 0:
00557 if ((bnd_box[c+4] <= coor) && (bnd_box[c+5] > coor)) select = 1;
00558 break;
00559 case 1:
00560 if ((bnd_box[c+2] <= coor) && (bnd_box[c+3] > coor)) select = 1;
00561 break;
00562 case 2:
00563 if ((bnd_box[c ] <= coor) && (bnd_box[c+1] > coor)) select = 1;
00564 break;
00565 }
00566
00567 if (select != 1)
00568 continue;
00569
00570 rebin_factor = 1;
00571 for (i = 0; i < lrefmax-lrefine[b]; i++) {
00572 rebin_factor *= 2;
00573 }
00574 N1r = N1b * rebin_factor;
00575 N2r = N2b * rebin_factor;
00576
00577 start[0] = b;
00578 h5dset->opID = H5DATASET_OP_READ;
00579 ret_val = h5ObjRequest(conn, h5dset, H5OBJECT_DATASET);
00580 if (ret_val < 0)
00581 PRINT_ERROR("H5DATASET_OP_READ failed.");
00582
00583 if (is_float)
00584 data_f = (float *)h5dset->value;
00585 else
00586 data_d = (double *)h5dset->value;
00587
00588 switch (pos) {
00589 case 0:
00590 i1 = (int)((bnd_box[c ] - xm1)/dx1);
00591 i2 = (int)((bnd_box[c+1] - xm1)/dx1);
00592 j1 = (int)((bnd_box[c+2] - xm2)/dx2);
00593 j2 = (int)((bnd_box[c+3] - xm2)/dx2);
00594 k = (int)((coor - bnd_box[c+4])/dz);
00595
00596 for (j = j1; j < j2; j++) {
00597 for (i = i1; i < i2; i++) {
00598 ii = (i-i1)/rebin_factor;
00599 jj = (j-j1)/rebin_factor;
00600 if (is_float)
00601 slice[i+j*N1] = data_f[ii+jj*N1b+k*N1b*N2b];
00602 else
00603 slice[i+j*N1] = (float)data_d[ii+jj*N1b+k*N1b*N2b];
00604 }
00605 }
00606 break;
00607 case 1:
00608 i1 = (int)((bnd_box[c ] - xm1)/dx1);
00609 i2 = (int)((bnd_box[c+1] - xm1)/dx1);
00610 j1 = (int)((bnd_box[c+4] - xm2)/dx2);
00611 j2 = (int)((bnd_box[c+5] - xm2)/dx2);
00612 k = (int)((coor - bnd_box[c+2])/dy);
00613 for (j = j1; j < j2; j++) {
00614 for (i = i1; i < i2; i++) {
00615 ii = (i-i1)/rebin_factor;
00616 jj = (j-j1)/rebin_factor;
00617 if (is_float)
00618 slice[i+j*N1] = data_f[ii+k*N1b+jj*N1b*N2b];
00619 else
00620 slice[i+j*N1] = (float)data_d[ii+k*N1b+jj*N1b*N2b];
00621 }
00622 }
00623 break;
00624 case 2:
00625 i1 = (int)((bnd_box[c+2] - xm1)/dx1);
00626 i2 = (int)((bnd_box[c+3] - xm1)/dx1);
00627 j1 = (int)((bnd_box[c+4] - xm2)/dx2);
00628 j2 = (int)((bnd_box[c+5] - xm2)/dx2);
00629 k = (int)((coor - bnd_box[c ])/dx);
00630 for (j = j1; j < j2; j++) {
00631 for (i = i1; i < i2; i++) {
00632 ii = (i-i1)/rebin_factor;
00633 jj = (j-j1)/rebin_factor;
00634 if (is_float)
00635 slice[i+j*N1] = data_f[k+ii*N1b+jj*N1b*N2b];
00636 else
00637 slice[i+j*N1] = (float)data_d[k+ii*N1b+jj*N1b*N2b];
00638 }
00639 }
00640 break;
00641 }
00642 }
00643
00644 sizes[0] = N1;
00645 sizes[1] = N2;
00646
00647 #ifdef DEBUG
00648 puts("slice data:");
00649 for (i=0; i<10; i++) {
00650 for (j=0; j<10; j++) {
00651 printf("%f, ", slice[N2*i+j]);
00652 }
00653 printf("\n");
00654 }
00655 puts("\n");
00656 #endif
00657
00658 done:
00659
00660
00661 if (h5file) {
00662 if (h5file->fid > 0) {
00663 h5file->opID = H5FILE_OP_CLOSE;
00664 h5ObjRequest(conn, h5file, H5OBJECT_FILE);
00665 }
00666 H5File_dtor(h5file);
00667 free(h5file);
00668 }
00669
00670 if (ret_val<0 && slice) {
00671 free (slice);
00672 slice = NULL;
00673 }
00674
00675 if (bnd_box_cp)
00676 free (bnd_box_cp);
00677
00678 rcDisconnect(conn);
00679
00680 return slice;
00681
00682 }
00683
00684
00685
00686
00687 static void
00688 write_float(const char *infile, char *outfile, int pos, const char *var,
00689 float coor, int dims[], float *slice)
00690 {
00691 FILE* outfilePtr;
00692
00693 if (strlen(outfile)<1)
00694 get_defualt_outfile_name (infile, outfile, pos, var, coor, dims, "out");
00695
00696 outfilePtr = fopen(outfile, "w");
00697 fwrite(slice, dims[0]*dims[1]*sizeof(float), 1, outfilePtr);
00698 fclose(outfilePtr);
00699 }
00700
00701 #ifdef JPEG
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718 static int
00719 read_palette(const char *palfile, unsigned char *pal)
00720 {
00721 int nentries=0, i, j, idx;
00722 float v, r, g, b, ratio, max_v, min_v, max_color, min_color;
00723 float tbl[NCOLOR][4];
00724 unsigned char tbl_char[NCOLOR][3];
00725 char line[LINE_LENGTH];
00726 FILE *pFile;
00727
00728 if (! palfile || !pal)
00729 return 0;
00730
00731 pFile = fopen (palfile,"r");
00732 if (pFile == NULL)
00733 return 0;
00734
00735
00736 idx = 0;
00737 while (fgets(line, LINE_LENGTH, pFile) != NULL) {
00738 if (sscanf(line, "%f %f %f %f", &v, &r, &g, &b) != 4)
00739 continue;
00740
00741 tbl[idx][0] = v;
00742 tbl[idx][1] = r;
00743 tbl[idx][2] = g;
00744 tbl[idx][3] = b;
00745
00746 if (idx==0) {
00747 max_v = min_v = v;
00748 max_color = min_color = r;
00749 }
00750
00751 max_v = fmaxf(max_v, v);
00752 min_v = fminf(min_v, v);
00753 max_color = fmaxf(max_color, r);
00754 max_color = fmaxf(max_color, g);
00755 max_color = fmaxf(max_color, b);
00756 min_color = fminf(min_color, r);
00757 min_color = fminf(min_color, g);
00758 min_color = fminf(min_color, b);
00759
00760 idx++;
00761 if (idx >= NCOLOR)
00762 break;
00763 }
00764
00765 fclose(pFile);
00766
00767 nentries = idx;
00768 if (max_color <= 1) {
00769 ratio = (min_color==max_color) ? 1.0f : (float)((NCOLOR-1.0)/(max_color-min_color));
00770
00771 for (i=0; i<nentries; i++) {
00772 for (j=1; j<4; j++)
00773 tbl[i][j] = (tbl[i][j]-min_color)*ratio;
00774 }
00775 }
00776
00777 for (i=0; i<NCOLOR; i++) {
00778 for (j=0; j<3; j++)
00779 tbl_char[i][j] = 0;
00780 }
00781
00782 idx = 0;
00783 memset(pal, 0, NCOLOR*3);
00784 ratio = (min_v==max_v) ? 1.0f : (float)((NCOLOR-1.0)/(max_v-min_v));
00785
00786 for (i=0; i<nentries; i++) {
00787 idx = (int) ((tbl[i][0]-min_v)*ratio);
00788 for (j=0; j<3; j++)
00789 tbl_char[idx][j] = (unsigned char) tbl[i][j+1];
00790 }
00791
00792 for (i=1; i<NCOLOR; i++) {
00793 if ( (tbl_char[i][0]+tbl_char[i][1]+tbl_char[i][2]) == 0 ) {
00794 j = i+1;
00795 while ( (tbl_char[j][0]+tbl_char[j][1]+tbl_char[j][2]) == 0 ) {
00796 j++;
00797 if (j>=NCOLOR)
00798 return 0;
00799 }
00800
00801 tbl_char[i][0] = (unsigned char) (tbl_char[i-1][0] + (float)(tbl_char[j][0]-tbl_char[i-1][0])/(j-i));
00802 tbl_char[i][1] = (unsigned char) (tbl_char[i-1][1] + (float)(tbl_char[j][1]-tbl_char[i-1][1])/(j-i));
00803 tbl_char[i][2] = (unsigned char) (tbl_char[i-1][2] + (float)(tbl_char[j][2]-tbl_char[i-1][2])/(j-i));
00804 }
00805
00806 pal[i*3] = tbl_char[i][0];
00807 pal[i*3+1] = tbl_char[i][1];
00808 pal[i*3+2] = tbl_char[i][2];
00809 }
00810
00811 return nentries;
00812 }
00813
00814
00815 static void
00816 write_jpeg(const char *infile, char *outfile, int pos, const char *var,
00817 float coor, int dims[], float *slice, const char *palfile)
00818 {
00819 float min, max, ratio;
00820 int i, idx, npoints;
00821 unsigned char pal[NCOLOR*3];
00822 unsigned char *idx_values=NULL;
00823 unsigned char *pixel_values=NULL;
00824 unsigned char is_log_scaling=0;
00825 unsigned char has_color_table=0;
00826
00827 if (dims == NULL || slice==NULL || (npoints = dims[0] * dims[1])<1 )
00828 goto done;
00829
00830 if (strlen(outfile)<1)
00831 get_defualt_outfile_name (infile, outfile, pos, var, coor, dims, "jpg");
00832
00833 if (read_palette(palfile, pal))
00834 has_color_table = 1;
00835
00836
00837 min = max = slice[0];
00838 for (i=0; i<npoints; i++) {
00839 min = fminf(min, slice[i]);
00840 max = fmaxf(max, slice[i]);
00841 }
00842
00843 if (max < min)
00844 max = (float) ((2.0 * fabs(min)) + 1.0);
00845
00846
00847 if (min>=0) {
00848 is_log_scaling = 1;
00849 min = log10f(min);
00850 max = log10f(max);
00851 }
00852
00853
00854 idx_values = (unsigned char*)malloc(npoints);
00855 pixel_values = (unsigned char*)malloc(npoints*3);
00856 ratio = (min==max) ? 1.0f : (float)((NCOLOR-1.0)/(max-min));
00857 for (i=0; i<npoints; i++) {
00858 if (is_log_scaling)
00859 idx_values[i] = idx = (unsigned char) ((log10f(slice[i])-min)*ratio);
00860 else
00861 idx_values[i] = idx = (unsigned char) ((slice[i]-min)*ratio);
00862
00863 if (has_color_table) {
00864 pixel_values[i*3] = pal[idx*3];
00865 pixel_values[i*3+1] = pal[idx*3+1];
00866 pixel_values[i*3+2] = pal[idx*3+2];
00867 } else {
00868
00869 pixel_values[i*3] = pixel_values[i*3+1] = pixel_values[i*3+2] = idx_values[i];
00870 }
00871 }
00872
00873 write_jpeg_file(outfile, (JSAMPLE *)pixel_values, dims[0], dims[1], 3, 100);
00874
00875 done:
00876 if (idx_values)
00877 free(idx_values);
00878
00879 if (pixel_values)
00880 free(pixel_values);
00881 }
00882 #endif
00883
00884 int main(int argc, char* argv[])
00885 {
00886
00887 char infile[NAME_LENGTH];
00888 char outfile[NAME_LENGTH];
00889 char jpgfile[NAME_LENGTH];
00890 char palfile[NAME_LENGTH];
00891 int i, pos=0, dims[2];
00892 char var[VAR_LENGTH];
00893 float coor=0.0, *slice=NULL;
00894 int ret_val = 0;
00895
00896
00897
00898
00899 if (argc < 2) {
00900 (void) fprintf(stderr, "Invalid command line arguments.\n");
00901 usage(argv[0]);
00902 exit(1);
00903 }
00904
00905 outfile[0] = infile[0] = jpgfile[0] = palfile[0] = '\0';
00906 strcpy(var, "dens");
00907 for (i=1; i<argc; i++) {
00908 if (strncmp(argv[i], "-h", 2) == 0) {
00909 help(argv[0]);
00910 exit(1);
00911 } else if (strncmp(argv[i], "-o", 2) == 0) {
00912 strcpy(outfile, argv[++i]);
00913 } else if (strncmp(argv[i], "-p", 2) == 0) {
00914 pos = atoi(argv[++i]);
00915 } else if (strncmp(argv[i], "-m", 2) == 0) {
00916 memset(var, 0, VAR_LENGTH);
00917 strcpy(var, argv[++i]);
00918 } else if (strncmp(argv[i], "-c", 2) == 0) {
00919 coor = (float)atof(argv[++i]);
00920 } else if (strncmp(argv[i], "-j", 2) == 0) {
00921 strcpy(jpgfile, argv[++i]);
00922 } else if (strncmp(argv[i], "-t", 2) == 0) {
00923 strcpy(palfile, argv[++i]);
00924 }
00925 }
00926 strcpy(infile, argv[--i]);
00927
00928 trim(infile, NAME_LENGTH);
00929 trim(var, 5);
00930
00931 if ( (slice=get_slice(infile, pos, var, coor, dims)) ) {
00932 write_float(infile, outfile, pos, var, coor, dims, slice);
00933
00934
00935 }
00936
00937 #ifdef JPEG
00938 if (slice)
00939 write_jpeg(infile, jpgfile, pos, var, coor, dims, slice, palfile);
00940 #endif
00941
00942
00943 if (slice)
00944 free (slice);
00945
00946 return ret_val;
00947 }
00948