IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 35852


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.

Location:
trunk/ppSkycell/src
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppSkycell/src/ppSkycell.h

    r28375 r35852  
    1111    psString imagesName;                // Filename with images
    1212    psString masksName;                 // Filename with masks
     13    psString wcsrefName;                // Filename with WCS references
    1314    psString outRoot;                   // Output root name
    1415    int numInputs;                      // Number of inputs

  • trunk/ppSkycell/src/ppSkycellArguments.c

    r23996 r35852  
    5151    psMetadataAddStr(arguments, PS_LIST_TAIL, "-images", 0, "Filename with input images", NULL);
    5252    psMetadataAddStr(arguments, PS_LIST_TAIL, "-masks", 0, "Filename with input masks", NULL);
     53    psMetadataAddStr(arguments, PS_LIST_TAIL, "-wcsref", 0, "Filename with WCS references", NULL);
    5354    if (argc == 1 || !psArgumentParse(arguments, &argc, argv) || argc != 2) {
    5455        usage(argv[0], arguments, data);
    … …  
    5758    data->imagesName = psMemIncrRefCounter(psMetadataLookupStr(NULL, arguments, "-images"));
    5859    data->masksName = psMemIncrRefCounter(psMetadataLookupStr(NULL, arguments, "-masks"));
     60    data->wcsrefName = psMemIncrRefCounter(psMetadataLookupStr(NULL, arguments, "-wcsref"));
     61   
    5962    data->outRoot = psStringCopy(argv[1]);
    6063

  • trunk/ppSkycell/src/ppSkycellCamera.c

    r34800 r35852  
    7373    }
    7474
     75    psArray *wcsref = NULL;             // Names of WCS reference images
     76    if (data->wcsrefName) {
     77      wcsref = fileList(data->wcsrefName);
     78      if (!wcsref) {
     79        psError(psErrorCodeLast(), false, "No WCSrefs provided.");
     80        psFree(images);
     81        psFree(masks);
     82        return(false);
     83      }
     84      if (wcsref->n != data->numInputs) {
     85        psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Number of images (%ld) and wcsrefs (%ld) do not match",
     86                images->n, wcsref->n);
     87        psFree(images);
     88        psFree(masks);
     89        psFree(wcsref);
     90        return(false);
     91      }
     92    }
     93   
    7594    psMetadataAddStr(data->config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "Output root", data->outRoot);
    7695
    … …  
    93112            }
    94113        }
     114
     115        if (data->wcsrefName) {
     116          fileArguments("WCSREF", wcsref->data[i], "Name of the WCS reference", data->config);
     117          pmFPAfile *wcsref = pmFPAfileDefineFromArgs(&status, data->config, "PPSKYCELL.WCSREF", "WCSREF");
     118          if (!status || !wcsref) {
     119            psError(PS_ERR_IO, false, "Failed to build file from PPSKYCELL.WCSREF");
     120            return(false);
     121          }
     122/*        if (!pmFPAfileBindFromArgs(&status, image, data->config, "PPSKYCELL.WCSREF", "WCSREF") || !status) { */
     123/*          psError(PS_ERR_IO, false, "Failed to build file from PPSKYCELL.WCSREF"); */
     124/*          return false; */
     125/*        } */
     126        }
     127       
    95128    }
    96129

  • trunk/ppSkycell/src/ppSkycellData.c

    r24073 r35852  
    1515    psFree(data->imagesName);
    1616    psFree(data->masksName);
     17    psFree(data->wcsrefName);
    1718    psFree(data->outRoot);
    1819    psFree(data->config);
    … …  
    2829    data->imagesName = NULL;
    2930    data->masksName = NULL;
     31    data->wcsrefName = NULL;
    3032    data->outRoot = NULL;
    3133    data->numInputs = 0;

  • 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.