IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Apr 7, 2013, 6:23:04 AM (13 years ago)
Author:
eugene
Message:

cleanup code to be more readable (add pswarpStatsFile structure, pswarpCleanup, move code into pswarpDefineBackground and pswarpLoopBackground

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20130307/pswarp/src/pswarpLoop.c

    r34800 r35373  
    11/** @file pswarpLoop.c
    22 *
    3  *  @brief
    4  *
     3 *  @brief mail processing loop for pswarp
    54 *  @ingroup pswarp
    65 *
    … …  
    1716#define WCS_NONLIN_TOL 0.001            // Non-linear tolerance for header WCS
    1817#define TESTING 0                       // Testing output?
    19 
    20 
    2118
    2219// Loop over the inputs, warp them to the output skycell and then write out the output.
    … …  
    155152            psError(psErrorCodeLast(), false, "Unable to read bilevel mosaic astrometry for input FPA.");
    156153            psFree(view);
    157             psFree(stats);
    158154            goto DONE;
    159155        }
    … …  
    185181                psError(psErrorCodeLast(), false, "Unable to read bilevel chip astrometry for input FPA.");
    186182                psFree(view);
    187                 psFree(stats);
    188183                goto DONE;
    189184            }
    … …  
    193188                psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for input FPA.");
    194189                psFree(view);
    195                 psFree(stats);
    196190                goto DONE;
    197191            }
    … …  
    301295    if (!pmConceptsAverageCells(outCell, cells, NULL, NULL, false)) {
    302296        psError(psErrorCodeLast(), false, "Unable to average cell concepts.");
    303         psFree(stats);
    304297        psFree(cells);
    305298        psFree(view);
    … …  
    313306    if (!psMetadataCopy(outFPA->concepts, input->fpa->concepts)) {
    314307        psError(psErrorCodeLast(), false, "Unable to copy FPA concepts from input to output.");
    315         psFree(stats);
    316308        psFree(view);
    317309        goto DONE;
    … …  
    345337    if (!pmAstromWriteWCS(hdu->header, outFPA, outChip, WCS_NONLIN_TOL)) {
    346338        psError(psErrorCodeLast(), false, "Unable to generate WCS header.");
    347         psFree(stats);
    348339        goto DONE;
    349340    }
    … …  
    452443    return true;
    453444}
    454 
    455 // Loop over the inputs, warp them to the output skycell and then write out the output.
    456 bool pswarpLoopBackground(pmConfig *config, psMetadata *stats)
    457 {
    458     bool status;
    459     bool mdok;                          // Status of MD lookup
    460     const char *skyCamera = psMetadataLookupStr(NULL, config->arguments,
    461                                                 "SKYCELL.CAMERA");  // Name of camera for skycell
    462     pmConfigCamerasCull(config, skyCamera);
    463     pmConfigRecipesCull(config, "PSWARP,PPSTATS,PSPHOT,PSASTRO,MASKS,JPEG");
    464 
    465     // load the recipe
    466     psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSWARP_RECIPE);
    467     if (!recipe) {
    468         psError(PSWARP_ERR_CONFIG, false, "missing recipe %s", PSWARP_RECIPE);
    469         return false;
    470     }
    471 
    472     if (!pswarpSetMaskBits(config)) {
    473         psError(psErrorCodeLast(), false, "failed to set mask bits");
    474         return NULL;
    475     }
    476 
    477     // select the input data sources
    478     pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PSWARP.BKGMODEL");
    479     if (!input) {
    480         psError(PSWARP_ERR_CONFIG, true, "Can't find input data!\n");
    481         return false;
    482     }
    483 
    484     // use the external astrometry source if supplied
    485     pmFPAfile *astrom = psMetadataLookupPtr(NULL, config->files, "PSWARP.ASTROM");
    486     if (!astrom) {
    487         astrom = input;
    488     }
    489 
    490     if (astrom->camera != input->camera) {
    491         psError(PSWARP_ERR_DATA, true, "Input camera and astrometry camera do not match.");
    492         return false;
    493     }
    494 
    495     // select the output readout
    496     pmFPAview *view = pmFPAviewAlloc(0);
    497     view->chip = 0;
    498     view->cell = 0;
    499     view->readout = 0;
    500     pmReadout *output = pmFPAfileThisReadout(config->files, view, "PSWARP.OUTPUT.BKGMODEL");
    501     if (!output) {
    502         psError(PSWARP_ERR_CONFIG, true, "Can't find output background data!\n");
    503         return false;
    504     }
    505     psFree (view);
    506     // Turn all skycell files on to generate them, and then turn them off for the loop over the input images
    507     // the input, which is in a different format.
    508     {
    509         pswarpFileActivation(config, detectorFiles, false);
    510         pswarpFileActivation(config, photFiles, false);
    511         pswarpFileActivation(config, independentFiles, false);
    512         pswarpFileActivation(config, skycellFiles, true);
    513         if (!pswarpIOChecksBefore(config)) {
    514             psError(psErrorCodeLast(), false, "Unable to read files.");
    515             goto DONE;
    516         }
    517         pswarpFileActivation(config, skycellFiles, false);
    518     }
    519     // Read the input astrometry
    520     // XXX rather than use the activations here, this should just explicitly loop over the desired filerule
    521     {
    522 
    523       pmFPAfileActivate(config->files, true, "PSWARP.ASTROM");
    524 
    525         pmChip *chip;
    526         pmFPAview *view = pmFPAviewAlloc(0);
    527         if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
    528             psError(psErrorCodeLast(), false, "Unable to read files.");
    529             goto DONE;
    530         }
    531 
    532         while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
    533 #if 0
    534           // This needs to be removed because otherwise it throws an error of duplicate PSPHOT.DETECTIONS.
    535             if (!chip->process || !chip->file_exists) { continue; }
    536             if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
    537                 psError(psErrorCodeLast(), false, "Unable to read files.");
    538                 goto DONE;
    539             }
    540 #endif
    541             pmCell *cell;
    542             while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
    543               psTrace ("pswarp", 4, "ACell %d: %x %x %d\n", view->cell, cell->file_exists, cell->process,psErrorCodeLast());
    544                 if (!cell->process || !cell->file_exists) { continue; }
    545                 if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE) ||
    546                     !pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) {
    547                     psError(psErrorCodeLast(), false, "Unable to read files.");
    548                     goto DONE;
    549                 }
    550             }
    551             if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) {
    552                 psError(psErrorCodeLast(), false, "Unable to write files.");
    553                 goto DONE;
    554             }
    555         }
    556         if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) {
    557             psError(psErrorCodeLast(), false, "Unable to write files.");
    558             goto DONE;
    559         }
    560         psFree(view);
    561         pswarpFileActivation(config, detectorFiles, true);
    562         pmFPAfileActivate(config->files, false, "PSWARP.ASTROM");
    563     }
    564 
    565     // Don't care about the skycell anymore --- we've read it, and that's all we need to do.
    566     pmFPAfileActivate(config->files, false, "PSWARP.SKYCELL");
    567     view = pmFPAviewAlloc(0);
    568 
    569     // find the FPA phu
    570     bool bilevelAstrometry = false;
    571     pmHDU *phu = pmFPAviewThisPHU(view, astrom->fpa);
    572 
    573     //    pmAstromWCS *WCSF = pmAstromWCSfromHeader(phu->header);
    574    
    575     if (phu) {
    576         char *ctype = psMetadataLookupStr(NULL, phu->header, "CTYPE1");
    577         if (ctype) {
    578             bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
    579         }
    580     }
    581     if (bilevelAstrometry) {
    582         if (!pmAstromReadBilevelMosaic(input->fpa, phu->header)) {
    583             psError(psErrorCodeLast(), false, "Unable to read bilevel mosaic astrometry for input FPA.");
    584             psFree(view);
    585             psFree(stats);
    586             goto DONE;
    587         }
    588     }
    589 
    590     psList *cells = psListAlloc(NULL);  // List of cells, for concepts averaging
    591 
    592     // files associated with the science image
    593     if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
    594         psError(psErrorCodeLast(), false, "Unable to read files.");
    595         goto DONE;
    596     }
    597     pmChip *chip;
    598     while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
    599         psTrace ("pswarp", 4, "DChip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
    600         if (!chip->process || !chip->file_exists) { continue; }
    601         if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
    602             psError(psErrorCodeLast(), false, "Unable to read files.");
    603             goto DONE;
    604         }
    605 
    606         // Pull information from the header of the background files so we can use it to set values.
    607         pmHDU *hdu = pmFPAviewThisHDU(view,input->fpa);
    608         psMetadata *header = hdu->header;
    609        
    610         int IMAXIS1 = psMetadataLookupS32(NULL,header,"IMNAXIS1");
    611         int IMAXIS2 = psMetadataLookupS32(NULL,header,"IMNAXIS2");
    612         int NAXIS1 = psMetadataLookupS32(NULL,header,"NAXIS1");
    613         int NAXIS2 = psMetadataLookupS32(NULL,header,"NAXIS2");
    614         char *CCDSUM = psMetadataLookupStr(NULL,header,"CCDSUM");
    615         int CCDSUM1 = atoi(strtok(CCDSUM," "));
    616         int CCDSUM2 = atoi(strtok(NULL," "));
    617        
    618         psMetadataAddF32(config->arguments,PS_LIST_TAIL,"BKG_WARP_XOFFSET", PS_META_REPLACE,
    619                          "xoffset for background model data", (NAXIS1 * CCDSUM1 - IMAXIS1) / (2.0 * CCDSUM1));
    620         psMetadataAddF32(config->arguments,PS_LIST_TAIL,"BKG_WARP_YOFFSET", PS_META_REPLACE,
    621                          "yoffset for background model data", (NAXIS2 * CCDSUM2 - IMAXIS2) / (2.0 * CCDSUM2));
    622         psTrace("pswarp",5,"%d %d %d %d %d %d %g %g %d %d",
    623                 psMetadataLookupS32(NULL,header,"IMNAXIS1"),
    624                 psMetadataLookupS32(NULL,header,"IMNAXIS2"),
    625                 psMetadataLookupS32(NULL,header,"NAXIS1"),
    626                 psMetadataLookupS32(NULL,header,"NAXIS2"),
    627                 CCDSUM1,CCDSUM2,
    628                 psMetadataLookupF32(NULL,config->arguments,"BKG_WARP_XOFFSET"),
    629                 psMetadataLookupF32(NULL,config->arguments,"BKG_WARP_YOFFSET"),
    630                 psMetadataLookupS32(NULL,config->arguments,"BKG.XGRID"),
    631                 psMetadataLookupS32(NULL,config->arguments,"BKG.YGRID"));
    632        
    633        
    634         // read WCS data from the corresponding header
    635         hdu = pmFPAviewThisHDU (view, astrom->fpa);
    636 
    637         pmAstromWCS *WCS = pmAstromWCSfromHeader(hdu->header);
    638 
    639         double cd1f = (1.0 * CCDSUM1);// * (1.0 * psMetadataLookupS32(NULL,config->arguments,"BKG.XGRID"));
    640         double cd2f = (1.0 * CCDSUM2);// * (1.0 * psMetadataLookupS32(NULL,config->arguments,"BKG.YGRID"));
    641 
    642         WCS->cdelt1 *= cd1f;
    643         WCS->cdelt2 *= cd2f;
    644         WCS->crpix1 = WCS->crpix1 / cd1f;
    645         WCS->crpix2 = WCS->crpix2 / cd2f;
    646 
    647         // WCS->trans->x->nX/nY
    648         for (int q = 0; q <= WCS->trans->x->nX; q++) {
    649           for (int r = 0; r <= WCS->trans->x->nY; r++) {
    650             WCS->trans->x->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r);
    651           }
    652         }
    653         for (int q = 0; q <= WCS->trans->y->nX; q++) {
    654           for (int r = 0; r <= WCS->trans->y->nY; r++) {
    655             WCS->trans->y->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r);
    656           }
    657         }
    658         pmAstromWCStoHeader (hdu->header,WCS);
    659        
    660         if (bilevelAstrometry) {
    661             if (!pmAstromReadBilevelChip (chip, hdu->header)) {
    662                 psError(psErrorCodeLast(), false, "Unable to read bilevel chip astrometry for input FPA.");
    663                 psFree(view);
    664                 psFree(stats);
    665                 goto DONE;
    666             }
    667         } else {
    668             // we use a default FPA pixel scale of 1.0
    669             if (!pmAstromReadWCS (astrom->fpa, chip, hdu->header, 1.0)) {
    670                 psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for input FPA.");
    671                 psFree(view);
    672                 psFree(stats);
    673                 goto DONE;
    674             }
    675         }
    676         // Modify structure here.
    677 
    678        
    679         pmCell *cell;
    680         while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
    681             psTrace ("pswarp", 4, "DCell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
    682             if (!cell->process || !cell->file_exists) { continue; }
    683             if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
    684                 psError(psErrorCodeLast(), false, "Unable to read files.");
    685                 goto DONE;
    686             }
    687 
    688             psListAdd(cells, PS_LIST_TAIL, cell);
    689 
    690             // process each of the readouts
    691             pmReadout *readout;
    692             while ((readout = pmFPAviewNextReadout(view, input->fpa, 1)) != NULL) {
    693                 if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
    694                     psError(psErrorCodeLast(), false, "Unable to read files.");
    695                     goto DONE;
    696                 }
    697                 if (!readout->data_exists) {
    698                     continue;
    699                 }
    700 
    701                 for (int x = 0; x < readout->image->numCols; x++) {
    702                   for (int y = 0; y < readout->image->numRows; y++) {
    703                     readout->image->data.F32[y][x] = readout->image->data.F32[y][x] * (cd1f * cd2f) /
    704                       (psMetadataLookupS32(&mdok,config->arguments,"BKG.XGRID") *
    705                        psMetadataLookupS32(&mdok,config->arguments,"BKG.YGRID"));
    706                   }
    707                 }
    708                
    709                 psMetadataAddS32(config->arguments,PS_LIST_TAIL, "INTERPOLATION.MODE", PS_META_REPLACE, "", 8);
    710 
    711                 psMetadataAddBool(config->arguments,PS_LIST_TAIL, "BACKGROUND_WARPING", PS_META_REPLACE, "", true);
    712                 pswarpTransformReadout(output, readout, config);
    713                 psMetadataAddBool(config->arguments,PS_LIST_TAIL, "BACKGROUND_WARPING", PS_META_REPLACE, "", false);
    714 
    715                 if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
    716                     psError(psErrorCodeLast(), false, "Unable to write files.");
    717                     goto DONE;
    718                 }
    719             }
    720             if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
    721                 psError(psErrorCodeLast(), false, "Unable to write files.");
    722                 goto DONE;
    723             }
    724         }
    725         if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
    726             psError(psErrorCodeLast(), false, "Unable to write files.");
    727             goto DONE;
    728         }
    729     }
    730     if (!output->data_exists) {
    731         psWarning("No overlap between input and skycell.");
    732         psphotFilesActivate(config, false);
    733         psFree(cells);
    734         psFree(view);
    735         goto DONE;
    736     }
    737     pmCell *outCell = output->parent;   ///< Output cell
    738     pmChip *outChip = outCell->parent;  ///< Output chip
    739     pmFPA *outFPA = outChip->parent;    ///< Output FP
    740 
    741 /*     if (!pswarpPixelsLit(output, stats, config)) { */
    742 /*         psError(psErrorCodeLast(), false, "Unable to calculate pixel regions."); */
    743 /*         psFree(cells); */
    744 /*         psFree(view); */
    745 /*         goto DONE; */
    746 /*     } */
    747     psRegion *trimsec = psMetadataLookupPtr(NULL, outCell->concepts, "CELL.TRIMSEC"); ///< Trim section
    748     trimsec->x0 = trimsec->x1 = trimsec->y0 = trimsec->y1 = 0; ///< All pixels
    749 
    750     if (!psMetadataCopy(outFPA->concepts, input->fpa->concepts)) {
    751         psError(psErrorCodeLast(), false, "Unable to copy FPA concepts from input to output.");
    752         psFree(stats);
    753         psFree(view);
    754         goto DONE;
    755     }
    756 
    757     pmHDU *hdu = outFPA->hdu;           ///< HDU for the output warped image
    758 
    759     // Copy header from target
    760     {
    761         pmFPAview *skyView = pmFPAviewAlloc(0); ///< View into skycell
    762         skyView->chip = skyView->cell = 0;
    763         pmCell *cell = pmFPAfileThisCell(config->files, skyView, "PSWARP.SKYCELL"); // Skycell cell
    764         psFree(skyView);
    765         pmHDU *skyHDU = pmHDUFromCell(cell); ///< HDU
    766         if (!skyHDU) {
    767             psError(PSWARP_ERR_DATA, false, "Unable to find skycell HDU.");
    768             psFree(view);
    769             goto DONE;
    770         }
    771         hdu->header = psMetadataCopy(hdu->header, skyHDU->header);
    772     }
    773     pswarpVersionHeader(hdu->header);
    774    
    775     if (!pmAstromWriteWCS(hdu->header, outFPA, outChip, WCS_NONLIN_TOL)) {
    776         psError(psErrorCodeLast(), false, "Unable to generate WCS header.");
    777         psFree(stats);
    778         goto DONE;
    779     }
    780     if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
    781         psError(psErrorCodeLast(), false, "Unable to write files.");
    782         goto DONE;
    783     }
    784     // Done with the detector side of things
    785     pswarpFileActivation(config, detectorFiles, false);
    786     pswarpFileActivation(config, independentFiles, false);
    787 
    788 
    789     // Add MD5 information for readout
    790     const char *chipName = psMetadataLookupStr(NULL, output->parent->parent->concepts, "CHIP.NAME");
    791     const char *cellName = psMetadataLookupStr(NULL, output->parent->concepts, "CELL.NAME");
    792     psString headerName = NULL; ///< Header name for MD5
    793     psStringAppend(&headerName, "MD5_%s_%s_%d", chipName, cellName, view->readout);
    794     psVector *md5 = psImageMD5(output->image); ///< md5 hash
    795     psString md5string = psMD5toString(md5); ///< String
    796     psFree(md5);
    797     psMetadataAddStr(hdu->header, PS_LIST_TAIL, headerName, PS_META_REPLACE,
    798                      "Image MD5", md5string);
    799     psFree(md5string);
    800     psFree(headerName);
    801     psFree(view);
    802  DONE:
    803 
    804     return true;
    805 }
Note: See TracChangeset for help on using the changeset viewer.