IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5575


Ignore:
Timestamp:
Nov 22, 2005, 10:38:21 AM (21 years ago)
Author:
eugene
Message:

fixed up WCS

Location:
trunk/psastro
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/psastro/doc/mktest.txt

    r5565 r5575  
    2222  within floating point errors.
    2323
     24examples:
     25
     26build a fake dataset
     27# psastro-mktest doc/psastro-mktest.conf test.cmp test.cat test.raw test.ref
     28
     29run psastrom on the fake data:
     30# psastro doc/psastro.conf test.cmp test.dat
     31

  • trunk/psastro/src/psastro-mktest.c

    r5565 r5575  
    8080        chip->toFPA->y->coeff[1][1] = 0;
    8181
    82         psMetadataAdd (header, PS_LIST_TAIL, "PC001001", PS_DATA_F32 | PS_META_REPLACE, "foo", chip->toFPA->x->coeff[1][0]);
    83         psMetadataAdd (header, PS_LIST_TAIL, "PC001002", PS_DATA_F32 | PS_META_REPLACE, "foo", chip->toFPA->x->coeff[0][1]);
    84         psMetadataAdd (header, PS_LIST_TAIL, "PC002001", PS_DATA_F32 | PS_META_REPLACE, "foo", chip->toFPA->y->coeff[1][0]);
    85         psMetadataAdd (header, PS_LIST_TAIL, "PC002002", PS_DATA_F32 | PS_META_REPLACE, "foo", chip->toFPA->y->coeff[0][1]);
    86 
    8782        chip->fromFPA = p_psPlaneTransformLinearInvert (chip->toFPA);
    8883    }
    … …  
    136131
    137132        // write output header WCS info
    138         psMetadataAdd (header, PS_LIST_TAIL, "CRVAL1", PS_DATA_F32    | PS_META_REPLACE, "foo", RA*DEG_RAD);
    139         psMetadataAdd (header, PS_LIST_TAIL, "CRVAL2", PS_DATA_F32    | PS_META_REPLACE, "foo", DEC*DEG_RAD);
    140         psMetadataAdd (header, PS_LIST_TAIL, "CDELT1", PS_DATA_F32    | PS_META_REPLACE, "foo", PS*DEG_RAD);
    141         psMetadataAdd (header, PS_LIST_TAIL, "CDELT2", PS_DATA_F32    | PS_META_REPLACE, "foo", PS*DEG_RAD);
    142         psMetadataAdd (header, PS_LIST_TAIL, "CTYPE1", PS_DATA_STRING | PS_META_REPLACE, "foo", "RA---SIN");
    143         psMetadataAdd (header, PS_LIST_TAIL, "CTYPE2", PS_DATA_STRING | PS_META_REPLACE, "foo", "DEC--SIN");
    144         psMetadataAdd (header, PS_LIST_TAIL, "CRPIX1", PS_DATA_F32    | PS_META_REPLACE, "foo", tmp1.x);
    145         psMetadataAdd (header, PS_LIST_TAIL, "CRPIX2", PS_DATA_F32    | PS_META_REPLACE, "foo", tmp1.y);
     133        psMetadataAdd (header, PS_LIST_TAIL, "CTYPE1",   PS_DATA_STRING | PS_META_REPLACE, "foo", "RA---SIN");
     134        psMetadataAdd (header, PS_LIST_TAIL, "CTYPE2",   PS_DATA_STRING | PS_META_REPLACE, "foo", "DEC--SIN");
     135        psMetadataAdd (header, PS_LIST_TAIL, "CRVAL1",   PS_DATA_F32    | PS_META_REPLACE, "foo", RA*DEG_RAD);
     136        psMetadataAdd (header, PS_LIST_TAIL, "CRVAL2",   PS_DATA_F32    | PS_META_REPLACE, "foo", DEC*DEG_RAD);
     137        psMetadataAdd (header, PS_LIST_TAIL, "CRPIX1",   PS_DATA_F32    | PS_META_REPLACE, "foo", tmp1.x);
     138        psMetadataAdd (header, PS_LIST_TAIL, "CRPIX2",   PS_DATA_F32    | PS_META_REPLACE, "foo", tmp1.y);
     139        psMetadataAdd (header, PS_LIST_TAIL, "CDELT1",   PS_DATA_F32    | PS_META_REPLACE, "foo", PS*DEG_RAD);
     140        psMetadataAdd (header, PS_LIST_TAIL, "CDELT2",   PS_DATA_F32    | PS_META_REPLACE, "foo", PS*DEG_RAD);
     141        psMetadataAdd (header, PS_LIST_TAIL, "PC001001", PS_DATA_F32    | PS_META_REPLACE, "foo", chip->toFPA->x->coeff[1][0]);
     142        psMetadataAdd (header, PS_LIST_TAIL, "PC001002", PS_DATA_F32    | PS_META_REPLACE, "foo", chip->toFPA->x->coeff[0][1]);
     143        psMetadataAdd (header, PS_LIST_TAIL, "PC002001", PS_DATA_F32    | PS_META_REPLACE, "foo", chip->toFPA->y->coeff[1][0]);
     144        psMetadataAdd (header, PS_LIST_TAIL, "PC002002", PS_DATA_F32    | PS_META_REPLACE, "foo", chip->toFPA->y->coeff[0][1]);
    146145    }
    147146

  • trunk/psastro/src/psastro.h

    r5565 r5575  
    1111# define toSky projection
    1212
    13 bool              pmCellInterpretWCS (pmCell *cell, psMetadata *header) ;
     13bool              pmAstromReadWCS (psPlaneTransform **toFPA, psProjection **toSky, psMetadata *header);
     14bool              pmAstromWriteWCS (psPlaneTransform *toFPA, psProjection *toSky, psMetadata *header);
     15
    1416pmFPA            *pmFPACopyAstrom (pmFPA *inFPA);
    1517psMetadata       *psastroArguments (int *argc, char **argv);
    … …  
    2123bool              psastroWriteCMP (pmFPA *fpa, char *filename);
    2224
    23 psMetadata *testArguments (int *argc, char **argv);
    2425
    25 psPlane *psCoordReadoutToCell (psPlane *cellpix, psPlane *readpix, pmReadout *readout);
    26 psPlane *psCoordCellToReadout (psPlane *readpix, psPlane *cellpix, pmReadout *readout);
    27 psPlane* psCoordChipToCell_EAM(psPlane* cellCoord, const psPlane* chipCoord, const pmCell* cell);
     26psPlane          *psCoordReadoutToCell (psPlane *cellpix, psPlane *readpix, pmReadout *readout);
     27psPlane          *psCoordCellToReadout (psPlane *readpix, psPlane *cellpix, pmReadout *readout);
     28psPlane          *psCoordChipToCell_EAM(psPlane* cellCoord, const psPlane* chipCoord, const pmCell* cell);
    2829
    29 bool testWriteCMP (psMetadata *header, char *filename, psArray *sources);
    30 bool testWriteRef (char *filename, psArray *sources);
    31 bool testWriteRaw (char *filename, psArray *sources);
     30bool              testWriteCMP (psMetadata *header, char *filename, psArray *sources);
     31bool              testWriteRef (char *filename, psArray *sources);
     32bool              testWriteRaw (char *filename, psArray *sources);
     33psMetadata       *testArguments (int *argc, char **argv);
    3234
    3335bool              psastroProjectFPA (pmFPA *fpa, char *starlist, bool toSky);
    … …  
    3537bool              psastroProjectRefstars (psArray *stars, pmReadout *readout);
    3638
    37 psPlane* psCoordChipToCell_EAM(
    38     psPlane* out,                      ///< a plane struct to recycle. If NULL, a new struct is created
    39     const psPlane* in,                 ///< the Chip coordinate
    40     const pmCell* cell                 ///< the cell of interest
    41 );
    42 
    4339psPlaneTransform *p_psPlaneTransformLinearInvert_EAM(psPlaneTransform *transform);
    44 pmFPA *pmFPACopyAstrom (pmFPA *inFPA);
     40pmFPA            *pmFPACopyAstrom (pmFPA *inFPA);

  • trunk/psastro/src/psastroBuildFPA.c

    r5560 r5575  
    2020    pmFPA *fpa = pmFPAAlloc (NULL, NULL);
    2121
     22    // identify matrix for fpa to tpa
     23    fpa->toTPA   = psPlaneDistortAlloc (1, 1, 0, 0);
     24    fpa->toTPA->x->coeff[1][0][0][0] = 1;
     25    fpa->toTPA->x->mask[1][1][0][0]  = 1;
     26
     27    fpa->toTPA->y->coeff[0][1][0][0] = 1;
     28    fpa->toTPA->y->mask[1][1][0][0]  = 1;
     29    fpa->fromTangentPlane = psPlaneDistortInvert (fpa->toTangentPlane);
     30
    2231    psArrayRealloc (fpa->chips, Nchips);
    2332    for (int i = 0; i < Nchips; i++) {
    … …  
    3039            // XXX EAM : extend this function to take more than one header
    3140            cell->header = header;     
    32             pmCellInterpretWCS (cell, header);
     41            pmAstromReadWCS (&chip->toFPA, &fpa->toSky, header);
     42            chip->fromFPA = p_psPlaneTransformLinearInvert (chip->toFPA);
     43
     44            // allocate identity matrix for cell to chip
     45            cell->toChip   = psPlaneTransformAlloc (1, 1);
     46            cell->toChip->x->coeff[1][0] = 1;
     47            cell->toChip->x->mask[1][1]  = 1;
     48            cell->toChip->y->coeff[0][1] = 1;
     49            cell->toChip->y->mask[1][1]  = 1;
    3350
    3451            psArrayRealloc (cell->readouts, Nreadouts);
    … …  
    5269}
    5370
    54 // interpret header WCS
    55 bool pmCellInterpretWCS (pmCell *cell, psMetadata *header) {
    56 
    57     pmChip *chip;
    58     pmFPA  *fpa;
     71// interpret header WCS (only handles traditional WCS for the moment)
     72bool pmAstromReadWCS (psPlaneTransform **toFPAOut, psProjection **toSkyOut, psMetadata *header) {
     73
     74    psPlaneTransform *toFPA;
     75    psProjection *toSky;
     76    psProjectionType type;
    5977    bool status;
    60     psProjectionType type;
    6178    float crval1, crval2, crpix1, crpix2, cdelt1, cdelt2;
    6279    float pc1_1, pc1_2, pc2_1, pc2_2;
    … …  
    136153got_matrix:
    137154
    138     // XXX EAM : these must already be set
    139     chip = cell->parent;
    140     fpa  = chip->parent;
    141 
    142     // set toChip to identity as default
    143     // XXX EAM : can we actually use higher order?
    144     int nX = psMetadataLookupS32 (&status, header, "PSASTRO.CHIP.NX");
    145     if (!status) nX = 1;
    146     int nY = psMetadataLookupS32 (&status, header, "PSASTRO.CHIP.NY");
    147     if (!status) nY = 1;
    148 
    149     cell->toChip   = psPlaneTransformAlloc (1, 1);
    150     cell->toChip->x->coeff[1][0] = 1;
    151     cell->toChip->x->mask[1][1]  = 1;
    152 
    153     cell->toChip->y->coeff[0][1] = 1;
    154     cell->toChip->y->mask[1][1]  = 1;
    155 
    156155    // XXX EAM : if fpa->toSky and fpa->toTPA are already defined, then the
    157156    //           toFPA must be modified to match the crval(i), scale(i) and crpix(i)
    … …  
    161160    // XXX EAM : psPlaneTransformAlloc uses nTerm not nOrder (bug 581)
    162161    // XXX EAM : I've fixed this in pslib eam_rel8_b2
    163     chip->toFPA  = psPlaneTransformAlloc (1, 1);
     162    toFPA = psPlaneTransformAlloc (1, 1);
    164163   
    165     chip->toFPA->x->coeff[0][0] = crpix1;
    166     chip->toFPA->x->coeff[1][0] = pc1_1;
    167     chip->toFPA->x->coeff[0][1] = pc1_2;
    168     chip->toFPA->x->mask[1][1]  = 1;
    169 
    170     chip->toFPA->y->coeff[0][0] = crpix2;
    171     chip->toFPA->y->coeff[1][0] = pc2_1;
    172     chip->toFPA->y->coeff[0][1] = pc2_2;
    173     chip->toFPA->y->mask[1][1]  = 1;
    174 
    175     chip->fromFPA = p_psPlaneTransformLinearInvert (chip->toFPA);
    176 
    177     // set toTPA to identity as default
    178     // XXX EAM : psPlaneDistortAlloc uses nTerm not nOrder (bug 581)
    179     if (fpa->toTPA == NULL) {
    180         fpa->toTPA   = psPlaneDistortAlloc (1, 1, 0, 0);
    181         fpa->toTPA->x->coeff[1][0][0][0] = 1;
    182         fpa->toTPA->x->mask[1][1][0][0]  = 1;
    183 
    184         fpa->toTPA->y->coeff[0][1][0][0] = 1;
    185         fpa->toTPA->y->mask[1][1][0][0]  = 1;
    186         fpa->fromTangentPlane = psPlaneDistortInvert (fpa->toTangentPlane);
    187     } else {
    188         psLogMsg ("psastro", 2, "warning: fpa distortion already defined\n");
    189     }
     164    toFPA->x->coeff[0][0] = crpix1;
     165    toFPA->x->coeff[1][0] = pc1_1;
     166    toFPA->x->coeff[0][1] = pc1_2;
     167    toFPA->x->mask[1][1]  = 1;
     168
     169    toFPA->y->coeff[0][0] = crpix2;
     170    toFPA->y->coeff[1][0] = pc2_1;
     171    toFPA->y->coeff[0][1] = pc2_2;
     172    toFPA->y->mask[1][1]  = 1;
    190173
    191174    // center of projection is (0,0) coordinate of TPA
    192     fpa->toSky = psProjectionAlloc (crval1*RAD_DEG, crval2*RAD_DEG, cdelt1*RAD_DEG, cdelt2*RAD_DEG, type);
     175    toSky = psProjectionAlloc (crval1*RAD_DEG, crval2*RAD_DEG, cdelt1*RAD_DEG, cdelt2*RAD_DEG, type);
     176
     177    *toFPAOut = toFPA;
     178    *toSkyOut = toSky;
     179
    193180    return true;
    194181}
     182
     183
     184// convert toFPA / toSky components to traditional WCS
     185bool pmAstromWriteWCS (psPlaneTransform *toFPA, psProjection *toSky, psMetadata *header) {
     186
     187    switch (toSky->type) {
     188      case PS_PROJ_SIN:
     189        psMetadataAdd (header, PS_LIST_TAIL, "CTYPE1", PS_DATA_STRING | PS_META_REPLACE, "", "RA---SIN");
     190        psMetadataAdd (header, PS_LIST_TAIL, "CTYPE2", PS_DATA_STRING | PS_META_REPLACE, "", "DEC--SIN");
     191        break;
     192      case PS_PROJ_TAN:
     193        psMetadataAdd (header, PS_LIST_TAIL, "CTYPE1", PS_DATA_STRING | PS_META_REPLACE, "", "RA---TAN");
     194        psMetadataAdd (header, PS_LIST_TAIL, "CTYPE2", PS_DATA_STRING | PS_META_REPLACE, "", "DEC--TAN");
     195        break;
     196      case PS_PROJ_AIT:
     197        psMetadataAdd (header, PS_LIST_TAIL, "CTYPE1", PS_DATA_STRING | PS_META_REPLACE, "", "RA---AIT");
     198        psMetadataAdd (header, PS_LIST_TAIL, "CTYPE2", PS_DATA_STRING | PS_META_REPLACE, "", "DEC--AIT");
     199        break;
     200      case PS_PROJ_PAR:
     201        psMetadataAdd (header, PS_LIST_TAIL, "CTYPE1", PS_DATA_STRING | PS_META_REPLACE, "", "RA---PAR");
     202        psMetadataAdd (header, PS_LIST_TAIL, "CTYPE2", PS_DATA_STRING | PS_META_REPLACE, "", "DEC--PAR");
     203        break;
     204      default:
     205        psLogMsg ("psastro", 2, "warning: unknown projection type %s\n", toSky->type);
     206        return false;
     207    }
     208
     209    psMetadataAdd (header, PS_LIST_TAIL, "CRVAL1", PS_DATA_F32 | PS_META_REPLACE, "", toSky->R*DEG_RAD);
     210    psMetadataAdd (header, PS_LIST_TAIL, "CRVAL2", PS_DATA_F32 | PS_META_REPLACE, "", toSky->D*DEG_RAD);
     211    psMetadataAdd (header, PS_LIST_TAIL, "CRPIX1", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->x->coeff[0][0]);
     212    psMetadataAdd (header, PS_LIST_TAIL, "CRPIX2", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->y->coeff[0][0]);
     213    psMetadataAdd (header, PS_LIST_TAIL, "CDELT1", PS_DATA_F32 | PS_META_REPLACE, "", toSky->Xs*DEG_RAD);
     214    psMetadataAdd (header, PS_LIST_TAIL, "CDELT2", PS_DATA_F32 | PS_META_REPLACE, "", toSky->Ys*DEG_RAD);
     215
     216    psMetadataAdd (header, PS_LIST_TAIL, "PC001001", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->x->coeff[1][0]);
     217    psMetadataAdd (header, PS_LIST_TAIL, "PC001002", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->x->coeff[0][1]);
     218    psMetadataAdd (header, PS_LIST_TAIL, "PC002001", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->y->coeff[1][0]);
     219    psMetadataAdd (header, PS_LIST_TAIL, "PC002002", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->y->coeff[0][1]);
     220
     221    // alternative representations use
     222    // CD1_1 = PC001001*CDELT1, etc
     223    // make these representations optional
     224
     225    return true;
     226}

  • trunk/psastro/src/psastroUtils.c

    r5565 r5575  
    102102                match = pmAstromRadiusMatch (rawstars, refstars, config);
    103103
    104                 // fit astrometric terms
     104                // improved fit for astrometric terms
    105105                pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, config);
     106                chip->fromFPA = p_psPlaneTransformLinearInvert (chip->toFPA);
    106107            }
     108            pmAstromWriteWCS (chip->toFPA, fpa->toSky, cell->header);
    107109        }
    108110    }
Note: See TracChangeset for help on using the changeset viewer.