IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jul 24, 2013, 4:56:48 PM (13 years ago)
Author:
watersc1
Message:

Addition of WCS-reference option to fix the cases where the WCS is missing. Fix bug in WCS calculation/propagation that shows up when corners aren't populated.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppSkycell/src/ppSkycellLoop.c

    r34800 r35852  
    1010
    1111#define BUFFER 16                       // Size of buffer for projections
     12
     13#if (0)
     14// This is all junk that is basically useless, because it calculates things we already know.
    1215
    1316static void regionMinMax(psRegion *base,// Base region; modified
    … …  
    2225}
    2326
     27
    2428static psRegion *skycellRegion(pmAstromWCS *wcs, // World Coordinate System
    2529                               int numCols, int numRows // Size of image
    … …  
    3438    fromCoords->y = (Y); \
    3539    psPlaneTransformApply(toCoords, wcs->trans, fromCoords); \
     40    printf("%f %f %f %f\n",toCoords->x,toCoords->y,fromCoords->x,fromCoords->y); \
    3641    toCoords->x /= wcs->cdelt1; \
    3742    toCoords->y /= wcs->cdelt2; \
    38     toCoords->x += wcs->crpix1; \
    39     toCoords->y += wcs->crpix2; \
     43    toCoords->x -= wcs->crpix1;         \
     44    toCoords->y += wcs->crpix2;         \
    4045    region->x0 = PS_MAX(region->x0, toCoords->x); \
    4146    region->x1 = PS_MIN(region->x1, toCoords->x); \
    4247    region->y0 = PS_MIN(region->y0, toCoords->y); \
    4348    region->y1 = PS_MAX(region->y1, toCoords->y);
    44 
     49    //
    4550    REGION_RANGE(0, 0);                 // Lower left
    4651    psTrace("ppSkycell", 9, "Start: %.0f %.0f", toCoords->x, toCoords->y);
    … …  
    4954    REGION_RANGE(numCols, numRows);     // Upper right
    5055    psTrace("ppSkycell", 9, "Stop: %.0f %.0f\n", toCoords->x, toCoords->y);
    51 
     56    psTrace("ppSkycell", 7, "%g %g %g %g\n",wcs->cdelt1,wcs->cdelt2,wcs->crpix1,wcs->crpix2);
    5257    return region;
    5358}
     59#endif
    5460
    5561// Activate/deactivate a single element for a list
    … …  
    150156        pmAstromWCS *wcs = pmAstromWCSfromHeader(hdu->header); // World Coordinate System
    151157        if (!wcs) {
     158          if (data->wcsrefName) {
     159            pmFPAfile *wcsref = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.WCSREF", i);
     160            wcs = pmAstromWCSfromHeader(wcsref->fpa->hdu->header);
     161          }
     162          if (!wcs) {     
    152163            psError(psErrorCodeLast(), false, "Unable to read WCS for image %d", i);
    153164            return false;
    154         }
    155 
    156         psRegion *region = imageRegions->data[i] = skycellRegion(wcs, numCols, numRows); // Region of image
     165          }
     166        }
     167
     168        psRegion *region = psRegionAlloc(-1.0 * wcs->crpix1,numCols - wcs->crpix1,
     169                                         -1.0 * wcs->crpix2,numRows - wcs->crpix2);
     170        imageRegions->data[i] = region;
    157171        psTrace("ppSkycell", 5, "Image region %d is: [%.0f:%.0f,%.0f:%.0f]\n",
    158172                i, region->x0, region->x1, region->y0, region->y1);
    … …  
    162176            if (wcs->crval1 == crval1->data.F64[j] && wcs->crval2 == crval2->data.F64[j] &&
    163177                wcs->cdelt1 == cdelt1->data.F64[j] && wcs->cdelt2 == cdelt1->data.F64[j]) {
    164                 regionMinMax(projRegions->data[j], region);
     178              //                regionMinMax(projRegions->data[j], region);
     179              psRegion *proj = projRegions->data[j];
     180              proj->x0 = PS_MIN(region->x0,proj->x0);
     181              proj->x1 = PS_MAX(region->x1,proj->x1);
     182              proj->y0 = PS_MIN(region->y0,proj->y0);
     183              proj->y1 = PS_MAX(region->y1,proj->y1);
    165184                target->data.S32[i] = j;
    166185                found = true;
    … …  
    188207    for (int i = 0; i < numProj; i++) {
    189208        psRegion *projRegion = projRegions->data[i]; // Region for skycell projection
     209        //      projRegion->x0 += 5760;
     210        //projRegion->x1 += 5760;
     211        //      projRegion->x0 +=
     212        //      projRegion->x1 +=
    190213        psTrace("ppSkycell", 2, "Projection %d: [%.0f:%.0f,%.0f:%.0f]\n",
    191214                i, projRegion->x0, projRegion->x1, projRegion->y0, projRegion->y1);
    192         int xSize = -projRegion->x1 + projRegion->x0 + 1; // Size of unbinned image
     215        int xSize = projRegion->x1 - projRegion->x0 + 1; // Size of unbinned image
    193216        int ySize = projRegion->y1 - projRegion->y0 + 1; // Size of unbinned image
    194217        // Size of binned image 1
    … …  
    209232            psImageInit(mask2, 0xFF);
    210233        }
     234        // HDU containing the WCS we plan on using
    211235        pmHDU *projhdu = NULL;
     236
     237        // Do we need to modify the WCS?  This flips to zero after we've done it once.
    212238        int modify_wcs1 = 1;
    213239        int modify_wcs2 = 1;
     240
     241        // Because we may have holes, we need to ensure that we set the CRPIX in the binned image correction.
     242        // Find the minimum/maximum, so we know where the zero is.
     243        float maxCRPIX1   = -99e99;
     244        float maxCRPIX2   = -99e99;
    214245        // Loop over inputs to this projection.
    215246        for (int j = 0; j < data->numInputs; j++) {
    … …  
    229260            }
    230261
     262            // Read the HDU/WCS information from the first entry, and use that as the reference.
    231263            pmFPAfile *file = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.IMAGE", j);
     264            pmAstromWCS *wcs = pmAstromWCSfromHeader(file->fpa->hdu->header); // World Coordinate System
     265            if (!wcs) {
     266              if (data->wcsrefName) {
     267                pmFPAfile *wcsref = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.WCSREF", j);
     268                wcs = pmAstromWCSfromHeader(wcsref->fpa->hdu->header);
     269                if (!projhdu) {
     270                  projhdu = wcsref->fpa->hdu;
     271                }
     272              }
     273              if (!wcs) {         
     274                psError(psErrorCodeLast(), false, "Unable to read WCS for image %d", j);
     275                return false;
     276              }
     277            }
    232278            if (!projhdu) {
    233279              projhdu = file->fpa->hdu;
    234               //              psMetadataCopy(projhdu->header,file->fpa->hdu->header);
    235             }
    236 #if 0
    237             else {
    238               if ((psMetadataLookupF32(NULL,file->fpa->hdu->header,"CRPIX1") >=
    239                    psMetadataLookupF32(NULL,projhdu->header,"CRPIX1"))&&
    240                   (psMetadataLookupF32(NULL,file->fpa->hdu->header,"CRPIX2") >=
    241                    psMetadataLookupF32(NULL,projhdu->header,"CRPIX1"))) {
    242                 projhdu = file->fpa->hdu;
    243                 psMetadataCopy(projhdu->header,file->fpa->hdu->header);
    244               }
    245             }
    246 #endif   
    247            
     280            }
     281
     282            // However, we need to check to see if we've found the maximum CRPIX
     283            if (wcs->crpix1 > maxCRPIX1) {
     284              maxCRPIX1 = wcs->crpix1;
     285            }
     286            if (wcs->crpix2 < maxCRPIX2) {
     287              maxCRPIX2 = wcs->crpix2;
     288            }
     289
     290            // Actuall do the binning.
    248291            pmReadout *inRO = pmFPAviewThisReadout(view, file->fpa); // Readout with input
    249292            psFree(view);
    250293
    251             //      data->maskVal = 0xffff;
    252294            pmReadout *bin1RO = pmReadoutAlloc(NULL), *bin2RO = pmReadoutAlloc(NULL); // Binned readouts
    253295            if (!pmReadoutRebin(bin1RO, inRO, data->maskVal, data->bin1, data->bin1)) {
    … …  
    262304            }
    263305
     306#if (0)
    264307            psRegion *imageRegion = imageRegions->data[j]; // Region for image
    265             // Offsets for image on skycell
    266             int xOffset1 = (-imageRegion->x0 + projRegion->x0) / (float)data->bin1;
    267             int yOffset1 = (-imageRegion->y1 + projRegion->y1) / (float)data->bin1;
     308            fprintf(stderr,"%d IM %f %f %f %f PR %f %f %f %f\n",
     309                    j,
     310                    imageRegion->x0,imageRegion->x1,
     311                    imageRegion->y0,imageRegion->y1,
     312                    projRegion->x0,projRegion->x1,
     313                    projRegion->y0,projRegion->y1);
     314#endif
     315            // Offsets for image on this projection cell are just differences in CRPIX positions.
     316            int xOffset1 = (-1 * wcs->crpix1 - projRegion->x0) / (float)(data->bin1);
     317            int yOffset1 = (-1 * wcs->crpix2 - projRegion->y0) / (float)(data->bin1);
     318           
    268319            int xOffset2 = xOffset1 / (float)data->bin2, yOffset2 = yOffset1 / (float)data->bin2;
     320            psTrace("ppSkycell",5,"Offsets: %d %d : %d %d",
     321                    xOffset1,yOffset1,xOffset2,yOffset2);
     322
     323            // Overlay the data onto the appropriate pixels in the final outputs
    269324            // XXX Completely neglecting rotations
    270325            // The skycells are divided up neatly with them all having the same orientation
    … …  
    276331            }
    277332
     333            // Cleanup on input loop.
    278334            psFree(bin1RO);
    279335            psFree(bin2RO);
    … …  
    345401            WCS->cdelt1 *= cd1f;
    346402            WCS->cdelt2 *= cd2f;
    347             WCS->crpix1 = WCS->crpix1 / cd1f;
    348             WCS->crpix2 = WCS->crpix2 / cd2f;
     403            // Fudge the CRPIX incase we have missing corners
     404            if (maxCRPIX1 > WCS->crpix1) {
     405              WCS->crpix1 = maxCRPIX1 / cd1f;
     406            }
     407            else {
     408              WCS->crpix1 = WCS->crpix1 / cd1f;
     409            }
     410            if (maxCRPIX2 > WCS->crpix2) {
     411              WCS->crpix2 = maxCRPIX2 / cd2f;
     412            }
     413            else {
     414              WCS->crpix2 = WCS->crpix2 / cd2f;
     415            }
    349416#if WCS_DEBUG
    350417            fprintf(stderr,">>> %d %d (%g %g) (%g %g) (%g %g) %d %d %d %d\n",data->bin1,data->bin2,
    … …  
    376443            }
    377444            pmAstromWCStoHeader (fits1->fpa->hdu->header,WCS);
     445#if WCS_DEBUG
    378446            WCS = pmAstromWCSfromHeader(fits1->fpa->hdu->header);
    379 #if WCS_DEBUG
    380447            fprintf(stderr,">>> %d %d (%g %g) (%g %g) (%g %g)\n",data->bin1,data->bin2,
    381448                    cd1f,cd2f,WCS->cdelt1,WCS->cdelt2,
    … …  
    448515            WCS->crpix2 = WCS->crpix2 / cd2f;
    449516           
    450             for (int q = 0; q < WCS->trans->x->nX; q++) {
    451               for (int r = 0; r < WCS->trans->x->nY; r++) {
     517            for (int q = 0; q <= WCS->trans->x->nX; q++) {
     518              for (int r = 0; r <= WCS->trans->x->nY; r++) {
    452519                WCS->trans->x->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r);
    453520              }
    454521            }
    455             for (int q = 0; q < WCS->trans->y->nX; q++) {
    456               for (int r = 0; r < WCS->trans->y->nY; r++) {
     522            for (int q = 0; q <= WCS->trans->y->nX; q++) {
     523              for (int r = 0; r <= WCS->trans->y->nY; r++) {
    457524                WCS->trans->y->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r);
    458525              }
Note: See TracChangeset for help on using the changeset viewer.