IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Apr 26, 2013, 9:07:23 AM (13 years ago)
Author:
eugene
Message:

major mods to allow multichip output

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r35424 r35434  
    2020bool pswarpLoop(pmConfig *config, psMetadata *stats)
    2121{
    22     bool status;
    23     bool mdok;                          // Status of MD lookup
    24 
    2522    const char *skyCamera = psMetadataLookupStr(NULL, config->arguments,
    2623                                                "SKYCELL.CAMERA");  // Name of camera for skycell
    … …  
    2926
    3027    // load the recipe
     28    bool status = false;
    3129    psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSWARP_RECIPE);
    3230    if (!recipe) {
    … …  
    3533    }
    3634
     35    // XXX this should probably happen elsewhere (logically it belongs to file I/O, right?)
    3736    if (!pswarpSetMaskBits(config)) {
    3837        psError(psErrorCodeLast(), false, "failed to set mask bits");
    3938        return NULL;
    4039    }
    41 
    42     // output mask bits
    43     psImageMaskType maskValue = psMetadataLookupImageMask(&status, recipe, "MASK.OUTPUT");
    44     psAssert (status, "MASK.OUTPUT was not defined");
    4540
    4641    // select the input data sources
    … …  
    6257    }
    6358
    64     // select the output readout
    65     pmFPAview *view = pmFPAviewAlloc(0);
    66     view->chip = 0;
    67     view->cell = 0;
    68     view->readout = 0;
    69     pmReadout *output = pmFPAfileThisReadout(config->files, view, "PSWARP.OUTPUT");
     59    // select the input data sources
     60    pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PSWARP.OUTPUT");
    7061    if (!output) {
    7162        psError(PSWARP_ERR_CONFIG, true, "Can't find output data!\n");
    7263        return false;
    7364    }
    74     psFree (view);
    75 
    76 
    77     // Turn all skycell files on to generate them, and then turn them off for the loop over the input images
    78     // the input, which is in a different format.
    79     // XXX why is this being done here, doesn't it duplicate the work in pswarpDefine.c??
     65
     66    // use the external astrometry source if supplied
     67    pmFPAfile *skycell = psMetadataLookupPtr(NULL, config->files, "PSWARP.SKYCELL");
     68    if (!skycell) {
     69        psError(PSWARP_ERR_DATA, true, "Cannot find output astrometry.");
     70        return false;
     71    }
     72
     73    // XXX NOTE : 'output' was a readout, now we need to loop over all output targets
     74
     75    // ensure everyone is off except what we need below
    8076    pswarpFileActivation(config, detectorFiles, false);
    8177    pswarpFileActivation(config, photFiles, false);
    … …  
    8379    pswarpFileActivation(config, skycellFiles, false);
    8480   
    85     // pswarpFileActivation(config, skycellFiles, true);
    86     // if (0) {
    87     //   if (!pswarpIOChecksBefore(config)) {
    88     //     psError(psErrorCodeLast(), false, "Unable to read files.");
    89     //     goto DONE;
    90     //   }
    91     // }
    92    
    93     // Read the input astrometry
    94     // XXX rather than use the activations here, this should just explicitly loop over the desired filerule
    95     {
    96         pmFPAfileActivate(config->files, true, "PSWARP.ASTROM");
    97 
    98         pmChip *chip;
    99         pmFPAview *view = pmFPAviewAlloc(0);
    100         if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
    101             psError(psErrorCodeLast(), false, "Unable to read files.");
    102             goto DONE;
    103         }
    104         while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
    105             psTrace ("pswarp", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
    106             if (!chip->process || !chip->file_exists) { continue; }
    107             if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
    108                 psError(psErrorCodeLast(), false, "Unable to read files.");
    109                 goto DONE;
    110             }
    111             pmCell *cell;
    112             while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
    113                 psTrace ("pswarp", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
    114                 if (!cell->process || !cell->file_exists) { continue; }
    115                 if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE) ||
    116                     !pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) {
    117                     psError(psErrorCodeLast(), false, "Unable to read files.");
    118                     goto DONE;
    119                 }
    120             }
    121             if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) {
    122                 psError(psErrorCodeLast(), false, "Unable to write files.");
    123                 goto DONE;
    124             }
    125         }
    126         if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) {
    127             psError(psErrorCodeLast(), false, "Unable to write files.");
    128             goto DONE;
    129         }
    130         psFree(view);
    131 
    132         pswarpFileActivation(config, detectorFiles, true);
    133         pmFPAfileActivate(config->files, false, "PSWARP.ASTROM");
    134     }
    135 
    13681    // Turn on the source output --- we need to get rid of these so that we can measure the PSF
    137     pmFPAfileActivate(config->files, true, "PSWARP.OUTPUT.SOURCES");
    138 
    139     // Don't care about the skycell anymore --- we've read it, and that's all we need to do.
    140     pmFPAfileActivate(config->files, false, "PSWARP.SKYCELL");
    141     view = pmFPAviewAlloc(0);
    142 
    143     // find the FPA phu
    144     bool bilevelAstrometry = false;
    145     pmHDU *phu = pmFPAviewThisPHU(view, astrom->fpa);
    146     if (phu) {
    147         char *ctype = psMetadataLookupStr(NULL, phu->header, "CTYPE1");
    148         if (ctype) {
    149             bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
    150         }
    151     }
    152     if (bilevelAstrometry) {
    153         if (!pmAstromReadBilevelMosaic(input->fpa, phu->header)) {
    154             psError(psErrorCodeLast(), false, "Unable to read bilevel mosaic astrometry for input FPA.");
    155             psFree(view);
    156             goto DONE;
    157         }
    158     }
    159 
    160     psList *cells = psListAlloc(NULL);  // List of cells, for concepts averaging
     82    // XXX NOTE : do not tackle source issues yet...
     83    // XXX pmFPAfileActivate(config->files, true, "PSWARP.OUTPUT.SOURCES");
     84
     85    // Done with the detector side of things
     86    pswarpFileActivation(config, detectorFiles, true);
     87    pswarpFileActivation(config, independentFiles, true);
     88
     89    pmFPAview *view = pmFPAviewAlloc(0);
     90
     91    // XXX : need to generate a list of input cells for each output cell
     92    // XXX : psList *cells = psListAlloc(NULL);  // List of cells, for concepts averaging
    16193
    16294    // files associated with the science image
    … …  
    16698    }
    16799
    168     // *** main transformation block
     100    // *** main transformation block
     101    // *** this section loops over the input chips/cells and reads them one at a time
     102    // *** the output chips/cells are filled where appropriate, but not yet written to disk
    169103    pmChip *chip;
    170104    while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
    … …  
    176110        }
    177111
    178         // read WCS data from the corresponding header
    179         pmHDU *hdu = pmFPAviewThisHDU (view, astrom->fpa);
    180 
    181        
    182         if (bilevelAstrometry) {
    183             if (!pmAstromReadBilevelChip (chip, hdu->header)) {
    184                 psError(psErrorCodeLast(), false, "Unable to read bilevel chip astrometry for input FPA.");
    185                 psFree(view);
    186                 goto DONE;
    187             }
    188         } else {
    189             // we use a default FPA pixel scale of 1.0
    190             if (!pmAstromReadWCS (input->fpa, chip, hdu->header, 1.0)) {
    191                 psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for input FPA.");
    192                 psFree(view);
    193                 goto DONE;
    194             }
    195         }
    196        
    197112        pmCell *cell;
    198113        while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
    … …  
    203118                goto DONE;
    204119            }
    205 
    206             psListAdd(cells, PS_LIST_TAIL, cell);
    207120
    208121            // process each of the readouts
    … …  
    217130                }
    218131
     132/// XXX ignore sources for now
     133# if (0)
    219134                // Copy the detections from the astrometry carrier to the input, so they can be accessed by
    220135                // pswarpTransformReadout
    … …  
    224139                    psMetadataAddPtr(readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_ARRAY, "Sources from input astrometry", detections);
    225140                }
    226 
    227                 pswarpTransformReadout(output, readout, config);
     141# endif
     142
     143                pswarpTransformToTarget(output->fpa, readout, config);
    228144               
    229145                if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
    … …  
    242158        }
    243159    }
    244 
    245     if (!output->data_exists) {
    246         psWarning("No overlap between input and skycell.");
    247         if (stats) {
    248             psMetadataAddS32(stats, PS_LIST_TAIL, "QUALITY", PS_META_REPLACE,
    249                              "No overlap between input and skycell", PSWARP_ERR_NO_OVERLAP);
    250         }
    251         psphotFilesActivate(config, false);
    252         psFree(cells);
    253         psFree(view);
    254         goto DONE;
    255     }
    256    
    257     pmCell *outCell = output->parent;   ///< Output cell
    258     pmChip *outChip = outCell->parent;  ///< Output chip
    259     pmFPA *outFPA = outChip->parent;    ///< Output FP
    260 
    261     if (!pswarpPixelsLit(output, stats, config)) {
    262         psError(psErrorCodeLast(), false, "Unable to calculate pixel regions.");
    263         psFree(cells);
    264         psFree(view);
    265         goto DONE;
    266     }
    267     bool doStats = psMetadataLookupBool(&mdok,recipe,"MASK.STATS");
    268     if (doStats) {
    269       if (!pswarpMaskStats(output, stats, config)) {
    270         psError(psErrorCodeLast(), false, "Unable to calculate mask stats.");
    271         psFree(cells);
    272         psFree(view);
    273         goto DONE;
    274       }
    275     }
    276     // Set covariance matrix for output
    277     {
    278         psList *covariances = psMetadataLookupPtr(&mdok, output->analysis,
    279                                                   PSWARP_ANALYSIS_COVARIANCES); // Covariance matrices
    280         psAssert(covariances, "Should be there");
    281         psArray *covars = psListToArray(covariances); // Array of covariance matrices
    282         psKernel *covar = psImageCovarianceAverage(covars);
    283         psFree(covars);
    284         psMetadataRemoveKey(output->analysis, PSWARP_ANALYSIS_COVARIANCES);
    285 
    286         // Correct covariance matrix scale for the mean (square root of the) Jacobian
    287         double jacobian = psMetadataLookupF64(NULL, output->analysis, PSWARP_ANALYSIS_JACOBIAN); // Jacobian
    288         int goodPixels = psMetadataLookupS32(NULL, output->analysis, PSWARP_ANALYSIS_GOODPIX);   // Good pixels
    289         jacobian /= goodPixels;
    290         output->covariance = psImageCovarianceScale(covar, jacobian);
    291         psFree(covar);
    292 
    293         if (output->variance) {
    294             psImageCovarianceTransfer(output->variance, output->covariance);
    295         }
    296     }
    297 
    298     if (!pmConceptsAverageCells(outCell, cells, NULL, NULL, false)) {
    299         psError(psErrorCodeLast(), false, "Unable to average cell concepts.");
    300         psFree(cells);
    301         psFree(view);
    302         goto DONE;
    303     }
    304     psFree(cells);
    305 
    306     psRegion *trimsec = psMetadataLookupPtr(NULL, outCell->concepts, "CELL.TRIMSEC"); ///< Trim section
    307     trimsec->x0 = trimsec->x1 = trimsec->y0 = trimsec->y1 = 0; ///< All pixels
    308 
    309     if (!psMetadataCopy(outFPA->concepts, input->fpa->concepts)) {
    310         psError(psErrorCodeLast(), false, "Unable to copy FPA concepts from input to output.");
    311         psFree(view);
    312         goto DONE;
    313     }
    314 
    315     // Update ZP from the astrometry
    316     {
    317         psMetadataItem *item = psMetadataLookup(outFPA->concepts, "FPA.ZP");
    318         item->data.F32 = psMetadataLookupF32(NULL, astrom->fpa->concepts, "FPA.ZP");
    319     }
    320 
    321     pmHDU *hdu = outFPA->hdu;           ///< HDU for the output warped image
    322 
    323     // Copy header from target
    324     {
    325         pmFPAview *skyView = pmFPAviewAlloc(0); ///< View into skycell
    326         skyView->chip = skyView->cell = 0;
    327         pmCell *cell = pmFPAfileThisCell(config->files, skyView, "PSWARP.SKYCELL"); // Skycell cell
    328         psFree(skyView);
    329         pmHDU *skyHDU = pmHDUFromCell(cell); ///< HDU
    330         if (!skyHDU) {
    331             psError(PSWARP_ERR_DATA, false, "Unable to find skycell HDU.");
    332             psFree(view);
    333             goto DONE;
    334         }
    335         hdu->header = psMetadataCopy(hdu->header, skyHDU->header);
    336     }
    337 
    338     pswarpVersionHeader(hdu->header);
    339    
    340     if (!pmAstromWriteWCS(hdu->header, outFPA, outChip, WCS_NONLIN_TOL)) {
    341         psError(psErrorCodeLast(), false, "Unable to generate WCS header.");
    342         goto DONE;
    343     }
    344 
    345160    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
    346161        psError(psErrorCodeLast(), false, "Unable to write files.");
    … …  
    352167    pswarpFileActivation(config, independentFiles, false);
    353168
    354 
    355     // We need a new PSF model for the warped frame.  It would be good to generate this analytically, but
    356     // that's going to be tricky.  We have a list of sources, so we use those to redetermine the PSF model.
     169    // XXXX turn these on later .. pswarpUpdateStatistics (output);
     170    // XXXX turn these on later .. pswarpMakePSF (output);
     171
     172    // XXX have not yet written the output to disk
     173
     174    psFree(view);
     175
     176 DONE:
     177    return true;
     178}
     179
     180# if (0)
     181// once the output fpa elements have been built, loop over the fpa and generate stats
     182// for each readout
     183bool pswarpUpdateStatistics (pmFPA *output, pmFPA *input, pmConfig *config)  {
     184
     185    // output mask bits
     186    psImageMaskType maskValue = psMetadataLookupImageMask(&status, recipe, "MASK.OUTPUT");
     187    psAssert (status, "MASK.OUTPUT was not defined");
     188
     189    pmFPAview *view = pmFPAviewAlloc(0);
     190   
     191    pmChip *chip;
     192    while ((chip = pmFPAviewNextChip (view, output, 1)) != NULL) {
     193        psTrace ("pswarp", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
     194        if (!chip->process || !chip->file_exists) { continue; }
     195
     196        pmCell *cell;
     197        while ((cell = pmFPAviewNextCell (view, output->fpa, 1)) != NULL) {
     198            psTrace ("pswarp", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
     199            if (!cell->process || !cell->file_exists) { continue; }
     200
     201            // process each of the readouts
     202            pmReadout *readout;
     203            while ((readout = pmFPAviewNextReadout(view, output->fpa, 1)) != NULL) {
     204                if (!readout->data_exists) {
     205                    psWarning("No overlap between input and skycell.");
     206                    if (stats) {
     207                        psMetadataAddS32(stats, PS_LIST_TAIL, "QUALITY", PS_META_REPLACE,
     208                                         "No overlap between input and skycell", PSWARP_ERR_NO_OVERLAP);
     209                    }
     210                    psphotFilesActivate(config, false);
     211                    psFree(cells);
     212                    psFree(view);
     213                    goto DONE;
     214                }
     215   
     216                if (!pswarpPixelsLit(readout, stats, config)) {
     217                    psError(psErrorCodeLast(), false, "Unable to calculate pixel regions.");
     218                    psFree(view);
     219                    goto DONE;
     220                }
     221                bool doStats = psMetadataLookupBool(&mdok,recipe,"MASK.STATS");
     222                if (doStats) {
     223                    if (!pswarpMaskStats(readout, stats, config)) {
     224                        psError(psErrorCodeLast(), false, "Unable to calculate mask stats.");
     225                        psFree(cells);
     226                        psFree(view);
     227                        goto DONE;
     228                    }
     229                }
     230
     231                // Set covariance matrix for output
     232                {
     233                    psList *covariances = psMetadataLookupPtr(&mdok, output->analysis,
     234                                                              PSWARP_ANALYSIS_COVARIANCES); // Covariance matrices
     235                    psAssert(covariances, "Should be there");
     236                    psArray *covars = psListToArray(covariances); // Array of covariance matrices
     237                    psKernel *covar = psImageCovarianceAverage(covars);
     238                    psFree(covars);
     239                    psMetadataRemoveKey(output->analysis, PSWARP_ANALYSIS_COVARIANCES);
     240
     241                    // Correct covariance matrix scale for the mean (square root of the) Jacobian
     242                    double jacobian = psMetadataLookupF64(NULL, output->analysis, PSWARP_ANALYSIS_JACOBIAN); // Jacobian
     243                    int goodPixels = psMetadataLookupS32(NULL, output->analysis, PSWARP_ANALYSIS_GOODPIX);   // Good pixels
     244                    jacobian /= goodPixels;
     245                    output->covariance = psImageCovarianceScale(covar, jacobian);
     246                    psFree(covar);
     247
     248                    if (output->variance) {
     249                        psImageCovarianceTransfer(output->variance, output->covariance);
     250                    }
     251                }
     252            }
     253
     254            // XXX need to get the input Cells for this output cell (save on the analysis somewhere
     255            // XXX skip this for now
     256# if (0)
     257            psArray *inputCells = psMetadataLookupPtr (NULL, cell->analysis, "something");
     258            if (!pmConceptsAverageCells(cell, inputCells, NULL, NULL, false)) {
     259                psError(psErrorCodeLast(), false, "Unable to average cell concepts.");
     260                psFree(view);
     261                goto DONE;
     262            }
     263# endif
     264            // XXX Is this too ad-hoc?
     265            psRegion *trimsec = psMetadataLookupPtr(NULL, cell->concepts, "CELL.TRIMSEC"); ///< Trim section
     266            trimsec->x0 = trimsec->x1 = trimsec->y0 = trimsec->y1 = 0; ///< All pixels
     267
     268        }
     269
     270        // Perform statistics on the output image
     271        if (stats) {
     272            // XXX this should be done at each level...
     273            if (!ppStatsFPA(stats, output->parent->parent->parent, view, maskValue, config)) {
     274                psWarning("Unable to perform statistics on warped image.");
     275            }
     276        }
     277    }
     278    if (!psMetadataCopy(output->concepts, input->concepts)) {
     279        psError(psErrorCodeLast(), false, "Unable to copy FPA concepts from input to output.");
     280        psFree(view);
     281        goto DONE;
     282    }
     283
     284    // Update ZP from the astrometry
     285    {
     286        psMetadataItem *item = psMetadataLookup(outFPA->concepts, "FPA.ZP");
     287        item->data.F32 = psMetadataLookupF32(NULL, astrom->fpa->concepts, "FPA.ZP");
     288    }
     289
     290    // need to update the chip and fpa level astrometry appropriate, see
     291    // code in psastro
     292   
     293# if (0)
     294    // Add MD5 information for readout
     295    const char *chipName = psMetadataLookupStr(NULL, output->parent->parent->concepts, "CHIP.NAME");
     296    const char *cellName = psMetadataLookupStr(NULL, output->parent->concepts, "CELL.NAME");
     297    psString headerName = NULL; ///< Header name for MD5
     298    psStringAppend(&headerName, "MD5_%s_%s_%d", chipName, cellName, view->readout);
     299    psVector *md5 = psImageMD5(output->image); ///< md5 hash
     300    psString md5string = psMD5toString(md5); ///< String
     301    psFree(md5);
     302    psMetadataAddStr(hdu->header, PS_LIST_TAIL, headerName, PS_META_REPLACE,
     303                     "Image MD5", md5string);
     304    psFree(md5string);
     305    psFree(headerName);
     306# endif
     307
     308# if (0)
     309    if (!pmAstromWriteWCS(hdu->header, outFPA, outChip, WCS_NONLIN_TOL)) {
     310        psError(psErrorCodeLast(), false, "Unable to generate WCS header.");
     311        goto DONE;
     312    }
     313# endif
     314    return true;
     315}
     316# endif
     317
     318# if (0)
     319// We need a new PSF model for the warped frame.  It would be good to generate this analytically, but
     320// that's going to be tricky.  We have a list of sources, so we use those to redetermine the PSF model.
     321bool pswarpMakePSF (pmConfig *config) {
     322
     323    // load the recipe
     324    psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSWARP_RECIPE);
     325    if (!recipe) {
     326        psError(PSWARP_ERR_CONFIG, false, "missing recipe %s", PSWARP_RECIPE);
     327        return false;
     328    }
    357329
    358330    if (psMetadataLookupBool(&mdok, recipe, "PSF")) {
    … …  
    419391        psFree(view);
    420392    }
    421 
    422     // Perform statistics on the output image
    423     if (stats) {
    424         if (!ppStatsFPA(stats, output->parent->parent->parent, view, maskValue, config)) {
    425             psWarning("Unable to perform statistics on warped image.");
    426         }
    427     }
    428    
    429 
    430     // Add MD5 information for readout
    431     const char *chipName = psMetadataLookupStr(NULL, output->parent->parent->concepts, "CHIP.NAME");
    432     const char *cellName = psMetadataLookupStr(NULL, output->parent->concepts, "CELL.NAME");
    433     psString headerName = NULL; ///< Header name for MD5
    434     psStringAppend(&headerName, "MD5_%s_%s_%d", chipName, cellName, view->readout);
    435     psVector *md5 = psImageMD5(output->image); ///< md5 hash
    436     psString md5string = psMD5toString(md5); ///< String
    437     psFree(md5);
    438     psMetadataAddStr(hdu->header, PS_LIST_TAIL, headerName, PS_META_REPLACE,
    439                      "Image MD5", md5string);
    440     psFree(md5string);
    441     psFree(headerName);
    442     psFree(view);
    443 
    444  DONE:
    445 
    446393    return true;
    447394}
     395# endif
     396
     397// once the output fpa elements have been built, loop over the fpa and generate stats
     398// for each readout
     399bool pswarpTransformToTarget (pmFPA *output, pmReadout *input, pmConfig *config)  {
     400
     401    pmFPAview *view = pmFPAviewAlloc(0);
     402   
     403    pmChip *chip;
     404    while ((chip = pmFPAviewNextChip (view, output, 1)) != NULL) {
     405        psTrace ("pswarp", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
     406        if (!chip->process || !chip->file_exists) { continue; }
     407
     408        pmCell *cell;
     409        while ((cell = pmFPAviewNextCell (view, output, 1)) != NULL) {
     410            psTrace ("pswarp", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
     411            if (!cell->process || !cell->file_exists) { continue; }
     412
     413            // process each of the readouts
     414            pmReadout *readout;
     415            while ((readout = pmFPAviewNextReadout(view, output, 1)) != NULL) {
     416                pswarpTransformReadout (readout, input, config);
     417            }
     418        }
     419    }
     420    return true;
     421}
     422
Note: See TracChangeset for help on using the changeset viewer.