IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
May 9, 2013, 12:26:48 PM (13 years ago)
Author:
eugene
Message:

major upgrades to pswarp to enable skycell warp -> chip

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/pswarp/src/pswarpLoop.c

    r34800 r35563  
    11/** @file pswarpLoop.c
    22 *
    3  *  @brief
    4  *
     3 *  @brief main processing loop for pswarp
    54 *  @ingroup pswarp
    65 *
    … …  
    1211
    1312#include "pswarp.h"
    14 #include <ppStats.h>
    15 #include "pswarpFileNames.h"            // Lists of file rules used at different stages
    1613
    17 #define WCS_NONLIN_TOL 0.001            // Non-linear tolerance for header WCS
    18 #define TESTING 0                       // Testing output?
    19 
    20 
    21 
    22 // Loop over the inputs, warp them to the output skycell and then write out the output.
     14// Loop over the inputs, warp them to the output skycell and then update metadata
    2315bool pswarpLoop(pmConfig *config, psMetadata *stats)
    2416{
    25     bool status;
    26     bool mdok;                          // Status of MD lookup
    27 
    28     const char *skyCamera = psMetadataLookupStr(NULL, config->arguments,
    29                                                 "SKYCELL.CAMERA");  // Name of camera for skycell
    30     pmConfigCamerasCull(config, skyCamera);
    31     pmConfigRecipesCull(config, "PSWARP,PPSTATS,PSPHOT,PSASTRO,MASKS,JPEG");
    3217
    3318    // load the recipe
     19    bool status = false;
    3420    psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSWARP_RECIPE);
    3521    if (!recipe) {
    … …  
    3824    }
    3925
    40     if (!pswarpSetMaskBits(config)) {
    41         psError(psErrorCodeLast(), false, "failed to set mask bits");
    42         return NULL;
    43     }
    44 
    45     // output mask bits
    46     psImageMaskType maskValue = psMetadataLookupImageMask(&status, recipe, "MASK.OUTPUT");
    47     psAssert (status, "MASK.OUTPUT was not defined");
    48 
    4926    // select the input data sources
    50     pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PSWARP.INPUT");
    51     if (!input) {
    52         psError(PSWARP_ERR_CONFIG, true, "Can't find input data!\n");
     27    pmFPAfile *output = psMetadataLookupPtr(&status, config->files, "PSWARP.OUTPUT");
     28    if (!output) {
     29        psError(PSWARP_ERR_CONFIG, true, "Can't find output data!\n");
    5330        return false;
    5431    }
    5532
    5633    // use the external astrometry source if supplied
    57     pmFPAfile *astrom = psMetadataLookupPtr(NULL, config->files, "PSWARP.ASTROM");
    58     if (!astrom) {
    59         astrom = input;
    60     }
    61 
    62     if (astrom->camera != input->camera) {
    63         psError(PSWARP_ERR_DATA, true, "Input camera and astrometry camera do not match.");
     34    pmFPAfile *skycell = psMetadataLookupPtr(&status, config->files, "PSWARP.SKYCELL");
     35    if (!skycell) {
     36        psError(PSWARP_ERR_DATA, true, "Cannot find output astrometry.");
    6437        return false;
    6538    }
    6639
    67     // select the output readout
    6840    pmFPAview *view = pmFPAviewAlloc(0);
    69     view->chip = 0;
    70     view->cell = 0;
    71     view->readout = 0;
    72     pmReadout *output = pmFPAfileThisReadout(config->files, view, "PSWARP.OUTPUT");
    73     if (!output) {
    74         psError(PSWARP_ERR_CONFIG, true, "Can't find output data!\n");
     41
     42    int nInputs = psMetadataLookupS32(&status, config->arguments, "NUM_INPUTS");
     43    if (!status) {
     44        psError(PSWARP_ERR_DATA, true, "number of inputs is not defined (programming error)");
    7545        return false;
    7646    }
    77     psFree (view);
    7847
     48    // load in the input pixel data (ex. background model)
     49    pmFPAfileActivate(config->files, false, NULL);
     50    pmFPAfileActivate(config->files, true, "PSWARP.INPUT");
     51    pmFPAfileActivate(config->files, true, "PSWARP.MASK");
     52    pmFPAfileActivate(config->files, true, "PSWARP.VARIANCE");
    7953
    80     // Turn all skycell files on to generate them, and then turn them off for the loop over the input images
    81     // the input, which is in a different format.
    82     {
    83         pswarpFileActivation(config, detectorFiles, false);
    84         pswarpFileActivation(config, photFiles, false);
    85         pswarpFileActivation(config, independentFiles, false);
    86         pswarpFileActivation(config, skycellFiles, true);
    87         if (!pswarpIOChecksBefore(config)) {
    88             psError(psErrorCodeLast(), false, "Unable to read files.");
    89             goto DONE;
    90         }
    91         pswarpFileActivation(config, skycellFiles, false);
     54    // We re-activate the CMF load so we can transform the sources as well as the pixels.
     55    // We only need to read in these if the astrometry source is CMF.
     56    pmFPAfileActivate(config->files, true, "PSWARP.ASTROM");
     57
     58    // loop over this section once per input group
     59    for (int i = 0; i < nInputs; i++) {
     60        // select the input data sources
     61        pmFPAfile *input = pmFPAfileSelectSingle(config->files, "PSWARP.INPUT", i);
     62        if (!input) {
     63            psError(PSWARP_ERR_CONFIG, true, "Can't find input data!\n");
     64            return false;
     65        }
     66
     67        // select the input data sources
     68        pmFPAfile *astrom = pmFPAfileSelectSingle(config->files, "PSWARP.ASTROM", i);
     69        if (!astrom) {
     70            astrom = input;
     71        }
     72
     73        pmFPAviewReset (view);
     74
     75        // files associated with the science image
     76        if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
     77            psError(psErrorCodeLast(), false, "Unable to read files.");
     78            goto FAIL;
     79        }
     80
     81        // *** main transformation block
     82        // *** this section loops over the input chips/cells and reads them one at a time
     83        // *** the output chips/cells are filled where appropriate, but not yet written to disk
     84        pmChip *chip;
     85        while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
     86            psTrace ("pswarp", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
     87            if (!chip->process || !chip->file_exists) { continue; }
     88            if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
     89                psError(psErrorCodeLast(), false, "Unable to read files.");
     90                goto FAIL;
     91            }
     92
     93            pmCell *cell;
     94            while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
     95                psTrace ("pswarp", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
     96                if (!cell->process || !cell->file_exists) { continue; }
     97                if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
     98                    psError(psErrorCodeLast(), false, "Unable to read files.");
     99                    goto FAIL;
     100                }
     101
     102                // process each of the readouts
     103                pmReadout *readout;
     104                while ((readout = pmFPAviewNextReadout(view, input->fpa, 1)) != NULL) {
     105                    if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
     106                        psError(psErrorCodeLast(), false, "Unable to read files.");
     107                        goto FAIL;
     108                    }
     109                    if (!readout->data_exists) {
     110                        continue;
     111                    }
     112
     113                    // Copy the detections from the astrometry carrier to the input, so they can be accessed by
     114                    // pswarpTransformReadout
     115                    if (astrom != input) {
     116                        pmReadout *astromRO = pmFPAviewThisReadout(view, astrom->fpa); // Readout for astrometry
     117                        pmDetections *detections = psMetadataLookupPtr(&status, astromRO->analysis, "PSPHOT.DETECTIONS"); // Sources from astrometry
     118                        if (detections) {
     119                            psMetadataAddPtr(readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_ARRAY, "Sources from input astrometry", detections);
     120                        }
     121                    }
     122
     123                    pswarpTransformToTarget(output->fpa, readout, config, false);
     124               
     125                    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
     126                        psError(psErrorCodeLast(), false, "Unable to write files.");
     127                        goto FAIL;
     128                    }
     129                }
     130                if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
     131                    psError(psErrorCodeLast(), false, "Unable to write files.");
     132                    goto FAIL;
     133                }
     134            }
     135            if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
     136                psError(psErrorCodeLast(), false, "Unable to write files.");
     137                goto FAIL;
     138            }
     139        }
     140        if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
     141            psError(psErrorCodeLast(), false, "Unable to write files.");
     142            goto FAIL;
     143        }
     144
     145        if (!pswarpUpdateStatistics (output->fpa, stats, input->fpa, astrom->fpa, config)) {
     146            psError(psErrorCodeLast(), false, "problem generating statistics.");
     147            goto FAIL;
     148        }
     149        if (!pswarpUpdateMetadata (output->fpa, skycell->fpa, input->fpa, astrom->fpa, config, true)) {
     150            psError(psErrorCodeLast(), false, "problem generating statistics.");
     151            goto FAIL;
     152        }
    92153    }
    93154
    94     // Read the input astrometry
    95     // XXX rather than use the activations here, this should just explicitly loop over the desired filerule
    96     {
    97         pmFPAfileActivate(config->files, true, "PSWARP.ASTROM");
    98 
    99         pmChip *chip;
    100         pmFPAview *view = pmFPAviewAlloc(0);
    101         if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
    102             psError(psErrorCodeLast(), false, "Unable to read files.");
    103             goto DONE;
    104         }
    105         while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
    106             psTrace ("pswarp", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
    107             if (!chip->process || !chip->file_exists) { continue; }
    108             if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
    109                 psError(psErrorCodeLast(), false, "Unable to read files.");
    110                 goto DONE;
    111             }
    112             pmCell *cell;
    113             while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
    114                 psTrace ("pswarp", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
    115                 if (!cell->process || !cell->file_exists) { continue; }
    116                 if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE) ||
    117                     !pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) {
    118                     psError(psErrorCodeLast(), false, "Unable to read files.");
    119                     goto DONE;
    120                 }
    121             }
    122             if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) {
    123                 psError(psErrorCodeLast(), false, "Unable to write files.");
    124                 goto DONE;
    125             }
    126         }
    127         if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) {
    128             psError(psErrorCodeLast(), false, "Unable to write files.");
    129             goto DONE;
    130         }
    131         psFree(view);
    132 
    133         pswarpFileActivation(config, detectorFiles, true);
    134         pmFPAfileActivate(config->files, false, "PSWARP.ASTROM");
     155    if (!pswarpMakePSF (config, output, stats)) {
     156      psError(psErrorCodeLast(), false, "problem generating PSF.");
     157      goto FAIL;
    135158    }
    136159
    137     // Turn on the source output --- we need to get rid of these so that we can measure the PSF
    138     pmFPAfileActivate(config->files, true, "PSWARP.OUTPUT.SOURCES");
     160    psFree(view);
     161    return true;
    139162
    140     // Don't care about the skycell anymore --- we've read it, and that's all we need to do.
    141     pmFPAfileActivate(config->files, false, "PSWARP.SKYCELL");
    142     view = pmFPAviewAlloc(0);
     163 FAIL:
     164    psFree (view);
     165    return false;
     166}
    143167
    144     // find the FPA phu
    145     bool bilevelAstrometry = false;
    146     pmHDU *phu = pmFPAviewThisPHU(view, astrom->fpa);
    147     if (phu) {
    148         char *ctype = psMetadataLookupStr(NULL, phu->header, "CTYPE1");
    149         if (ctype) {
    150             bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
    151         }
    152     }
    153     if (bilevelAstrometry) {
    154         if (!pmAstromReadBilevelMosaic(input->fpa, phu->header)) {
    155             psError(psErrorCodeLast(), false, "Unable to read bilevel mosaic astrometry for input FPA.");
    156             psFree(view);
    157             psFree(stats);
    158             goto DONE;
    159         }
    160     }
     168// once the output fpa elements have been built, loop over the fpa and generate stats
     169// for each readout
     170bool pswarpTransformToTarget (pmFPA *output, pmReadout *input, pmConfig *config, bool backgroundWarp)  {
    161171
    162     psList *cells = psListAlloc(NULL);  // List of cells, for concepts averaging
    163 
    164     // files associated with the science image
    165     if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
    166         psError(psErrorCodeLast(), false, "Unable to read files.");
    167         goto DONE;
    168     }
    169 
     172    pmFPAview *view = pmFPAviewAlloc(0);
     173   
    170174    pmChip *chip;
    171     while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
     175    while ((chip = pmFPAviewNextChip (view, output, 1)) != NULL) {
    172176        psTrace ("pswarp", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
    173177        if (!chip->process || !chip->file_exists) { continue; }
    174         if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
    175             psError(psErrorCodeLast(), false, "Unable to read files.");
    176             goto DONE;
    177         }
    178178
    179         // read WCS data from the corresponding header
    180         pmHDU *hdu = pmFPAviewThisHDU (view, astrom->fpa);
    181 
    182        
    183         if (bilevelAstrometry) {
    184             if (!pmAstromReadBilevelChip (chip, hdu->header)) {
    185                 psError(psErrorCodeLast(), false, "Unable to read bilevel chip astrometry for input FPA.");
    186                 psFree(view);
    187                 psFree(stats);
    188                 goto DONE;
    189             }
    190         } else {
    191             // we use a default FPA pixel scale of 1.0
    192             if (!pmAstromReadWCS (input->fpa, chip, hdu->header, 1.0)) {
    193                 psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for input FPA.");
    194                 psFree(view);
    195                 psFree(stats);
    196                 goto DONE;
    197             }
    198         }
    199        
    200179        pmCell *cell;
    201         while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
     180        while ((cell = pmFPAviewNextCell (view, output, 1)) != NULL) {
    202181            psTrace ("pswarp", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
    203182            if (!cell->process || !cell->file_exists) { continue; }
    204             if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
    205                 psError(psErrorCodeLast(), false, "Unable to read files.");
    206                 goto DONE;
    207             }
    208 
    209             psListAdd(cells, PS_LIST_TAIL, cell);
    210183
    211184            // process each of the readouts
    212185            pmReadout *readout;
    213             while ((readout = pmFPAviewNextReadout(view, input->fpa, 1)) != NULL) {
    214                 if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
    215                     psError(psErrorCodeLast(), false, "Unable to read files.");
    216                     goto DONE;
    217                 }
    218                 if (!readout->data_exists) {
    219                     continue;
    220                 }
    221 
    222                 // Copy the detections from the astrometry carrier to the input, so they can be accessed by
    223                 // pswarpTransformReadout
    224                 pmReadout *astromRO = pmFPAviewThisReadout(view, astrom->fpa); // Readout for astrometry
    225                 pmDetections *detections = psMetadataLookupPtr(&mdok, astromRO->analysis, "PSPHOT.DETECTIONS"); // Sources from astrometry
    226                 if (detections) {
    227                     psMetadataAddPtr(readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_ARRAY, "Sources from input astrometry", detections);
    228                 }
    229 
    230                 pswarpTransformReadout(output, readout, config);
    231                
    232                 if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
    233                     psError(psErrorCodeLast(), false, "Unable to write files.");
    234                     goto DONE;
    235                 }
    236             }
    237             if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
    238                 psError(psErrorCodeLast(), false, "Unable to write files.");
    239                 goto DONE;
    240             }
    241         }
    242         if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
    243             psError(psErrorCodeLast(), false, "Unable to write files.");
    244             goto DONE;
    245         }
     186            while ((readout = pmFPAviewNextReadout(view, output, 1)) != NULL) {
     187                pswarpTransformReadout (readout, input, config, backgroundWarp);
     188            }
     189        }
    246190    }
    247 
    248     if (!output->data_exists) {
    249         psWarning("No overlap between input and skycell.");
    250         if (stats) {
    251             psMetadataAddS32(stats, PS_LIST_TAIL, "QUALITY", PS_META_REPLACE,
    252                              "No overlap between input and skycell", PSWARP_ERR_NO_OVERLAP);
    253         }
    254         psphotFilesActivate(config, false);
    255         psFree(cells);
    256         psFree(view);
    257         goto DONE;
    258     }
    259    
    260     pmCell *outCell = output->parent;   ///< Output cell
    261     pmChip *outChip = outCell->parent;  ///< Output chip
    262     pmFPA *outFPA = outChip->parent;    ///< Output FP
    263 
    264     if (!pswarpPixelsLit(output, stats, config)) {
    265         psError(psErrorCodeLast(), false, "Unable to calculate pixel regions.");
    266         psFree(cells);
    267         psFree(view);
    268         goto DONE;
    269     }
    270     bool doStats = psMetadataLookupBool(&mdok,recipe,"MASK.STATS");
    271     if (doStats) {
    272       if (!pswarpMaskStats(output, stats, config)) {
    273         psError(psErrorCodeLast(), false, "Unable to calculate mask stats.");
    274         psFree(cells);
    275         psFree(view);
    276         goto DONE;
    277       }
    278     }
    279     // Set covariance matrix for output
    280     {
    281         psList *covariances = psMetadataLookupPtr(&mdok, output->analysis,
    282                                                   PSWARP_ANALYSIS_COVARIANCES); // Covariance matrices
    283         psAssert(covariances, "Should be there");
    284         psArray *covars = psListToArray(covariances); // Array of covariance matrices
    285         psKernel *covar = psImageCovarianceAverage(covars);
    286         psFree(covars);
    287         psMetadataRemoveKey(output->analysis, PSWARP_ANALYSIS_COVARIANCES);
    288 
    289         // Correct covariance matrix scale for the mean (square root of the) Jacobian
    290         double jacobian = psMetadataLookupF64(NULL, output->analysis, PSWARP_ANALYSIS_JACOBIAN); // Jacobian
    291         int goodPixels = psMetadataLookupS32(NULL, output->analysis, PSWARP_ANALYSIS_GOODPIX);   // Good pixels
    292         jacobian /= goodPixels;
    293         output->covariance = psImageCovarianceScale(covar, jacobian);
    294         psFree(covar);
    295 
    296         if (output->variance) {
    297             psImageCovarianceTransfer(output->variance, output->covariance);
    298         }
    299     }
    300 
    301     if (!pmConceptsAverageCells(outCell, cells, NULL, NULL, false)) {
    302         psError(psErrorCodeLast(), false, "Unable to average cell concepts.");
    303         psFree(stats);
    304         psFree(cells);
    305         psFree(view);
    306         goto DONE;
    307     }
    308     psFree(cells);
    309 
    310     psRegion *trimsec = psMetadataLookupPtr(NULL, outCell->concepts, "CELL.TRIMSEC"); ///< Trim section
    311     trimsec->x0 = trimsec->x1 = trimsec->y0 = trimsec->y1 = 0; ///< All pixels
    312 
    313     if (!psMetadataCopy(outFPA->concepts, input->fpa->concepts)) {
    314         psError(psErrorCodeLast(), false, "Unable to copy FPA concepts from input to output.");
    315         psFree(stats);
    316         psFree(view);
    317         goto DONE;
    318     }
    319 
    320     // Update ZP from the astrometry
    321     {
    322         psMetadataItem *item = psMetadataLookup(outFPA->concepts, "FPA.ZP");
    323         item->data.F32 = psMetadataLookupF32(NULL, astrom->fpa->concepts, "FPA.ZP");
    324     }
    325 
    326     pmHDU *hdu = outFPA->hdu;           ///< HDU for the output warped image
    327 
    328     // Copy header from target
    329     {
    330         pmFPAview *skyView = pmFPAviewAlloc(0); ///< View into skycell
    331         skyView->chip = skyView->cell = 0;
    332         pmCell *cell = pmFPAfileThisCell(config->files, skyView, "PSWARP.SKYCELL"); // Skycell cell
    333         psFree(skyView);
    334         pmHDU *skyHDU = pmHDUFromCell(cell); ///< HDU
    335         if (!skyHDU) {
    336             psError(PSWARP_ERR_DATA, false, "Unable to find skycell HDU.");
    337             psFree(view);
    338             goto DONE;
    339         }
    340         hdu->header = psMetadataCopy(hdu->header, skyHDU->header);
    341     }
    342 
    343     pswarpVersionHeader(hdu->header);
    344    
    345     if (!pmAstromWriteWCS(hdu->header, outFPA, outChip, WCS_NONLIN_TOL)) {
    346         psError(psErrorCodeLast(), false, "Unable to generate WCS header.");
    347         psFree(stats);
    348         goto DONE;
    349     }
    350 
    351     if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
    352         psError(psErrorCodeLast(), false, "Unable to write files.");
    353         goto DONE;
    354     }
    355 
    356     // Done with the detector side of things
    357     pswarpFileActivation(config, detectorFiles, false);
    358     pswarpFileActivation(config, independentFiles, false);
    359 
    360 
    361     // We need a new PSF model for the warped frame.  It would be good to generate this analytically, but
    362     // that's going to be tricky.  We have a list of sources, so we use those to redetermine the PSF model.
    363 
    364     if (psMetadataLookupBool(&mdok, recipe, "PSF")) {
    365         pswarpFileActivation(config, photFiles, true);
    366         if (!pswarpIOChecksBefore(config)) {
    367             psError(psErrorCodeLast(), false, "Unable to read files.");
    368             goto DONE;
    369         }
    370 
    371         // supply the readout and fpa of interest to psphot
    372         pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT");
    373         pmFPACopy(photFile->fpa, outFPA);
    374 
    375         pmFPAview *view = pmFPAviewAlloc(0); ///< View into skycell
    376         view->chip = view->cell = view->readout = 0;
    377 
    378         // grab the sources of interest from the storage location (pmFPAfile PSPHOT.INPUT.CMF)
    379         psArray *sources = psphotLoadPSFSources (config, view);
    380         if (!sources) {
    381             psError(psErrorCodeLast(), false, "No sources supplied to measure PSF");
    382             goto DONE;
    383         }
    384 
    385         pmModelClassSetLimits(PM_MODEL_LIMITS_STRICT);
    386 
    387         // measure the PSF using these sources
    388         if (!psphotReadoutFindPSF(config, view, "PSPHOT.INPUT", sources)) {
    389             // This is likely a data quality issue
    390             // XXX Split into multiple cases using error codes?
    391             psErrorStackPrint(stderr, "Unable to determine PSF");
    392             psWarning("Unable to determine PSF --- suspect bad data quality.");
    393             if (stats && psMetadataLookupS32(NULL, stats, "QUALITY") == 0) {
    394                 psMetadataAddS32(stats, PS_LIST_TAIL, "QUALITY", PS_META_REPLACE,
    395                                  "Unable to determine PSF", psErrorCodeLast());
    396             }
    397             psErrorClear();
    398             psphotFilesActivate(config, false);
    399         }
    400 
    401         // Ensure seeing is carried over
    402         pmChip *photChip = pmFPAviewThisChip(view, photFile->fpa);                 // Chip with seeing
    403         psMetadataItem *item = psMetadataLookup(outChip->concepts, "CHIP.SEEING"); // Concept with seeing
    404         item->data.F32 = psMetadataLookupF32(NULL, photChip->concepts, "CHIP.SEEING");
    405 
    406 // XXX EAM : put this in a visualization function
    407 #if (TESTING)
    408         {
    409             #define PSF_SIZE 20         ///< Half-size of PSF
    410             #define PSF_FLUX 10000      ///< Central flux for PSF
    411             pmChip *photChip = pmFPAviewThisChip(view, photFile->fpa);
    412             pmPSF *psf = psMetadataLookupPtr(NULL, photChip->analysis, "PSPHOT.PSF");
    413             psImage *image = psImageAlloc(2 * PSF_SIZE + 1, 2 * PSF_SIZE + 1, PS_TYPE_F32);
    414             psImageInit(image, 0);
    415             pmModel *model = pmModelFromPSFforXY(psf, PSF_SIZE, PSF_SIZE, PSF_FLUX);
    416             pmModelAdd(image, NULL, model, PM_MODEL_OP_FULL, 0);
    417             psFree(model);
    418             psFits *fits = psFitsOpen("psf.fits", "w");
    419             psFitsWriteImage(fits, NULL, image, 0, NULL);
    420             psFitsClose(fits);
    421             psFree(image);
    422         }
    423 #endif
    424 
    425         psFree(view);
    426     }
    427 
    428     // Perform statistics on the output image
    429     if (stats) {
    430         if (!ppStatsFPA(stats, output->parent->parent->parent, view, maskValue, config)) {
    431             psWarning("Unable to perform statistics on warped image.");
    432         }
    433     }
    434    
    435 
    436     // Add MD5 information for readout
    437     const char *chipName = psMetadataLookupStr(NULL, output->parent->parent->concepts, "CHIP.NAME");
    438     const char *cellName = psMetadataLookupStr(NULL, output->parent->concepts, "CELL.NAME");
    439     psString headerName = NULL; ///< Header name for MD5
    440     psStringAppend(&headerName, "MD5_%s_%s_%d", chipName, cellName, view->readout);
    441     psVector *md5 = psImageMD5(output->image); ///< md5 hash
    442     psString md5string = psMD5toString(md5); ///< String
    443     psFree(md5);
    444     psMetadataAddStr(hdu->header, PS_LIST_TAIL, headerName, PS_META_REPLACE,
    445                      "Image MD5", md5string);
    446     psFree(md5string);
    447     psFree(headerName);
    448     psFree(view);
    449 
    450  DONE:
    451 
    452191    return true;
    453192}
    454193
    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.