- Timestamp:
- Dec 10, 2013, 2:09:59 PM (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/eam_branches/ipp-20130904/ppSkycell/src/ppSkycellLoop.c
r35852 r36368 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 … … 139 92 psArray *imageRegions = psArrayAlloc(data->numInputs); // Region for image 140 93 psArray *regionHDUs = psArrayAlloc(data->numInputs); 94 95 psVector *exptimes = psVectorAlloc(data->numInputs, PS_TYPE_F64); 96 141 97 // Determine which projection cells we have to deal with. 142 98 for (int i = 0; i < data->numInputs; i++) { 143 99 pmFPAfile *file = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.IMAGE", i); // File to examine 144 100 // Header in the FPA should have been read as a part of defining the file... 145 #if 0146 bool pmFPAReadHeaderSet(pmFPA *fpa, // FPA to read into147 psFits *fits, // FITS file from which to read148 pmConfig *config // Configuration149 );150 #endif151 101 pmHDU *hdu = file->fpa->hdu; // Header of interest 152 102 … … 187 137 } 188 138 } 189 190 if (!found) { 139 140 if (!found) { // Add new projection cell if we didn't find one. 191 141 psVectorAppend(crval1, wcs->crval1); 192 142 psVectorAppend(crval2, wcs->crval2); 193 143 psVectorAppend(cdelt1, wcs->cdelt1); 194 144 psVectorAppend(cdelt2, wcs->cdelt2); 145 195 146 psRegion *projRegion = psRegionAlloc(region->x0, region->x1, region->y0, region->y1); 196 147 psArrayAdd(projRegions, projRegions->n, projRegion); … … 200 151 psArrayAdd(regionHDUs,1,hdu); 201 152 numProj++; 153 202 154 } 203 } 204 155 156 // Save exptime 157 exptimes->data.F64[i] = psMetadataLookupF64(NULL, hdu->header, "EXPTIME"); 158 } 159 160 float exptime_target = 1.0; 161 162 // psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); 163 // psVectorStats(stats,exptimes,NULL, NULL, 0); 164 // exptime_target = stats->robustMedian; 165 205 166 pmFPAfileActivate(data->config->files, false, NULL); 206 167 // Loop over projections 207 168 for (int i = 0; i < numProj; i++) { 208 169 psRegion *projRegion = projRegions->data[i]; // Region for skycell projection 209 // projRegion->x0 += 5760;210 //projRegion->x1 += 5760;211 // projRegion->x0 +=212 // projRegion->x1 +=213 170 psTrace("ppSkycell", 2, "Projection %d: [%.0f:%.0f,%.0f:%.0f]\n", 214 171 i, projRegion->x0, projRegion->x1, projRegion->y0, projRegion->y1); 215 int xSize = projRegion->x1 - projRegion->x0 + 1; // Size of unbinned image 216 int ySize = projRegion->y1 - projRegion->y0 + 1; // Size of unbinned image 217 // Size of binned image 1 172 173 // Size of unbinned image 174 int xSize = projRegion->x1 - projRegion->x0 + 1; 175 int ySize = projRegion->y1 - projRegion->y0 + 1; 176 177 // Size of binned images 218 178 int numCols1 = xSize / (float)data->bin1 + 1.5, numRows1 = ySize / (float)data->bin1 + 1.5; 219 // Size of binned image 2220 179 int numCols2 = numCols1 / (float)data->bin2 + 1.5, numRows2 = numRows1 / (float)data->bin2 + 1.5; 221 180 222 psImage *image1 = psImageAlloc(numCols1, numRows1, PS_TYPE_F32); // Binned image 223 psImage *image2 = psImageAlloc(numCols2, numRows2, PS_TYPE_F32); // Binned image 181 // Binned image containers 182 psImage *image1 = psImageAlloc(numCols1, numRows1, PS_TYPE_F32); 183 psImage *image2 = psImageAlloc(numCols2, numRows2, PS_TYPE_F32); 224 184 psImageInit(image1,NAN); 225 185 psImageInit(image2,NAN); 186 187 // Binned image radius values. Used to determine primacy. 188 psImage *radius1 = psImageAlloc(numCols1, numRows1, PS_TYPE_F32); 189 psImage *radius2 = psImageAlloc(numCols2, numRows2, PS_TYPE_F32); 190 psImageInit(radius1,pow(numCols1 + numRows1,2)); // These values can be anything, just need to be larger than a binned radius. 191 psImageInit(radius2,pow(numCols1 + numRows1,2)); 226 192 227 psImage *mask1 = NULL, *mask2 = NULL; // Binned masks228 if (data->masksName) {229 mask1 = psImageAlloc(numCols1, numRows1, PS_TYPE_IMAGE_MASK);230 mask2 = psImageAlloc(numCols2, numRows2, PS_TYPE_IMAGE_MASK);231 psImageInit(mask1, 0xFF);232 psImageInit(mask2, 0xFF);233 }234 193 // HDU containing the WCS we plan on using 235 194 pmHDU *projhdu = NULL; … … 249 208 } 250 209 pmFPAfileActivateSingle(data->config->files, true, "PPSKYCELL.IMAGE", j); 251 if (data->masksName) {252 pmFPAfileActivateSingle(data->config->files, true, "PPSKYCELL.MASK", j);253 }254 210 255 211 pmFPAview *view = filesIterateDown(data->config); // View to readout … … 304 260 } 305 261 306 #if (0) 307 psRegion *imageRegion = imageRegions->data[j]; // Region for image 308 fprintf(stderr,"%d IM %f %f %f %f PR %f %f %f %f\n", 309 j, 310 imageRegion->x0,imageRegion->x1, 311 imageRegion->y0,imageRegion->y1, 312 projRegion->x0,projRegion->x1, 313 projRegion->y0,projRegion->y1); 314 #endif 262 // Scale by exposure time 263 double exptime_factor = pow(exptimes->data.F64[j]/exptime_target,data->exptimeOrder); 264 if ((isfinite(exptime_factor))&&(exptime_factor != 0.0)) { 265 psBinaryOp(bin1RO->image,bin1RO->image,"/", 266 psScalarAlloc(exptime_factor, PS_TYPE_F32)); 267 psBinaryOp(bin2RO->image,bin2RO->image,"/", 268 psScalarAlloc(exptime_factor, PS_TYPE_F32)); 269 } 270 // End scaling 271 315 272 // Offsets for image on this projection cell are just differences in CRPIX positions. 316 273 int xOffset1 = (-1 * wcs->crpix1 - projRegion->x0) / (float)(data->bin1); 317 274 int yOffset1 = (-1 * wcs->crpix2 - projRegion->y0) / (float)(data->bin1); 318 275 319 int xOffset2 = xOffset1 / (float)data->bin2, yOffset2 = yOffset1 / (float)data->bin2; 276 int xOffset2 = xOffset1 / (float)data->bin2; 277 int yOffset2 = yOffset1 / (float)data->bin2; 320 278 psTrace("ppSkycell",5,"Offsets: %d %d : %d %d", 321 279 xOffset1,yOffset1,xOffset2,yOffset2); 322 280 281 // Check each pixel for primacy. A pixel is primary 282 // if it is closer to the central pixel of its skycell 283 // than any other pixel is to theirs. 284 285 // Let's just do the overlay here, instead of doing math, then handing it 286 // off to another function. That seems silly. 287 288 int u,v,x,y; 289 290 double u0 = xOffset1 + (bin1RO->image->numCols) / 2.0; 291 double v0 = yOffset1 + (bin1RO->image->numRows) / 2.0; 292 for (x = 0; x < bin1RO->image->numCols; x++) { 293 for (y = 0; y < bin1RO->image->numRows; y++) { 294 if (!isfinite(bin1RO->image->data.F32[y][x])) { 295 continue; 296 } 297 u = x + xOffset1; 298 v = y + yOffset1; 299 double R2 = pow(u - u0,2) + pow(v - v0,2); 300 301 if (R2 < radius1->data.F32[v][u]) { 302 radius1->data.F32[v][u] = R2; 303 image1->data.F32[v][u] = bin1RO->image->data.F32[y][x]; 304 } 305 } 306 } 307 308 u0 = xOffset2 + (bin2RO->image->numCols) / 2.0; 309 v0 = yOffset2 + (bin2RO->image->numRows) / 2.0; 310 for (x = 0; x < bin2RO->image->numCols; x++) { 311 for (y = 0; y < bin2RO->image->numRows; y++) { 312 if (!isfinite(bin2RO->image->data.F32[y][x])) { 313 continue; 314 } 315 u = x + xOffset2; 316 v = y + yOffset2; 317 double R2 = pow(u - u0,2) + pow(v - v0,2); 318 if (R2 < radius2->data.F32[v][u]) { 319 radius2->data.F32[v][u] = R2; 320 image2->data.F32[v][u] = bin2RO->image->data.F32[y][x]; 321 } 322 } 323 } 324 325 326 327 323 328 // Overlay the data onto the appropriate pixels in the final outputs 324 329 // XXX Completely neglecting rotations 325 330 // The skycells are divided up neatly with them all having the same orientation 326 psImageOverlaySection(image1, bin1RO->image, xOffset1, yOffset1, "E"); 327 psImageOverlaySection(image2, bin2RO->image, xOffset2, yOffset2, "E"); 328 if (data->masksName) { 329 psImageOverlaySection(mask1, bin1RO->mask, xOffset1, yOffset1, "M"); 330 psImageOverlaySection(mask2, bin2RO->mask, xOffset2, yOffset2, "M"); 331 } 331 // psImageOverlaySection(image1, bin1RO->image, xOffset1, yOffset1, "E"); 332 // psImageOverlaySection(image2, bin2RO->image, xOffset2, yOffset2, "E"); 332 333 333 334 // Cleanup on input loop. … … 337 338 psFree(file->fpa); 338 339 file->fpa = NULL; 339 if (data->masksName) {340 pmFPAfile *file = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.MASK", j);341 psFree(file->fpa);342 file->fpa = NULL;343 }344 340 pmFPAfileActivate(data->config->files, false, NULL); 345 341 } … … 350 346 pmFPAfileActivate(data->config->files, true, "PPSKYCELL.BIN1"); 351 347 pmFPAfileActivate(data->config->files, true, "PPSKYCELL.BIN2"); 352 if (data->masksName) {353 pmFPAfileActivate(data->config->files, true, "PPSKYCELL.BIN1.MASK");354 pmFPAfileActivate(data->config->files, true, "PPSKYCELL.BIN2.MASK");355 }356 348 } 357 349 … … 366 358 ro1->image = image1; 367 359 ro2->image = image2; 368 ro1->mask = mask1;369 ro2->mask = mask2;370 360 371 361 ro1->data_exists = cell1->data_exists = cell1->parent->data_exists = true; … … 386 376 fits1->fpa->hdu->header = psMetadataAlloc(); 387 377 psMetadataCopy(fits1->fpa->hdu->header,projhdu->header); 388 378 389 379 // Change wcs here 390 #define WCS_DEBUG 0391 380 if (modify_wcs1) { 392 381 pmAstromWCS *WCS = pmAstromWCSfromHeader(fits1->fpa->hdu->header); 393 382 double cd1f = 1.0 * data->bin1; 394 383 double cd2f = 1.0 * data->bin1; 395 #if WCS_DEBUG 396 fprintf(stderr,">>> %d %d (%g %g) (%g %g) (%g %g)\n",data->bin1,data->bin2, 397 cd1f,cd2f,WCS->cdelt1,WCS->cdelt2, 398 WCS->crpix1,WCS->crpix2 399 ); 400 #endif 384 401 385 WCS->cdelt1 *= cd1f; 402 386 WCS->cdelt2 *= cd2f; … … 414 398 WCS->crpix2 = WCS->crpix2 / cd2f; 415 399 } 416 #if WCS_DEBUG 417 fprintf(stderr,">>> %d %d (%g %g) (%g %g) (%g %g) %d %d %d %d\n",data->bin1,data->bin2, 418 cd1f,cd2f,WCS->cdelt1,WCS->cdelt2, 419 WCS->crpix1,WCS->crpix2, 420 WCS->trans->x->nX, WCS->trans->x->nY, 421 WCS->trans->y->nX, WCS->trans->y->nY 422 ); 423 #endif 400 424 401 for (int q = 0; q <= WCS->trans->x->nX; q++) { 425 402 for (int r = 0; r <= WCS->trans->x->nY; r++) { 426 403 WCS->trans->x->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r); 427 #if WCS_DEBUG428 fprintf(stderr,"PC: %d %d %g %g\n",429 q,r,pow(cd1f,q) * pow(cd2f,r),430 WCS->trans->x->coeff[q][r]);431 #endif432 404 } 433 405 } … … 435 407 for (int r = 0; r <= WCS->trans->y->nY; r++) { 436 408 WCS->trans->y->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r); 437 #if WCS_DEBUG438 fprintf(stderr,"PC: %d %d %g %g\n",439 q,r,pow(cd1f,q) * pow(cd2f,r),440 WCS->trans->y->coeff[q][r]);441 #endif442 409 } 443 410 } 444 411 pmAstromWCStoHeader (fits1->fpa->hdu->header,WCS); 445 #if WCS_DEBUG446 WCS = pmAstromWCSfromHeader(fits1->fpa->hdu->header);447 fprintf(stderr,">>> %d %d (%g %g) (%g %g) (%g %g)\n",data->bin1,data->bin2,448 cd1f,cd2f,WCS->cdelt1,WCS->cdelt2,449 WCS->crpix1,WCS->crpix2450 );451 #endif452 412 modify_wcs1 = 0; 453 413 } … … 459 419 460 420 pmCell *Fcell1 = pmFPAfileThisCell(data->config->files, view, "PPSKYCELL.BIN1"); // Rebinned cell 1 461 /* // This is a hack to get a functioning header created so the fits images can be written out. */ 421 // This is a hack to get a functioning header created so the fits images can be written out. 462 422 psMetadataAddS32(Fcell1->concepts,PS_LIST_TAIL,"CELL.XPARITY", PS_META_REPLACE,"",1); 463 423 psMetadataAddS32(Fcell1->concepts,PS_LIST_TAIL,"CELL.YPARITY", PS_META_REPLACE,"",1); 464 424 psMetadataAddS32(Fcell1->concepts,PS_LIST_TAIL,"CELL.READDIR", PS_META_REPLACE,"",1); 465 425 426 // I am baffled that this is the way to get the exposure time updated correctly. 427 psMetadataItem *item1; 428 item1 = psMetadataLookup(Fcell1->concepts,"CELL.EXPOSURE"); 429 item1->data.F32 = exptime_target; 430 466 431 pmReadout *Fro1 = pmReadoutAlloc(Fcell1); 467 432 Fro1->image = image1; 468 Fro1->mask = mask1;469 433 Fro1->data_exists = Fcell1->data_exists = Fcell1->parent->data_exists = true; 470 434 … … 472 436 fits1->fileIndex = i; 473 437 474 // Do the mask if we need to475 if (data->masksName) {476 /* pmFPAfile *Mask1 = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.BIN1.MASK", 0); */477 /* // Mask1->fpa->hdu = pmHDUAlloc(NULL); */478 /* // Mask1->fpa->hdu->header = psMetadataAlloc(); */479 /* // psMetadataCopy(Mask1->fpa->hdu->header,fits1->fpa->hdu->header); */480 /* pmChip *Mchip1 = pmFPAfileThisChip(data->config->files, view, "PPSKYCELL.BIN1.MASK"); */481 /* psMetadataAddS32(Mchip1->concepts,PS_LIST_TAIL,"CHIP.XPARITY", PS_META_REPLACE,"",1); */482 /* psMetadataAddS32(Mchip1->concepts,PS_LIST_TAIL,"CHIP.YPARITY", PS_META_REPLACE,"",1); */483 484 /* pmCell *Mcell1 = pmFPAfileThisCell(data->config->files, view, "PPSKYCELL.BIN1.MASK"); // Rebinned cell 1 */485 /* /\* // This is a hack to get a functioning header created so the fits images can be written out. *\/ */486 /* psMetadataAddS32(Mcell1->concepts,PS_LIST_TAIL,"CELL.XPARITY", PS_META_REPLACE,"",1); */487 /* psMetadataAddS32(Mcell1->concepts,PS_LIST_TAIL,"CELL.YPARITY", PS_META_REPLACE,"",1); */488 /* psMetadataAddS32(Mcell1->concepts,PS_LIST_TAIL,"CELL.READDIR", PS_META_REPLACE,"",1); */489 490 /* pmReadout *Mro1 = pmReadoutAlloc(Mcell1); */491 /* Mro1->image = image1; */492 /* Mro1->mask = mask1; */493 /* Mro1->data_exists = Mcell1->data_exists = Mcell1->parent->data_exists = true; */494 495 /* Mask1->save = true; */496 /* Mask1->fileIndex = i; */497 }498 499 500 501 438 // Repeat with second binned image 502 439 pmFPAfile *fits2 = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.BIN2", 0); … … 504 441 fits2->fpa->hdu->header = psMetadataAlloc(); 505 442 psMetadataCopy(fits2->fpa->hdu->header,projhdu->header); 506 443 507 444 if (modify_wcs2) { 508 445 pmAstromWCS *WCS = pmAstromWCSfromHeader(fits2->fpa->hdu->header); … … 538 475 psMetadataAddS32(Fcell2->concepts,PS_LIST_TAIL,"CELL.YPARITY", PS_META_REPLACE,"",1); 539 476 psMetadataAddS32(Fcell2->concepts,PS_LIST_TAIL,"CELL.READDIR", PS_META_REPLACE,"",1); 477 478 psMetadataItem *item2 = psMetadataLookup(Fcell2->concepts,"CELL.EXPOSURE"); 479 item2->data.F32 = exptime_target; 540 480 541 481 pmReadout *Fro2 = pmReadoutAlloc(Fcell2); 542 482 Fro2->image = image2; 543 Fro2->mask = mask2;544 483 Fro2->data_exists = Fcell2->data_exists = Fcell2->parent->data_exists = true; 545 484 … … 551 490 552 491 psFree(view); 553 #if 0554 {555 psString filename = NULL; // Filename for image556 psStringAppend(&filename, "skycell_%d_1.fits", i);557 psFits *fits = psFitsOpen(filename, "w");558 psFree(filename);559 psFitsWriteImage(fits, NULL, image1, 0, NULL);560 psFitsClose(fits);561 }562 563 {564 psString filename = NULL; // Filename for image565 psStringAppend(&filename, "skycell_%d_2.fits", i);566 psFits *fits = psFitsOpen(filename, "w");567 psFree(filename);568 psFitsWriteImage(fits, NULL, image2, 0, NULL);569 psFitsClose(fits);570 }571 #endif572 492 573 493 filesIterateUp(data->config);
Note:
See TracChangeset
for help on using the changeset viewer.
