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/pswarpDefine.c

    r34800 r35563  
    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
    … …  
    4238    }
    4339   
    44     // open the full skycell file; no need to defer different depths. only load the header data
    45     pmFPAview *view = pmFPAviewAlloc (0);
    46     pmFPAfileOpen (skycell, view, config);
     40    // Read the output astrometry
     41    // XXX this section loops over the input skycell to load the header info
     42    {
     43        pmFPAview *view = pmFPAviewAlloc(0);
    4744
    48     // Read header and create target
    49     {
    50         if (!pmFPAReadHeaderSet(skycell->fpa, skycell->fits, config)) {
    51             psError(psErrorCodeLast(), false, "Unable to read headers for skycell.");
    52             psFree(view);
    53             return false;
    54         }
    55         view->chip = 0;
    56         view->cell = 0;
    57         view->readout = 0;
    58         pmCell *source = pmFPAfileThisCell(config->files, view, "PSWARP.SKYCELL"); ///< Source cell
    59         pmHDU *hdu = pmHDUFromCell(source); ///< HDU for source
    60         if (!hdu || !hdu->header) {
    61             psError(PM_ERR_PROG, false, "Unable to find header for sky cell.");
    62             psFree(view);
    63             return false;
    64         }
    65         int numCols = psMetadataLookupS32(NULL, hdu->header, "NAXIS1"); ///< Number of columns
    66         int numRows = psMetadataLookupS32(NULL, hdu->header, "NAXIS2"); ///< Number of rows
    67         if ((numCols == 0) || (numRows == 0)) {
    68             psError(PSWARP_ERR_DATA, false, "skycell has invalid dimensions %d x %d", numCols, numRows);
    69             psFree(view);
    70             return false;
    71         }
     45        pmChip *chip;
     46        if (!pmFPAfileRead (skycell, view, config)) {
     47            psError(PS_ERR_IO, false, "failed READ at FPA %s", skycell->name);
     48            psFree(view);
     49            return false;
     50        }
     51        while ((chip = pmFPAviewNextChip (view, skycell->fpa, 1)) != NULL) {
     52            psTrace ("pswarp", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
     53            if (!chip->process || !chip->file_exists) { continue; }
     54            if (!pmFPAfileRead (skycell, view, config)) {
     55                psError(psErrorCodeLast(), false, "failed READ at CHIP %s", skycell->name);
     56                psFree(view);
     57                return false;
     58            }
     59            pmCell *cell;
     60            while ((cell = pmFPAviewNextCell (view, skycell->fpa, 1)) != NULL) {
     61                psTrace ("pswarp", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
     62                if (!cell->process || !cell->file_exists) { continue; }
     63                if (!pmFPAfileRead (skycell, view, config)) {
     64                    psError(psErrorCodeLast(), false, "failed READ at CELL %s", skycell->name);
     65                    psFree(view);
     66                    return false;
     67                }
     68                if (!pmFPAfileClose(skycell, view)) {
     69                    psError(psErrorCodeLast(), false, "failed CLOSE at CELL %s", skycell->name);
     70                    psFree (view);
     71                    return false;
     72                }
    7273
    73         pmCell *target = pmFPAviewThisCell(view, output->fpa); ///< Target cell
    74         pmReadout *readout = pmReadoutAlloc(target); ///< Target readout
    75         readout->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    76         psImageInit(readout->image, NAN);
    77         psFree(readout);                // Drop reference
    78 
    79         bool status = false;
    80         psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.XBIN");
    81         psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.YBIN");
    82         psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.XSIZE");
    83         psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.YSIZE");
    84         psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.XPARITY");
    85         psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.YPARITY");
    86         psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.X0");
    87         psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.Y0");
    88 
    89        
     74                // we've got the output astrom header
     75                pmHDU *hdu = pmHDUFromCell(cell); ///< HDU for source
     76                if (!hdu || !hdu->header) {
     77                    psError(PM_ERR_PROG, false, "Unable to find header for sky cell.");
     78                    psFree(view);
     79                    return false;
     80                }
     81                int numCols = psMetadataLookupS32(NULL, hdu->header, "NAXIS1"); ///< Number of columns
     82                int numRows = psMetadataLookupS32(NULL, hdu->header, "NAXIS2"); ///< Number of rows
     83                if ((numCols == 0) || (numRows == 0)) {
     84                    psError(PSWARP_ERR_DATA, false, "skycell has invalid dimensions %d x %d", numCols, numRows);
     85                    psFree(view);
     86                    return false;
     87                }
     88               
     89                // XXX generate the output pixels
     90                pmCell *target = pmFPAviewThisCell(view, output->fpa); ///< Target cell
     91                pmReadout *readout = pmReadoutAlloc(target); ///< Target readout
     92                readout->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     93                psImageInit(readout->image, NAN);
     94                psFree(readout);                // Drop reference
     95               
     96                bool status = false;
     97                psMetadataItemSupplement(&status, target->concepts, cell->concepts, "CELL.XBIN");
     98                psMetadataItemSupplement(&status, target->concepts, cell->concepts, "CELL.YBIN");
     99                psMetadataItemSupplement(&status, target->concepts, cell->concepts, "CELL.XSIZE");
     100                psMetadataItemSupplement(&status, target->concepts, cell->concepts, "CELL.YSIZE");
     101                psMetadataItemSupplement(&status, target->concepts, cell->concepts, "CELL.XPARITY");
     102                psMetadataItemSupplement(&status, target->concepts, cell->concepts, "CELL.YPARITY");
     103                psMetadataItemSupplement(&status, target->concepts, cell->concepts, "CELL.X0");
     104                psMetadataItemSupplement(&status, target->concepts, cell->concepts, "CELL.Y0");
     105            }
     106            if (!pmFPAfileClose(skycell, view)) {
     107                psError(psErrorCodeLast(), false, "failed CLOSE at CHIP %s", skycell->name);
     108                psFree (view);
     109                return false;
     110            }
     111        }
     112        if (!pmFPAfileClose(skycell, view)) {
     113            psError(psErrorCodeLast(), false, "failed CLOSE at FPA %s", skycell->name);
     114            psFree (view);
     115            return false;
     116        }
     117        psFree(view);
    90118    }
    91119
    92     // XXX this is not a sufficient test
     120    // the section below converts the header astrometry info (in skycell) to the pmFPA
     121    // astrometry structures, saving them on the output fpa structure.
     122
     123    // XXX this block should be modified to loop over all chips
     124
     125    pmFPAview *view = pmFPAviewAlloc (0);
    93126    view->chip = 0;
    94127    view->cell = 0;
    … …  
    134167}
    135168
    136 bool pswarpDefineBackground (pmConfig *config) {
    137 
    138     // load the PSWARP recipe
    139     psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSWARP_RECIPE);
    140     if (!recipe) {
    141         psError(PSWARP_ERR_CONFIG, true, "Can't find PSWARP recipe!\n");
    142         return false;
    143     }
    144 
    145     // select the input data sources
    146     pmFPAfile *skycell = psMetadataLookupPtr (NULL, config->files, "PSWARP.SKYCELL");
    147     if (!skycell) {
    148         psError(PSWARP_ERR_CONFIG, false, "Can't find skycell data!\n");
    149         return false;
    150     }
    151     pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PSWARP.OUTPUT.BKGMODEL");
    152     if (!output) {
    153         psError(PSWARP_ERR_CONFIG, false, "Can't find output data!\n");
    154         return false;
    155     }
    156     pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PSWARP.INPUT");
    157     if (!input) {
    158         psError(PSWARP_ERR_CONFIG, false, "Can't find input data!\n");
    159         return false;
    160     }
    161    
    162     // open the full skycell file; no need to defer different depths. only load the header data
    163     pmFPAview *view = pmFPAviewAlloc (0);
    164     pmFPAfileOpen (skycell, view, config);
    165 
    166     // Read header and create target
    167     {
    168         if (!pmFPAReadHeaderSet(skycell->fpa, skycell->fits, config)) {
    169             psError(psErrorCodeLast(), false, "Unable to read headers for skycell.");
    170             psFree(view);
    171             return false;
    172         }
    173         view->chip = 0;
    174         view->cell = 0;
    175         view->readout = 0;
    176         pmCell *source = pmFPAfileThisCell(config->files, view, "PSWARP.SKYCELL"); ///< Source cell
    177         pmHDU *hdu = pmHDUFromCell(source); ///< HDU for source
    178         if (!hdu || !hdu->header) {
    179             psError(PM_ERR_PROG, false, "Unable to find header for sky cell.");
    180             psFree(view);
    181             return false;
    182         }
    183         int numCols = psMetadataLookupS32(NULL, hdu->header, "NAXIS1"); ///< Number of columns
    184         int numRows = psMetadataLookupS32(NULL, hdu->header, "NAXIS2"); ///< Number of rows
    185         if ((numCols == 0) || (numRows == 0)) {
    186             psError(PSWARP_ERR_DATA, false, "skycell has invalid dimensions %d x %d", numCols, numRows);
    187             psFree(view);
    188             return false;
    189         }
    190 
    191         pmCell *target = pmFPAviewThisCell(view, output->fpa); ///< Target cell
    192         pmReadout *readout = pmReadoutAlloc(target); ///< Target readout
    193         readout->image = psImageAlloc(numCols / output->xBin, numRows / output->yBin, PS_TYPE_F32);
    194         psImageInit(readout->image, NAN);
    195         psFree(readout);                // Drop reference
    196 
    197         bool status = false;
    198         psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.XBIN");
    199         psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.YBIN");
    200 /*      psMetadataAddS32(target->concepts,PS_LIST_TAIL,"CELL.XBIN",PS_META_REPLACE,"",output->xBin); */
    201 /*      psMetadataAddS32(target->concepts,PS_LIST_TAIL,"CELL.YBIN",PS_META_REPLACE,"",output->yBin); */
    202         psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.XSIZE");
    203         psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.YSIZE");
    204 /*      psMetadataAddS32(target->concepts,PS_LIST_TAIL,"CELL.XSIZE",PS_META_REPLACE,"",numCols / output->xBin); */
    205 /*      psMetadataAddS32(target->concepts,PS_LIST_TAIL,"CELL.YSIZE",PS_META_REPLACE,"",numRows / output->yBin); */
    206         psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.XPARITY");
    207         psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.YPARITY");
    208         psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.X0");
    209         psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.Y0");
    210 
    211        
    212     }
    213 
    214     // XXX this is not a sufficient test
    215     view->chip = 0;
    216     view->cell = 0;
    217     view->readout = -1;
    218     pmHDU *phu = pmFPAviewThisPHU(view, skycell->fpa); ///< Skycell PHU
    219     pmHDU *hdu = pmFPAviewThisHDU(view, skycell->fpa); ///< Skycell header
    220 
    221     pmAstromWCS *WCS = pmAstromWCSfromHeader(hdu->header);
    222 
    223     double cd1f = 1.0 * output->xBin;
    224     double cd2f = 1.0 * output->yBin;
    225    
    226     WCS->crpix1 = WCS->crpix1 / cd1f;
    227     WCS->crpix2 = WCS->crpix2 / cd2f;
    228    
    229     WCS->cdelt1 *= cd1f;
    230     WCS->cdelt2 *= cd2f;
    231 
    232     WCS->trans->x->coeff[1][0] *= cd1f;
    233     WCS->trans->x->coeff[0][1] *= cd2f;
    234     WCS->trans->y->coeff[1][0] *= cd1f;
    235     WCS->trans->y->coeff[0][1] *= cd2f;
    236    
    237    
    238     pmAstromWCStoHeader (hdu->header,WCS);
    239 
    240     bool bilevelAstrometry = false;
    241     if (phu) {
    242         char *ctype = psMetadataLookupStr(NULL, phu->header, "CTYPE1");
    243         if (ctype) {
    244             bilevelAstrometry = !strcmp(&ctype[4], "-DIS");
    245         }
    246     }
    247 
    248     // We read from the skycell into the output.  i.e., the output receives the desired astrometry.
    249     pmChip *outputChip = pmFPAviewThisChip(view, output->fpa); ///< Chip in the output
    250     if (bilevelAstrometry) {
    251         if (!pmAstromReadBilevelMosaic(output->fpa, phu->header)) {
    252             psError(psErrorCodeLast(), false, "Unable to read bilevel mosaic astrometry for skycell.");
    253             psFree(view);
    254             return false;
    255         }
    256         if (!pmAstromReadBilevelChip(outputChip, hdu->header)) {
    257             psError(psErrorCodeLast(), false, "Unable to read bilevel chip astrometry for skycell.");
    258             psFree(view);
    259             return false;
    260         }
    261     } else {
    262         // we use a default FPA pixel scale of 1.0
    263       if (!pmAstromReadWCS(output->fpa, outputChip, hdu->header, 1.0)) {
    264             psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for skycell.");
    265             psFree(view);
    266             return false;
    267         }
    268     }
    269 
    270    
    271     view->chip = view->cell = view->readout = -1;
    272     pmFPAAddSourceFromView(output->fpa, view, output->format);
    273 
    274 
    275     psFree (view);
    276     return true;
    277 }
Note: See TracChangeset for help on using the changeset viewer.