Changeset 36357 for trunk/ppSkycell/src/ppSkycellLoop.c
- Timestamp:
- Dec 5, 2013, 7:02:38 PM (13 years ago)
- File:
-
- 1 edited
-
trunk/ppSkycell/src/ppSkycellLoop.c (modified) (17 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/ppSkycell/src/ppSkycellLoop.c
r36356 r36357 11 11 #define BUFFER 16 // Size of buffer for projections 12 12 13 #if (0)14 // This is all junk that is basically useless, because it calculates things we already know.15 16 static void regionMinMax(psRegion *base,// Base region; modified17 const psRegion *compare // Comparison region18 )19 {20 base->x0 = PS_MAX(base->x0, compare->x0);21 base->x1 = PS_MIN(base->x1, compare->x1);22 base->y0 = PS_MIN(base->y0, compare->y0);23 base->y1 = PS_MAX(base->y1, compare->y1);24 return;25 }26 27 28 static psRegion *skycellRegion(pmAstromWCS *wcs, // World Coordinate System29 int numCols, int numRows // Size of image30 )31 {32 psPlane *fromCoords = psPlaneAlloc(), *toCoords = psPlaneAlloc(); // Coordinates for transforms33 psRegion *region = psRegionAlloc(-INFINITY, INFINITY, INFINITY, -INFINITY); // Region for skycell34 35 // Limit the region using the nominated point (X,Y)36 #define REGION_RANGE(X,Y) \37 fromCoords->x = (X); \38 fromCoords->y = (Y); \39 psPlaneTransformApply(toCoords, wcs->trans, fromCoords); \40 printf("%f %f %f %f\n",toCoords->x,toCoords->y,fromCoords->x,fromCoords->y); \41 toCoords->x /= wcs->cdelt1; \42 toCoords->y /= wcs->cdelt2; \43 toCoords->x -= wcs->crpix1; \44 toCoords->y += wcs->crpix2; \45 region->x0 = PS_MAX(region->x0, toCoords->x); \46 region->x1 = PS_MIN(region->x1, toCoords->x); \47 region->y0 = PS_MIN(region->y0, toCoords->y); \48 region->y1 = PS_MAX(region->y1, toCoords->y);49 //50 REGION_RANGE(0, 0); // Lower left51 psTrace("ppSkycell", 9, "Start: %.0f %.0f", toCoords->x, toCoords->y);52 REGION_RANGE(0, numRows); // Upper left53 REGION_RANGE(numCols, 0); // Lower right54 REGION_RANGE(numCols, numRows); // Upper right55 psTrace("ppSkycell", 9, "Stop: %.0f %.0f\n", toCoords->x, toCoords->y);56 psTrace("ppSkycell", 7, "%g %g %g %g\n",wcs->cdelt1,wcs->cdelt2,wcs->crpix1,wcs->crpix2);57 return region;58 }59 #endif60 13 61 14 // Activate/deactivate a single element for a list … … 140 93 psArray *regionHDUs = psArrayAlloc(data->numInputs); 141 94 142 psVector *exptimes = psVectorAlloc Empty(data->numInputs, PS_TYPE_F64);95 psVector *exptimes = psVectorAlloc(data->numInputs, PS_TYPE_F64); 143 96 144 97 // Determine which projection cells we have to deal with. … … 146 99 pmFPAfile *file = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.IMAGE", i); // File to examine 147 100 // Header in the FPA should have been read as a part of defining the file... 148 #if 0149 bool pmFPAReadHeaderSet(pmFPA *fpa, // FPA to read into150 psFits *fits, // FITS file from which to read151 pmConfig *config // Configuration152 );153 #endif154 101 pmHDU *hdu = file->fpa->hdu; // Header of interest 155 102 … … 218 165 for (int i = 0; i < numProj; i++) { 219 166 psRegion *projRegion = projRegions->data[i]; // Region for skycell projection 220 // projRegion->x0 += 5760;221 //projRegion->x1 += 5760;222 // projRegion->x0 +=223 // projRegion->x1 +=224 167 psTrace("ppSkycell", 2, "Projection %d: [%.0f:%.0f,%.0f:%.0f]\n", 225 168 i, projRegion->x0, projRegion->x1, projRegion->y0, projRegion->y1); 226 int xSize = projRegion->x1 - projRegion->x0 + 1; // Size of unbinned image 227 int ySize = projRegion->y1 - projRegion->y0 + 1; // Size of unbinned image 228 // Size of binned image 1 169 170 // Size of unbinned image 171 int xSize = projRegion->x1 - projRegion->x0 + 1; 172 int ySize = projRegion->y1 - projRegion->y0 + 1; 173 174 // Size of binned images 229 175 int numCols1 = xSize / (float)data->bin1 + 1.5, numRows1 = ySize / (float)data->bin1 + 1.5; 230 // Size of binned image 2231 176 int numCols2 = numCols1 / (float)data->bin2 + 1.5, numRows2 = numRows1 / (float)data->bin2 + 1.5; 232 177 233 psImage *image1 = psImageAlloc(numCols1, numRows1, PS_TYPE_F32); // Binned image 234 psImage *image2 = psImageAlloc(numCols2, numRows2, PS_TYPE_F32); // Binned image 178 // Binned image containers 179 psImage *image1 = psImageAlloc(numCols1, numRows1, PS_TYPE_F32); 180 psImage *image2 = psImageAlloc(numCols2, numRows2, PS_TYPE_F32); 235 181 psImageInit(image1,NAN); 236 182 psImageInit(image2,NAN); 183 184 // Binned image radius values. Used to determine primacy. 185 psImage *radius1 = psImageAlloc(numCols1, numRows1, PS_TYPE_F32); 186 psImage *radius2 = psImageAlloc(numCols2, numRows2, PS_TYPE_F32); 187 psImageInit(radius1,pow(numCols1 + numRows1,2)); // These values can be anything, just need to be larger than a binned radius. 188 psImageInit(radius2,pow(numCols1 + numRows1,2)); 237 189 238 psImage *mask1 = NULL, *mask2 = NULL; // Binned masks239 if (data->masksName) {240 mask1 = psImageAlloc(numCols1, numRows1, PS_TYPE_IMAGE_MASK);241 mask2 = psImageAlloc(numCols2, numRows2, PS_TYPE_IMAGE_MASK);242 psImageInit(mask1, 0xFF);243 psImageInit(mask2, 0xFF);244 }245 190 // HDU containing the WCS we plan on using 246 191 pmHDU *projhdu = NULL; … … 260 205 } 261 206 pmFPAfileActivateSingle(data->config->files, true, "PPSKYCELL.IMAGE", j); 262 if (data->masksName) {263 pmFPAfileActivateSingle(data->config->files, true, "PPSKYCELL.MASK", j);264 }265 207 266 208 pmFPAview *view = filesIterateDown(data->config); // View to readout … … 317 259 // Scale by exposure time 318 260 double exptime_factor = pow(exptimes->data.F64[j]/stats->robustMedian,data->exptimeOrder); 319 psBinaryOp(bin1RO->image,bin1RO->image,"/", 320 psScalarAlloc(exptime_factor, PS_TYPE_F64)); 321 psBinaryOp(bin2RO->image,bin2RO->image,"/", 322 psScalarAlloc(exptime_factor, PS_TYPE_F32)); 261 if ((isfinite(exptime_factor))&&(exptime_factor != 0.0)) { 262 psBinaryOp(bin1RO->image,bin1RO->image,"/", 263 psScalarAlloc(exptime_factor, PS_TYPE_F32)); 264 psBinaryOp(bin2RO->image,bin2RO->image,"/", 265 psScalarAlloc(exptime_factor, PS_TYPE_F32)); 266 } 323 267 // End scaling 324 #if (0) 325 psRegion *imageRegion = imageRegions->data[j]; // Region for image 326 fprintf(stderr,"%d IM %f %f %f %f PR %f %f %f %f\n", 327 j, 328 imageRegion->x0,imageRegion->x1, 329 imageRegion->y0,imageRegion->y1, 330 projRegion->x0,projRegion->x1, 331 projRegion->y0,projRegion->y1); 332 #endif 268 333 269 // Offsets for image on this projection cell are just differences in CRPIX positions. 334 270 int xOffset1 = (-1 * wcs->crpix1 - projRegion->x0) / (float)(data->bin1); 335 271 int yOffset1 = (-1 * wcs->crpix2 - projRegion->y0) / (float)(data->bin1); 336 272 337 int xOffset2 = xOffset1 / (float)data->bin2, yOffset2 = yOffset1 / (float)data->bin2; 273 int xOffset2 = xOffset1 / (float)data->bin2; 274 int yOffset2 = yOffset1 / (float)data->bin2; 338 275 psTrace("ppSkycell",5,"Offsets: %d %d : %d %d", 339 276 xOffset1,yOffset1,xOffset2,yOffset2); 340 277 341 // This is where any primary code would go. I'm unclear on how it would work. 342 278 // Check each pixel for primacy. A pixel is primary 279 // if it is closer to the central pixel of its skycell 280 // than any other pixel is to theirs. 281 282 // Let's just do the overlay here, instead of doing math, then handing it 283 // off to another function. That seems silly. 284 285 int u,v,x,y; 286 287 double u0 = xOffset1 + (bin1RO->image->numCols) / 2.0; 288 double v0 = yOffset1 + (bin1RO->image->numRows) / 2.0; 289 for (x = 0; x < bin1RO->image->numCols; x++) { 290 for (y = 0; y < bin1RO->image->numRows; y++) { 291 if (!isfinite(bin1RO->image->data.F32[y][x])) { 292 continue; 293 } 294 u = x + xOffset1; 295 v = y + yOffset1; 296 double R2 = pow(u - u0,2) + pow(v - v0,2); 297 298 if (R2 < radius1->data.F32[v][u]) { 299 radius1->data.F32[v][u] = R2; 300 image1->data.F32[v][u] = bin1RO->image->data.F32[y][x]; 301 } 302 } 303 } 304 305 u0 = xOffset2 + (bin2RO->image->numCols) / 2.0; 306 v0 = yOffset2 + (bin2RO->image->numRows) / 2.0; 307 for (x = 0; x < bin2RO->image->numCols; x++) { 308 for (y = 0; y < bin2RO->image->numRows; y++) { 309 if (!isfinite(bin2RO->image->data.F32[y][x])) { 310 continue; 311 } 312 u = x + xOffset2; 313 v = y + yOffset2; 314 double R2 = pow(u - u0,2) + pow(v - v0,2); 315 if (R2 < radius2->data.F32[v][u]) { 316 radius2->data.F32[v][u] = R2; 317 image2->data.F32[v][u] = bin2RO->image->data.F32[y][x]; 318 } 319 } 320 } 343 321 344 322 323 345 324 346 325 // Overlay the data onto the appropriate pixels in the final outputs 347 326 // XXX Completely neglecting rotations 348 327 // The skycells are divided up neatly with them all having the same orientation 349 psImageOverlaySection(image1, bin1RO->image, xOffset1, yOffset1, "E"); 350 psImageOverlaySection(image2, bin2RO->image, xOffset2, yOffset2, "E"); 351 if (data->masksName) { 352 psImageOverlaySection(mask1, bin1RO->mask, xOffset1, yOffset1, "M"); 353 psImageOverlaySection(mask2, bin2RO->mask, xOffset2, yOffset2, "M"); 354 } 328 // psImageOverlaySection(image1, bin1RO->image, xOffset1, yOffset1, "E"); 329 // psImageOverlaySection(image2, bin2RO->image, xOffset2, yOffset2, "E"); 355 330 356 331 // Cleanup on input loop. … … 360 335 psFree(file->fpa); 361 336 file->fpa = NULL; 362 if (data->masksName) {363 pmFPAfile *file = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.MASK", j);364 psFree(file->fpa);365 file->fpa = NULL;366 }367 337 pmFPAfileActivate(data->config->files, false, NULL); 368 338 } … … 373 343 pmFPAfileActivate(data->config->files, true, "PPSKYCELL.BIN1"); 374 344 pmFPAfileActivate(data->config->files, true, "PPSKYCELL.BIN2"); 375 if (data->masksName) {376 pmFPAfileActivate(data->config->files, true, "PPSKYCELL.BIN1.MASK");377 pmFPAfileActivate(data->config->files, true, "PPSKYCELL.BIN2.MASK");378 }379 345 } 380 346 … … 389 355 ro1->image = image1; 390 356 ro2->image = image2; 391 ro1->mask = mask1;392 ro2->mask = mask2;393 357 394 358 ro1->data_exists = cell1->data_exists = cell1->parent->data_exists = true; … … 409 373 fits1->fpa->hdu->header = psMetadataAlloc(); 410 374 psMetadataCopy(fits1->fpa->hdu->header,projhdu->header); 411 375 412 376 // Change wcs here 413 #define WCS_DEBUG 0414 377 if (modify_wcs1) { 415 378 pmAstromWCS *WCS = pmAstromWCSfromHeader(fits1->fpa->hdu->header); 416 379 double cd1f = 1.0 * data->bin1; 417 380 double cd2f = 1.0 * data->bin1; 418 #if WCS_DEBUG 419 fprintf(stderr,">>> %d %d (%g %g) (%g %g) (%g %g)\n",data->bin1,data->bin2, 420 cd1f,cd2f,WCS->cdelt1,WCS->cdelt2, 421 WCS->crpix1,WCS->crpix2 422 ); 423 #endif 381 424 382 WCS->cdelt1 *= cd1f; 425 383 WCS->cdelt2 *= cd2f; … … 437 395 WCS->crpix2 = WCS->crpix2 / cd2f; 438 396 } 439 #if WCS_DEBUG 440 fprintf(stderr,">>> %d %d (%g %g) (%g %g) (%g %g) %d %d %d %d\n",data->bin1,data->bin2, 441 cd1f,cd2f,WCS->cdelt1,WCS->cdelt2, 442 WCS->crpix1,WCS->crpix2, 443 WCS->trans->x->nX, WCS->trans->x->nY, 444 WCS->trans->y->nX, WCS->trans->y->nY 445 ); 446 #endif 397 447 398 for (int q = 0; q <= WCS->trans->x->nX; q++) { 448 399 for (int r = 0; r <= WCS->trans->x->nY; r++) { 449 400 WCS->trans->x->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r); 450 #if WCS_DEBUG451 fprintf(stderr,"PC: %d %d %g %g\n",452 q,r,pow(cd1f,q) * pow(cd2f,r),453 WCS->trans->x->coeff[q][r]);454 #endif455 401 } 456 402 } … … 458 404 for (int r = 0; r <= WCS->trans->y->nY; r++) { 459 405 WCS->trans->y->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r); 460 #if WCS_DEBUG461 fprintf(stderr,"PC: %d %d %g %g\n",462 q,r,pow(cd1f,q) * pow(cd2f,r),463 WCS->trans->y->coeff[q][r]);464 #endif465 406 } 466 407 } 467 408 pmAstromWCStoHeader (fits1->fpa->hdu->header,WCS); 468 #if WCS_DEBUG469 WCS = pmAstromWCSfromHeader(fits1->fpa->hdu->header);470 fprintf(stderr,">>> %d %d (%g %g) (%g %g) (%g %g)\n",data->bin1,data->bin2,471 cd1f,cd2f,WCS->cdelt1,WCS->cdelt2,472 WCS->crpix1,WCS->crpix2473 );474 #endif475 409 modify_wcs1 = 0; 476 410 } … … 482 416 483 417 pmCell *Fcell1 = pmFPAfileThisCell(data->config->files, view, "PPSKYCELL.BIN1"); // Rebinned cell 1 484 /* // This is a hack to get a functioning header created so the fits images can be written out. */ 418 // This is a hack to get a functioning header created so the fits images can be written out. 485 419 psMetadataAddS32(Fcell1->concepts,PS_LIST_TAIL,"CELL.XPARITY", PS_META_REPLACE,"",1); 486 420 psMetadataAddS32(Fcell1->concepts,PS_LIST_TAIL,"CELL.YPARITY", PS_META_REPLACE,"",1); 487 421 psMetadataAddS32(Fcell1->concepts,PS_LIST_TAIL,"CELL.READDIR", PS_META_REPLACE,"",1); 488 psMetadataAddF64(Fcell1->concepts,PS_LIST_TAIL,"EXPTIME", PS_META_REPLACE, "", stats->robustMedian); 422 423 // I am baffled that this is the way to get the exposure time updated correctly. 424 psMetadataItem *item1; 425 item1 = psMetadataLookup(Fcell1->concepts,"CELL.EXPOSURE"); 426 item1->data.F32 = stats->robustMedian; 427 489 428 pmReadout *Fro1 = pmReadoutAlloc(Fcell1); 490 429 Fro1->image = image1; 491 Fro1->mask = mask1;492 430 Fro1->data_exists = Fcell1->data_exists = Fcell1->parent->data_exists = true; 493 431 … … 495 433 fits1->fileIndex = i; 496 434 497 // Do the mask if we need to498 if (data->masksName) {499 /* pmFPAfile *Mask1 = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.BIN1.MASK", 0); */500 /* // Mask1->fpa->hdu = pmHDUAlloc(NULL); */501 /* // Mask1->fpa->hdu->header = psMetadataAlloc(); */502 /* // psMetadataCopy(Mask1->fpa->hdu->header,fits1->fpa->hdu->header); */503 /* pmChip *Mchip1 = pmFPAfileThisChip(data->config->files, view, "PPSKYCELL.BIN1.MASK"); */504 /* psMetadataAddS32(Mchip1->concepts,PS_LIST_TAIL,"CHIP.XPARITY", PS_META_REPLACE,"",1); */505 /* psMetadataAddS32(Mchip1->concepts,PS_LIST_TAIL,"CHIP.YPARITY", PS_META_REPLACE,"",1); */506 507 /* pmCell *Mcell1 = pmFPAfileThisCell(data->config->files, view, "PPSKYCELL.BIN1.MASK"); // Rebinned cell 1 */508 /* /\* // This is a hack to get a functioning header created so the fits images can be written out. *\/ */509 /* psMetadataAddS32(Mcell1->concepts,PS_LIST_TAIL,"CELL.XPARITY", PS_META_REPLACE,"",1); */510 /* psMetadataAddS32(Mcell1->concepts,PS_LIST_TAIL,"CELL.YPARITY", PS_META_REPLACE,"",1); */511 /* psMetadataAddS32(Mcell1->concepts,PS_LIST_TAIL,"CELL.READDIR", PS_META_REPLACE,"",1); */512 513 /* pmReadout *Mro1 = pmReadoutAlloc(Mcell1); */514 /* Mro1->image = image1; */515 /* Mro1->mask = mask1; */516 /* Mro1->data_exists = Mcell1->data_exists = Mcell1->parent->data_exists = true; */517 518 /* Mask1->save = true; */519 /* Mask1->fileIndex = i; */520 }521 522 523 524 435 // Repeat with second binned image 525 436 pmFPAfile *fits2 = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.BIN2", 0); … … 527 438 fits2->fpa->hdu->header = psMetadataAlloc(); 528 439 psMetadataCopy(fits2->fpa->hdu->header,projhdu->header); 529 440 530 441 if (modify_wcs2) { 531 442 pmAstromWCS *WCS = pmAstromWCSfromHeader(fits2->fpa->hdu->header); … … 561 472 psMetadataAddS32(Fcell2->concepts,PS_LIST_TAIL,"CELL.YPARITY", PS_META_REPLACE,"",1); 562 473 psMetadataAddS32(Fcell2->concepts,PS_LIST_TAIL,"CELL.READDIR", PS_META_REPLACE,"",1); 563 psMetadataAddF64(Fcell2->concepts,PS_LIST_TAIL,"EXPTIME", PS_META_REPLACE, "", stats->robustMedian); 474 475 psMetadataItem *item2 = psMetadataLookup(Fcell2->concepts,"CELL.EXPOSURE"); 476 item2->data.F32 = stats->robustMedian; 564 477 565 478 pmReadout *Fro2 = pmReadoutAlloc(Fcell2); 566 479 Fro2->image = image2; 567 Fro2->mask = mask2;568 480 Fro2->data_exists = Fcell2->data_exists = Fcell2->parent->data_exists = true; 569 481 … … 575 487 576 488 psFree(view); 577 #if 0578 {579 psString filename = NULL; // Filename for image580 psStringAppend(&filename, "skycell_%d_1.fits", i);581 psFits *fits = psFitsOpen(filename, "w");582 psFree(filename);583 psFitsWriteImage(fits, NULL, image1, 0, NULL);584 psFitsClose(fits);585 }586 587 {588 psString filename = NULL; // Filename for image589 psStringAppend(&filename, "skycell_%d_2.fits", i);590 psFits *fits = psFitsOpen(filename, "w");591 psFree(filename);592 psFitsWriteImage(fits, NULL, image2, 0, NULL);593 psFitsClose(fits);594 }595 #endif596 489 597 490 filesIterateUp(data->config);
Note:
See TracChangeset
for help on using the changeset viewer.
