00001
00002
00003
00004
00142 #include "mypng.h"
00143 #include "jsoc_main.h"
00144
00145 #define PNGDIE(msg) {fprintf(stderr,"%s\n",msg); \
00146 if (png_ptr) png_destroy_write_struct(&png_ptr, &info_ptr); \
00147 if (row_pointers) free(row_pointers); \
00148 if (fp) fclose(fp); \
00149 return(1);}
00150
00151 #define DIE(msg) {fflush(stdout);fprintf(stderr,"%s, status=%d\n",msg,status); return(status);}
00152
00153 #define MINMAX 1
00154 #define MAXMIN 2
00155 #define HISTEQ 3
00156 #define MINMAX99 4
00157 #define MAXMIN99 5
00158 #define MINMAXGIVEN 6
00159 #define MAG 7
00160 #define MAGMINMAXGIVEN 8
00161 #define LOG 9
00162 #define SQRT 10
00163
00164 #define NEAR_VERTICAL (0.1) // threshold to use image_cutout instead of upNcenter
00165
00166 char *module_name = "render_image";
00167
00168 ModuleArgs_t module_args[] =
00169 {
00170 {ARG_STRING, "in", "NOT SPECIFIED", "Input data series."},
00171 {ARG_INT, "n", "0", "Limit of number of records, 0 for no limin; >0 count from start; <0 count from end"},
00172 {ARG_STRING, "out", "NOT SPECIFIED", "Output data series or full path of directory or pipe via ppm program."},
00173 {ARG_STRING, "outname", "NOT SPECIFIED", "Output quantity name for filenames"},
00174 {ARG_STRING, "scaling", "MINMAX", "Color table type, minmax, minmax99, minmaxgiven, histeq, log, sqrt, ..."},
00175 {ARG_STRING, "pallette", "none", "Color table, GREY or filename"},
00176 {ARG_STRING, "palette", "GREY", "Color table, GREY or filename"},
00177 {ARG_STRING, "type", "png", "Image protocol, png, ppm, ..."},
00178 {ARG_STRING, "outid", "#", "output identifier, #=record counter, time=time as yyyymmdd_hhmmss"},
00179 {ARG_STRING, "tkey", "T_REC", "Time keyword name used if outid==time, defaults to T_REC"},
00180 {ARG_FLOAT, "min", "DRMS_MISSING_FLOAT", "Min value for scaling"},
00181 {ARG_FLOAT, "max", "DRMS_MISSING_FLOAT", "Max value for scaling"},
00182 {ARG_FLOAT, "bias", "30.0", "Max value for mag scaling"},
00183 {ARG_FLOAT, "prescale", "1.0", "Multiplication scaling prior to color scaling"},
00184 {ARG_INTS, "scale", "1", "Reduction factors"},
00185 {ARG_FLAG, "c", "0", "Crop flag, causes crop to RSUN_OBS"},
00186 {ARG_FLAG, "r", "1", "Rotate to solar north up and reposition to image center"},
00187 {ARG_FLAG, "u", "0", "use unchanged, is-is for rotation and centering"},
00188 {ARG_FLAG, "w", "0", "use white for missing pixels, instead of black"},
00189 {ARG_FLAG, "x", "0", "High-quality flag, sets bytespercolor to 2 instead of 1"},
00190 {ARG_END}
00191 };
00192
00193 #define Deg2Rad (M_PI/180.0)
00194 #define Rad2arcsec (3600.0/Deg2Rad)
00195 #define arcsec2Rad (Deg2Rad/3600.0)
00196 #define Rad2Deg (180.0/M_PI)
00197
00198 struct ObsInfo_struct
00199 {
00200
00201 TIME t_obs;
00202 double rsun_obs, obs_vr, obs_vw, obs_vn;
00203 double crpix1, crpix2, cdelt1, cdelt2, crota2;
00204 double crval1, crval2;
00205 double cosa, sina;
00206 double obs_b0;
00207
00208 int i,j;
00209
00210 double x,y,r;
00211 double rho;
00212 double lon;
00213 double lat;
00214 double sinlat, coslat;
00215 double sig;
00216 double mu;
00217 double chi;
00218 double obs_v;
00219 };
00220
00221 typedef struct ObsInfo_struct ObsInfo_t;
00222
00223 void rebinArraySF(DRMS_Array_t *out, DRMS_Array_t *in);
00224
00225 int add_png(char *imgPath, DRMS_Array_t *imgArray, int scaletype, int usewhite, char *palette, int colors,
00226 int bytespercolor, double *minp, double *maxp, double prescale);
00227
00228 int add_ppm(char *imgPath, DRMS_Array_t *imgArray, int scaletype, int usewhite, char *palette, int colors,
00229 int bytespercolor, double *minp, double *maxp, double prescale);
00230
00231 int make_png(char *imgPath, unsigned char *data, int height, int width, char *palette, int bytepercolor, int colors);
00232
00233 char *set_scaling(DRMS_Array_t *in, double *minp, double *maxp, int *nmissp, char *palette, int colors,
00234 int bytepercolor, int scaling, int usewhite, double prescale);
00235
00236 int upNcenter(DRMS_Array_t *arr, ObsInfo_t *ObsLoc);
00237 int reposition(DRMS_Array_t *arr, ObsInfo_t *ObsLoc);
00238 int crop_image(DRMS_Array_t *arr, ObsInfo_t *ObsLoc);
00239
00240 ObsInfo_t *GetObsInfo(DRMS_Segment_t *seg, ObsInfo_t *pObsLoc, int *rstatus);
00241
00242 int image_cutout(void *inarr, int nx, int ny, int in_type, float crota2, float mag, float dx, float dy,
00243 void **outarr, int wide, int high, float xc, float yc, int regridtype, int do_stretch);
00244
00245 #define PNG 1
00246 #define PPM 2
00247 #define PPMPIPE 3
00248
00249 int DoIt(void)
00250 {
00251 int status = DRMS_SUCCESS;
00252 DRMS_RecordSet_t *inRS, *outRS = NULL;
00253 DRMS_Record_t *inRec, *outRec;
00254 int irec, nrecs;
00255 int fileonly=0;
00256 int imgtype;
00257 int quality;
00258
00259 char *inQuery = (char *)params_get_str(&cmdparams, "in");
00260 char *outQuery = (char *)params_get_str(&cmdparams, "out");
00261 char *outName = (char *)params_get_str(&cmdparams, "outname");
00262 char *outid = (char *)params_get_str(&cmdparams, "outid");
00263 char *tkey = (char *)params_get_str(&cmdparams, "tkey");
00264 char *palette = (char *)params_get_str(&cmdparams, "palette");
00265 if (strcmp((char *)params_get_str(&cmdparams, "pallette"), "none"))
00266 palette = (char *)params_get_str(&cmdparams, "pallette");
00267 int limit = params_get_int(&cmdparams, "n");
00268 int bytespercolor = (params_isflagset(&cmdparams, "x") ? 2 : 1);
00269 int missingwhite = params_isflagset(&cmdparams, "w");
00270 int crop = params_isflagset(&cmdparams, "c");
00271 int rotate = params_isflagset(&cmdparams, "r");
00272 int asis = params_isflagset(&cmdparams, "u");
00273 char *type = (char *)params_get_str(&cmdparams, "type");
00274 char *scaling = (char *)params_get_str(&cmdparams, "scaling");
00275 double called_min = params_get_double(&cmdparams, "min");
00276 double called_max = params_get_double(&cmdparams, "max");
00277 double prescale = params_get_double(&cmdparams, "prescale");
00278 double min, max;
00279 int *scales = NULL;
00280 int iscale, nscales;
00281 int colors;
00282 int scaletype;
00283 int srcNx = 0, srcNy = 0;
00284 char *inSegName = NULL, *segnamep;
00285
00286 if (strcasecmp(palette, "grey") == 0 || strcasecmp(palette, "gray") == 0)
00287 {
00288 palette = "GREY";
00289 colors = 1;
00290 }
00291 else
00292 colors = 3;
00293
00294 if (strcasecmp(type, "png")==0) imgtype = PNG;
00295 else if (strcasecmp(type, "ppm")==0) imgtype = PPM;
00296 else imgtype = PPMPIPE;
00297
00298 nscales = cmdparams_get_intarr(&cmdparams, "scale", &scales, &status);
00299 if (status)
00300 DIE("scale parameter not found");
00301
00302
00303
00304
00305
00306 if (strcasecmp(scaling, "minmax") == 0) scaletype = MINMAX;
00307 else if (strcasecmp(scaling, "maxmin") == 0) scaletype = MAXMIN;
00308 else if (strcasecmp(scaling, "minmax") == 0) scaletype = MINMAX99;
00309 else if (strcasecmp(scaling, "maxmin99") == 0) scaletype = MAXMIN99;
00310 else if (strcasecmp(scaling, "histeq") == 0) scaletype = HISTEQ;
00311 else if (strcasecmp(scaling, "minmaxgiven") == 0) scaletype = MINMAXGIVEN;
00312 else if (strncasecmp(scaling, "mag", 3) == 0) scaletype = MAG;
00313 else if (strncasecmp(scaling, "log", 3) == 0) scaletype = LOG;
00314 else if (strncasecmp(scaling, "sqrt", 4) == 0) scaletype = SQRT;
00315
00316 else scaletype = 0;
00317 fprintf(stdout, " using scaling %d\n", scaletype);
00318 if (isnan(called_min) || isnan(called_max))
00319 {
00320 if (scaletype == MINMAXGIVEN)
00321 DIE("min and max must be specified for type MINMAXGIVEN\n");
00322 }
00323 else
00324 {
00325 if (scaletype == MINMAX)
00326 scaletype = MINMAXGIVEN;
00327 else if (scaletype == MAG)
00328 scaletype = MAGMINMAXGIVEN;
00329 }
00330
00331 if (limit == 0)
00332 inRS = drms_open_records(drms_env, inQuery, &status);
00333 else
00334 inRS = drms_open_nrecords(drms_env, inQuery, limit, &status);
00335 if (status || inRS->n == 0)
00336 DIE("No input data found");
00337 drms_stage_records(inRS, 1, 1);
00338 nrecs = inRS->n;
00339 fprintf(stderr,"nrecs=%d\n",nrecs);
00340 segnamep = index(inQuery, '{');
00341 if (segnamep)
00342 {
00343 segnamep++;
00344 inSegName = strdup(strtok(segnamep, ",}"));
00345 }
00346 else
00347 {
00348 DRMS_Segment_t *seg0 = drms_segment_lookupindex(inRS->records[0], 0, 0);
00349
00350 if (seg0)
00351 {
00352 inSegName = strdup(seg0->info->name);
00353 }
00354 else
00355 {
00356 DIE("Series has no segments.");
00357 }
00358 }
00359
00360 if (!outQuery[0])
00361 DIE("No valid output place specified");
00362 if (outQuery[0] == '/' || (outQuery[0] == '.' && outQuery[1] && outQuery[1] == '/') || outQuery[0] == '|')
00363 {
00364
00365 fileonly = 1;
00366 }
00367 else
00368 {
00369
00370 fileonly = 0;
00371 if (strcmp(inRS->records[0]->seriesinfo->seriesname, outQuery) == 0)
00372 outRS = drms_clone_records(inRS, DRMS_PERMANENT, DRMS_SHARE_SEGMENTS, &status);
00373 else
00374 outRS = drms_create_records(drms_env, nrecs, outQuery, DRMS_PERMANENT, &status);
00375 if (status)
00376 DIE("Output recordset not created");
00377 }
00378
00379 for (irec=0; irec<nrecs; irec++)
00380 {
00381 int status;
00382 DRMS_Array_t *srcArray;
00383 DRMS_Segment_t *srcSeg;
00384 ObsInfo_t *ObsLoc;
00385
00386 char imageID[100];
00387 inRec = inRS->records[irec];
00388 quality = drms_getkey_int(inRec, "QUALITY", &status);
00389 if (!status && quality < 0)
00390 {
00391 fprintf(stderr,"Quality bad for rec %d\n", irec);
00392 continue;
00393 }
00394 else
00395 {
00396 if (inSegName)
00397 {
00398 srcSeg = drms_segment_lookup(inRec, inSegName);
00399 if (srcSeg)
00400 {
00401 srcArray = drms_segment_read(srcSeg, DRMS_TYPE_FLOAT, &status);
00402 }
00403 else
00404 {
00405 DIE("Failure looking up segment by name.")
00406 }
00407 }
00408 else
00409 {
00410 DIE("Could not obtain segment name.");
00411 }
00412
00413 if (!srcArray || status)
00414 {
00415 fprintf(stderr," No data file found rec %d, status=%d\n", irec,status);
00416 if (srcArray)
00417 {
00418 drms_free_array(srcArray);
00419 srcArray = NULL;
00420 }
00421 continue;
00422 }
00423 ObsLoc = GetObsInfo(srcSeg, NULL, &status);
00424 if (ObsLoc)
00425 {
00426 if (!asis)
00427 {
00428 if (rotate)
00429 reposition(srcArray, ObsLoc);
00430 else
00431 upNcenter(srcArray, ObsLoc);
00432 }
00433 if (crop)
00434 crop_image(srcArray, ObsLoc);
00435 }
00436 srcNx = srcArray->axis[0];
00437 srcNy = srcArray->axis[1];
00438 }
00439 if (strcmp(outid,"#") == 0)
00440 sprintf(imageID, "%04d", irec);
00441 else if (strncasecmp(outid,"none",4) == 0)
00442 imageID[0] = '\0';
00443 else if (strncasecmp(outid,"time",4) == 0)
00444 {
00445 int timestatus;
00446 char timebuf[100];
00447 static int t[] = {0,1,2,3,5,6,8,9,10,11,12,14,15,17,18};
00448 int i;
00449 int maxtplace = 15;
00450 char *tplace = index(outid,':');
00451 TIME now = drms_getkey_time(inRec, tkey, ×tatus);
00452 if (timestatus)
00453 {
00454 fprintf(stderr,"Time key given in 'tkey' param, %s, is not found in the input series.\n",tkey);
00455 return(timestatus);
00456 }
00457 sprint_time(timebuf, now, "TAI", 0);
00458 if (tplace)
00459 sscanf(tplace+1,"%d",&maxtplace);
00460 for (i=0; i<maxtplace; i++)
00461 imageID[i] = timebuf[t[i]];
00462 imageID[i] = '\0';
00463 }
00464 else
00465 {
00466 fprintf(stderr, "Invalid outid param: |%s|\n", outid);
00467 return(1);
00468 }
00469 for (iscale=0; iscale<nscales; iscale++)
00470 {
00471 int dimxy;
00472 char dimxytxt[50];
00473 DRMS_Array_t *imageArray = NULL;
00474 char imgPath[DRMS_MAXPATHLEN];
00475 int scale = scales[iscale];
00476 int imageDims[2];
00477 int tmparray = 0;
00478 char imageFileName[DRMS_MAXSEGFILENAME];
00479 min = called_min;
00480 max = called_max;
00481 imageDims[0] = srcNx/scale;
00482 imageDims[1] = srcNy/scale;
00483 dimxy = imageDims[0];
00484 sprintf(dimxytxt,"%d",dimxy);
00485
00486
00487 if (scale != 1)
00488 {
00489 imageArray = drms_array_create(DRMS_TYPE_FLOAT, 2, imageDims, NULL, &status);
00490 rebinArraySF(imageArray, srcArray);
00491 tmparray = 1;
00492 }
00493 else
00494 {
00495 imageArray = srcArray;
00496 tmparray = 0;
00497 }
00498
00499 if (imageID[0] == '\0')
00500 sprintf(imageFileName, "%s", outName);
00501 else
00502 sprintf(imageFileName, "%s%s%s_%s.%s", imageID, (strcmp(outid,"#")==0 ? "." : "_"), outName,
00503 (dimxy == 4096 ? "4k" :
00504 (dimxy == 2048 ? "2k" :
00505 (dimxy == 1024 ? "1k" :
00506 (dimxy == 512 ? "512" :
00507 (dimxy == 256 ? "256" :
00508 (dimxy == 128 ? "128" : dimxytxt)))))), type);
00509
00510 if (fileonly)
00511 {
00512 if (outQuery[0] == '|')
00513 {
00514 char buf[4096];
00515 char *pbuf = buf, *c = outQuery;
00516 while (*c && (pbuf-buf) < sizeof(buf))
00517 {
00518 if (*c == '{')
00519 {
00520 c++;
00521 if (*c == '%')
00522 {
00523 float x;
00524 int ix, minx;
00525 c++;
00526 x = strtof(c, &c);
00527 ix = round(0.01*x*imageDims[1]);
00528 if (*c == ':')
00529 {
00530 c++;
00531 minx = (int)strtol(c, &c, 10);
00532 }
00533 else minx = 0;
00534 pbuf += sprintf(pbuf,"%d",(ix < minx ? minx : ix));
00535 }
00536 else if (strncmp(c, "ID", 2) == 0)
00537 {
00538 pbuf += sprintf(pbuf, "%s", imageID);
00539 c += 2;
00540 }
00541 if (*c++ != '}')
00542 DIE("Unrecognized token in pipe string");
00543 }
00544 else
00545 *pbuf++ = *c++;
00546 }
00547 *pbuf = '\0';
00548 sprintf(imgPath, "%s > %s", buf, imageFileName);
00549 }
00550 else
00551 sprintf(imgPath, "%s/%s", outQuery, imageFileName);
00552 }
00553 else
00554 {
00555 DRMS_Segment_t *imageSeg;
00556 char recdir[DRMS_MAXPATHLEN];
00557 outRec = outRS->records[irec];
00558 drms_copykeys(outRec, inRec, 0, kDRMS_KeyClass_Explicit);
00559 drms_record_directory(outRec, recdir, 0);
00560 imageSeg = drms_segment_lookup(outRec, "image");
00561 strcpy(imageSeg->filename, imageFileName);
00562 drms_keyword_setdate(outRec);
00563 sprintf(imgPath, "%s/%s", recdir, imageFileName);
00564 }
00565
00566 if (imgtype == PNG)
00567 add_png(imgPath, imageArray, scaletype, missingwhite, palette, colors, bytespercolor, &min, &max, prescale);
00568 else if (imgtype == PPM || imgtype == PPMPIPE)
00569 add_ppm(imgPath, imageArray, scaletype, missingwhite, palette, colors, bytespercolor, &min, &max, prescale);
00570
00571 if (tmparray)
00572 drms_free_array(imageArray);
00573 }
00574 drms_free_array(srcArray);
00575 }
00576 fprintf(stderr,"render_image done, irec=%d\n",irec);
00577 drms_close_records(inRS, DRMS_FREE_RECORD);
00578 fprintf(stderr,"render_image close in records done\n");
00579 if (!fileonly && outRS)
00580 drms_close_records(outRS, DRMS_INSERT_RECORD);
00581 fprintf(stderr,"leaving module\n");
00582 return (DRMS_SUCCESS);
00583 }
00584
00585
00586
00587
00588
00589
00590 int add_ppm(char *imgPath, DRMS_Array_t *imageArray, int scaletype, int usewhite, char *palette, int colors, int bytespercolor, double *minp, double *maxp, double prescale)
00591 {
00592 int i, Nx, Ny;
00593 unsigned char *imgData = NULL;
00594 int nmissing;
00595 int maxpixval = (bytespercolor == 1 ? 255 : 65535);
00596 FILE *out;
00597 char *ppmcode;
00598
00599 imgData = (unsigned char*)set_scaling(imageArray, minp, maxp, &nmissing, palette, colors, bytespercolor, scaletype, usewhite, prescale);
00600 if (!imgData)
00601 {
00602 fprintf(stderr,"scaling failed\n");
00603 return(1);
00604 }
00605
00606 ppmcode = (colors == 1 ? "P5" : "P6");
00607
00608 Nx = imageArray->axis[0];
00609 Ny = imageArray->axis[1];
00610
00611 if (imgPath[0] == '|')
00612 out = popen(imgPath+1, "w");
00613 else
00614 out = fopen(imgPath, "w");
00615
00616 fprintf(out, "%s\n%d %d %d\n", ppmcode, Nx, Ny, maxpixval);
00617 for (i=Ny-1; i >=0; i--)
00618 fwrite(imgData + i*Nx*bytespercolor*colors, bytespercolor, colors*Nx, out);
00619 if (imgPath[0] == '|')
00620 pclose(out);
00621 else
00622 fclose(out);
00623
00624 free(imgData);
00625 return(0);
00626 }
00627
00628 int add_png(char *imgPath, DRMS_Array_t *imageArray, int scaletype, int usewhite, char *palette, int colors, int bytespercolor, double *minp, double *maxp, double prescale)
00629 {
00630 int srcNx, srcNy;
00631 int status=0;
00632 unsigned char *imgData = NULL;
00633 int nmissing;
00634
00635 imgData = (unsigned char*)set_scaling(imageArray, minp, maxp, &nmissing, palette, colors, bytespercolor, scaletype, usewhite, prescale);
00636
00637 if (!imgData)
00638 {
00639 fprintf(stderr,"scaling failed, status=%d\n",status);
00640 return(status);
00641 }
00642
00643
00644
00645 if (make_png(imgPath, imgData, imageArray->axis[0], imageArray->axis[1], palette, bytespercolor, colors) != 0)
00646 {
00647 fprintf(stderr,"png failed, status=%d\n",status);
00648 free(imgData);
00649 return(status);
00650 }
00651 free(imgData);
00652 return(0);
00653 }
00654
00655
00656 void rebinArraySF(DRMS_Array_t *out, DRMS_Array_t *in)
00657 {
00658 float *inData = in->data;
00659 float *outData = out->data;
00660 float *inp, *outp;
00661 int inNx=in->axis[0], inNy=in->axis[1];
00662 int outNx=out->axis[0], outNy=out->axis[1];
00663 int inI, inJ, outI, outJ, i;
00664 int binsize = inNx/outNx;
00665 double inRow[inNx], outRow[outNx];
00666 int inRowN[inNx], outRowN[outNx];
00667
00668 for (outJ=0; outJ < outNy; outJ++)
00669 {
00670 int inRow0 = outJ * binsize;
00671 for (outI=0; outI<outNx; outI++)
00672 {
00673 outRowN[outI] = 0;
00674 outRow[outI] = 0.0;
00675 }
00676 for (inJ = inRow0; inJ < inRow0+binsize; inJ++)
00677 {
00678 inp = inData + inJ*inNx;
00679 for (outI=0; outI<outNx; outI++)
00680 for (i=0; i<binsize; i++, inp++)
00681 if (*inp != DRMS_MISSING_FLOAT)
00682 {
00683 outRow[outI] += *inp;
00684 outRowN[outI]++;
00685 }
00686 }
00687 for (outI=0; outI < outNx; outI++)
00688 outData[outJ*outNx + outI] = (outRowN[outI] ?
00689 (outRow[outI]/outRowN[outI]) : DRMS_MISSING_FLOAT);
00690 }
00691 }
00692
00693
00694
00695 unsigned short *init_color_table(char * palette, int colors, int bytepercolor)
00696 {
00697 int c, i;
00698 int tablesize = (bytepercolor == 1 ? 256 : 65536);
00699 int maxcolor = tablesize - 2;
00700 int lastcolor = tablesize - 1;
00701 int mincolor = 1;
00702 unsigned short *table = (unsigned short *)malloc(colors*tablesize*sizeof(unsigned short));
00703
00704 if (!table) return NULL;
00705 if (strcmp(palette,"GREY") == 0)
00706 for (i=0; i<tablesize; i++)
00707 for (c = 0; c<colors; c++)
00708 table[i*colors + c] = i;
00709
00710 else if (colors == 3 && access(palette, R_OK) == 0)
00711 {
00712 FILE *ct = fopen(palette,"r");
00713 void proc_color(FILE *ct, unsigned short *t, int max);
00714 char buf[256];
00715 char *dot = rindex(palette, '.');
00716 if (dot && strcmp(dot, ".lut")==0)
00717 {
00718 float fr, fg, fb;
00719 table[0] = table[1] = table[2] = 0;
00720 table[colors*lastcolor] = table[colors*lastcolor+1] = table[colors*lastcolor+2] = lastcolor;
00721 for (i=1; i<lastcolor; i++)
00722 {
00723 fscanf(ct, "%f%f%f", &fr, &fg, &fb);
00724 c = 0;
00725 table[i*colors + c++] = maxcolor*fr;
00726 table[i*colors + c++] = maxcolor*fg;
00727 table[i*colors + c++] = maxcolor*fb;
00728 }
00729 }
00730 else
00731 {
00732
00733 unsigned short r[tablesize], g[tablesize], b[tablesize];
00734 for (i=0; i<=tablesize; i++)
00735 r[i] = g[i] = b[i] = 0;
00736 while (fgets(buf,256,ct))
00737 {
00738 if (strncmp(buf,"RED:",4)==0)
00739 proc_color(ct,r,maxcolor);
00740 if (strncmp(buf,"GREEN:",6)==0)
00741 proc_color(ct,g,maxcolor);
00742 if (strncmp(buf,"BLUE:",5)==0)
00743 proc_color(ct,b,maxcolor);
00744 }
00745 table[0] = table[1] = table[2] = 0;
00746 table[colors*lastcolor] = table[colors*lastcolor+1] = table[colors*lastcolor+2] = lastcolor;
00747 for (i=1; i<lastcolor; i++)
00748 {
00749 c = 0;
00750 table[i*colors + c++] = r[i];
00751 table[i*colors + c++] = g[i];
00752 table[i*colors + c++] = b[i];
00753 }
00754 }
00755 fclose(ct);
00756 }
00757 return(table);
00758 }
00759
00760 void proc_color(FILE *ct, unsigned short *t, int max)
00761 {
00762 float cut, val;
00763 int pos = -1;
00764 double prev = 0.0;
00765 int prevloc = 0;
00766 double slope;
00767 pos = -1;
00768 prev = 0.0;
00769 prevloc = 0;
00770 while (fscanf(ct,"(%f,%f)", &cut, &val) == 2)
00771 {
00772 int loc = max*cut + 0.5;
00773 val *= max;
00774
00775 if (pos < 0)
00776 {
00777 pos = 0;
00778 prev = val;
00779 prevloc = 0;
00780 }
00781 if (loc == prevloc)
00782 slope = 0.0;
00783 else
00784 slope = (val - prev)/(loc - prevloc);
00785 while (pos <= loc)
00786 {
00787 t[pos] = slope * (pos-prevloc) + prev;
00788 pos++;
00789 }
00790 prev = t[loc];
00791 prevloc = loc;
00792 }
00793 while (pos <= max)
00794 t[pos++] = prev;
00795 }
00796
00797
00798 char *set_scaling(DRMS_Array_t *in, double *minp, double *maxp,
00799 int *nmissp, char *palette, int colors, int bytepercolor, int scaling, int usewhite, double prescale)
00800 {
00801 int idata, ndata, nmiss = 0;
00802 float *inData = in->data;
00803 ndata = in->axis[0] * in->axis[1];
00804 int tablesize = (bytepercolor == 1 ? 256 : 65536);
00805 int color;
00806 int grey = strcmp(palette, "GREY") == 0;
00807 static unsigned short *table = NULL;
00808 static unsigned short *greytable = NULL;
00809 static unsigned short *colortable = NULL;
00810 char *out;
00811 out = (char *)malloc(ndata * colors * bytepercolor * sizeof(char));
00812 int maxcolor = (bytepercolor == 1 ? 254 : 65534);
00813 int missingcolor = (usewhite ? maxcolor + 1 : 0);
00814
00815 fprintf(stderr,"set_scaling called, scaling=%d, palette=%s, min=%f, max=%f, prescale=%f",scaling,palette,*minp,*maxp,prescale);
00816 fprintf(stderr,"\n missingcolor %d maxcolor %d \n", missingcolor, maxcolor);
00817 if (!out)
00818 return(NULL);
00819
00820
00821
00822 if (grey && scaling == MINMAX && colors == 1)
00823 {
00824 if (!greytable)
00825 {
00826 int i;
00827 greytable = (unsigned short *)malloc(65536 * sizeof(unsigned short));
00828 for (i=0; i<65536; i++)
00829 greytable[i] = i;
00830 colortable = greytable;
00831 }
00832 }
00833 else
00834 {
00835 if (!table)
00836 table = init_color_table(palette, colors, bytepercolor);
00837 colortable = table;
00838 }
00839
00840 switch(scaling)
00841 {
00842 case MAG:
00843 case MAGMINMAXGIVEN:
00844 {
00845 float bias = params_get_double(&cmdparams, "bias");
00846 for (idata=0; idata<ndata; idata++)
00847 {
00848 float val = inData[idata];
00849 if (!isnan(val))
00850 val *= prescale;
00851
00852 inData[idata] = 400.0*val/pow((fabs(val)+bias), 0.75);
00853 }
00854 if (isnan(*minp)) *minp = -1500.0;
00855 if (isnan(*maxp)) *maxp = -*minp;
00856
00857 }
00858 case MINMAX:
00859 case MAXMIN:
00860 case MINMAXGIVEN:
00861 case LOG:
00862 case SQRT:
00863 {
00864 double min;
00865 double max;
00866 if (bytepercolor != 1 && bytepercolor != 2)
00867 {fprintf(stderr,"wrong bytes per color\n"); return(NULL);}
00868 if (scaling == MINMAXGIVEN || scaling == MAG || scaling == MAGMINMAXGIVEN)
00869 {
00870 min = *minp;
00871 max = *maxp;
00872
00873 }
00874 else
00875 {
00876 if (isnan(*minp) || isnan(*maxp))
00877 {
00878 int n=0;
00879 min = 1.0e20;
00880 max = -min;
00881
00882 for (idata=0; idata<ndata; idata++)
00883 {
00884 float val = inData[idata];
00885 if (!isnan(val))
00886 {
00887 val *= prescale;
00888 n++;
00889 if (val > max) max = val;
00890 if (val < min) min = val;
00891 }
00892 }
00893 if (n==0)
00894 {
00895 fprintf(stderr,"set_scaling, min or max is missing and no data values found.\n");
00896 return(NULL);
00897 }
00898 }
00899 if (isnan(*minp)) *minp = min;
00900 if (isnan(*maxp)) *maxp = max;
00901 min = *minp;
00902 max = *maxp;
00903 }
00904 double use_min, use_max;
00905 if (scaling == LOG)
00906 {
00907 if (min <= 0)
00908 min = 1.0e-10;
00909 if (max <= 0)
00910 max = 1.0e-10;
00911 use_min = log10(min);
00912 use_max = log10(max);
00913 }
00914 else if (scaling == SQRT)
00915 {
00916 if (min < 0)
00917 min = 0;
00918 if (max < 0)
00919 max = 0;
00920 use_min = sqrt(min);
00921 use_max = sqrt(max);
00922 }
00923 for (idata=0; idata<ndata; idata++)
00924 {
00925
00926 float val = inData[idata];
00927 int newval;
00928 if (!isnan(val))
00929 {
00930 val *= prescale;
00931 if (val <= min) val = min;
00932 if (val >= max) val = max;
00933 if (scaling == LOG)
00934 newval = (maxcolor-1)*((log10(val)-use_min)/(use_max-use_min)) + 0.5;
00935 else if (scaling == SQRT )
00936 newval = (maxcolor-1)*((sqrt(val)-use_min)/(use_max-use_min)) + 0.5;
00937 else
00938 {
00939 newval = (maxcolor-1)*(val-min)/(max-min) + 0.5;
00940 if (scaling == MAXMIN)
00941 newval = maxcolor - newval;
00942 }
00943 if (newval >= maxcolor) newval = maxcolor;
00944 if (newval < 1) newval = 1;
00945 }
00946 else
00947 {
00948 newval = missingcolor;
00949 nmiss += 1;
00950 }
00951 for (color=0; color<colors; color++)
00952 {
00953 if (bytepercolor == 1)
00954 *((unsigned char *)out + colors*idata + color) = colortable[newval*colors + color];
00955 else
00956 *((unsigned short *)out + colors*idata + color) = colortable[newval*colors + color];
00957 }
00958 }
00959 *maxp = max;
00960 *minp = min;
00961 *nmissp = nmiss;
00962 fprintf(stderr," min %f max %f nmissp %d ndata %d\n", *minp, *maxp, nmiss, ndata);
00963 break;
00964 }
00965
00966 case HISTEQ:
00967 {
00968 int ihist, nhist=65536;
00969 int *hist = (int *)malloc(nhist * sizeof(int));
00970 int icolor, ncolor = (bytepercolor == 1 ? 255 : 65535);
00971 int total, want;
00972 unsigned int newval;
00973 double valspercolor;
00974 char *outa;
00975 outa = set_scaling(in, minp, maxp, nmissp, "GREY", 1, 2, MINMAX, 0, prescale);
00976 valspercolor = (ndata - *nmissp)/(double)ncolor;
00977 for (ihist = 0; ihist < nhist; ihist++)
00978 hist[ihist] = 0;
00979 for (idata=0; idata<ndata; idata++)
00980 hist[*((unsigned short *)outa+idata)] += 1;
00981 for (icolor=0, total=0, ihist=0; icolor<=ncolor; icolor++)
00982 {
00983 want = icolor * valspercolor;
00984 for ( ; ihist<nhist; ihist++)
00985 {
00986 if (total >= want)
00987 break;
00988 total += hist[ihist];
00989 hist[ihist] = icolor;
00990 }
00991 }
00992 for ( ; ihist<nhist; ihist++)
00993 hist[ihist] = ncolor;
00994 for (idata=0; idata<ndata; idata++)
00995 {
00996 if (!isnan(inData[idata]))
00997 {
00998 newval = hist[*((unsigned short *)outa+idata)];
00999 if (newval > maxcolor) newval = maxcolor;
01000 if (newval < 1) newval = 1;
01001 }
01002 else
01003 newval = missingcolor;
01004 for (color=0; color<colors; color++)
01005 {
01006 if (bytepercolor == 1)
01007 *((unsigned char *)out + colors*idata + color) = colortable[newval*colors + color];
01008 else
01009 *((unsigned short *)out + colors*idata + color) = colortable[newval*colors + color];
01010 }
01011 }
01012
01013 free(outa);
01014 free(hist);
01015 break;
01016 }
01017
01018 case MINMAX99:
01019 case MAXMIN99:
01020 {
01021 fprintf(stderr,"DOING set_scaling case MAXMIN99 \n");
01022 int ihist, nhist=65536;
01023 int *hist = (int *)malloc(nhist * sizeof(int));
01024 int icolor, ncolor = (bytepercolor == 1 ? 255 : 65535);
01025 int total, want;
01026 unsigned int newval;
01027 double newmin, newmax;
01028 char *outa;
01029 fprintf(stderr," CALLING set_scaling case MINMAX %d \n", MINMAX);
01030 *minp = 0;
01031 *maxp = 0;
01032 outa = set_scaling(in, minp, maxp, nmissp, "GREY", 1, 2, MINMAX, 0, prescale);
01033
01034 fprintf(stderr,"CONTINUING: set_scaling, case MAXMIN99 %f %f \n", *minp, *maxp);
01035
01036 for (ihist=0; ihist<nhist; ihist++)
01037 hist[ihist] = 0;
01038 for (idata=0; idata<ndata; idata++)
01039 {
01040 hist[*((unsigned short *)outa+idata)] += 1;
01041 }
01042
01043 want = ndata/600;
01044 fprintf(stderr,"ndata %d want %d nhist %d \n", ndata, want, nhist);
01045 total = 0;
01046 for (ihist=0; ihist<nhist; total += hist[ihist++])
01047 if (total >= want)
01048 break;
01049 newmin = *minp + ((double)ihist/(double)nhist) * (*maxp - *minp);
01050
01051 total = 0;
01052 for (ihist=nhist - 1; ihist>=0; total += hist[ihist--])
01053 if (total >= want)
01054 break;
01055 newmax = *minp + ((double)ihist/(double)nhist) * (*maxp - *minp);
01056 fprintf(stderr,"ihist %d newmax=%f newmin=%f \n", ihist, newmax, newmin);
01057 if (newmax <= newmin)
01058 {
01059 newmax = *maxp;
01060 newmin = *minp;
01061 }
01062
01063 fprintf(stderr,"newmax %f minp %f \n", newmax, *minp);
01064 for (idata=0; idata<ndata; idata++)
01065 {
01066 float val = inData[idata];
01067 unsigned short newval;
01068 if (!isnan(val)){
01069 if (val <= newmin) val=newmin;
01070 if (val >= newmax) val=newmax;
01071 newval = maxcolor*(val-newmin)/(newmax-newmin);
01072 if (scaling == MAXMIN99)
01073 newval = maxcolor - newval;
01074 if (newval > maxcolor) newval = maxcolor;
01075 if (newval < 1) newval = 1;
01076 }
01077
01078 else
01079 newval = missingcolor;
01080 for (color=0; color<colors; color++)
01081 {
01082 if (bytepercolor == 1)
01083 *((unsigned char *)out + colors*idata + color) = colortable[newval*colors + color];
01084 else
01085 *((unsigned short *)out + colors*idata + color) = colortable[newval*colors + color];
01086 }
01087 }
01088 free(outa);
01089 free(hist);
01090 break;
01091 }
01092 }
01093 return(out);
01094 }
01095
01096
01097 int make_png(char *filename, unsigned char *data,
01098 int height, int width, char *palette, int bytepercolor, int colors)
01099 {
01100 int row;
01101 FILE *fp;
01102 png_byte **row_pointers = NULL;
01103 png_structp png_ptr = NULL;
01104 png_infop info_ptr = NULL;
01105
01106 fp = fopen(filename, "w");
01107 if (!fp)
01108 PNGDIE("cant open output file");
01109
01110 png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL,
01111 NULL, NULL);
01112 if (!png_ptr)
01113 PNGDIE("cant create png_struct");
01114
01115 info_ptr = png_create_info_struct(png_ptr);
01116 if (!info_ptr)
01117 PNGDIE("cant create png_infop");
01118
01119 if (setjmp(png_jmpbuf(png_ptr)))
01120 PNGDIE("png fatal error");
01121
01122 png_init_io(png_ptr, fp);
01123
01124 png_set_IHDR(png_ptr, info_ptr, width, height, 8*bytepercolor,
01125 (colors == 1 ? PNG_COLOR_TYPE_GRAY : PNG_COLOR_TYPE_RGB), PNG_INTERLACE_ADAM7,
01126 PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT);
01127 png_write_info(png_ptr, info_ptr);
01128 if (bytepercolor > 1)
01129 png_set_swap(png_ptr);
01130 png_set_gAMA(png_ptr, info_ptr, 2.2);
01131 png_set_sRGB(png_ptr, info_ptr, PNG_sRGB_INTENT_ABSOLUTE);
01132 row_pointers = (png_byte **)malloc(height * sizeof(png_byte *));
01133
01134 for (row=0; row<height ; row++)
01135 {
01136 if (bytepercolor == 1)
01137 row_pointers[height - row - 1] = (png_byte *)((char *)data + colors*row*width);
01138 else
01139 row_pointers[height - row - 1] = (png_byte *)((short *)data + colors*row*width);
01140 }
01141 png_set_rows(png_ptr, info_ptr, row_pointers);
01142 png_write_png(png_ptr, info_ptr, PNG_TRANSFORM_SWAP_ENDIAN, NULL);
01143 fclose(fp);
01144
01145 png_destroy_write_struct(&png_ptr, &info_ptr);
01146 free(row_pointers);
01147 return(0);
01148 }
01149
01150
01151
01152
01153
01154
01155
01156
01157
01158 int reposition(DRMS_Array_t *arr, ObsInfo_t *ObsLoc)
01159 {
01160 int status;
01161 int regrid = 1;
01162 int stretch = 0;
01163 int type_float = 3;
01164 int nx = arr->axis[0];
01165 int ny = arr->axis[1];
01166 float xc = 0.0;
01167 float yc = 0.0;
01168 float mag = 1.0;
01169 float dx = (nx + 1.0)/2.0 - ObsLoc->crpix1;
01170 float dy = (ny + 1.0)/2.0 - ObsLoc->crpix2;
01171 float *outarray;
01172 status = image_cutout(arr->data, nx, ny, type_float, ObsLoc->crota2,
01173 mag, dx, dy, &(void *)outarray, nx, ny, xc, yc, regrid, stretch);
01174 if (status) {fprintf(stderr,"image_cutout failed, staus=%d\n", status);}
01175 free(arr->data);
01176 arr->data = outarray;
01177 ObsLoc->crota2 = 0.0;
01178 ObsLoc->crpix1 += dx;
01179 ObsLoc->crpix2 += dy;
01180 return(status);
01181 }
01182
01183
01184
01185
01186
01187 int upNcenter(DRMS_Array_t *arr, ObsInfo_t *ObsLoc)
01188 {
01189 int nx, ny, ix, iy, i, j, xoff, yoff, max_off;
01190 double rot, x0, y0, midx, midy;
01191 float *data;
01192 float *data2;
01193 if (!arr || !ObsLoc)
01194 return(1);
01195 data = arr->data;
01196 nx = arr->axis[0];
01197 ny = arr->axis[1];
01198 x0 = ObsLoc->crpix1 - 1;
01199 y0 = ObsLoc->crpix2 - 1;
01200 midx = (nx-1.0)/2.0;
01201 midy = (ny-1.0)/2.0;
01202 if ((rot = fabs(ObsLoc->crota2)) > 179 && rot < 181)
01203 {
01204
01205 float val;
01206 int half = ny / 2;
01207 int odd = ny & 1;
01208 for (iy=0; iy<half; iy++)
01209 {
01210 for (ix=0; ix<nx; ix++)
01211 {
01212 i = iy*nx + ix;
01213 j = (ny - 1 - iy)*nx + (nx - 1 - ix);
01214 val = data[i];
01215 data[i] = data[j];
01216 data[j] = val;
01217 }
01218 }
01219
01220 if (odd) {
01221 for (ix=0; ix<nx/2; ++ix) {
01222 i = nx*half + ix;
01223 j = nx*half + nx - ix - 1;
01224 val = data[i];
01225 data[i] = data[j];
01226 data[j] = val;
01227 }
01228 }
01229 x0 = nx - 1 - x0;
01230 y0 = ny - 1 - y0;
01231 rot = ObsLoc->crota2 - 180.0;
01232 if (rot < -90.0) rot += 360.0;
01233 ObsLoc->crota2 = rot;
01234 }
01235
01236 xoff = round(x0 - midx);
01237 yoff = round(y0 - midy);
01238 max_off = 20.0 / ObsLoc->cdelt1;
01239 if (arr->parent_segment &&
01240 arr->parent_segment->record &&
01241 arr->parent_segment->record->seriesinfo &&
01242 arr->parent_segment->record->seriesinfo->seriesname &&
01243 strncasecmp(arr->parent_segment->record->seriesinfo->seriesname, "aia", 3) &&
01244 abs(xoff) < max_off && abs(yoff) < max_off)
01245 {
01246 if (abs(xoff) >= 1 || abs(yoff) >= 1)
01247 {
01248 data2 = malloc(4*nx*ny);
01249 for (iy=0; iy<ny; iy++)
01250 {
01251 int jy = iy + yoff;
01252 for (ix=0; ix<nx; ix++)
01253 {
01254 int jx = ix + xoff;
01255 int idx = jy*nx + jx;
01256 int idx2 = iy*nx + ix;
01257 if (jx<0 || jx>=nx || jy<0 || jy>=ny)
01258 data2[idx2] = DRMS_MISSING_FLOAT;
01259 else
01260 data2[idx2] = data[idx];
01261 }
01262 }
01263 x0 -= xoff;
01264 y0 -= yoff;
01265 free(data);
01266 arr->data = data2;
01267 }
01268 }
01269
01270 ObsLoc->crpix1 = x0 + 1;
01271 ObsLoc->crpix2 = y0 + 1;
01272 return(0);
01273 }
01274
01275
01276
01277 int crop_image(DRMS_Array_t *arr, ObsInfo_t *ObsLoc)
01278 {
01279 int nx, ny, ix, iy, i, j, xoff, yoff;
01280 double x0, y0;
01281 double rsun = ObsLoc->rsun_obs/ObsLoc->cdelt1;
01282 double scale, crop_limit;
01283 float *data;
01284 if (!arr || !ObsLoc)
01285 return(1);
01286 data = arr->data;
01287 nx = arr->axis[0];
01288 ny = arr->axis[1];
01289 x0 = ObsLoc->crpix1 - 1;
01290 y0 = ObsLoc->crpix2 - 1;
01291 scale = 1.0/rsun;
01292 crop_limit = 0.99975;
01293 for (iy=0; iy<ny; iy++)
01294 for (ix=0; ix<nx; ix++)
01295 {
01296 double x, y, R2;
01297 float *Ip = data + iy*nx + ix;
01298 if (drms_ismissing_float(*Ip))
01299 continue;
01300 x = ((double)ix - x0) * scale;
01301 y = ((double)iy - y0) * scale;
01302 R2 = x*x + y*y;
01303 if (R2 > crop_limit)
01304 *Ip = DRMS_MISSING_FLOAT;
01305 }
01306 return(0);
01307 }
01308
01309 #define CHECK(keyname) {if (status) {fprintf(stderr,"Keyword failure to find: %s, status=%d\n",keyname,status); *rstatus=status; return(NULL);}}
01310
01311 ObsInfo_t *GetObsInfo(DRMS_Segment_t *seg, ObsInfo_t *pObsLoc, int *rstatus)
01312 {
01313 TIME t_prev;
01314 DRMS_Record_t *rec;
01315 TIME t_obs;
01316 double dv;
01317 ObsInfo_t *ObsLoc;
01318 int status;
01319 char t_key[100];
01320
01321 if (!seg || !(rec = seg->record))
01322 { *rstatus = 1; return(NULL); }
01323
01324 ObsLoc = (pObsLoc ? pObsLoc : (ObsInfo_t *)malloc(sizeof(ObsInfo_t)));
01325 if (!pObsLoc)
01326 memset(ObsLoc, 0, sizeof(ObsInfo_t));
01327
01328 t_prev = ObsLoc->t_obs;
01329 if (drms_keyword_lookup(rec, "T_OBS", 0))
01330 strcpy(t_key, "T_OBS");
01331 else
01332 strcpy(t_key, "T_REC");
01333 t_obs = drms_getkey_time(rec, t_key, &status); CHECK(t_key);
01334
01335 if (t_obs <= 0.0)
01336 { *rstatus = 2; return(NULL); }
01337
01338 if (t_obs != t_prev)
01339 {
01340 ObsLoc->crpix1 = drms_getkey_double(rec, "CRPIX1", &status); CHECK("CRPIX1");
01341 ObsLoc->crpix2 = drms_getkey_double(rec, "CRPIX2", &status); CHECK("CRPIX2");
01342 ObsLoc->crval1 = drms_getkey_double(rec, "CRVAL1", &status); CHECK("CRVAL1");
01343 ObsLoc->crval2 = drms_getkey_double(rec, "CRVAL2", &status); CHECK("CRVAL2");
01344 ObsLoc->cdelt1 = drms_getkey_double(rec, "CDELT1", &status); CHECK("CDELT1");
01345 ObsLoc->cdelt2 = drms_getkey_double(rec, "CDELT2", &status); CHECK("CDELT2");
01346 ObsLoc->crota2 = drms_getkey_double(rec, "CROTA2", &status); if (status) ObsLoc->crota2 = 0.0;
01347 ObsLoc->sina = sin(ObsLoc->crota2*Deg2Rad);
01348 ObsLoc->cosa = sqrt (1.0 - ObsLoc->sina*ObsLoc->sina);
01349 ObsLoc->rsun_obs = drms_getkey_double(rec, "RSUN_OBS", &status);
01350 if (status)
01351 {
01352 double dsun_obs = drms_getkey_double(rec, "DSUN_OBS", &status); CHECK("DSUN_OBS");
01353 ObsLoc->rsun_obs = asin(696000000.0/dsun_obs)/arcsec2Rad;
01354 }
01355 ObsLoc->obs_vr = drms_getkey_double(rec, "OBS_VR", &status); CHECK("OBS_VR");
01356 ObsLoc->obs_vw = drms_getkey_double(rec, "OBS_VW", &status); CHECK("OBS_VW");
01357 ObsLoc->obs_vn = drms_getkey_double(rec, "OBS_VN", &status); CHECK("OBS_VN");
01358 ObsLoc->obs_b0 = drms_getkey_double(rec, "CRLT_OBS", &status); CHECK("CRLT_OBS");
01359 ObsLoc->t_obs = t_obs;
01360 }
01361 *rstatus = 0;
01362 return(ObsLoc);
01363 }
01364
01365
01366
