IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 27, 2021, 8:16:29 AM (5 years ago)
Author:
eugene
Message:

adding stdev of median residual in bins as a more-robust measure of astrometry success

Location:
branches/eam_branches/ipp-dev-20210817
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-dev-20210817/psModules/src/astrom/pmAstrometryObjects.c

    r41813 r41814  
    340340    psVector *yResGood = psVectorAllocEmpty (match->n, PS_TYPE_F32);
    341341
     342    // we measure the stdev of the median residual in NxN bins.
     343    // use only valid (not NAN) measurements
     344    psVector *xPosValid = psVectorAllocEmpty (match->n, PS_TYPE_F32);
     345    psVector *yPosValid = psVectorAllocEmpty (match->n, PS_TYPE_F32);
     346    psVector *xResValid = psVectorAllocEmpty (match->n, PS_TYPE_F32);
     347    psVector *yResValid = psVectorAllocEmpty (match->n, PS_TYPE_F32);
     348
    342349    for (int i = 0; i < match->n; i++) {
    343350        if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue;
    … …  
    345352        pmAstromObj *rawStar = raw->data[pair->raw];
    346353        if (!isfinite(rawStar->dMag)) continue;
     354
     355        bool isValid = true;
     356        isValid = isValid & isfinite (xFit->data.F32[i]);
     357        isValid = isValid & isfinite (yFit->data.F32[i]);
     358        isValid = isValid & isfinite (xRes->data.F32[i]);
     359        isValid = isValid & isfinite (yRes->data.F32[i]);
     360        if (isValid) {
     361          psVectorAppend (xPosValid, xFit->data.F32[i]);
     362          psVectorAppend (yPosValid, yFit->data.F32[i]);
     363          psVectorAppend (xResValid, xRes->data.F32[i]);
     364          psVectorAppend (yResValid, yRes->data.F32[i]);
     365        }
     366
     367        // for the systematic error, use only high S/N stars
    347368        if (rawStar->dMag > 0.02) continue;
    348369
    … …  
    377398      results->dYrange = pmAstromVectorRange (yResGood, 0.1, 0.9, results->yStats->clippedStdev);
    378399
    379       results->dXstdev = pmAstrom2DSystematics (xFit, yFit, xRes);
    380       results->dYstdev = pmAstrom2DSystematics (xFit, yFit, yRes);
    381     }
    382 
    383     psTrace ("psModules.astrom", 3, "dXsys: %f, dXrange: %f\n", results->dXsys, results->dXrange);
    384     psTrace ("psModules.astrom", 3, "dYsys: %f, dYrange: %f\n", results->dYsys, results->dYrange);
     400      results->dXstdev = pmAstrom2DSystematics (xPosValid, yPosValid, xResValid);
     401      results->dYstdev = pmAstrom2DSystematics (xPosValid, yPosValid, yResValid);
     402    }
     403
     404    psTrace ("psModules.astrom", 3, "dXsys: %f, dXrange: %f, dXstdev: %f\n", results->dXsys, results->dXrange, results->dXstdev);
     405    psTrace ("psModules.astrom", 3, "dYsys: %f, dYrange: %f, dYstdev: %f\n", results->dYsys, results->dYrange, results->dYstdev);
    385406
    386407    psFree (xErr);
    … …  
    390411    psFree (xResGood);
    391412    psFree (yResGood);
     413    psFree (xPosValid);
     414    psFree (yPosValid);
     415    psFree (xResValid);
     416    psFree (yResValid);
    392417
    393418    psFree (x);

  • branches/eam_branches/ipp-dev-20210817/psModules/src/astrom/pmAstrometryObjects.h

    r39926 r41814  
    102102    double  dXrange;                    ///< 10% - 90% range X residuals (unmasked, high S/N)
    103103    double  dYrange;                    ///< 10% - 90% range Y residuals (unmasked, high S/N)
     104    double  dXstdev;                    ///< stdev of median residual in NxN bins
     105    double  dYstdev;                    ///< stdev of median residual in NxN bins
    104106}
    105107pmAstromFitResults;

  • branches/eam_branches/ipp-dev-20210817/psastro/src/psastroMosaicOneChip.c

    r39926 r41814  
    134134    int astNstar = results->yStats->clippedNvalues;
    135135
     136    // astStdev is the average 1D stdev of median residuals in arcsec ('results' are in FPA units = microns)
     137    // the median residuals are calculated in a grid of N x N bins
     138    float astStdev = 0.5*(results->xStdev + results->yStdev) * plateScale;
     139
    136140    // if we clip away too many stars, the order may be invalid
    137141    if (order == 3) { minNstar = PS_MAX (15, minNstar); }
    … …  
    142146
    143147    // 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);
     148    // psLogMsg ("psastro", PS_LOG_INFO, "astrometry solution: error: %f arcsec, Nstars: %d", astError, astNstar);
     149    psLogMsg ("psastro", PS_LOG_INFO, "astrometry solution: error: %f arcsec, Nstars: %d, stdev: %f arcsec", astError, astNstar, astStdev);
    145150    if ((maxError > 0) && (astError > maxError)) {
    146151        psLogMsg("psastro", PS_LOG_INFO, "residual error is too large, failed to find a solution: %f > %f", astError, maxError);
    … …  
    155160    psMetadataAddF32 (updates, PS_LIST_TAIL, "PERROR",   PS_META_REPLACE, "astrometry error (pixels)", pixError);
    156161    psMetadataAddF32 (updates, PS_LIST_TAIL, "CERROR",   PS_META_REPLACE, "astrometry error (arcsec)", astError);
     162    psMetadataAddF32 (updates, PS_LIST_TAIL, "CERSTD",   PS_META_REPLACE, "astrometry stdev (arcsec)", astStdev);
    157163    if (validSolution) {
    158164        psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", astError/sqrt(astNstar));

  • branches/eam_branches/ipp-dev-20210817/psastro/src/psastroOneChipFit.c

    r41285 r41814  
    192192    float astError = 0.5*(rawXstdev + rawYstdev) * plateScale;
    193193
     194    // astStdev is the average 1D stdev of median residuals in arcsec ('results' are in FPA units = microns)
     195    // the median residuals are calculated in a grid of N x N bins
     196    float astStdev = 0.5*(results->xStdev + results->yStdev) * plateScale;
     197
    194198    // x and y are forced to use the same subset of values:
    195199    int astNstar = results->yStats->clippedNvalues;
    … …  
    198202
    199203    // 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);
     204    psLogMsg ("psastro", PS_LOG_INFO, "astrometry solution: error: %f arcsec, Nstars: %d, stdev: %f arcsec", astError, astNstar, astStdev);
    201205    if (astError > maxError) {
    202206        psLogMsg("psastro", PS_LOG_INFO, "residual error is too large, failed to find a solution: %f > %f", astError, maxError);
    … …  
    211215    psMetadataAddF32 (updates, PS_LIST_TAIL, "PERROR",   PS_META_REPLACE, "astrometry error (pixels)", pixError);
    212216    psMetadataAddF32 (updates, PS_LIST_TAIL, "CERROR",   PS_META_REPLACE, "astrometry error (arcsec)", astError);
     217    psMetadataAddF32 (updates, PS_LIST_TAIL, "CERSTD",   PS_META_REPLACE, "astrometry stdev (arcsec)", astStdev);
    213218    if (validSolution) {
    214219        psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", astError/sqrt(astNstar));
Note: See TracChangeset for help on using the changeset viewer.