IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 18, 2012, 6:10:35 PM (14 years ago)
Author:
watersc1
Message:

Update to ppSkycell to add a proper header and a proper WCS. Finally.

Fix WCS bug in pswarpLoop that was revealed when working on ppSkycell. Finally.

Changed psImagePixelManip to only copy finite pixel values. This is likely "slow" relative to the old version, and it looks like I'm getting some random default values. In any case, this fixes the NAN striping that was present in the ppSkycell outputs.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/czw_branch/20120906/ppSkycell/src/ppSkycellLoop.c

    r28375 r34449  
    1515                         )
    1616{
    17     base->x0 = PS_MIN(base->x0, compare->x0);
    18     base->x1 = PS_MAX(base->x1, compare->x1);
     17    base->x0 = PS_MAX(base->x0, compare->x0);
     18    base->x1 = PS_MIN(base->x1, compare->x1);
    1919    base->y0 = PS_MIN(base->y0, compare->y0);
    2020    base->y1 = PS_MAX(base->y1, compare->y1);
    … …  
    2727{
    2828    psPlane *fromCoords = psPlaneAlloc(), *toCoords = psPlaneAlloc(); // Coordinates for transforms
    29     psRegion *region = psRegionAlloc(INFINITY, -INFINITY, INFINITY, -INFINITY); // Region for skycell
     29    psRegion *region = psRegionAlloc(-INFINITY, INFINITY, INFINITY, -INFINITY); // Region for skycell
    3030
    3131// Limit the region using the nominated point (X,Y)
    … …  
    3838    toCoords->x += wcs->crpix1; \
    3939    toCoords->y += wcs->crpix2; \
    40     region->x0 = PS_MIN(region->x0, toCoords->x); \
    41     region->x1 = PS_MAX(region->x1, toCoords->x); \
     40    region->x0 = PS_MAX(region->x0, toCoords->x); \
     41    region->x1 = PS_MIN(region->x1, toCoords->x); \
    4242    region->y0 = PS_MIN(region->y0, toCoords->y); \
    4343    region->y1 = PS_MAX(region->y1, toCoords->y);
    … …  
    132132    psVector *target = psVectorAlloc(data->numInputs, PS_TYPE_S32); // Target for each input
    133133    psArray *imageRegions = psArrayAlloc(data->numInputs); // Region for image
    134 
     134    psArray *regionHDUs = psArrayAlloc(data->numInputs);
    135135    for (int i = 0; i < data->numInputs; i++) {
    136136        pmFPAfile *file = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.IMAGE", i); // File to examine
    … …  
    178178            target->data.S32[i] = numProj;
    179179            psTrace("ppSkycell", 3, "Image %d uses new projection\n", i);
    180 
     180            psArrayAdd(regionHDUs,1,hdu);
    181181            numProj++;
    182182        }
    … …  
    189189        psTrace("ppSkycell", 2, "Projection %d: [%.0f:%.0f,%.0f:%.0f]\n",
    190190                i, projRegion->x0, projRegion->x1, projRegion->y0, projRegion->y1);
    191         int xSize = projRegion->x1 - projRegion->x0 + 1; // Size of unbinned image
     191        int xSize = -projRegion->x1 + projRegion->x0 + 1; // Size of unbinned image
    192192        int ySize = projRegion->y1 - projRegion->y0 + 1; // Size of unbinned image
    193193        // Size of binned image 1
    … …  
    208208            psImageInit(mask2, 0xFF);
    209209        }
    210 
     210        pmHDU *projhdu = NULL;
     211        int modify_wcs1 = 1;
     212        int modify_wcs2 = 1;
    211213        for (int j = 0; j < data->numInputs; j++) {
    212214            if (target->data.S32[j] != i) {
    … …  
    227229
    228230            pmFPAfile *file = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.IMAGE", j);
     231            if (!projhdu) {
     232              projhdu = file->fpa->hdu;
     233            }
    229234            pmReadout *inRO = pmFPAviewThisReadout(view, file->fpa); // Readout with input
    230235            psFree(view);
    231236
    232237            // Flip images; no idea why this has to be done, but apparently it does
    233             {
    234                 psImage *rot = psImageRotate(NULL, inRO->image, M_PI, NAN, PS_INTERPOLATE_BILINEAR);
    235                 psFree(inRO->image);
    236                 inRO->image = rot;
    237             }
    238             if (inRO->mask) {
    239                 psImage *rot = psImageRotate(NULL, inRO->mask, M_PI, 0, PS_INTERPOLATE_FLAT);
    240                 psFree(inRO->mask);
    241                 inRO->mask = rot;
    242             }
     238/*             { */
     239/*                 psImage *rot = psImageRotate(NULL, inRO->image, M_PI, NAN, PS_INTERPOLATE_BILINEAR); */
     240/*                 psFree(inRO->image); */
     241/*                 inRO->image = rot; */
     242/*             } */
     243/*             if (inRO->mask) { */
     244/*                 psImage *rot = psImageRotate(NULL, inRO->mask, M_PI, 0, PS_INTERPOLATE_FLAT); */
     245/*                 psFree(inRO->mask); */
     246/*                 inRO->mask = rot; */
     247/*             } */
    243248
    244249            pmReadout *bin1RO = pmReadoutAlloc(NULL), *bin2RO = pmReadoutAlloc(NULL); // Binned readouts
    … …  
    256261            psRegion *imageRegion = imageRegions->data[j]; // Region for image
    257262            // Offsets for image on skycell
    258             int xOffset1 = (imageRegion->x0 - projRegion->x0) / (float)data->bin1;
     263            int xOffset1 = (-imageRegion->x0 + projRegion->x0) / (float)data->bin1;
    259264            int yOffset1 = (imageRegion->y0 - projRegion->y0) / (float)data->bin1;
    260265            int xOffset2 = xOffset1 / (float)data->bin2, yOffset2 = yOffset1 / (float)data->bin2;
    … …  
    313318
    314319        if (data->doFits) {
     320
     321          pmFPAfile *fits1 = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.BIN1", 0);
     322          // Copy header from projection root hdu
     323          fits1->fpa->hdu = pmHDUAlloc(NULL);
     324          fits1->fpa->hdu->header = psMetadataAlloc();
     325          psMetadataCopy(fits1->fpa->hdu->header,projhdu->header);
     326
     327          // Change wcs here
     328          if (modify_wcs1) {
     329            pmAstromWCS *WCS = pmAstromWCSfromHeader(fits1->fpa->hdu->header);
     330            double cd1f = 1.0 * data->bin1;
     331            double cd2f = 1.0 * data->bin1;
     332            fprintf(stderr,">>> %d %d (%g %g) (%g %g) (%g %g)\n",data->bin1,data->bin2,
     333                    cd1f,cd2f,WCS->cdelt1,WCS->cdelt2,
     334                    WCS->crpix1,WCS->crpix2
     335                    );
     336
     337            WCS->cdelt1 *= cd1f;
     338            WCS->cdelt2 *= cd2f;
     339            WCS->crpix1 = WCS->crpix1 / cd1f;
     340            WCS->crpix2 = WCS->crpix2 / cd2f;
     341            fprintf(stderr,">>> %d %d (%g %g) (%g %g) (%g %g) %d %d %d %d\n",data->bin1,data->bin2,
     342                    cd1f,cd2f,WCS->cdelt1,WCS->cdelt2,
     343                    WCS->crpix1,WCS->crpix2,
     344                    WCS->trans->x->nX,              WCS->trans->x->nY,
     345                    WCS->trans->y->nX,              WCS->trans->y->nY
     346                    );
     347           
     348            for (int q = 0; q <= WCS->trans->x->nX; q++) {
     349              for (int r = 0; r <= WCS->trans->x->nY; r++) {
     350                WCS->trans->x->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r);
     351                fprintf(stderr,"PC: %d %d %g %g\n",
     352                        q,r,pow(cd1f,q) * pow(cd2f,r),
     353                        WCS->trans->x->coeff[q][r]);
     354              }
     355            }
     356            for (int q = 0; q <= WCS->trans->y->nX; q++) {
     357              for (int r = 0; r <= WCS->trans->y->nY; r++) {
     358                WCS->trans->y->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r);
     359                fprintf(stderr,"PC: %d %d %g %g\n",
     360                         q,r,pow(cd1f,q) * pow(cd2f,r),
     361                         WCS->trans->y->coeff[q][r]);
     362
     363              }
     364            }
     365            pmAstromWCStoHeader (fits1->fpa->hdu->header,WCS);
     366            WCS = pmAstromWCSfromHeader(fits1->fpa->hdu->header);
     367            fprintf(stderr,">>> %d %d (%g %g) (%g %g) (%g %g)\n",data->bin1,data->bin2,
     368                    cd1f,cd2f,WCS->cdelt1,WCS->cdelt2,
     369                    WCS->crpix1,WCS->crpix2
     370                    );
     371           
     372            modify_wcs1 = 0;
     373          }
     374
     375         
     376          pmChip *Fchip1 = pmFPAfileThisChip(data->config->files, view, "PPSKYCELL.BIN1");
     377          psMetadataAddS32(Fchip1->concepts,PS_LIST_TAIL,"CHIP.XPARITY", PS_META_REPLACE,"",1);
     378          psMetadataAddS32(Fchip1->concepts,PS_LIST_TAIL,"CHIP.YPARITY", PS_META_REPLACE,"",1);
     379
    315380          pmCell *Fcell1 = pmFPAfileThisCell(data->config->files, view, "PPSKYCELL.BIN1"); // Rebinned cell 1
    316           pmCell *Fcell2 = pmFPAfileThisCell(data->config->files, view, "PPSKYCELL.BIN2"); // Rebinned cell 2
    317 
    318           // This is a hack to get a functioning header created so the fits images can be written out.
     381/*        // This is a hack to get a functioning header created so the fits images can be written out. */
    319382          psMetadataAddS32(Fcell1->concepts,PS_LIST_TAIL,"CELL.XPARITY", PS_META_REPLACE,"",1);
    320383          psMetadataAddS32(Fcell1->concepts,PS_LIST_TAIL,"CELL.YPARITY", PS_META_REPLACE,"",1);
    321384          psMetadataAddS32(Fcell1->concepts,PS_LIST_TAIL,"CELL.READDIR", PS_META_REPLACE,"",1);
     385
     386          pmReadout *Fro1 = pmReadoutAlloc(Fcell1);
     387          Fro1->image = image1;
     388          Fro1->mask = mask1;
     389          Fro1->data_exists = Fcell1->data_exists = Fcell1->parent->data_exists = true;
     390         
     391          fits1->save = true;
     392          fits1->fileIndex = i;
     393
     394          // Repeat with second binned image
     395          pmFPAfile *fits2 = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.BIN2", 0);
     396          fits2->fpa->hdu = pmHDUAlloc(NULL);
     397          fits2->fpa->hdu->header = psMetadataAlloc();
     398          psMetadataCopy(fits2->fpa->hdu->header,projhdu->header);
     399
     400          if (modify_wcs2) {
     401            pmAstromWCS *WCS = pmAstromWCSfromHeader(fits2->fpa->hdu->header);
     402            double cd1f = 1.0 * data->bin2 * data->bin1;
     403            double cd2f = 1.0 * data->bin2 * data->bin1;
     404           
     405            WCS->cdelt1 *= cd1f;
     406            WCS->cdelt2 *= cd2f;
     407            WCS->crpix1 = WCS->crpix1 / cd1f;
     408            WCS->crpix2 = WCS->crpix2 / cd2f;
     409           
     410            for (int q = 0; q < WCS->trans->x->nX; q++) {
     411              for (int r = 0; r < WCS->trans->x->nY; r++) {
     412                WCS->trans->x->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r);
     413              }
     414            }
     415            for (int q = 0; q < WCS->trans->y->nX; q++) {
     416              for (int r = 0; r < WCS->trans->y->nY; r++) {
     417                WCS->trans->y->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r);
     418              }
     419            }
     420            pmAstromWCStoHeader (fits2->fpa->hdu->header,WCS);
     421            modify_wcs2 = 0;
     422          }
     423
     424
     425         
     426          pmChip *Fchip2 = pmFPAfileThisChip(data->config->files, view, "PPSKYCELL.BIN2");
     427          psMetadataAddS32(Fchip2->concepts,PS_LIST_TAIL,"CHIP.XPARITY", PS_META_REPLACE,"",1);
     428          psMetadataAddS32(Fchip2->concepts,PS_LIST_TAIL,"CHIP.YPARITY", PS_META_REPLACE,"",1);
     429
     430         
     431          pmCell *Fcell2 = pmFPAfileThisCell(data->config->files, view, "PPSKYCELL.BIN2"); // Rebinned cell 2
    322432          psMetadataAddS32(Fcell2->concepts,PS_LIST_TAIL,"CELL.XPARITY", PS_META_REPLACE,"",1);
    323433          psMetadataAddS32(Fcell2->concepts,PS_LIST_TAIL,"CELL.YPARITY", PS_META_REPLACE,"",1);
    324434          psMetadataAddS32(Fcell2->concepts,PS_LIST_TAIL,"CELL.READDIR", PS_META_REPLACE,"",1);
    325 
    326           psMetadataAddS32(Fcell1->parent->concepts,PS_LIST_TAIL,"CHIP.XPARITY", PS_META_REPLACE,"",1);
    327           psMetadataAddS32(Fcell1->parent->concepts,PS_LIST_TAIL,"CHIP.YPARITY", PS_META_REPLACE,"",1);
    328           psMetadataAddS32(Fcell2->parent->concepts,PS_LIST_TAIL,"CHIP.XPARITY", PS_META_REPLACE,"",1);
    329           psMetadataAddS32(Fcell2->parent->concepts,PS_LIST_TAIL,"CHIP.YPARITY", PS_META_REPLACE,"",1);
    330 
    331           pmReadout *Fro1 = pmReadoutAlloc(Fcell1), *Fro2 = pmReadoutAlloc(Fcell2); // Binned readouts
    332435         
    333           Fro1->image = image1;
     436          pmReadout *Fro2 = pmReadoutAlloc(Fcell2);
    334437          Fro2->image = image2;
    335           Fro1->mask = mask1;
    336438          Fro2->mask = mask2;
    337          
    338           Fro1->data_exists = Fcell1->data_exists = Fcell1->parent->data_exists = true;
    339439          Fro2->data_exists = Fcell2->data_exists = Fcell2->parent->data_exists = true;
    340          
    341           pmFPAfile *fits1 = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.BIN1", 0);
    342           fits1->save = true;
    343           fits1->fileIndex = i;
    344           pmFPAfile *fits2 = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.BIN2", 0);
     440
    345441          fits2->save = true;
    346442          fits2->fileIndex = i;
Note: See TracChangeset for help on using the changeset viewer.