IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 35373


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

Location:
branches/eam_branches/ipp-20130307/pswarp/src
Files:
4 added
5 edited
1 moved

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20130307/pswarp/src/Makefile.am

    r28043 r35373  
    2828        pswarp.c                        \
    2929        pswarpArguments.c               \
    30         pswarpExit.c                    \
     30        pswarpExitCode.c                \
    3131        pswarpLoop.c                    \
     32        pswarpLoopBackground.c          \
     33        pswarpCleanup.c                 \
     34        pswarpStatsFile.c               \
    3235        pswarpDefine.c                  \
     36        pswarpDefineBackground.c        \
    3337        pswarpDefineSkycell.c           \
    3438        pswarpErrorCodes.c              \

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

    r34800 r35373  
    1212
    1313#include "pswarp.h"
    14 #include "pswarpFileNames.h"
    1514
    1615int main (int argc, char **argv)
    … …  
    2221    psphotInit();
    2322
    24     const char *statsName = NULL;       // Filename for statistics
    25     psMetadata *stats = NULL;           // Container for statistics
    26     FILE *statsFile = NULL;             // File stream for statistics
    27 
    2823    pmConfig *config = pswarpArguments(argc, argv);
    2924    if (!config) {
    30         goto DIE;
     25        pswarpCleanup(config, NULL);
    3126    }
     27
     28    pswarpStatsFile *statsFile = pswarpStatsFileOpen (config);
    3229
    3330    pswarpVersionPrint();
    … …  
    3532    // load identify the data sources
    3633    if (!pswarpParseCamera(config)) {
    37         goto DIE;
     34        pswarpCleanup(config, statsFile);
    3835    }
    3936
    4037    if (!pswarpOptions(config)) {
    41         goto DIE;
     38        pswarpCleanup(config, statsFile);
    4239    }
    4340
    4441    // load the skycell layout information
    4542    if (!pswarpDefine(config)) {
    46         goto DIE;
     43        pswarpCleanup(config, statsFile);
    4744    }
    48     if (psMetadataLookupBool(NULL,config->arguments,"BACKGROUND.MODEL")) {
    49       if (!pswarpDefineBackground(config)) {
    50         goto DIE;
    51       }
    52     }
    53 
    54     // Open the statistics file
    55     bool mdok;                          ///< Status of MD lookup
    56     statsName = psMetadataLookupStr(&mdok, config->arguments, "STATS"); ///< Filename for statistics
    57     if (mdok && statsName && strlen(statsName) > 0) {
    58         psString resolved = pmConfigConvertFilename(statsName, config, true, true);
    59         statsFile = fopen(resolved, "w");
    60         if (!statsFile) {
    61             psError(PS_ERR_IO, true, "Unable to open statistics file %s for writing.\n", resolved);
    62             psFree(resolved);
    63             goto DIE;
    64         }
    65         psFree(resolved);
    66         stats = psMetadataAlloc();
    67         psMetadataAddS32(stats, PS_LIST_TAIL, "QUALITY", 0, "No problems", 0);
     45    if (!pswarpDefineBackground(config)) {
     46        pswarpCleanup(config, statsFile);
    6847    }
    6948
    7049    // load and warp
    71     if (!pswarpLoop(config, stats)) {
    72         goto DIE;
     50    if (!pswarpLoop(config, statsFile->md)) {
     51        pswarpCleanup(config, statsFile);
    7352    }
    74     if (psMetadataLookupBool(NULL,config->arguments,"BACKGROUND.MODEL")) {
    75       if (!pswarpLoopBackground(config, stats)) {
     53    if (!pswarpLoopBackground(config, statsFile->md)) {
    7654        fprintf(stderr,"Dying!\n");
    77         goto DIE;
    78       }
     55        pswarpCleanup(config, statsFile);
    7956    }
    8057
    8158    psLogMsg("pswarp", PS_LOG_INFO, "complete pswarp run: %f sec\n", psTimerMark("pswarp"));
     59    pswarpCleanup(config, statsFile);
     60}
    8261
    83 DIE:
    84     {
    85         psExit exitValue = pswarpExitCode(PS_EXIT_SUCCESS); // Exit code
    86 
    87         // Ensure everything is written out, at every level
    88         pswarpFileActivation(config, detectorFiles, true);
    89         pswarpFileActivation(config, skycellFiles, true);
    90         pswarpFileActivation(config, photFiles, true);
    91         pswarpFileActivation(config, independentFiles, true);
    92         if (!pswarpIOChecksAfter(config)) {
    93             psError(psErrorCodeLast(), false, "Unable to write files.");
    94             exitValue = pswarpExitCode(exitValue);
    95             pmFPAfileFreeSetStrict(false);
    96         }
    97 
    98         // Write out summary statistics
    99         if (stats && statsFile) {
    100             psMetadataAddF32(stats, PS_LIST_TAIL, "DT_WARP", 0, "Time for warp completion",
    101                              psTimerMark("pswarp"));
    102 
    103             const char *statsMDC = psMetadataConfigFormat(stats);
    104             if (!statsMDC) {
    105                 psError(psErrorCodeLast(), false, "Unable to get statistics file.");
    106                 exitValue = pswarpExitCode(exitValue);
    107             }
    108             psFree(stats);
    109             if (fprintf(statsFile, "%s", statsMDC) != strlen(statsMDC)) {
    110                 psError(PSWARP_ERR_IO, true, "Unable to write statistics file.");
    111                 exitValue = pswarpExitCode(exitValue);
    112             }
    113             psFree(statsMDC);
    114             if (fclose(statsFile) == EOF) {
    115                 psError(PSWARP_ERR_IO, true, "Unable to close statistics file.");
    116                 exitValue = pswarpExitCode(exitValue);
    117             }
    118             pmConfigRunFilenameAddWrite(config, "STATS", statsName);
    119         }
    120 
    121         // Dump configuration
    122         psString dump_file = psMetadataLookupStr(&mdok, config->arguments, "DUMP_CONFIG");
    123         if (dump_file) {
    124             if (!pmConfigDump(config, dump_file)) {
    125                 psError(psErrorCodeLast(), false, "Unable to dump configuration");
    126                 exitValue = pswarpExitCode(exitValue);
    127             }
    128         }
    129 
    130         psThreadPoolFinalize();
    131         psMemCheckCorruption(stderr, true);
    132 
    133         psFree(config);
    134 
    135         psTimerStop();
    136         pmVisualClose();
    137         pmModelClassCleanup();
    138         pmConceptsDone();
    139         pmConfigDone();
    140         psLibFinalize();
    141 
    142         exitValue = pswarpExitCode(exitValue);
    143         return exitValue;
    144     }
    145 }
     62// NOTE: pswarpCleanup calls exit

  • branches/eam_branches/ipp-20130307/pswarp/src/pswarp.h

    r34800 r35373  
    7777    double jacobian;                    ///< (Square root of) local Jacobian
    7878} pswarpTransformTileArgs;
     79
     80typedef struct {
     81    FILE *f;                            // File stream for statistics
     82    char *name;                         // Filename for statistics
     83    psMetadata *md;                     // Container for statistics
     84} pswarpStatsFile;
    7985
    8086pswarpTransformTileArgs *pswarpTransformTileArgsAlloc(void);
    … …  
    160166    );
    161167
     168pswarpStatsFile *pswarpStatsFileAlloc ();
     169pswarpStatsFile *pswarpStatsFileOpen (pmConfig *config);
     170bool pswarpStatsFileSave (pmConfig *config, pswarpStatsFile *statsFile);
     171
     172// cleanup memory and exit
     173void pswarpCleanup (pmConfig *config, pswarpStatsFile *statsFile) PS_ATTR_NORETURN;

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

    r35372 r35373  
    11/** @file pswarpDefine.c
    22 *
    3  *  @brief
    4  *
     3 *  @brief generates the output images
    54 *  @ingroup pswarp
    65 *
    … …  
    1312# include "pswarp.h"
    1413
    15 /**
    16  * loads the skycell metadata
    17  */
    1814bool pswarpDefine (pmConfig *config) {
    1915
    … …  
    132128}
    133129
    134 bool pswarpDefineBackground (pmConfig *config) {
    135 
    136     // load the PSWARP recipe
    137     psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSWARP_RECIPE);
    138     if (!recipe) {
    139         psError(PSWARP_ERR_CONFIG, true, "Can't find PSWARP recipe!\n");
    140         return false;
    141     }
    142 
    143     // select the input data sources
    144     pmFPAfile *skycell = psMetadataLookupPtr (NULL, config->files, "PSWARP.SKYCELL");
    145     if (!skycell) {
    146         psError(PSWARP_ERR_CONFIG, false, "Can't find skycell data!\n");
    147         return false;
    148     }
    149     pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PSWARP.OUTPUT.BKGMODEL");
    150     if (!output) {
    151         psError(PSWARP_ERR_CONFIG, false, "Can't find output data!\n");
    152         return false;
    153     }
    154     pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PSWARP.INPUT");
    155     if (!input) {
    156         psError(PSWARP_ERR_CONFIG, false, "Can't find input data!\n");
    157         return false;
    158     }
    159    
    160     // open the full skycell file; no need to defer different depths. only load the header data
    161     pmFPAview *view = pmFPAviewAlloc (0);
    162     pmFPAfileOpen (skycell, view, config);
    163 
    164     // Read header and create target
    165     {
    166         if (!pmFPAReadHeaderSet(skycell->fpa, skycell->fits, config)) {
    167             psError(psErrorCodeLast(), false, "Unable to read headers for skycell.");
    168             psFree(view);
    169             return false;
    170         }
    171         view->chip = 0;
    172         view->cell = 0;
    173         view->readout = 0;
    174         pmCell *source = pmFPAfileThisCell(config->files, view, "PSWARP.SKYCELL"); ///< Source cell
    175         pmHDU *hdu = pmHDUFromCell(source); ///< HDU for source
    176         if (!hdu || !hdu->header) {
    177             psError(PM_ERR_PROG, false, "Unable to find header for sky cell.");
    178             psFree(view);
    179             return false;
    180         }
    181         int numCols = psMetadataLookupS32(NULL, hdu->header, "NAXIS1"); ///< Number of columns
    182         int numRows = psMetadataLookupS32(NULL, hdu->header, "NAXIS2"); ///< Number of rows
    183         if ((numCols == 0) || (numRows == 0)) {
    184             psError(PSWARP_ERR_DATA, false, "skycell has invalid dimensions %d x %d", numCols, numRows);
    185             psFree(view);
    186             return false;
    187         }
    188 
    189         pmCell *target = pmFPAviewThisCell(view, output->fpa); ///< Target cell
    190         pmReadout *readout = pmReadoutAlloc(target); ///< Target readout
    191         readout->image = psImageAlloc(numCols / output->xBin, numRows / output->yBin, PS_TYPE_F32);
    192         psImageInit(readout->image, NAN);
    193         psFree(readout);                // Drop reference
    194 
    195         bool status = false;
    196         psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.XBIN");
    197         psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.YBIN");
    198 /*      psMetadataAddS32(target->concepts,PS_LIST_TAIL,"CELL.XBIN",PS_META_REPLACE,"",output->xBin); */
    199 /*      psMetadataAddS32(target->concepts,PS_LIST_TAIL,"CELL.YBIN",PS_META_REPLACE,"",output->yBin); */
    200         psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.XSIZE");
    201         psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.YSIZE");
    202 /*      psMetadataAddS32(target->concepts,PS_LIST_TAIL,"CELL.XSIZE",PS_META_REPLACE,"",numCols / output->xBin); */
    203 /*      psMetadataAddS32(target->concepts,PS_LIST_TAIL,"CELL.YSIZE",PS_META_REPLACE,"",numRows / output->yBin); */
    204         psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.XPARITY");
    205         psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.YPARITY");
    206         psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.X0");
    207         psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.Y0");
    208 
    209        
    210     }
    211 
    212     // XXX this is not a sufficient test
    213     view->chip = 0;
    214     view->cell = 0;
    215     view->readout = -1;
    216     pmHDU *phu = pmFPAviewThisPHU(view, skycell->fpa); ///< Skycell PHU
    217     pmHDU *hdu = pmFPAviewThisHDU(view, skycell->fpa); ///< Skycell header
    218 
    219     pmAstromWCS *WCS = pmAstromWCSfromHeader(hdu->header);
    220 
    221     double cd1f = 1.0 * output->xBin;
    222     double cd2f = 1.0 * output->yBin;
    223    
    224     WCS->crpix1 = WCS->crpix1 / cd1f;
    225     WCS->crpix2 = WCS->crpix2 / cd2f;
    226    
    227     WCS->cdelt1 *= cd1f;
    228     WCS->cdelt2 *= cd2f;
    229 
    230     WCS->trans->x->coeff[1][0] *= cd1f;
    231     WCS->trans->x->coeff[0][1] *= cd2f;
    232     WCS->trans->y->coeff[1][0] *= cd1f;
    233     WCS->trans->y->coeff[0][1] *= cd2f;
    234    
    235    
    236     pmAstromWCStoHeader (hdu->header,WCS);
    237 
    238     bool bilevelAstrometry = false;
    239     if (phu) {
    240         char *ctype = psMetadataLookupStr(NULL, phu->header, "CTYPE1");
    241         if (ctype) {
    242             bilevelAstrometry = !strcmp(&ctype[4], "-DIS");
    243         }
    244     }
    245 
    246     // We read from the skycell into the output.  i.e., the output receives the desired astrometry.
    247     pmChip *outputChip = pmFPAviewThisChip(view, output->fpa); ///< Chip in the output
    248     if (bilevelAstrometry) {
    249         if (!pmAstromReadBilevelMosaic(output->fpa, phu->header)) {
    250             psError(psErrorCodeLast(), false, "Unable to read bilevel mosaic astrometry for skycell.");
    251             psFree(view);
    252             return false;
    253         }
    254         if (!pmAstromReadBilevelChip(outputChip, hdu->header)) {
    255             psError(psErrorCodeLast(), false, "Unable to read bilevel chip astrometry for skycell.");
    256             psFree(view);
    257             return false;
    258         }
    259     } else {
    260         // we use a default FPA pixel scale of 1.0
    261       if (!pmAstromReadWCS(output->fpa, outputChip, hdu->header, 1.0)) {
    262             psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for skycell.");
    263             psFree(view);
    264             return false;
    265         }
    266     }
    267 
    268    
    269     view->chip = view->cell = view->readout = -1;
    270     pmFPAAddSourceFromView(output->fpa, view, output->format);
    271 
    272 
    273     psFree (view);
    274     return true;
    275 }

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

    r35372 r35373  
    1212
    1313# include "pswarp.h"
    14 
    1514
    1615psExit pswarpExitCode(psExit exitValue)

  • 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.