IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
May 8, 2013, 4:41:42 PM (13 years ago)
Author:
eugene
Message:

warping of background model appears to work

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r35421 r35535  
    1111
    1212#include "pswarp.h"
    13 #include <ppStats.h>
    1413#include "pswarpFileNames.h"            // Lists of file rules used at different stages
    1514
    … …  
    2221    bool status;
    2322    bool mdok;                          // Status of MD lookup
    24     const char *skyCamera = psMetadataLookupStr(NULL, config->arguments,
    25                                                 "SKYCELL.CAMERA");  // Name of camera for skycell
     23
    2624    if (!psMetadataLookupBool(NULL,config->arguments,"BACKGROUND.MODEL")) return true;
    2725
    28     pmConfigCamerasCull(config, skyCamera);
    29     pmConfigRecipesCull(config, "PSWARP,PPSTATS,PSPHOT,PSASTRO,MASKS,JPEG");
     26    int xGrid = psMetadataLookupS32(&mdok,config->arguments,"BKG.XGRID");
     27    int yGrid = psMetadataLookupS32(&mdok,config->arguments,"BKG.YGRID");
     28
     29    // WHAT IS THIS?? should it be disabled at the end?
     30    psMetadataAddS32(config->arguments, PS_LIST_TAIL, "INTERPOLATION.MODE", PS_META_REPLACE, "", 8);
    3031
    3132    // load the recipe
    … …  
    3637    }
    3738
    38     if (!pswarpSetMaskBits(config)) {
    39         psError(psErrorCodeLast(), false, "failed to set mask bits");
    40         return NULL;
    41     }
    42 
    4339    // select the input data sources
    44     pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PSWARP.BKGMODEL");
    45     if (!input) {
    46         psError(PSWARP_ERR_CONFIG, true, "Can't find input data!\n");
     40    pmFPAfile *output = psMetadataLookupPtr(&status, config->files, "PSWARP.OUTPUT.BKGMODEL");
     41    if (!output) {
     42        psError(PSWARP_ERR_CONFIG, true, "Can't find output data!\n");
    4743        return false;
    4844    }
    4945
    5046    // use the external astrometry source if supplied
    51     pmFPAfile *astrom = psMetadataLookupPtr(NULL, config->files, "PSWARP.ASTROM");
    52     if (!astrom) {
    53         astrom = input;
    54     }
    55 
    56     if (astrom->camera != input->camera) {
    57         psError(PSWARP_ERR_DATA, true, "Input camera and astrometry camera do not match.");
    58         return false;
    59     }
    60 
    61     // select the output readout
     47    pmFPAfile *skycell = psMetadataLookupPtr(&status, config->files, "PSWARP.SKYCELL");
     48    if (!skycell) {
     49        psError(PSWARP_ERR_DATA, true, "Cannot find output astrometry.");
     50        return false;
     51    }
     52
    6253    pmFPAview *view = pmFPAviewAlloc(0);
    63     view->chip = 0;
    64     view->cell = 0;
    65     view->readout = 0;
    66     pmReadout *output = pmFPAfileThisReadout(config->files, view, "PSWARP.OUTPUT.BKGMODEL");
    67     if (!output) {
    68         psError(PSWARP_ERR_CONFIG, true, "Can't find output background data!\n");
    69         return false;
    70     }
     54
     55    int nInputs = psMetadataLookupS32(&status, config->arguments, "NUM_INPUTS");
     56    if (!status) {
     57        psError(PSWARP_ERR_DATA, true, "number of inputs is not defined (programming error)");
     58        return false;
     59    }
     60
     61    // loop over this section once per input group
     62    for (int i = 0; i < nInputs; i++) {
     63
     64        // select the input data source : we are reading & transforming pixels for the background model
     65        pmFPAfile *input = pmFPAfileSelectSingle(config->files, "PSWARP.BKGMODEL", i);
     66        if (!input) {
     67            psError(PSWARP_ERR_CONFIG, true, "Can't find input data!\n");
     68            return false;
     69        }
     70
     71        // select the input data sources
     72        pmFPAfile *astrom = pmFPAfileSelectSingle(config->files, "PSWARP.ASTROM", i);
     73        if (!astrom) {
     74            astrom = input;
     75        }
     76
     77        // check input astrometry mode
     78        bool bilevelAstrometry = psMetadataLookupBool (NULL, astrom->fpa->analysis, "ASTROMETRY.BILEVEL");
     79        if (bilevelAstrometry) {
     80            // top-level elements of the input model astrometry
     81            input->fpa->toTPA = psMemIncrRefCounter (astrom->fpa->toTPA);
     82            input->fpa->fromTPA = psMemIncrRefCounter (astrom->fpa->fromTPA);
     83            input->fpa->toSky = psMemIncrRefCounter (astrom->fpa->toSky);
     84        }
     85
     86        // we only need to load the background models here
     87        pswarpFileActivation(config, detectorFiles, false);
     88        pswarpFileActivation(config, photFiles, false);
     89        pswarpFileActivation(config, independentFiles, false);
     90        pswarpFileActivation(config, skycellFiles, false);
     91
     92        // we are only reading the background models, not writing anything
     93        pmFPAfileActivate(config->files, true, "PSWARP.BKGMODEL");
     94
     95        pmFPAviewReset (view);
     96
     97        pmFPAview *view = pmFPAviewAlloc(0);
     98        if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
     99            psError(psErrorCodeLast(), false, "Unable to read files.");
     100            goto FAIL;
     101        }
     102
     103        pmChip *chip;
     104        while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
     105            if (!chip->process || !chip->file_exists) { continue; }
     106            if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
     107                psError(psErrorCodeLast(), false, "Unable to read files.");
     108                goto FAIL;
     109            }
     110
     111            // adjust output astrometry
     112            double xBin = NAN;
     113            double yBin = NAN;
     114            pswarpGetInputScales (&xBin, &yBin, config, view, chip);
     115            double fluxRatio = xBin * yBin / (float) (xGrid * yGrid);
     116           
     117            if (!pswarpModifyChipAstrom (config, view, chip, astrom, bilevelAstrometry, xBin, yBin)) {
     118                psError(psErrorCodeLast(), false, "failed to set BKGMODEL astrometry.");
     119                goto FAIL;
     120            }
     121
     122            pmCell *cell;
     123            while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
     124                psTrace ("pswarp", 4, "ACell %d: %x %x %d\n", view->cell, cell->file_exists, cell->process,psErrorCodeLast());
     125                if (!cell->process || !cell->file_exists) { continue; }
     126                if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
     127                    psError(psErrorCodeLast(), false, "Unable to read files.");
     128                    goto FAIL;
     129                }
     130
     131                // process each of the readouts
     132                pmReadout *readout;
     133                while ((readout = pmFPAviewNextReadout(view, input->fpa, 1)) != NULL) {
     134                    if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
     135                        psError(psErrorCodeLast(), false, "Unable to read files.");
     136                        goto FAIL;
     137                    }
     138                    if (!readout->data_exists) {
     139                        continue;
     140                    }
     141
     142                    // re-normalize the BKGMODEL pixels by modified astrometry
     143                    for (int x = 0; x < readout->image->numCols; x++) {
     144                        for (int y = 0; y < readout->image->numRows; y++) {
     145                            readout->image->data.F32[y][x] *= fluxRatio;
     146                        }
     147                    }
     148               
     149                    // transform the actual BKGMODEL pixels from readout to output
     150                    // XXX replace the use of this temporary metadata item with an argument?
     151                    // XXX move the source transform out of pswarpTransformReadout?
     152                    psMetadataAddBool(config->arguments,PS_LIST_TAIL, "BACKGROUND_WARPING", PS_META_REPLACE, "", true);
     153                    pswarpTransformToTarget (output->fpa, readout, config);
     154                    psMetadataAddBool(config->arguments,PS_LIST_TAIL, "BACKGROUND_WARPING", PS_META_REPLACE, "", false);
     155
     156                    if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) { // Readout
     157                        psError(psErrorCodeLast(), false, "Unable to close files.");
     158                        goto FAIL;
     159                    }
     160                }
     161
     162                if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) { // Cell
     163                    psError(psErrorCodeLast(), false, "Unable to close files.");
     164                    goto FAIL;
     165                }
     166            }
     167
     168            if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) { // Chip
     169                psError(psErrorCodeLast(), false, "Unable to close files.");
     170                goto FAIL;
     171            }
     172        }
     173
     174        if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) { // FPA
     175            psError(psErrorCodeLast(), false, "Unable to close files.");
     176            goto FAIL;
     177        }
     178
     179        if (false && !pswarpUpdateMetadata (output->fpa, skycell->fpa, input->fpa, astrom->fpa, config, false)) {
     180            psError(psErrorCodeLast(), false, "problem generating statistics.");
     181            goto FAIL;
     182        }
     183    }
     184
    71185    psFree (view);
    72     // Turn all skycell files on to generate them, and then turn them off for the loop over the input images
    73     // the input, which is in a different format.
    74     {
    75         pswarpFileActivation(config, detectorFiles, false);
    76         pswarpFileActivation(config, photFiles, false);
    77         pswarpFileActivation(config, independentFiles, false);
    78         pswarpFileActivation(config, skycellFiles, true);
    79         if (!pswarpIOChecksBefore(config)) {
    80             psError(psErrorCodeLast(), false, "Unable to read files.");
    81             goto DONE;
    82         }
    83         pswarpFileActivation(config, skycellFiles, false);
    84     }
    85     // Read the input astrometry
    86     // XXX rather than use the activations here, this should just explicitly loop over the desired filerule
    87     {
    88 
    89       pmFPAfileActivate(config->files, true, "PSWARP.ASTROM");
    90 
    91         pmChip *chip;
    92         pmFPAview *view = pmFPAviewAlloc(0);
    93         if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
    94             psError(psErrorCodeLast(), false, "Unable to read files.");
    95             goto DONE;
    96         }
    97 
    98         while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
    99 #if 0
    100           // This needs to be removed because otherwise it throws an error of duplicate PSPHOT.DETECTIONS.
    101             if (!chip->process || !chip->file_exists) { continue; }
    102             if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
    103                 psError(psErrorCodeLast(), false, "Unable to read files.");
    104                 goto DONE;
    105             }
    106 #endif
    107             pmCell *cell;
    108             while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
    109               psTrace ("pswarp", 4, "ACell %d: %x %x %d\n", view->cell, cell->file_exists, cell->process,psErrorCodeLast());
    110                 if (!cell->process || !cell->file_exists) { continue; }
    111                 if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE) ||
    112                     !pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) {
    113                     psError(psErrorCodeLast(), false, "Unable to read files.");
    114                     goto DONE;
    115                 }
    116             }
    117             if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) {
    118                 psError(psErrorCodeLast(), false, "Unable to write 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         psFree(view);
    127         pswarpFileActivation(config, detectorFiles, true);
    128         pmFPAfileActivate(config->files, false, "PSWARP.ASTROM");
    129     }
    130 
    131     // Don't care about the skycell anymore --- we've read it, and that's all we need to do.
    132     pmFPAfileActivate(config->files, false, "PSWARP.SKYCELL");
    133     view = pmFPAviewAlloc(0);
    134 
    135     // find the FPA phu
    136     bool bilevelAstrometry = false;
    137     pmHDU *phu = pmFPAviewThisPHU(view, astrom->fpa);
    138 
    139     //    pmAstromWCS *WCSF = pmAstromWCSfromHeader(phu->header);
    140    
    141     if (phu) {
    142         char *ctype = psMetadataLookupStr(NULL, phu->header, "CTYPE1");
    143         if (ctype) {
    144             bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
    145         }
    146     }
    147     if (bilevelAstrometry) {
    148         if (!pmAstromReadBilevelMosaic(input->fpa, phu->header)) {
    149             psError(psErrorCodeLast(), false, "Unable to read bilevel mosaic astrometry for input FPA.");
    150             psFree(view);
    151             psFree(stats);
    152             goto DONE;
    153         }
    154     }
    155 
    156     psList *cells = psListAlloc(NULL);  // List of cells, for concepts averaging
    157 
    158     // files associated with the science image
    159     if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
    160         psError(psErrorCodeLast(), false, "Unable to read files.");
    161         goto DONE;
    162     }
    163     pmChip *chip;
    164     while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
    165         psTrace ("pswarp", 4, "DChip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
    166         if (!chip->process || !chip->file_exists) { continue; }
    167         if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
    168             psError(psErrorCodeLast(), false, "Unable to read files.");
    169             goto DONE;
    170         }
    171 
    172         // Pull information from the header of the background files so we can use it to set values.
    173         pmHDU *hdu = pmFPAviewThisHDU(view,input->fpa);
    174         psMetadata *header = hdu->header;
     186    return true;
     187
     188FAIL:
     189    psFree (view);
     190    return false;
     191}
     192
     193// **** input == (small) BKGMODEL, astrom == (fullsize) ASTROM
     194// use astrometry from input ref astrometry and BINNING info from input model to modify the astrometry
     195// elements of input chip.  scale is the area of an input pixel is raw image pixels
     196bool pswarpGetInputScales (double *xBin, double *yBin, pmConfig *config, pmFPAview *view, pmChip *chip) {
     197
     198    // Pull information from the header of the background files so we can use it to set values.
     199    pmHDU *hdu = pmFPAviewThisHDU(view, chip->parent);
     200    psMetadata *header = hdu->header;
    175201       
    176         int IMAXIS1 = psMetadataLookupS32(NULL,header,"IMNAXIS1");
    177         int IMAXIS2 = psMetadataLookupS32(NULL,header,"IMNAXIS2");
    178         int NAXIS1 = psMetadataLookupS32(NULL,header,"NAXIS1");
    179         int NAXIS2 = psMetadataLookupS32(NULL,header,"NAXIS2");
    180         char *CCDSUM = psMetadataLookupStr(NULL,header,"CCDSUM");
    181         int CCDSUM1 = atoi(strtok(CCDSUM," "));
    182         int CCDSUM2 = atoi(strtok(NULL," "));
     202    int IMAXIS1 = psMetadataLookupS32(NULL,header,"IMNAXIS1");
     203    int IMAXIS2 = psMetadataLookupS32(NULL,header,"IMNAXIS2");
     204    int NAXIS1 = psMetadataLookupS32(NULL,header,"NAXIS1");
     205    int NAXIS2 = psMetadataLookupS32(NULL,header,"NAXIS2");
     206    char *CCDSUM = psMetadataLookupStr(NULL,header,"CCDSUM");
     207    int CCDSUM1 = atoi(strtok(CCDSUM," "));
     208    int CCDSUM2 = atoi(strtok(NULL," "));
    183209       
    184         psMetadataAddF32(config->arguments,PS_LIST_TAIL,"BKG_WARP_XOFFSET", PS_META_REPLACE,
    185                          "xoffset for background model data", (NAXIS1 * CCDSUM1 - IMAXIS1) / (2.0 * CCDSUM1));
    186         psMetadataAddF32(config->arguments,PS_LIST_TAIL,"BKG_WARP_YOFFSET", PS_META_REPLACE,
    187                          "yoffset for background model data", (NAXIS2 * CCDSUM2 - IMAXIS2) / (2.0 * CCDSUM2));
    188         psTrace("pswarp",5,"%d %d %d %d %d %d %g %g %d %d",
    189                 psMetadataLookupS32(NULL,header,"IMNAXIS1"),
    190                 psMetadataLookupS32(NULL,header,"IMNAXIS2"),
    191                 psMetadataLookupS32(NULL,header,"NAXIS1"),
    192                 psMetadataLookupS32(NULL,header,"NAXIS2"),
    193                 CCDSUM1,CCDSUM2,
    194                 psMetadataLookupF32(NULL,config->arguments,"BKG_WARP_XOFFSET"),
    195                 psMetadataLookupF32(NULL,config->arguments,"BKG_WARP_YOFFSET"),
    196                 psMetadataLookupS32(NULL,config->arguments,"BKG.XGRID"),
    197                 psMetadataLookupS32(NULL,config->arguments,"BKG.YGRID"));
    198        
    199        
    200         // read WCS data from the corresponding header
    201         hdu = pmFPAviewThisHDU (view, astrom->fpa);
    202 
    203         pmAstromWCS *WCS = pmAstromWCSfromHeader(hdu->header);
    204 
    205         double cd1f = (1.0 * CCDSUM1);// * (1.0 * psMetadataLookupS32(NULL,config->arguments,"BKG.XGRID"));
    206         double cd2f = (1.0 * CCDSUM2);// * (1.0 * psMetadataLookupS32(NULL,config->arguments,"BKG.YGRID"));
    207 
    208         WCS->cdelt1 *= cd1f;
    209         WCS->cdelt2 *= cd2f;
    210         WCS->crpix1 = WCS->crpix1 / cd1f;
    211         WCS->crpix2 = WCS->crpix2 / cd2f;
    212 
    213         // WCS->trans->x->nX/nY
    214         for (int q = 0; q <= WCS->trans->x->nX; q++) {
    215           for (int r = 0; r <= WCS->trans->x->nY; r++) {
    216             WCS->trans->x->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r);
    217           }
    218         }
    219         for (int q = 0; q <= WCS->trans->y->nX; q++) {
    220           for (int r = 0; r <= WCS->trans->y->nY; r++) {
    221             WCS->trans->y->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r);
    222           }
    223         }
    224         pmAstromWCStoHeader (hdu->header,WCS);
    225        
    226         if (bilevelAstrometry) {
    227             if (!pmAstromReadBilevelChip (chip, hdu->header)) {
    228                 psError(psErrorCodeLast(), false, "Unable to read bilevel chip astrometry for input FPA.");
    229                 psFree(view);
    230                 psFree(stats);
    231                 goto DONE;
    232             }
    233         } else {
    234             // we use a default FPA pixel scale of 1.0
    235             if (!pmAstromReadWCS (astrom->fpa, chip, hdu->header, 1.0)) {
    236                 psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for input FPA.");
    237                 psFree(view);
    238                 psFree(stats);
    239                 goto DONE;
    240             }
    241         }
    242         // Modify structure here.
    243 
    244        
    245         pmCell *cell;
    246         while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
    247             psTrace ("pswarp", 4, "DCell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
    248             if (!cell->process || !cell->file_exists) { continue; }
    249             if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
    250                 psError(psErrorCodeLast(), false, "Unable to read files.");
    251                 goto DONE;
    252             }
    253 
    254             psListAdd(cells, PS_LIST_TAIL, cell);
    255 
    256             // process each of the readouts
    257             pmReadout *readout;
    258             while ((readout = pmFPAviewNextReadout(view, input->fpa, 1)) != NULL) {
    259                 if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
    260                     psError(psErrorCodeLast(), false, "Unable to read files.");
    261                     goto DONE;
    262                 }
    263                 if (!readout->data_exists) {
    264                     continue;
    265                 }
    266 
    267                 for (int x = 0; x < readout->image->numCols; x++) {
    268                   for (int y = 0; y < readout->image->numRows; y++) {
    269                     readout->image->data.F32[y][x] = readout->image->data.F32[y][x] * (cd1f * cd2f) /
    270                       (psMetadataLookupS32(&mdok,config->arguments,"BKG.XGRID") *
    271                        psMetadataLookupS32(&mdok,config->arguments,"BKG.YGRID"));
    272                   }
    273                 }
    274                
    275                 psMetadataAddS32(config->arguments,PS_LIST_TAIL, "INTERPOLATION.MODE", PS_META_REPLACE, "", 8);
    276 
    277                 psMetadataAddBool(config->arguments,PS_LIST_TAIL, "BACKGROUND_WARPING", PS_META_REPLACE, "", true);
    278                 pswarpTransformReadout(output, readout, config);
    279                 psMetadataAddBool(config->arguments,PS_LIST_TAIL, "BACKGROUND_WARPING", PS_META_REPLACE, "", false);
    280 
    281                 if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
    282                     psError(psErrorCodeLast(), false, "Unable to write files.");
    283                     goto DONE;
    284                 }
    285             }
    286             if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
    287                 psError(psErrorCodeLast(), false, "Unable to write files.");
    288                 goto DONE;
    289             }
    290         }
    291         if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
    292             psError(psErrorCodeLast(), false, "Unable to write files.");
    293             goto DONE;
    294         }
    295     }
    296     if (!output->data_exists) {
    297         psWarning("No overlap between input and skycell.");
    298         psphotFilesActivate(config, false);
    299         psFree(cells);
    300         psFree(view);
    301         goto DONE;
    302     }
    303     pmCell *outCell = output->parent;   ///< Output cell
    304     pmChip *outChip = outCell->parent;  ///< Output chip
    305     pmFPA *outFPA = outChip->parent;    ///< Output FP
    306 
    307 /*     if (!pswarpPixelsLit(output, stats, config)) { */
    308 /*         psError(psErrorCodeLast(), false, "Unable to calculate pixel regions."); */
    309 /*         psFree(cells); */
    310 /*         psFree(view); */
    311 /*         goto DONE; */
    312 /*     } */
    313     psRegion *trimsec = psMetadataLookupPtr(NULL, outCell->concepts, "CELL.TRIMSEC"); ///< Trim section
    314     trimsec->x0 = trimsec->x1 = trimsec->y0 = trimsec->y1 = 0; ///< All pixels
    315 
    316     if (!psMetadataCopy(outFPA->concepts, input->fpa->concepts)) {
    317         psError(psErrorCodeLast(), false, "Unable to copy FPA concepts from input to output.");
    318         psFree(stats);
    319         psFree(view);
    320         goto DONE;
    321     }
    322 
    323     pmHDU *hdu = outFPA->hdu;           ///< HDU for the output warped image
    324 
    325     // Copy header from target
    326     {
    327         pmFPAview *skyView = pmFPAviewAlloc(0); ///< View into skycell
    328         skyView->chip = skyView->cell = 0;
    329         pmCell *cell = pmFPAfileThisCell(config->files, skyView, "PSWARP.SKYCELL"); // Skycell cell
    330         psFree(skyView);
    331         pmHDU *skyHDU = pmHDUFromCell(cell); ///< HDU
    332         if (!skyHDU) {
    333             psError(PSWARP_ERR_DATA, false, "Unable to find skycell HDU.");
    334             psFree(view);
    335             goto DONE;
    336         }
    337         hdu->header = psMetadataCopy(hdu->header, skyHDU->header);
    338     }
    339     pswarpVersionHeader(hdu->header);
    340    
    341     if (!pmAstromWriteWCS(hdu->header, outFPA, outChip, WCS_NONLIN_TOL)) {
    342         psError(psErrorCodeLast(), false, "Unable to generate WCS header.");
    343         psFree(stats);
    344         goto DONE;
    345     }
    346     if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
    347         psError(psErrorCodeLast(), false, "Unable to write files.");
    348         goto DONE;
    349     }
    350     // Done with the detector side of things
    351     pswarpFileActivation(config, detectorFiles, false);
    352     pswarpFileActivation(config, independentFiles, false);
    353 
    354     // Add MD5 information for readout
    355     const char *chipName = psMetadataLookupStr(NULL, output->parent->parent->concepts, "CHIP.NAME");
    356     const char *cellName = psMetadataLookupStr(NULL, output->parent->concepts, "CELL.NAME");
    357     psString headerName = NULL; ///< Header name for MD5
    358     psStringAppend(&headerName, "MD5_%s_%s_%d", chipName, cellName, view->readout);
    359     psVector *md5 = psImageMD5(output->image); ///< md5 hash
    360     psString md5string = psMD5toString(md5); ///< String
    361     psFree(md5);
    362     psMetadataAddStr(hdu->header, PS_LIST_TAIL, headerName, PS_META_REPLACE,
    363                      "Image MD5", md5string);
    364     psFree(md5string);
    365     psFree(headerName);
    366     psFree(view);
    367  DONE:
     210    psMetadataAddF32(config->arguments,PS_LIST_TAIL,"BKG_WARP_XOFFSET", PS_META_REPLACE, "xoffset for background model data", (NAXIS1 * CCDSUM1 - IMAXIS1) / (2.0 * CCDSUM1));
     211    psMetadataAddF32(config->arguments,PS_LIST_TAIL,"BKG_WARP_YOFFSET", PS_META_REPLACE, "yoffset for background model data", (NAXIS2 * CCDSUM2 - IMAXIS2) / (2.0 * CCDSUM2));
     212
     213    psTrace("pswarp",5,"%d %d %d %d %d %d %g %g %d %d",
     214            psMetadataLookupS32(NULL,header,"IMNAXIS1"),
     215            psMetadataLookupS32(NULL,header,"IMNAXIS2"),
     216            psMetadataLookupS32(NULL,header,"NAXIS1"),
     217            psMetadataLookupS32(NULL,header,"NAXIS2"),
     218            CCDSUM1,CCDSUM2,
     219            psMetadataLookupF32(NULL,config->arguments,"BKG_WARP_XOFFSET"),
     220            psMetadataLookupF32(NULL,config->arguments,"BKG_WARP_YOFFSET"),
     221            psMetadataLookupS32(NULL,config->arguments,"BKG.XGRID"),
     222            psMetadataLookupS32(NULL,config->arguments,"BKG.YGRID"));
     223
     224    *xBin = CCDSUM1;
     225    *yBin = CCDSUM2;
    368226
    369227    return true;
    370228}
     229
     230// **** input == (small) BKGMODEL, astrom == (fullsize) ASTROM
     231// use astrometry from input ref astrometry and BINNING info from input model to modify the astrometry
     232// elements of input chip.  scale is the area of an input pixel is raw image pixels
     233bool pswarpModifyChipAstrom (pmConfig *config, pmFPAview *view, pmChip *chip, pmFPAfile *astrom, bool bilevelAstrometry, double xBin, double yBin) {
     234
     235    // read WCS data from the corresponding header
     236    pmHDU *hdu = pmFPAviewThisHDU (view, astrom->fpa);
     237
     238    // generate a WCS structure from the ASTROM header keywords
     239    pmAstromWCS *WCS = pmAstromWCSfromHeader(hdu->header);
     240
     241    // re-scale the terms of the WCS to match the BKGMODEL scale
     242    double cd1f = xBin;
     243    double cd2f = yBin;
     244
     245    WCS->cdelt1 *= cd1f;
     246    WCS->cdelt2 *= cd2f;
     247    WCS->crpix1 = WCS->crpix1 / cd1f;
     248    WCS->crpix2 = WCS->crpix2 / cd2f;
     249
     250    // WCS->trans->x->nX/nY
     251    // adjust the scale of the WCS coeffs
     252    for (int q = 0; q <= WCS->trans->x->nX; q++) {
     253        for (int r = 0; r <= WCS->trans->x->nY; r++) {
     254            WCS->trans->x->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r);
     255        }
     256    }
     257    for (int q = 0; q <= WCS->trans->y->nX; q++) {
     258        for (int r = 0; r <= WCS->trans->y->nY; r++) {
     259            WCS->trans->y->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r);
     260        }
     261    }
     262
     263    // this is a bit crude : use a locally temporary header to store the WCS information
     264    // which is then applied to the input chip (why not modify the input chip header?)
     265    psMetadata *tmpHeader = psMetadataAlloc();
     266
     267    // write the modified WCS to the header
     268    // XXX this should probably not modify the astrom header
     269    pmAstromWCStoHeader (tmpHeader, WCS);
     270       
     271    // use the modified header to update the WCS elements of the BKGMODEL fpa/chip structures
     272    if (bilevelAstrometry) {
     273        if (!pmAstromReadBilevelChip (chip, tmpHeader)) {
     274            psError(psErrorCodeLast(), false, "Unable to read bilevel chip astrometry for input FPA.");
     275            return false;
     276        }
     277    } else {
     278        // we use a default FPA pixel scale of 1.0
     279        if (!pmAstromReadWCS (chip->parent, chip, tmpHeader, 1.0)) {
     280            psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for input FPA.");
     281            return false;
     282        }
     283    }
     284
     285    return true;
     286}
Note: See TracChangeset for help on using the changeset viewer.