|
ocssw
1.0
|
00001 /* 00002 Inputs: 00003 see usage... 00004 Output: 00005 Eight-bit PGM data stream. 00006 00007 Code draws heavily on swplatecar, smigen, and l2bin 00008 00009 Original Development: 00010 Sean Bailey -- 30 June 2006 00011 00012 Modifications: 00013 SWB, 17 Jul 2008, fixed the threshold parameter - long/float discrepancy 00014 00015 */ 00016 00017 #include <stdlib.h> 00018 #include <math.h> 00019 #include <string.h> 00020 #include <stdint.h> 00021 #include <png.h> 00022 #include "miscstruct.h" 00023 #include "miscanfill.h" 00024 #include "mipoly.h" 00025 #include "mfhdf.h" 00026 #include <clo.h> 00027 #include "l2mapgen.h" 00028 #include "l2mapgen_input.h" 00029 #include <genutils.h> 00030 #include <readL2scan.h> 00031 00032 #include <stdio.h> 00033 #include <X11/X.h> 00034 #include <X11/Xlib.h> 00035 #include <unistd.h> 00036 00037 #include <geotiffio.h> 00038 #include <geo_normalize.h> 00039 #include <geo_tiffp.h> 00040 #include <geo_keyp.h> 00041 #include <xtiffio.h> 00042 #include <geokeys.h> 00043 00044 00045 #define ROOT2 1.4142135623730950488016887242096981 00046 #define PI 3.14159265358979323846 00047 #define BINBELOWTHRESH 110 00048 00049 #define CALLOC(ptr,typ,num) { \ 00050 (ptr) = (typ *)calloc((num) , sizeof(typ)); \ 00051 if((ptr) == NULL){ \ 00052 fprintf(stderr,"-E- %s line %d: Memory allocation failure.\n", \ 00053 __FILE__,__LINE__); \ 00054 exit(EXIT_FAILURE); \ 00055 } \ 00056 } 00057 00058 /* function prototypes */ 00059 int scan_convert(XPoint *ptsIn); 00060 int collect_bins(int number_of_initial_points, XPoint *initial_point, 00061 int *span_width); 00062 int miCreateETandAET(); 00063 int miInsertionSort(); 00064 00065 static int *binsoverlain = NULL; 00066 static int binsperpixel = 0; 00067 static int numoverlain; 00068 static int width, height; 00069 static double scale; 00070 00071 static int32 sd_id; 00072 00073 int32 l3m_params = 0; 00074 char **parmname_list; 00075 char **parmname_short; 00076 char **unit_list; 00077 char **scaling_list; 00078 float32 *maximum_list; 00079 float32 *minimum_list; 00080 char **palette_list; 00081 char **precision_list; 00082 00083 int main(int argc, char *argv[]) { 00084 00085 char title[512]; 00086 static instr input; 00087 static l2_prod l2_str[MAXNFILES]; 00088 static meta_l2Type meta_l2; 00089 static char buf[65535]; 00090 char qual[9]; 00091 00092 int32 npix, nscans; 00093 int32 nflag; 00094 int32 nfiles; 00095 int status; 00096 int32 l2_flags; 00097 byte quality; 00098 byte goodpix; 00099 int prodidx = -1; 00100 int qualidx = -1; 00101 double lat, lon; 00102 double north, south, west, east; 00103 double latmax = -90.0; 00104 double latmin = 90.0; 00105 double lonmin = 180.0; 00106 double lonmax = -180.0; 00107 00108 double *sum; 00109 unsigned char *mean; 00110 unsigned char *mask; 00111 int *count; 00112 00113 uint32_t numbins, bin, numoutpix; 00114 XPoint corners[8]; /* actually, 4 pixel corners + 4 side midpoints */ 00115 float64 threshold, outpixpercent; 00116 float32 nrad, srad, wrad, erad; 00117 int i, j, k, m, n, ii; 00118 int progress; 00119 uint32 mask_252, mask_253, all_masks; 00120 uint32 mask_glint = 0; 00121 uint32 required_mask; 00122 00123 char scale_type[50]; 00124 uint8 default_palfile = 0; 00125 char *cptr; 00126 00127 static uint16 maxval = 250; 00128 float32 val; 00129 float32 slope = 1.0; 00130 float32 intercept = 0.0; 00131 int32 prod_num = -1; 00132 byte *rgb; 00133 short *r, *g, *b; 00134 00135 FILE *fp = NULL; 00136 FILE *outfp = NULL; 00137 00138 00139 /* hold all of the command line options */ 00140 clo_optionList_t* list = clo_createList(); 00141 00142 if (argc == 1 || strcmp(argv[1], "-version") == 0 || strcmp(argv[1], "--version") == 0) { 00143 want_verbose = 0; 00144 l2mapgen_init_options(list); 00145 l2mapgen_read_options(list, argc, argv); 00146 clo_printUsage(list); 00147 exit(EXIT_FAILURE); 00148 } 00149 00150 if (l2mapgen_input(argc, argv, list, &input) != 0) { 00151 clo_printUsage(list); 00152 exit(EXIT_FAILURE); 00153 } 00154 00155 if (strcasecmp(input.palfile, "default") == 0) { 00156 default_palfile = 1; 00157 } 00158 00159 /* Read product table */ 00160 fp = fopen(input.product_table, "r"); 00161 if (fp == 0x0) { 00162 fprintf(stderr, "SMIGEN product table \"%s\" not found.\n", input.product_table); 00163 exit(EXIT_FAILURE); 00164 } 00165 00166 l3m_params = 0; 00167 while (fgets(buf, 128, fp) != NULL) { 00168 if ((buf[0] >= 0x41) && (buf[0] <= 0x5a)) l3m_params++; 00169 } 00170 fseek(fp, 0, SEEK_SET); 00171 00172 parmname_list = (char**) calloc(l3m_params, sizeof (char *)); 00173 parmname_short = (char**) calloc(l3m_params, sizeof (char *)); 00174 unit_list = (char**) calloc(l3m_params, sizeof (char *)); 00175 scaling_list = (char**) calloc(l3m_params, sizeof (char *)); 00176 palette_list = (char**) calloc(l3m_params, sizeof (char *)); 00177 maximum_list = (float32 *) calloc(l3m_params, sizeof (float32)); 00178 minimum_list = (float32 *) calloc(l3m_params, sizeof (float32)); 00179 precision_list = (char **) calloc(l3m_params, sizeof (char *)); 00180 00181 i = 0; 00182 while (fgets(buf, 128, fp) != NULL) { 00183 if ((buf[0] >= 0x41) && (buf[0] <= 0x5a)) { 00184 00185 cptr = strtok(buf, ":"); 00186 parmname_list[i] = (char*) malloc(strlen(cptr) + 1); 00187 strcpy(parmname_list[i], cptr); 00188 00189 cptr = strtok(NULL, ":"); 00190 parmname_short[i] = (char*) malloc(strlen(cptr) + 1); 00191 strcpy(parmname_short[i], cptr); 00192 00193 cptr = strtok(NULL, ":"); 00194 unit_list[i] = (char*) malloc(strlen(cptr) + 1); 00195 strcpy(unit_list[i], cptr); 00196 00197 cptr = strtok(NULL, ":"); 00198 scaling_list[i] = (char*) malloc(strlen(cptr) + 1); 00199 strcpy(scaling_list[i], cptr); 00200 00201 cptr = strtok(NULL, ":"); 00202 minimum_list[i] = (float32) atof(cptr); 00203 00204 cptr = strtok(NULL, ":"); 00205 maximum_list[i] = (float32) atof(cptr); 00206 00207 cptr = strtok(NULL, ":"); 00208 precision_list[i] = (char*) malloc(strlen(cptr) + 1); 00209 strcpy(precision_list[i], cptr); 00210 00211 cptr = strtok(NULL, "\n"); 00212 palette_list[i] = (char*) malloc(strlen(cptr) + 1); 00213 strcpy(palette_list[i], cptr); 00214 00215 i++; 00216 } 00217 } 00218 fclose(fp); 00219 00220 00221 for (i = 0; i < l3m_params; i++) { 00222 00223 if (strcmp(parmname_short[i], input.prod) == 0) { 00224 00225 prod_num = i; 00226 00227 /* define scaling */ 00228 strcpy(scale_type, scaling_list[i]); 00229 if (input.stype == 1) strcpy(scale_type, "linear"); 00230 if (input.stype == 2) strcpy(scale_type, "logarithmic"); 00231 00232 if (input.datamin == 0.0) input.datamin = minimum_list[i]; 00233 if (input.datamax == 0.0) input.datamax = maximum_list[i]; 00234 00235 00236 if (strcmp(scale_type, "linear") == 0) { 00237 00238 strcpy(scale_type, "LINEAR"); 00239 intercept = input.datamin; 00240 slope = (input.datamax - intercept) / maxval; 00241 } else if (strcmp(scale_type, "logarithmic") == 0) { 00242 00243 strcpy(scale_type, "LOG"); 00244 00245 intercept = log10(input.datamin); 00246 slope = (log10(input.datamax) - intercept) / maxval; 00247 } 00248 00249 /* Read palette file */ 00250 if (default_palfile) { 00251 strcpy(input.palfile, input.palette_dir); 00252 strcat(input.palfile, "/"); 00253 strcat(input.palfile, palette_list[i]); 00254 strcat(input.palfile, ".pal"); 00255 } 00256 00257 if (!(r = (short *) calloc(256, sizeof (short)))) { 00258 fprintf(stderr, "smigen: Error allocating space for red.\n"); 00259 return -1; 00260 }; 00261 if (!(g = (short *) calloc(256, sizeof (short)))) { 00262 fprintf(stderr, "smigen: Error allocating space for green.\n"); 00263 return -1; 00264 }; 00265 if (!(b = (short *) calloc(256, sizeof (short)))) { 00266 fprintf(stderr, "smigen: Error allocating space for blue.\n"); 00267 return -1; 00268 }; 00269 00270 if (getlut_file(input.palfile, r, g, b)) { 00271 fprintf(stderr, "Error reading palette file %s\n", input.palfile); 00272 free(r); 00273 free(g); 00274 free(b); 00275 return -1; 00276 } 00277 if (input.mask) { 00278 r[252] = 128; 00279 g[252] = 128; 00280 b[252] = 128; 00281 r[253] = 160; 00282 g[253] = 82; 00283 b[253] = 45; 00284 r[254] = 255; 00285 g[254] = 255; 00286 b[254] = 255; 00287 r[255] = 0; 00288 g[255] = 0; 00289 b[255] = 0; 00290 } 00291 for (i = 0; i < 256; i++) { 00292 input.palette[i * 3] = r[i]; 00293 input.palette[i * 3 + 1] = g[i]; 00294 input.palette[i * 3 + 2] = b[i]; 00295 } 00296 00297 break; 00298 } // if prod matches 00299 } // for i 00300 00301 if (prod_num == -1) { 00302 /* define scaling */ 00303 strcpy(scale_type, "LINEAR"); 00304 if (input.stype == 1) strcpy(scale_type, "linear"); 00305 if (input.stype == 2) strcpy(scale_type, "logarithmic"); 00306 00307 if (input.datamin == 0.0) input.datamin = 0.001; 00308 if (input.datamax == 0.0) input.datamax = 1.0; 00309 00310 00311 if (strcmp(scale_type, "linear") == 0) { 00312 00313 strcpy(scale_type, "LINEAR"); 00314 intercept = input.datamin; 00315 slope = (input.datamax - intercept) / maxval; 00316 } 00317 00318 if (strcmp(scale_type, "logarithmic") == 0) { 00319 00320 strcpy(scale_type, "LOG"); 00321 00322 intercept = log10(input.datamin); 00323 slope = (log10(input.datamax) - intercept) / maxval; 00324 } 00325 00326 /* Read default.pal palette file */ 00327 strcpy(input.palfile, input.palette_dir); 00328 strcat(input.palfile, "/default.pal"); 00329 00330 if (!(r = (short *) calloc(256, sizeof (short)))) { 00331 fprintf(stderr, "smigen: Error allocating space for red.\n"); 00332 return -1; 00333 }; 00334 if (!(g = (short *) calloc(256, sizeof (short)))) { 00335 fprintf(stderr, "smigen: Error allocating space for green.\n"); 00336 return -1; 00337 }; 00338 if (!(b = (short *) calloc(256, sizeof (short)))) { 00339 fprintf(stderr, "smigen: Error allocating space for blue.\n"); 00340 return -1; 00341 }; 00342 00343 if (getlut_file(input.palfile, r, g, b)) { 00344 fprintf(stderr, "Error reading palette file %s\n", input.palfile); 00345 free(r); 00346 free(g); 00347 free(b); 00348 return -1; 00349 } 00350 if (input.mask) { 00351 r[252] = 128; 00352 g[252] = 128; 00353 b[252] = 128; 00354 r[253] = 160; 00355 g[253] = 82; 00356 b[253] = 45; 00357 r[254] = 255; 00358 g[254] = 255; 00359 b[254] = 255; 00360 r[255] = 0; 00361 g[255] = 0; 00362 b[255] = 0; 00363 } 00364 for (i = 0; i < 256; i++) { 00365 input.palette[i * 3] = r[i]; 00366 input.palette[i * 3 + 1] = g[i]; 00367 input.palette[i * 3 + 2] = b[i]; 00368 } 00369 } 00370 00371 00372 /* Single HDF input */ 00373 /* ---------------- */ 00374 if (Hishdf(input.ifile) == TRUE) { 00375 nfiles = 1; 00376 status = openL2(input.ifile, 0x0, &l2_str[0]); 00377 00378 status = readL2meta(&meta_l2, 0); 00379 if (meta_l2.northlat >= latmax) latmax = meta_l2.northlat; 00380 if (meta_l2.southlat <= latmin) latmin = meta_l2.southlat; 00381 if (meta_l2.westlon <= lonmin) lonmin = meta_l2.westlon; 00382 if (meta_l2.eastlon >= lonmax) lonmax = meta_l2.eastlon; 00383 00384 closeL2(&l2_str[0], 0); 00385 fprintf(stderr, "Single HDF input\n"); 00386 } else { 00387 00388 /* Filelist input - Determine number of input files */ 00389 /* ------------------------------------------------ */ 00390 nfiles = 0; 00391 fp = fopen(input.ifile, "r"); 00392 if (fp == NULL) { 00393 fprintf(stderr, "Input listing file: \"%s\" not found.\n", input.ifile); 00394 return -1; 00395 } 00396 while (fgets(buf, 256, fp) != NULL) nfiles++; 00397 fclose(fp); 00398 fprintf(stderr, "%d input files\n", nfiles); 00399 00400 00401 /* Open L2 input files */ 00402 /* ------------------- */ 00403 fp = fopen(input.ifile, "r"); 00404 for (i = 0; i < nfiles; i++) { 00405 fgets(buf, 256, fp); 00406 buf[strlen(buf) - 1] = 0; 00407 status = openL2(buf, 0x0, &l2_str[i]); 00408 00409 status = readL2meta(&meta_l2, i); 00410 if (meta_l2.northlat >= latmax) latmax = meta_l2.northlat; 00411 if (meta_l2.southlat <= latmin) latmin = meta_l2.southlat; 00412 if (meta_l2.westlon <= lonmin) lonmin = meta_l2.westlon; 00413 if (meta_l2.eastlon >= lonmax) lonmax = meta_l2.eastlon; 00414 00415 closeL2(&l2_str[i], i); 00416 00417 } /* ifile loop */ 00418 fclose(fp); 00419 } 00420 /* Setup flag masking */ 00421 /* --------------- */ 00422 if (input.flaguse[0] == 0) { 00423 strcpy(input.flaguse, DEF_FLAG); 00424 } else { 00425 input.mask = 1; 00426 } 00427 strcpy(buf, l2_str[0].flagnames); 00428 setupflags(buf, "HIGLINT", &mask_252, &required_mask, &status); 00429 setupflags(buf, "LAND", &mask_253, &required_mask, &status); 00430 setupflags(buf, input.flaguse, &all_masks, &required_mask, &status); 00431 if ((all_masks & mask_252) != 0) 00432 mask_glint = 1; 00433 00434 /* Make sure L2 product exists for every input L2 file */ 00435 /* ---------------------------------------------------------------- */ 00436 status = 0; 00437 for (j = 0; j < nfiles; j++) { 00438 for (i = 0; i < l2_str[j].nprod; i++) { 00439 if (strcmp(l2_str[j].prodname[i], input.prod) == 0) { 00440 status++; 00441 } 00442 } 00443 } 00444 00445 if (status != nfiles) { 00446 fprintf(stderr, "Product %s not found in all L2 file(s)\n", input.prod); 00447 exit(EXIT_FAILURE); 00448 } 00449 00450 /* Get the box boundaries. */ 00451 if (input.north == input.south && input.west == input.east) { 00452 input.north = latmax; 00453 input.south = latmin; 00454 input.west = lonmin; 00455 input.east = lonmax; 00456 } 00457 00458 nrad = input.north * PI / 180; 00459 srad = input.south * PI / 180; 00460 wrad = input.west * PI / 180; 00461 erad = input.east * PI / 180; 00462 00463 fprintf(stderr, "Mapping data to:\n N: %8.5f\n S: %8.5f\n W: %8.5f\n E: %8.5f\n", 00464 input.north, input.south, input.west, input.east); 00465 00466 fprintf(stderr, "Scale Type : %s\n", scale_type); 00467 fprintf(stderr, "Data Min (abs) : %8.4f\n", input.datamin); 00468 fprintf(stderr, "Data Max (abs) : %8.4f\n", input.datamax); 00469 fprintf(stderr, "Scale Slope : %8.4f\n", slope); 00470 fprintf(stderr, "Scale Intercept : %8.4f\n", intercept); 00471 if (input.apply_pal >= 1 || default_palfile == 0) 00472 fprintf(stderr, "Applying palette: %s\n", input.palfile); 00473 if (input.mask) 00474 fprintf(stderr, "Applying masking to flagged pixels\n"); 00475 00476 if (nrad <= srad) { 00477 fprintf(stderr, "The northernmost boundary must be greater than the "); 00478 fprintf(stderr, "southernmost boundary.\n"); 00479 exit(EXIT_FAILURE); 00480 } 00481 00482 /* Get the size of the output image. */ 00483 width = input.width; 00484 if (width < 32) { 00485 fprintf(stderr, "Width (%d) is too small to produce a useful image.\n", 00486 input.width); 00487 exit(EXIT_FAILURE); 00488 } 00489 if (wrad < erad) 00490 height = rint((nrad - srad) * width / (erad - wrad)); 00491 else 00492 height = rint((nrad - srad) * width / (erad - wrad + 2 * PI)); 00493 00494 numbins = width * height; 00495 scale = height / (nrad - srad); 00496 00497 /* Allocate memory for the output data. */ 00498 CALLOC(sum, double, numbins); 00499 CALLOC(count, int, numbins); 00500 CALLOC(mean, unsigned char, numbins); 00501 CALLOC(mask, unsigned char, numbins); 00502 CALLOC(rgb, byte, numbins * 3); 00503 /* Get the threshold percentage. */ 00504 threshold = (double) input.threshold; 00505 00506 /* For each input image... */ 00507 for (k = 0; k < nfiles; k++) { 00508 fprintf(stderr, "Opening HDF file, %s ...\n", l2_str[k].filename); 00509 status = reopenL2(k, &l2_str[k]); 00510 00511 nscans = l2_str[k].nrec; 00512 npix = l2_str[k].nsamp; 00513 fprintf(stderr, "pix: %d scan: %d\n", npix, nscans); 00514 00515 if (strcmp(input.prod, "sst") == 0 || strcmp(input.prod, "sst4") == 0) { 00516 strcpy(qual, "qual_"); 00517 strcat(qual, input.prod); 00518 for (i = 0; i < l2_str[k].nprod; i++) { 00519 if (strcmp(qual, l2_str[k].prodname[i]) == 0) { 00520 qualidx = i; 00521 break; 00522 } 00523 } 00524 } 00525 for (i = 0; i < l2_str[k].nprod; i++) 00526 if (strcmp(input.prod, l2_str[k].prodname[i]) == 0) { 00527 prodidx = i; 00528 break; 00529 } 00530 00531 if (npix <= 0 || nscans <= 0) { 00532 fprintf(stderr, 00533 "-E- %s line %d: Bad scene dimension: npix=%d nscans=%d in file, %s.\n", 00534 __FILE__, __LINE__, npix, nscans, l2_str[k].filename); 00535 exit(EXIT_FAILURE); 00536 } 00537 00538 /* Use the following variable to show this program's progress. */ 00539 progress = (int) ceil(((double) nscans / 78)); 00540 00541 00542 /* For each scan line ... */ 00543 fprintf(stderr, "Mapping swath pixels to Plate Carree projection...\n"); 00544 for (i = 0; i < nscans; i++) { 00545 00546 /* Read swath record from L2 */ 00547 /* ------------------------- */ 00548 status = readL2(&l2_str[k], k, i, -1, NULL); 00549 00550 /* For each pixel ... */ 00551 for (j = 0; j < npix; j++) { 00552 00553 double plat, plon; /* pixel center */ 00554 double sinplat, cosplat; 00555 double c[2]; /* edge/corner distance */ 00556 double sinc[2], cosc[2]; 00557 double sinaz[8], cosaz[8]; 00558 int out = 0; /* Number of pixel "corners" */ 00559 /* that fall outside the image */ 00560 00561 plat = l2_str[k].latitude[j] * PI / 180.0; 00562 plon = l2_str[k].longitude[j] * PI / 180.0; 00563 00564 sinplat = sin(plat); 00565 cosplat = cos(plat); 00566 00567 if (j < npix - 1) { 00568 00569 double nlat, nlon; /* next pixel center */ 00570 double dlat, dlon; /* delta lat and lon */ 00571 double sindlat2, sindlon2, cosnlat; 00572 double az; /* azimuth to next pixel */ 00573 double phw; /* pixel half width */ 00574 00575 00576 nlat = l2_str[k].latitude[j + 1] * PI / 180.0; 00577 nlon = l2_str[k].longitude[j + 1] * PI / 180.0; 00578 00579 cosnlat = cos(nlat); 00580 00581 dlat = nlat - plat; 00582 dlon = nlon - plon; 00583 00584 sindlat2 = sin(dlat / 2); 00585 sindlon2 = sin(dlon / 2); 00586 00587 phw = asin(sqrt(sindlat2 * sindlat2 + /* Page 30, */ 00588 sindlon2 * sindlon2 * /* Equation */ 00589 cosplat * cosnlat)); /* 5-3a */ 00590 00591 /* 00592 Make the pixel's coverage slightly broader to avoid 00593 black pin holes in the output image where the estimated 00594 pixel boundaries do not quite line up. Don't increase 00595 the size too much, however, because that will cause 00596 unwanted image blurring. 00597 */ 00598 phw *= 1.3; 00599 00600 c[0] = phw; /* center-to-edge distance */ 00601 c[1] = phw * ROOT2; /* center-to-corner distance */ 00602 00603 sinc[0] = sin(c[0]); 00604 sinc[1] = sin(c[1]); 00605 cosc[0] = cos(c[0]); 00606 cosc[1] = cos(c[1]); 00607 00608 az = /* Page 30, Equation 5-4b */ 00609 atan2( 00610 cosnlat * sin(dlon), 00611 cosplat * sin(nlat) - sinplat * cosnlat * cos(dlon) 00612 ); 00613 00614 sinaz[0] = sin(az); 00615 cosaz[0] = cos(az); 00616 for (ii = 1; ii < 8; ii++) { 00617 az += PI / 4; 00618 sinaz[ii] = sin(az); 00619 cosaz[ii] = cos(az); 00620 } 00621 } 00622 00623 for (ii = 0; ii < 8; ii++) { 00624 double lat, lon; 00625 int ec; /* edge = 0, corner = 1 */ 00626 00627 ec = ii & 1; 00628 00629 /* Page 31, Equation 5-5 */ 00630 lat = asin(sinplat * cosc[ec] 00631 + cosplat * sinc[ec] * cosaz[ii]); 00632 00633 lon = plon /* Page 31, Equation 5-6 */ 00634 + atan2(sinc[ec] * sinaz[ii], 00635 cosplat * cosc[ec] 00636 - sinplat * sinc[ec] * cosaz[ii]); 00637 00638 if (wrad > erad) { /* 180-degree meridian included */ 00639 if (lon >= wrad && lon > erad) { 00640 corners[ii].x = (lon - wrad) * scale; 00641 } else if (lon < wrad && lon <= erad) { 00642 corners[ii].x = (lon - wrad + 2 * PI) * scale; 00643 } else if (lon - erad < wrad - lon) { 00644 corners[ii].x = (lon - wrad + 2 * PI) * scale; 00645 } else { 00646 corners[ii].x = (lon - wrad) * scale; 00647 } 00648 } else { /* 180-degree meridian not included */ 00649 corners[ii].x = (lon - wrad) * scale; 00650 } 00651 corners[ii].y = (nrad - lat) * scale; 00652 00653 if (corners[ii].x < 0 || corners[ii].x >= width 00654 || corners[ii].y < 0 || corners[ii].y >= height) out++; 00655 } 00656 00657 /* 00658 If out == 8, then the entire pixel is outside the 00659 image area, so I skip to the next one. 00660 */ 00661 if (out == 8) continue; 00662 00663 l2_flags = l2_str[k].l2_flags[j]; 00664 if (qualidx >= 0) { 00665 quality = l2_str[k].l2_data[qualidx * l2_str[k].nsamp + j]; 00666 } 00667 00668 00669 /* Use scan conversion to determine 00670 * which bins are overlain by this pixel. */ 00671 numoverlain = 0; 00672 scan_convert(corners); 00673 00674 for (m = 0; m < numoverlain; m++) { 00675 n = binsoverlain[m]; 00676 00677 /* 00678 Keep track of the masks if so requested. 00679 I am basing the mask logic on the code in put_l2brs.c 00680 which is part of the level-2 browse file generator (l2brsgen). 00681 At the time I write this, pixel values 251 and 255 are unused 00682 in the SeaWiFS level-2 browse files. Norman Kuring 10-Jul-2001 00683 */ 00684 00685 if (input.mask && qualidx < 0) { 00686 if (l2_flags & mask_252 && mask_glint && 00687 strcmp(input.prod, "sst") != 0 && 00688 strcmp(input.prod, "sst4") != 0) { 00689 mask[n] = 252; 00690 } else if (l2_flags & mask_253) { 00691 if (mask[n] != 252) mask[n] = 253; 00692 } else if (l2_flags & all_masks) { 00693 if (mask[n] != 252 && mask[n] != 253) mask[n] = 254; 00694 } 00695 } else if (input.mask && qualidx >= 0) { 00696 if (l2_flags & mask_253) { 00697 mask[n] = 253; 00698 } else if (quality > input.quality) { 00699 if (mask[n] != 253) mask[n] = 254; 00700 } 00701 } 00702 00703 /* Accumulate the sums for each bin. Only count unmasked pixels. */ 00704 goodpix = 1; 00705 if (qualidx >= 0 && quality > input.quality) goodpix = 0; 00706 else if (qualidx < 0 && (l2_flags & all_masks) != 0) goodpix = 0; 00707 00708 if (goodpix) { 00709 val = l2_str[k].l2_data[prodidx * l2_str[k].nsamp + j]; 00710 if (val > input.datamax) val = input.datamax; 00711 if (val < input.datamin) val = input.datamin; 00712 sum[n] += val; 00713 count[n]++; 00714 } 00715 } 00716 } 00717 if (i != 0 && i % progress == 0) fprintf(stderr, "."); 00718 } 00719 fprintf(stderr, "\n"); 00720 } 00721 fprintf(stderr, "Finished Plate Carree projection mapping\n"); 00722 00723 00724 /* Calculate the means for each bin. */ 00725 fprintf(stderr, "Computing means...\n"); 00726 numoutpix = 0; 00727 for (n = 0; n < numbins; n++) { 00728 if (count[n] == 0) { 00729 if (input.mask) { 00730 if (mask[n] > 0) { 00731 mean[n] = mask[n]; 00732 numoutpix++; 00733 } else { 00734 mean[n] = 255; 00735 } 00736 } 00737 } else { 00738 if (strcmp(scale_type, "LINEAR") == 0) { 00739 mean[n] = (sum[n] / count[n] - intercept) / slope; 00740 } else if (strcmp(scale_type, "LOG") == 0) { 00741 mean[n] = (log10(sum[n] / count[n]) - intercept) / slope; 00742 } 00743 numoutpix++; 00744 } 00745 } 00746 00747 if ((double) 100 * numoutpix / numbins < threshold) { 00748 fprintf(stderr, "Number of output pixels (%u of a possible %u) ", 00749 numoutpix, numbins); 00750 fprintf(stderr, "< threshold (%.2f%%). Output image not generated.\n", 00751 threshold); 00752 exit(BINBELOWTHRESH); 00753 } 00754 00755 if (input.outmode == 1) { // ppm 00756 00757 // setup output file pointer 00758 if (strlen(input.ofile) > 0) { 00759 fprintf(stderr, "Writing out file: %s ...\n", input.ofile); 00760 outfp = fopen(input.ofile, "wb"); 00761 } else { 00762 fprintf(stderr, "Writing out data to STDOUT...\n"); 00763 outfp = stdout; 00764 } 00765 00766 if (input.apply_pal == 0 && default_palfile == 1) { 00767 fprintf(outfp, "P5\n%d %d\n255\n", width, height); 00768 for (n = 0; n < numbins; n++) { 00769 fputc((int) mean[n], outfp); 00770 } 00771 } else { 00772 /* 00773 * Convert from greyscale and palette to rgb array 00774 */ 00775 for (i = 0; i < numbins; i++) { 00776 rgb[3 * i] = r[mean[i]]; 00777 rgb[3 * i + 1] = g[mean[i]]; 00778 rgb[3 * i + 2] = b[mean[i]]; 00779 } 00780 00781 fprintf(outfp, "P6\n%d %d\n255\n", width, height); 00782 fwrite(&rgb[0], 1, width * height * 3, outfp); 00783 } 00784 00785 // close file if not using stdout 00786 if (stdout != outfp) 00787 fclose(outfp); 00788 00789 // end ppm 00790 00791 } else if (input.outmode == 2) { // png 00792 00793 // setup output file pointer 00794 if (strlen(input.ofile) > 0) { 00795 fprintf(stderr, "Writing out file: %s ...\n", input.ofile); 00796 outfp = fopen(input.ofile, "wb"); 00797 } else { 00798 fprintf(stderr, "Writing out data to STDOUT...\n"); 00799 outfp = stdout; 00800 } 00801 00802 png_structp png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, 00803 NULL, NULL, NULL); 00804 if (!png_ptr) { 00805 fprintf(stderr, "l2mapgen: Error: Unable to create PNG write structure.\n"); 00806 exit(EXIT_FAILURE); 00807 } 00808 00809 png_infop info_ptr = png_create_info_struct(png_ptr); 00810 if (!info_ptr) { 00811 png_destroy_write_struct(&png_ptr, (png_infopp) NULL); 00812 fprintf(stderr, "l2mapgen: Error: Unable to create PNG info structure.\n"); 00813 exit(EXIT_FAILURE); 00814 } 00815 if (setjmp(png_jmpbuf(png_ptr))) { 00816 png_destroy_write_struct(&png_ptr, (png_infopp) NULL); 00817 fprintf(stderr, "l2mapgen: Error: Unable to call PNG setjmp().\n"); 00818 exit(EXIT_FAILURE); 00819 } 00820 png_init_io(png_ptr, outfp); 00821 00822 uint8 * row_pointers[height]; 00823 for (i = 0; i < height; i++) 00824 row_pointers[i] = mean + (i * width); 00825 png_set_rows(png_ptr, info_ptr, row_pointers); 00826 00827 if (input.apply_pal == 0 && default_palfile == 1) { 00828 00829 // Grayscale 00830 png_set_IHDR(png_ptr, info_ptr, width, height, 00831 8, PNG_COLOR_TYPE_GRAY, PNG_INTERLACE_NONE, 00832 PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE); 00833 } else { 00834 00835 // color palette 00836 png_set_IHDR(png_ptr, info_ptr, width, height, 00837 8, PNG_COLOR_TYPE_PALETTE, PNG_INTERLACE_NONE, 00838 PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE); 00839 png_set_PLTE(png_ptr, info_ptr, (png_const_colorp) input.palette, 256); 00840 } 00841 00842 png_write_png(png_ptr, info_ptr, PNG_TRANSFORM_IDENTITY, NULL); 00843 00844 /* clean up after the write, and free any memory allocated */ 00845 png_destroy_write_struct(&png_ptr, (png_infopp) NULL); 00846 00847 // close file if not using stdout 00848 if (stdout != outfp) 00849 fclose(outfp); 00850 00851 // end png 00852 00853 } else if (input.outmode == 3) { // geotiff 00854 uint16 *rr, *gg, *bb; 00855 TIFF *tiff; 00856 GTIF *gtif; 00857 float *fmean; 00858 00859 if (strlen(input.ofile) == 0) { 00860 fprintf(stderr, "Can not write TIFF file to stdout.\n"); 00861 exit(EXIT_FAILURE); 00862 } 00863 00864 tiff = XTIFFOpen(input.ofile, "w"); 00865 if (tiff == NULL) { 00866 fprintf(stderr, "Could not open outgoing image\n"); 00867 exit(EXIT_FAILURE); 00868 } 00869 gtif = GTIFNew(tiff); 00870 if (gtif == NULL) { 00871 fprintf(stderr, "Could not create geoTIFF structure\n"); 00872 exit(EXIT_FAILURE); 00873 } 00874 00875 // calc geo TIFF tags 00876 double tiepoints[6] = {0, 0, 0, 0, 0, 0}; 00877 double pixscale[3] = {0, 0, 0}; 00878 00879 // pixel width 00880 pixscale[0] = (input.east - input.west) / width; 00881 00882 // pixel height 00883 pixscale[1] = (input.north - input.south) / height; 00884 00885 // set top left corner pixel lat, lon 00886 tiepoints[3] = input.west; 00887 tiepoints[4] = input.north; 00888 00889 TIFFSetField(tiff, TIFFTAG_IMAGEWIDTH, width); 00890 TIFFSetField(tiff, TIFFTAG_IMAGELENGTH, height); 00891 TIFFSetField(tiff, GTIFF_PIXELSCALE, 3, pixscale); 00892 TIFFSetField(tiff, GTIFF_TIEPOINTS, 6, tiepoints); 00893 00894 if (input.apply_pal == 0 && default_palfile == 1) { 00895 // Grayscale 00896 00897 TIFFSetField(tiff, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG); 00898 TIFFSetField(tiff, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_MINISBLACK); 00899 TIFFSetField(tiff, TIFFTAG_BITSPERSAMPLE, 32); 00900 TIFFSetField(tiff, TIFFTAG_SAMPLESPERPIXEL, 1); 00901 TIFFSetField(tiff, TIFFTAG_SAMPLEFORMAT, SAMPLEFORMAT_IEEEFP); 00902 00903 // write geo TIFF keys 00904 GTIFKeySet(gtif, GTModelTypeGeoKey, TYPE_SHORT, 1, ModelGeographic); 00905 GTIFKeySet(gtif, GTRasterTypeGeoKey, TYPE_SHORT, 1, RasterPixelIsArea); 00906 GTIFKeySet(gtif, GeographicTypeGeoKey, TYPE_SHORT, 1, GCS_WGS_84); 00907 00908 GTIFWriteKeys(gtif); 00909 00910 fmean = (float*) malloc(numbins * sizeof (float)); 00911 if (fmean == NULL) { 00912 fprintf(stderr, "Could not allocate memory for TIFF grayscale mean\n"); 00913 exit(EXIT_FAILURE); 00914 } 00915 00916 // calc the mean as a float 00917 for (i = 0; i < numbins; i++) { 00918 if (count[i] > 0) { 00919 fmean[i] = sum[i] / count[i]; 00920 } else { 00921 fmean[i] = -32767.0; 00922 } 00923 } 00924 00925 // Actually write the image 00926 if (TIFFWriteEncodedStrip(tiff, 0, fmean, numbins * sizeof (float)) == 0) { 00927 fprintf(stderr, "Could not write TIFF image\n"); 00928 exit(EXIT_FAILURE); 00929 } 00930 00931 free(fmean); 00932 00933 } else { 00934 // Colormap 00935 rr = (uint16*) malloc(256 * sizeof (uint16)); 00936 gg = (uint16*) malloc(256 * sizeof (uint16)); 00937 bb = (uint16*) malloc(256 * sizeof (uint16)); 00938 if (rr == NULL || gg == NULL || bb == NULL) { 00939 fprintf(stderr, "Could not allocate memory for TIFF color map\n"); 00940 exit(EXIT_FAILURE); 00941 } 00942 00943 // need a colormap of shorts not bytes 00944 for (i = 0; i < 256; i++) { 00945 rr[i] = r[i] << 8; 00946 gg[i] = g[i] << 8; 00947 bb[i] = b[i] << 8; 00948 } 00949 00950 TIFFSetField(tiff, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG); 00951 TIFFSetField(tiff, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_PALETTE); 00952 TIFFSetField(tiff, TIFFTAG_BITSPERSAMPLE, 8); 00953 TIFFSetField(tiff, TIFFTAG_SAMPLESPERPIXEL, 1); 00954 TIFFSetField(tiff, TIFFTAG_COLORMAP, rr, gg, bb); 00955 00956 // write geo TIFF keys 00957 GTIFKeySet(gtif, GTModelTypeGeoKey, TYPE_SHORT, 1, ModelGeographic); 00958 GTIFKeySet(gtif, GTRasterTypeGeoKey, TYPE_SHORT, 1, RasterPixelIsArea); 00959 GTIFKeySet(gtif, GeographicTypeGeoKey, TYPE_SHORT, 1, GCS_WGS_84); 00960 00961 GTIFWriteKeys(gtif); 00962 00963 // Actually write the image 00964 if (TIFFWriteEncodedStrip(tiff, 0, mean, width * height) == 0) { 00965 fprintf(stderr, "Could not write TIFF image\n"); 00966 exit(EXIT_FAILURE); 00967 } 00968 00969 free(rr); 00970 free(gg); 00971 free(bb); 00972 00973 } // color map 00974 00975 GTIFFree(gtif); 00976 XTIFFClose(tiff); 00977 00978 // end geotiff 00979 00980 } else { 00981 00982 fprintf(stderr, "l2mapgen: Error: invalid value for outmode - %d.\n", input.outmode); 00983 exit(EXIT_FAILURE); 00984 } 00985 00986 for (i = 0; i < l3m_params; i++) { 00987 free(parmname_list[i]); 00988 free(parmname_short[i]); 00989 free(unit_list[i]); 00990 free(scaling_list[i]); 00991 free(palette_list[i]); 00992 free(precision_list[i]); 00993 } 00994 free(parmname_list); 00995 free(parmname_short); 00996 free(unit_list); 00997 free(scaling_list); 00998 free(palette_list); 00999 free(maximum_list); 01000 free(minimum_list); 01001 free(precision_list); 01002 01003 free(sum); 01004 free(count); 01005 free(mean); 01006 free(mask); 01007 free(rgb); 01008 free(r); 01009 free(g); 01010 free(b); 01011 fprintf(stderr, "Done.\n"); 01012 01013 return (EXIT_SUCCESS); 01014 01015 } 01016 01017 01018 01019 /*------------------------------------------------------------------------------ 01020 01021 Function: scan_convert 01022 01023 Returns type: int 01024 01025 Description: This function is derived from the scan-conversion 01026 code used by an X-server for filling in polygons. 01027 01028 Parameters: (in calling order) 01029 Type Name I/O Description 01030 ---- ---- --- ----------- 01031 struct _DDXPoint* ptsIn In vertex specifications 01032 01033 Modification history: 01034 Programmer Date Description of change 01035 ---------- ---- --------------------- 01036 Norman Kuring 15-Sep-1992 Adaptation from X-server code 01037 for my own nefarious purposes. 01038 01039 ------------------------------------------------------------------------------*/ 01040 01041 /*********************************************************** 01042 Copyright 1987 by Digital Equipment Corporation, Maynard, Massachusetts, 01043 and the Massachusetts Institute of Technology, Cambridge, Massachusetts. 01044 01045 All Rights Reserved 01046 01047 Permission to use, copy, modify, and distribute this software and its 01048 documentation for any purpose and without fee is hereby granted, 01049 provided that the above copyright notice appear in all copies and that 01050 both that copyright notice and this permission notice appear in 01051 supporting documentation, and that the names of Digital or MIT not be 01052 used in advertising or publicity pertaining to distribution of the 01053 software without specific, written prior permission. 01054 01055 DIGITAL DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING 01056 ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, IN NO EVENT SHALL 01057 DIGITAL BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR 01058 ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, 01059 WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, 01060 ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS 01061 SOFTWARE. 01062 01063 ******************************************************************/ 01064 /* $XConsortium: mipolygen.c,v 5.0 89/06/09 15:08:35 keith Exp $ */ 01065 01066 extern void miloadAET(), micomputeWAET(), miFreeStorage(); 01067 01068 /* 01069 * 01070 * Written by Brian Kelleher; Oct. 1985 01071 * 01072 * Routine to fill a polygon. Two fill rules are 01073 * supported: frWINDING and frEVENODD. 01074 * 01075 * See fillpoly.h for a complete description of the algorithm. 01076 * 01077 * Modified by Norman Kuring to fill a bin-number array instead of 01078 * filling shapes on the server. 01079 */ 01080 01081 int scan_convert(XPoint *ptsIn) { 01082 int count = 8; 01083 register EdgeTableEntry *pAET; /* the Active Edge Table */ 01084 register int y; /* the current scanline */ 01085 register int nPts = 0; /* number of pts in buffer */ 01086 register EdgeTableEntry *pWETE; /* Winding Edge Table */ 01087 register ScanLineList *pSLL; /* Current ScanLineList */ 01088 register XPoint *ptsOut; /* ptr to output buffers */ 01089 int *wdth; 01090 XPoint FirstPoint[NUMPTSTOBUFFER]; /* the output buffers */ 01091 int FirstWidth[NUMPTSTOBUFFER]; 01092 EdgeTableEntry *pPrevAET; /* previous AET entry */ 01093 EdgeTable ET; /* Edge Table header node */ 01094 EdgeTableEntry AET; /* Active ET header node */ 01095 EdgeTableEntry *pETEs; /* Edge Table Entries buff */ 01096 ScanLineListBlock SLLBlock; /* header for ScanLineList */ 01097 int fixWAET = 0; 01098 01099 if (!(pETEs = (EdgeTableEntry *) malloc(sizeof (EdgeTableEntry) * count))) { 01100 return (0); 01101 } 01102 ptsOut = FirstPoint; 01103 wdth = FirstWidth; 01104 if (!miCreateETandAET(count, ptsIn, &ET, &AET, pETEs, &SLLBlock)) { 01105 free(pETEs); 01106 return (0); 01107 } 01108 pSLL = ET.scanlines.next; 01109 01110 /* 01111 * for each scanline 01112 */ 01113 for (y = ET.ymin; y < ET.ymax; y++) { 01114 /* 01115 * Add a new edge to the active edge table when we 01116 * get to the next edge. 01117 */ 01118 if (pSLL && y == pSLL->scanline) { 01119 miloadAET(&AET, pSLL->edgelist); 01120 micomputeWAET(&AET); 01121 pSLL = pSLL->next; 01122 } 01123 pPrevAET = &AET; 01124 pAET = AET.next; 01125 pWETE = pAET; 01126 01127 /* 01128 * for each active edge 01129 */ 01130 while (pAET) { 01131 /* 01132 * if the next edge in the active edge table is 01133 * also the next edge in the winding active edge 01134 * table. 01135 */ 01136 if (pWETE == pAET) { 01137 ptsOut->x = pAET->bres.minor; 01138 ptsOut++->y = y; 01139 *wdth++ = pAET->nextWETE->bres.minor - pAET->bres.minor; 01140 nPts++; 01141 01142 /* 01143 * send out the buffer 01144 */ 01145 if (nPts == NUMPTSTOBUFFER) { 01146 collect_bins(nPts, FirstPoint, FirstWidth); 01147 ptsOut = FirstPoint; 01148 wdth = FirstWidth; 01149 nPts = 0; 01150 } 01151 01152 pWETE = pWETE->nextWETE; 01153 while (pWETE != pAET) 01154 EVALUATEEDGEWINDING(pAET, pPrevAET, y, fixWAET); 01155 pWETE = pWETE->nextWETE; 01156 } 01157 EVALUATEEDGEWINDING(pAET, pPrevAET, y, fixWAET); 01158 } 01159 01160 /* 01161 * reevaluate the Winding active edge table if we 01162 * just had to resort it or if we just exited an edge. 01163 */ 01164 if (miInsertionSort(&AET) || fixWAET) { 01165 micomputeWAET(&AET); 01166 fixWAET = 0; 01167 } 01168 } 01169 01170 /* 01171 * Get any spans that we missed by buffering 01172 */ 01173 collect_bins(nPts, FirstPoint, FirstWidth); 01174 free(pETEs); 01175 miFreeStorage(SLLBlock.next); 01176 return (1); 01177 } 01178 01179 int collect_bins( 01180 int number_of_initial_points, 01181 XPoint *initial_point, 01182 int *span_width 01183 ) { 01184 01185 int i; 01186 short x, y; 01187 int end; 01188 01189 /* The variables, width, height, binsoverlain, and numoverlain, are 01190 * static global varibles. */ 01191 01192 for (i = 0; i < number_of_initial_points; i++) { 01193 01194 y = initial_point[i].y; 01195 if (y >= height || y < 0) continue; 01196 01197 for ( 01198 x = initial_point[i].x, 01199 end = initial_point[i].x + span_width[i]; 01200 x < end; 01201 x++ 01202 ) { 01203 if (x >= width || x < 0) continue; 01204 if (numoverlain >= binsperpixel) { 01205 binsperpixel += 1024; 01206 binsoverlain = (int *) realloc(binsoverlain, binsperpixel * sizeof (int)); 01207 if (binsoverlain == NULL) { 01208 fprintf(stderr, "-E- %s line %d: ", __FILE__, __LINE__); 01209 fprintf(stderr, 01210 "Memory allocation error (numoverlain=%d binsperpixel=%d\n", 01211 numoverlain, binsperpixel); 01212 exit(EXIT_FAILURE); 01213 } 01214 return (0); 01215 } 01216 binsoverlain[numoverlain++] = y * width + x; 01217 } 01218 } 01219 return (1); 01220 }
1.7.6.1