IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 41895 for trunk/psastro


Ignore:
Timestamp:
Nov 4, 2021, 6:08:21 PM (5 years ago)
Author:
eugene
Message:

merge changes from eam_branches/ipp-dev-20210817 (add CERSTD, NASTUSED to header, allow rejection no CERSTD, add psastroTIO)

Location:
trunk/psastro
Files:
14 edited
5 copied

Legend:

Unmodified
Added
Removed

  • trunk/psastro

  • trunk/psastro/src/Makefile.am

    r41705 r41895  
    1010libpsastro_la_LDFLAGS = $(PSASTRO_LIBS) $(PPSTATS_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
    1111
    12 bin_PROGRAMS = psastro psastroExtract psastroModel psastroModelFit gpcModel
     12bin_PROGRAMS = psastro psastroExtract psastroModel psastroModelFit gpcModel psastroTIO
    1313
    1414psastro_CFLAGS = $(PSASTRO_CFLAGS) $(PPSTATS_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
    … …  
    3131gpcModel_LDFLAGS = $(PSASTRO_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
    3232gpcModel_LDADD = libpsastro.la
     33
     34psastroTIO_CFLAGS = $(PSASTRO_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
     35psastroTIO_LDFLAGS = $(PSASTRO_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
     36psastroTIO_LDADD = libpsastro.la
    3337
    3438psastro_SOURCES = \
    … …  
    6266        psastroModelAdjust.c        \
    6367        psastroModelDataSave.c      \
     68        psastroCleanup.c
     69
     70psastroTIO_SOURCES = \
     71        psastroTIO.c                \
     72        psastroTIOArguments.c       \
     73        psastroTIOParseCamera.c   \
     74        psastroTIODataLoad.c      \
     75        psastroTIODataSave.c      \
    6476        psastroCleanup.c
    6577

  • trunk/psastro/src/psastroAstromGuess.c

    r39926 r41895  
    342342        if (!updates) goto skip_chip;
    343343
     344        // in psastroOneChipFit & psastroMosaicOneCihp, astrometry failures are marked with NASTRO = 0
     345        // these should be ignored when checking the overall solution
    344346        int nAstro = psMetadataLookupS32 (&status, updates, "NASTRO");
    345347        if (!nAstro) goto skip_chip;
    346348
     349        // it is not clear when astError = 0.0
    347350        float astError = psMetadataLookupF32 (&status, updates, "CERROR");
    348351        if (fabs(astError) < 1e-6) goto skip_chip;
    … …  
    397400    }
    398401
    399     psPlaneTransform *map = psPlaneTransformAlloc (1, 1);
     402    psPlaneTransform *map = psPlaneTransformAlloc (1, 1, PS_POLYNOMIAL_ORD);
    400403    map->x->coeffMask[1][1] = PS_POLY_MASK_SET;
    401404    map->y->coeffMask[1][1] = PS_POLY_MASK_SET;

  • trunk/psastro/src/psastroConvert.c

    r41435 r41895  
    172172    calStars->n = j;
    173173    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.CALSTARS", PS_DATA_ARRAY, "astrometry objects for masking", calStars);
     174    psFree (calStars);
    174175
    175176    psFree (index);

  • trunk/psastro/src/psastroFixChips.c

    r41536 r41895  
    170170    if (DEBUG) fclose (f);
    171171
    172     psPlaneTransform *map = psPlaneTransformAlloc (1, 1);
     172    psPlaneTransform *map = psPlaneTransformAlloc (1, 1, PS_POLYNOMIAL_ORD);
    173173
    174174    psVector *mask = psVectorAlloc (nPts, PS_TYPE_VECTOR_MASK);
    … …  
    289289        // apply the exiting fromTPA transformation to make the new toFPA consistent with the toTPA layter
    290290        // XXX this only works if toTPA is at most a linear transformation
    291         psPlaneTransform *toFPA = psPlaneTransformAlloc(refChip->toFPA->x->nX, refChip->toFPA->x->nY);
     291        psPlaneTransform *toFPA = psPlaneTransformAlloc(refChip->toFPA->x->nX, refChip->toFPA->x->nY, PS_POLYNOMIAL_ORD);
    292292        for (int i = 0; i <= refChip->toFPA->x->nX; i++) {
    293293            for (int j = 0; j <= refChip->toFPA->x->nY; j++) {

  • trunk/psastro/src/psastroMetadataStats.c

    r41071 r41895  
    3737        psMetadataItemSupplement(&status, stats, header, "ZPT_ERR");
    3838        psMetadataItemSupplement(&status, stats, header, "CERROR");
     39        psMetadataItemSupplement(&status, stats, header, "CERSTD");
    3940        psMetadataItemSupplement(&status, stats, header, "NASTRO");
    4041        psMetadataItemSupplement(&status, stats, header, "AST_R0");

  • trunk/psastro/src/psastroMosaicAstrom.c

    r41536 r41895  
    277277    }
    278278
    279     // the original transforms are (1, 1), but we will need higher order to fit the distortions
     279    // the original transforms are (1, 1), but we will need higher order to fit the distortions.
     280    // forward transformations (chip->fpa->tpa->sky) use Ordinary polynomials
    280281    psFree (fpa->toTPA);
    281 
    282     // the original transforms are (1, 1), but we will need higher order to fit the distortions
    283     fpa->toTPA = psPlaneTransformAlloc (order, order);
     282    fpa->toTPA = psPlaneTransformAlloc (order, order, PS_POLYNOMIAL_ORD);
    284283    for (int i = 0; i <= fpa->toTPA->x->nX; i++) {
    285284        for (int j = 0; j <= fpa->toTPA->x->nY; j++) {

  • trunk/psastro/src/psastroMosaicChipAstrom.c

    r41500 r41895  
    9191    return true;
    9292}
    93 
    94 // XXX count success and failures
    95 // XXX set bad quality if fraction of successes < threshold
    96 // XXX e.g., 5% success would mean 3 of 60 GPC1 chips are goo.
    97 // XXX or, 10% would be 6 of 60.
    98 // XXX seems like a good indicator of failure.
    99 // XXX also, test CERROR against a minimum value?

  • trunk/psastro/src/psastroMosaicCorrectDistortion.c

    r41536 r41895  
    2222    }
    2323
    24     // store the new coeffs in a new structure
    25     psPlaneTransform *toTPAnew = psPlaneTransformAlloc(fpa->toTPA->x->nX, fpa->toTPA->x->nY);
     24    // Store the new coeffs in a new structure.  Forward transformations (chip->fpa->tpa->sky)
     25    // use Ordinary polynomials
     26    psPlaneTransform *toTPAnew = psPlaneTransformAlloc(fpa->toTPA->x->nX, fpa->toTPA->x->nY, PS_POLYNOMIAL_ORD);
    2627
    2728    // set the new coeffs, or set the mask

  • trunk/psastro/src/psastroMosaicFPtoTP.c

    r26259 r41895  
    8080
    8181    // linear fit without xy cross term
    82     psPlaneTransform *map = psPlaneTransformAlloc (1, 1);
     82    psPlaneTransform *map = psPlaneTransformAlloc (1, 1, PS_POLYNOMIAL_ORD);
    8383    map->x->coeffMask[1][1] = PS_POLY_MASK_SET;
    8484    map->y->coeffMask[1][1] = PS_POLY_MASK_SET;

  • trunk/psastro/src/psastroMosaicGradients.c

    r39926 r41895  
    8484        return false;
    8585    }
     86
     87    // forward transformations (chip->fpa->tpa->sky) use Ordinary polynomials
    8688    psFree (fpa->toTPA);
    87     fpa->toTPA = psPlaneTransformAlloc (order, order);
     89    fpa->toTPA = psPlaneTransformAlloc (order, order, PS_POLYNOMIAL_ORD);
    8890    for (int i = 0; i <= fpa->toTPA->x->nX; i++) {
    8991        for (int j = 0; j <= fpa->toTPA->x->nY; j++) {

  • trunk/psastro/src/psastroMosaicOneChip.c

    r39926 r41895  
    2323    bool status;
    2424    char errorWord[64];
     25    char stdevWord[64];
    2526    char orderWord[64];
    2627
    … …  
    4445
    4546    // allowed limits for valid solutions
     47    // CERROR is currently tested during the iterations, but not CERSTD
    4648    snprintf (errorWord, 64, "PSASTRO.MOSAIC.MAX.ERROR.N%d", iteration);
    47     REQUIRED_RECIPE_VALUE (float maxError, errorWord, F32, "failed to find single-chip max allowed error\n");
    48     REQUIRED_RECIPE_VALUE (int minNstar, "PSASTRO.MOSAIC.MIN.NSTAR", S32, "failed to find single-chip min allowed stars\n");
     49    snprintf (stdevWord, 64, "PSASTRO.MOSAIC.MAX.STDEV.N%d", iteration);
     50    REQUIRED_RECIPE_VALUE (float maxError,                  errorWord, F32, "failed to find single-chip max allowed error\n");
     51    REQUIRED_RECIPE_VALUE (float maxStdev,                  stdevWord, F32, "failed to find single-chip max allowed stdev\n");
     52    REQUIRED_RECIPE_VALUE (int   minNstar, "PSASTRO.MOSAIC.MIN.NSTAR", S32, "failed to find single-chip min allowed stars\n");
    4953
    5054    // set the order of the per-chip fit (higher order only if iteration > 0)
    … …  
    8892        }
    8993    } else {
     94        // Forward transformations (chip->fpa->tpa->sky) use Ordinary polynomials
    9095        psFree (chip->toFPA);
    91         chip->toFPA = psPlaneTransformAlloc (order, order);
     96        chip->toFPA = psPlaneTransformAlloc (order, order, PS_POLYNOMIAL_ORD);
    9297        for (int i = 0; i <= chip->toFPA->x->nX; i++) {
    9398            for (int j = 0; j <= chip->toFPA->x->nY; j++) {
    … …  
    134139    int astNstar = results->yStats->clippedNvalues;
    135140
     141    // astStdev is the average 1D stdev of median residuals in arcsec ('results' are in FPA units = microns)
     142    // the median residuals are calculated in a grid of N x N bins
     143    float astStdev = 0.5*(results->dXstdev + results->dYstdev) * plateScale;
     144
    136145    // if we clip away too many stars, the order may be invalid
    137146    if (order == 3) { minNstar = PS_MAX (15, minNstar); }
    … …  
    141150    bool validSolution = true;
    142151
    143     // XXX should these result in errors or be handled another way?
    144     psLogMsg ("psastro", PS_LOG_INFO, "astrometry solution: error: %f arcsec, Nstars: %d", astError, astNstar);
     152    // We have options to exclude chips on the basis of NASTRO, CERROR, CERSTD
     153    psLogMsg ("psastro", PS_LOG_INFO, "astrometry solution: error: %f arcsec, Nstars: %d, stdev: %f arcsec", astError, astNstar, astStdev);
    145154    if ((maxError > 0) && (astError > maxError)) {
    146155        psLogMsg("psastro", PS_LOG_INFO, "residual error is too large, failed to find a solution: %f > %f", astError, maxError);
    147156        validSolution = false;
    148157    }
     158    if ((maxStdev > 0) && (astStdev > maxStdev)) {
     159        psLogMsg("psastro", PS_LOG_INFO, "residual stdev is too large, failed to find a solution: %f > %f", astStdev, maxStdev);
     160        validSolution = false;
     161    }
    149162    if (astNstar < minNstar) {
    150163        psLogMsg("psastro", PS_LOG_INFO, "solution uses too few stars: %d < %d", astNstar, minNstar);
    … …  
    152165    }
    153166
    154     // DVO expects NASTRO = 0 if we fail to find a solution
     167    // DVO expects NASTRO = 0 if we fail to find a solution, NASTUSED is the true number
    155168    psMetadataAddF32 (updates, PS_LIST_TAIL, "PERROR",   PS_META_REPLACE, "astrometry error (pixels)", pixError);
    156169    psMetadataAddF32 (updates, PS_LIST_TAIL, "CERROR",   PS_META_REPLACE, "astrometry error (arcsec)", astError);
     170    psMetadataAddF32 (updates, PS_LIST_TAIL, "CERSTD",   PS_META_REPLACE, "astrometry stdev (arcsec)", astStdev);
    157171    if (validSolution) {
    158172        psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", astError/sqrt(astNstar));
    159173        psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "number of astrometry stars", astNstar);
     174        psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTUSED", PS_META_REPLACE, "number of astrometry stars", astNstar);
    160175    } else {
    161176        psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", 0.0);
    162177        psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "number of astrometry stars", 0);
     178        psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTUSED", PS_META_REPLACE, "number of astrometry stars", astNstar);
    163179    }
    164180    psMetadataAddF32 (updates, PS_LIST_TAIL, "EQUINOX",  PS_META_REPLACE, "", 2000.0); // XXX this is bogus: should be defined based on equinox of refstars

  • trunk/psastro/src/psastroOneChipFit.c

    r41285 r41895  
    4040    // allowed limits for valid solutions
    4141    REQUIRED_RECIPE_VALUE (float maxError, "PSASTRO.MAX.ERROR", F32);
    42     REQUIRED_RECIPE_VALUE (int minNstar, "PSASTRO.MIN.NSTAR", S32);
     42    REQUIRED_RECIPE_VALUE (float maxStdev, "PSASTRO.MAX.STDEV", F32);
     43    REQUIRED_RECIPE_VALUE (int   minNstar, "PSASTRO.MIN.NSTAR", S32);
    4344
    4445    psArray *match = NULL;
    … …  
    121122        }
    122123
    123         // create output toFPA; set masks appropriate to the Elixir DVO astrometry format
     124        // Create output toFPA; set masks appropriate to the Elixir DVO astrometry format.
     125        // Forward transformations (chip->fpa->tpa->sky) use Ordinary polynomials
    124126        psFree (chip->toFPA);
    125         chip->toFPA = psPlaneTransformAlloc (order, order);
     127        chip->toFPA = psPlaneTransformAlloc (order, order, PS_POLYNOMIAL_ORD);
    126128        for (int i = 0; i <= chip->toFPA->x->nX; i++) {
    127129            for (int j = 0; j <= chip->toFPA->x->nY; j++) {
    … …  
    192194    float astError = 0.5*(rawXstdev + rawYstdev) * plateScale;
    193195
     196    // astStdev is the average 1D stdev of median residuals in arcsec ('results' are in FPA units = microns)
     197    // the median residuals are calculated in a grid of N x N bins
     198    float astStdev = 0.5*(results->dXstdev + results->dYstdev) * plateScale;
     199
    194200    // x and y are forced to use the same subset of values:
    195201    int astNstar = results->yStats->clippedNvalues;
    … …  
    198204
    199205    // XXX should these result in errors or be handled another way?
    200     psLogMsg ("psastro", PS_LOG_INFO, "astrometry solution: error: %f arcsec, Nstars: %d", astError, astNstar);
     206    psLogMsg ("psastro", PS_LOG_INFO, "astrometry solution: error: %f arcsec, Nstars: %d, stdev: %f arcsec", astError, astNstar, astStdev);
    201207    if (astError > maxError) {
    202208        psLogMsg("psastro", PS_LOG_INFO, "residual error is too large, failed to find a solution: %f > %f", astError, maxError);
    203209        validSolution = false;
    204210    }
     211    if (astStdev > maxStdev) {
     212        psLogMsg("psastro", PS_LOG_INFO, "residual stdev is too large, failed to find a solution: %f > %f", astStdev, maxStdev);
     213        validSolution = false;
     214    }
    205215    if (astNstar < minNstar) {
    206216        psLogMsg("psastro", PS_LOG_INFO, "solution uses too few stars: %d < %d", astNstar, minNstar);
    … …  
    211221    psMetadataAddF32 (updates, PS_LIST_TAIL, "PERROR",   PS_META_REPLACE, "astrometry error (pixels)", pixError);
    212222    psMetadataAddF32 (updates, PS_LIST_TAIL, "CERROR",   PS_META_REPLACE, "astrometry error (arcsec)", astError);
     223    psMetadataAddF32 (updates, PS_LIST_TAIL, "CERSTD",   PS_META_REPLACE, "astrometry stdev (arcsec)", astStdev);
    213224    if (validSolution) {
    214225        psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", astError/sqrt(astNstar));
    215226        psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "number of astrometry stars", astNstar);
     227        psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTUSED", PS_META_REPLACE, "number of astrometry stars", astNstar);
    216228    } else {
    217229        psastroChipFailureHeader (updates);
     230        psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTUSED", PS_META_REPLACE, "number of astrometry stars", astNstar);
    218231    }
    219232    psMetadataAddF32 (updates, PS_LIST_TAIL, "EQUINOX",  PS_META_REPLACE, "equinox of ref catalog", 2000.0); // XXX this is bogus: should be defined based on equinox of refstars

  • trunk/psastro/src/psastroUtils.c

    r41536 r41895  
    215215psPlaneTransform *psPlaneTransformCopy (psPlaneTransform *input) {
    216216
    217     psPlaneTransform *output = psPlaneTransformAlloc (input->x->nX, input->x->nY);
     217    psPlaneTransform *output = psPlaneTransformAlloc (input->x->nX, input->x->nY, input->x->type);
    218218
    219219    for (int i = 0; i < input->x->nX; i++) {
Note: See TracChangeset for help on using the changeset viewer.