IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 26260


Ignore:
Timestamp:
Nov 22, 2009, 3:00:47 PM (16 years ago)
Author:
eugene
Message:

various updates from eam_branches/20091113

1) psf model order limits now consistent in poly and map modes
2) added subpix visualization
3) added model fit in MATCHED_REFS
4) handle input smf files with empty extensions (earlier failure)
5) function for generating unique reference matches

Location:
trunk/psModules/src
Files:
17 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/astrom/pmAstrometryModel.c

    r19595 r26260  
    725725    double X = Xo + RX*cos(POS - To)*cos(Po) + RY*sin(POS - To)*sin(Po);
    726726    double Y = Yo + RY*sin(POS - To)*cos(Po) - RX*cos(POS - To)*sin(Po);
    727     psLogMsg ("psModules.astrom", 4, "Boresite coords on reference chip: %f, %f\n", X, Y);
     727    psLogMsg ("psModules.astrom", 4, "Boresite coords on reference chip: %f, %f pix = %f, %f sky\n", X, Y, PM_DEG_RAD*RA, PM_DEG_RAD*DEC);
    728728
    729729    // get reference chip from name

  • trunk/psModules/src/astrom/pmAstrometryObjects.c

    r24034 r26260  
    3737#include "pmKapaPlots.h"
    3838#include "pmAstrometryVisual.h"
     39
     40// XXX this is defined in pmPSFtry.h, which makes no sense
     41float psVectorSystematicError (psVector *residuals, psVector *errors, float clipFraction);
    3942
    4043#define PM_ASTROMETRYOBJECTS_DEBUG 1
    … …  
    283286            return results;
    284287        }
    285         psTrace ("psModules.astrom", 3, "x resid: %f +/- %f (%ld of %ld)\n", results->xStats->clippedMean, results->xStats->clippedStdev, results->xStats->clippedNvalues, x->n);
     288        // psTrace ("psModules.astrom", 3, "x resid: %f +/- %f (%ld of %ld)\n", results->xStats->clippedMean, results->xStats->clippedStdev, results->xStats->clippedNvalues, x->n);
     289        psTrace ("psModules.astrom", 3, "x resid: %f +/- %f (%ld of %ld)\n", results->xStats->robustMedian, results->xStats->robustStdev, results->xStats->clippedNvalues, x->n);
    286290
    287291        if (!psVectorClipFitPolynomial2D (map->y, results->yStats, mask, 0xff, y, wt, X, Y)) {
    … …  
    296300            return results;
    297301        }
    298         psTrace ("psModules.astrom", 3, "y resid: %f +/- %f (%ld of %ld)\n", results->yStats->clippedMean, results->yStats->clippedStdev, results->yStats->clippedNvalues, y->n);
     302        // psTrace ("psModules.astrom", 3, "y resid: %f +/- %f (%ld of %ld)\n", results->yStats->clippedMean, results->yStats->clippedStdev, results->yStats->clippedNvalues, y->n);
     303        psTrace ("psModules.astrom", 3, "y resid: %f +/- %f (%ld of %ld)\n", results->yStats->robustMedian, results->yStats->robustStdev, results->yStats->clippedNvalues, y->n);
    299304    }
    300305    results->xStats->clipIter = stats->clipIter;
    301306    results->yStats->clipIter = stats->clipIter;
     307
     308    // *** calculate the 90%-ile and the systematic scatter for each direction.
     309
     310    // generate the X residual vector
     311    psVector *xFit = psPolynomial2DEvalVector (map->x, X, Y);
     312    if (!xFit) abort();
     313    psVector *xRes = (psVector *) psBinaryOp (NULL, x, "-", xFit);
     314    if (!xRes) abort();
     315    psFree (xFit);
     316   
     317    psVector *yFit = psPolynomial2DEvalVector (map->y, X, Y);
     318    if (!yFit) abort();
     319    psVector *yRes = (psVector *) psBinaryOp (NULL, y, "-", yFit);
     320    if (!yRes) abort();
     321    psFree (yFit);
     322
     323    // extract a high-quality subset (unmasked, S/N > XXX) and position errors
     324    // XXX for now, generate a position error based on the magnitude error
     325    psVector *xErr     = psVectorAllocEmpty (match->n, PS_TYPE_F32);
     326    psVector *yErr     = psVectorAllocEmpty (match->n, PS_TYPE_F32);
     327    psVector *xResGood = psVectorAllocEmpty (match->n, PS_TYPE_F32);
     328    psVector *yResGood = psVectorAllocEmpty (match->n, PS_TYPE_F32);
     329
     330    for (int i = 0; i < match->n; i++) {
     331        if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue;
     332        pmAstromMatch *pair = match->data[i];
     333        pmAstromObj *rawStar = raw->data[pair->raw];
     334        if (!isfinite(rawStar->dMag)) continue;
     335        if (rawStar->dMag > 0.02) continue;
     336
     337        // two likely failure values: NAN or 0.0 --> use dMag in this case
     338        float xErrValue, yErrValue;
     339        if (isfinite(rawStar->chip->xErr) && (rawStar->chip->xErr > 0.0)) {
     340            xErrValue = rawStar->chip->xErr;
     341        } else {
     342            xErrValue = PS_MAX(0.005, rawStar->dMag);
     343        }
     344        if (isfinite(rawStar->chip->yErr) && (rawStar->chip->yErr > 0.0)) {
     345            yErrValue = rawStar->chip->yErr;
     346        } else {
     347            yErrValue = PS_MAX(0.005, rawStar->dMag);
     348        }
     349
     350        psVectorAppend (xErr,     xErrValue);
     351        psVectorAppend (yErr,     yErrValue);
     352        psVectorAppend (xResGood, xRes->data.F32[i]);
     353        psVectorAppend (yResGood, yRes->data.F32[i]);
     354    }
     355
     356    results->dXsys = psVectorSystematicError (xResGood, xErr, 0.05);
     357    results->dYsys = psVectorSystematicError (yResGood, yErr, 0.05);
     358
     359    results->dXrange = pmAstromVectorRange (xResGood, 0.1, 0.9, results->xStats->clippedStdev);
     360    results->dYrange = pmAstromVectorRange (yResGood, 0.1, 0.9, results->yStats->clippedStdev);
     361
     362    psTrace ("psModules.astrom", 3, "dXsys: %f, dXrange: %f\n", results->dXsys, results->dXrange);
     363    psTrace ("psModules.astrom", 3, "dYsys: %f, dYrange: %f\n", results->dYsys, results->dYrange);
     364
     365    psFree (xErr);
     366    psFree (yErr);
     367    psFree (xRes);
     368    psFree (yRes);
     369    psFree (xResGood);
     370    psFree (yResGood);
    302371
    303372    psFree (x);
    … …  
    311380}
    312381
     382// set the bin closest to the corresponding value.  if USE_END is +/- 1,
     383// out-of-range saturates on lower/upper bin REGARDLESS of actual value
     384#define PS_BIN_FOR_VALUE(RESULT, VECTOR, VALUE, USE_END) { \
     385        psVectorBinaryDisectResult result; \
     386        psScalar tmpScalar; \
     387        tmpScalar.type.type = PS_TYPE_F32; \
     388        tmpScalar.data.F32 = (VALUE); \
     389        RESULT = psVectorBinaryDisect (&result, VECTOR, &tmpScalar); \
     390        switch (result) { \
     391          case PS_BINARY_DISECT_PASS: \
     392            break; \
     393          case PS_BINARY_DISECT_OUTSIDE_RANGE: \
     394            psTrace("psModules.astrom", 6, "selected bin outside range"); \
     395            if (USE_END == -1) { RESULT = 0; } \
     396            if (USE_END == +1) { RESULT = VECTOR->n - 1; } \
     397            break; \
     398          case PS_BINARY_DISECT_INVALID_INPUT: \
     399          case PS_BINARY_DISECT_INVALID_TYPE: \
     400            psAbort ("programming error"); \
     401            break; \
     402        } }
     403
     404# define PS_BIN_INTERPOLATE(RESULT, VECTOR, BOUNDS, BIN, VALUE) { \
     405        float dX, dY, Xo, Yo, Xt; \
     406        if (BIN == BOUNDS->n - 1) { \
     407            dX = 0.5*(BOUNDS->data.F32[BIN+1] - BOUNDS->data.F32[BIN-1]); \
     408            dY = VECTOR->data.F32[BIN] - VECTOR->data.F32[BIN-1]; \
     409            Xo = 0.5*(BOUNDS->data.F32[BIN+1] + BOUNDS->data.F32[BIN]); \
     410            Yo = VECTOR->data.F32[BIN]; \
     411        } else { \
     412            dX = 0.5*(BOUNDS->data.F32[BIN+2] - BOUNDS->data.F32[BIN]); \
     413            dY = VECTOR->data.F32[BIN+1] - VECTOR->data.F32[BIN]; \
     414            Xo = 0.5*(BOUNDS->data.F32[BIN+1] + BOUNDS->data.F32[BIN]); \
     415            Yo = VECTOR->data.F32[BIN]; \
     416        } \
     417        if (dY != 0) { \
     418            Xt = (VALUE - Yo)*dX/dY + Xo; \
     419        } else { \
     420            Xt = Xo; \
     421        } \
     422        Xt = PS_MIN (BOUNDS->data.F32[BIN+1], PS_MAX(BOUNDS->data.F32[BIN], Xt)); \
     423        psTrace("psModules.astrom", 6, "(Xo, Yo, dX, dY, Xt, Yt) is (%.2f %.2f %.2f %.2f %.2f %.2f)\n", \
     424                Xo, Yo, dX, dY, Xt, VALUE); \
     425        RESULT = Xt; }
     426
     427float pmAstromVectorRange (psVector *myVector, float minFrac, float maxFrac, float stdevGuess) {
     428
     429    psStats *stats = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics for min and max
     430    psHistogram *histogram = NULL;      // Histogram of the data
     431    psHistogram *cumulative = NULL;     // Cumulative histogram of the data
     432    float min = NAN, max = NAN;         // Mimimum and maximum values
     433
     434    // Get the minimum and maximum values
     435    if (!psVectorStats(stats, myVector, NULL, NULL, 0)) {
     436        psFree(stats);
     437        return NAN;
     438    }
     439    min = stats->min;
     440    max = stats->max;
     441    if (isnan(min) || isnan(max)) {
     442        psFree(stats);
     443        return NAN;
     444    }
     445
     446    psTrace("psModules.astrom", 5, "Data min/max is (%.2f, %.2f)\n", min, max);
     447
     448    // If all data points have the same value, then we set the appropriate members of stats and return.
     449    if (fabs(max - min) <= FLT_EPSILON) {
     450        psFree (stats);
     451        return 0.0;
     452    }
     453   
     454    // Define the histogram bin size. 
     455    float binSize = 0.001;
     456    long numBins = PS_MIN(100000, (max - min) / binSize); // Number of bins
     457    psTrace("psModules.astrom", 5, "Numbins is %ld\n", numBins);
     458    psTrace("psModules.astrom", 5, "Creating a robust histogram from data range (%.2f - %.2f)\n", min, max);
     459
     460    // allocate the histogram containers
     461    histogram = psHistogramAlloc(min, max, numBins);
     462    cumulative = psHistogramAlloc(min, max, numBins);
     463
     464    if (!psVectorHistogram(histogram, myVector, NULL, NULL, 0)) {
     465        // if psVectorHistogram returns false, we have a programming error
     466        psAbort ("Unable to generate histogram");
     467    }
     468    if (psTraceGetLevel("psModules.astrom") >= 8) {
     469        PS_VECTOR_PRINT_F32(histogram->bounds);
     470        PS_VECTOR_PRINT_F32(histogram->nums);
     471    }
     472
     473    // Convert the specific histogram to a cumulative histogram
     474    // The cumulative histogram data points correspond to the UPPER bound value (N < Bin[i+1])
     475    cumulative->nums->data.F32[0] = histogram->nums->data.F32[0];
     476    for (long i = 1; i < histogram->nums->n; i++) {
     477        cumulative->nums->data.F32[i] = cumulative->nums->data.F32[i-1] + histogram->nums->data.F32[i];
     478        cumulative->bounds->data.F32[i-1] = histogram->bounds->data.F32[i];
     479    }
     480    if (psTraceGetLevel("psModules.astrom") >= 8) {
     481        PS_VECTOR_PRINT_F32(cumulative->bounds);
     482        PS_VECTOR_PRINT_F32(cumulative->nums);
     483    }
     484
     485    // Find the bin which contains the first data point above the limit
     486    long totalDataPoints = cumulative->nums->data.F32[numBins - 1];
     487    psTrace("psModules.astrom", 6, "Total data points is %ld\n", totalDataPoints);
     488
     489    // find bin which is the lower bound of the limit value (value[bin] < f < value[bin+1]
     490    long binMin;
     491    PS_BIN_FOR_VALUE(binMin, cumulative->nums, minFrac * totalDataPoints, 0);
     492    psTrace("psModules.astrom", 6, "The bin is %ld (%.4f to %.4f)\n", binMin, cumulative->bounds->data.F32[binMin], cumulative->bounds->data.F32[binMin+1]);
     493
     494    // Linear interpolation to the limit value in bin units
     495    float valueMin;
     496    PS_BIN_INTERPOLATE (valueMin, cumulative->nums, cumulative->bounds, binMin, totalDataPoints * minFrac);
     497    psTrace("psModules.astrom", 6, "limit value is %f\n", valueMin);
     498
     499    // find bin which is the lower bound of the limit value (value[bin] < f < value[bin+1]
     500    long binMax;
     501    PS_BIN_FOR_VALUE(binMax, cumulative->nums, maxFrac * totalDataPoints, 0);
     502    psTrace("psModules.astrom", 6, "The bin is %ld (%.4f to %.4f)\n", binMax, cumulative->bounds->data.F32[binMax], cumulative->bounds->data.F32[binMax+1]);
     503
     504    // Linear interpolation to the limit value in bin units
     505    float valueMax;
     506    PS_BIN_INTERPOLATE (valueMax, cumulative->nums, cumulative->bounds, binMax, totalDataPoints * maxFrac);
     507    psTrace("psModules.astrom", 6, "limit value is %f\n", valueMax);
     508
     509    // Clean up
     510    psFree(histogram);
     511    psFree(cumulative);
     512    psFree(stats);
     513
     514    return (valueMax - valueMin);
     515}
    313516
    314517/******************************************************************************
    … …  
    634837        }
    635838
    636 # if 0
    637         char line[16];
    638         psFits *fits = psFitsOpen ("grid.image.fits", "w");
    639         psFitsWriteImage (fits, NULL, gridNP, 0, NULL);
    640         psFitsClose (fits);
    641         fprintf (stderr, "wrote grid image, press return to continue\n");
    642         fgets (line, 15, stdin);
    643 # endif
     839        if (psTraceGetLevel("psModules.astrom") >= 5) {
     840            char line[16];
     841            psFits *fits = psFitsOpen ("grid.image.fits", "w");
     842            psFitsWriteImage (fits, NULL, gridNP, 0, NULL);
     843            psFitsClose (fits);
     844            fprintf (stderr, "wrote grid image, press return to continue\n");
     845            fgets (line, 15, stdin);
     846        }
    644847
    645848        // only check bins with at least 1/2 of max bin
    646849        // XXX requiring at least 3 matches in bin
    647850        int minNpts = PS_MAX (0.5*imStats->max, 5);
    648         psTrace("psModule.astrom", 5, "minNpts: %d, min: %d, max: %d, median: %f, stdev: %f", minNpts, (int)(imStats->min), (int)(imStats->max), imStats->sampleMedian, imStats->sampleStdev);
     851        psTrace("psModule.astrom", 4, "minNpts: %d, min: %d, max: %d, median: %f, stdev: %f", minNpts, (int)(imStats->min), (int)(imStats->max), imStats->sampleMedian, imStats->sampleStdev);
    649852
    650853        // find the 'best' bin
    … …  
    687890
    688891        // XXX this function is crashing
    689         // pmAstromVisualPlotGridMatch(raw, ref, gridNP, stats->offset.x, stats->offset.y, maxOffpix, Scale, Offset);
     892        pmAstromVisualPlotGridMatch(raw, ref, gridNP, stats->offset.x, stats->offset.y, maxOffpix, Scale, Offset);
     893        pmAstromVisualPlotGridMatchOverlay(raw, ref);
    690894
    691895        psFree (imStats);
    … …  
    9621166*/
    9631167
     1168/*****************************************************************************/
     1169static void pmAstromMatchInfoFree (pmAstromMatchInfo *info)
     1170{
     1171    if (info == NULL) return;
     1172    return;
     1173}
     1174
     1175
     1176/*****************************************************************************/
     1177pmAstromMatchInfo *pmAstromMatchInfoAlloc()
     1178{
     1179    pmAstromMatchInfo *info = psAlloc (sizeof(pmAstromMatchInfo));
     1180    psMemSetDeallocator(info, (psFreeFunc) pmAstromMatchInfoFree);
     1181
     1182    info->match = NULL;
     1183    info->radius = NAN;
     1184
     1185    return (info);
     1186}
     1187
     1188// generate a unique set of matches (choose closest match)
     1189psArray *pmAstromRadiusMatchUniq (psArray *rawstars, psArray *refstars, psArray *matches) {
     1190
     1191    // I have the matches between the refstars and the rawstars. 
     1192    // For each refstar, find the single match which has the smallest radius and reject
     1193    // all others. 
     1194
     1195    // create an array of refstars->n arrays, each containing all of the matches for the
     1196    // given refstar. 
     1197
     1198    psArray *refstarMatches = psArrayAlloc (refstars->n);
     1199
     1200    for (int i = 0; i < matches->n; i++) {
     1201
     1202        pmAstromMatch *match = matches->data[i];
     1203       
     1204        // accumulate this refstar match on the array for this refstar (create if needed)
     1205        psArray *refSet = refstarMatches->data[match->ref];
     1206        if (!refSet) {
     1207            refstarMatches->data[match->ref] = psArrayAllocEmpty (8);
     1208            refSet = refstarMatches->data[match->ref];
     1209        }
     1210
     1211        pmAstromMatchInfo *matchInfo = pmAstromMatchInfoAlloc();
     1212
     1213        pmAstromObj *refStar = refstars->data[match->ref];
     1214        pmAstromObj *rawStar = rawstars->data[match->raw];
     1215       
     1216        matchInfo->match = match; // reference to the match of interest
     1217        matchInfo->radius = hypot (refStar->FP->x - rawStar->FP->x, refStar->FP->y - rawStar->FP->y);
     1218
     1219        psArrayAdd (refSet, 8, matchInfo); // matchInfo->match is just a reference
     1220        psFree (matchInfo);
     1221    }
     1222
     1223    // we now have a set of matches for each refstar and their distances; create a new set
     1224    // keeping only the closest entry for each match
     1225
     1226    psArray *unique = psArrayAllocEmpty (PS_MAX(16, matches->n / 2));
     1227    for (int i = 0; i < refstars->n; i++) {
     1228
     1229        psArray *refSet = refstarMatches->data[i];
     1230        if (!refSet) continue;
     1231        if (refSet->n == 0) continue; // not certain how this can happen...
     1232
     1233        if (refSet->n == 1) {
     1234            pmAstromMatchInfo *matchInfo = refSet->data[0];
     1235            psArrayAdd (unique, 32, matchInfo->match);
     1236            continue;
     1237        }
     1238
     1239        pmAstromMatchInfo *matchInfo = refSet->data[0];
     1240        float minRadius = matchInfo->radius;
     1241        pmAstromMatch *minMatch = matchInfo->match;
     1242        for (int j = 1; j < refSet->n; j++) {
     1243            pmAstromMatchInfo *matchInfo = refSet->data[j];
     1244            if (minRadius < matchInfo->radius) continue;
     1245            minMatch = matchInfo->match;
     1246            minRadius = matchInfo->radius;
     1247        }
     1248       
     1249        psArrayAdd (unique, 32, minMatch); // minMatch is just a reference to a match on matches,
     1250    }
     1251
     1252    psLogMsg ("psModules.astrom", 3, "generate unique matches to reference stars: %ld matches -> %ld matches\n", matches->n, unique->n);
     1253    psFree (refstarMatches);
     1254
     1255    return unique;
     1256}

  • trunk/psModules/src/astrom/pmAstrometryObjects.h

    r24021 r26260  
    5454typedef struct
    5555{
    56     int raw;                             ///< What is this?
    57     int ref;                             ///< What is this?
     56    int raw;                             ///< reference to the rawstar entry
     57    int ref;                             ///< reference to the refstar entry
    5858}
    5959pmAstromMatch;
     60
     61
     62/*
     63 * The pmAstromMatchInfo structure is used to generate a unique set of matches
     64 */
     65typedef struct
     66{
     67    pmAstromMatch *match;               ///< reference to the match
     68    float radius;                       ///< distance between the object
     69}
     70pmAstromMatchInfo;
    6071
    6172
    … …  
    8596    int     nMatch;                     ///<
    8697    double  nSigma;                     ///<
     98    double  dXsys;                      ///< systematic error in X
     99    double  dYsys;                      ///< systematic error in Y
     100    double  dXrange;                    ///< 10% - 90% range X residuals (unmasked, high S/N)
     101    double  dYrange;                    ///< 10% - 90% range Y residuals (unmasked, high S/N)
    87102}
    88103pmAstromFitResults;
    … …  
    121136);
    122137
     138psArray *pmAstromRadiusMatchUniq (psArray *refstars, psArray *rawstars, psArray *matches);
    123139
    124140pmAstromStats *pmAstromStatsAlloc(void);
    … …  
    343359);
    344360
     361float pmAstromVectorRange (psVector *myVector, float minFrac, float maxFrac, float stdevGuess);
     362
    345363/// @}
    346364#endif // PM_ASTROMETRY_OBJECTS_H

  • trunk/psModules/src/astrom/pmAstrometryVisual.c

    r25729 r26260  
    450450
    451451    graphdata.color = KapaColorByName ("red");
    452     graphdata.style = 1;
     452    graphdata.style = 0;
    453453
    454454    //overplot clumpy regions excluded from analysis
    … …  
    905905    KapaPlotVector (kapa, gridNP->numCols, horizontalIndices, "x");
    906906    KapaPlotVector (kapa, gridNP->numCols, horizHistSlice, "y");
     907
    907908    float xslice[2] = {offsetX - Scale / 2., offsetX - Scale / 2.};
    908909    float yslice[2] = {-5, 100};
     910    graphdata.style = 0;
    909911    graphdata.color = KapaColorByName("red");
    910912    KapaPrepPlot(kapa, 2, &graphdata);
    … …  
    927929    KapaPlotVector (kapa, gridNP->numRows, vertHistSlice, "x");
    928930    KapaPlotVector (kapa, gridNP->numRows, verticalIndices, "y");
     931
    929932    yslice[0] = yslice[1] = offsetY - Scale / 2.;
    930933    xslice[0] = -5; xslice[1] = 100;
     934    graphdata.style = 0;
    931935    graphdata.color = KapaColorByName("red");
    932936    KapaPrepPlot(kapa, 2, &graphdata);
    … …  
    940944} // end of pmAstromVisualPlotGridMatch
    941945
     946
     947bool pmAstromVisualPlotGridMatchOverlay (const psArray *raw,
     948                                         const psArray *ref)
     949{
     950    //make sure we want to plot this
     951    if (!pmVisualIsVisual() || !plotGridMatch) return true;
     952    if (!pmVisualInitWindow(&kapa2, "psastro:plots")){
     953        return false;
     954    }
     955
     956    Graphdata graphdata;
     957    psVector *xPlot = psVectorAlloc (PS_MAX(raw->n, ref->n), PS_TYPE_F32); // x data points
     958    psVector *yPlot = psVectorAlloc (PS_MAX(raw->n, ref->n), PS_TYPE_F32); // y data points
     959    psVector *zPlot = psVectorAlloc (PS_MAX(raw->n, ref->n), PS_TYPE_F32); // y data points
     960
     961    // set up plot information
     962    KapaClearPlots(kapa2);
     963    KapaInitGraph(&graphdata);
     964
     965    KapaSetFont(kapa2, "helvetica", 14);
     966    KapaBox(kapa2, &graphdata);
     967    KapaSendLabel (kapa2, "X (FP)", KAPA_LABEL_XM);
     968    KapaSendLabel (kapa2, "Y (FP)", KAPA_LABEL_YM);
     969    KapaSendLabel (kapa2, "pmAstromGridAngle residuals. Box: Correlation Peak.", KAPA_LABEL_XP);
     970
     971    // plot the REF data.  (also calculate the plot ranges, accumulate the plot vectors)
     972    graphdata.xmin = +INT_MAX;
     973    graphdata.xmax = -INT_MAX;
     974    graphdata.ymin = +INT_MAX;
     975    graphdata.ymax = -INT_MAX;
     976    for (int i = 0; i < ref->n; i++) {
     977        pmAstromObj *obj = ref->data[i];
     978        graphdata.xmin = PS_MIN(graphdata.xmin, obj->FP->x);
     979        graphdata.xmax = PS_MAX(graphdata.xmax, obj->FP->x);
     980        graphdata.ymin = PS_MIN(graphdata.ymin, obj->FP->y);
     981        graphdata.ymax = PS_MAX(graphdata.ymax, obj->FP->y);
     982        xPlot->data.F32[i] = obj->FP->x;
     983        yPlot->data.F32[i] = obj->FP->y;
     984        zPlot->data.F32[i] = obj->Mag;
     985    }
     986    xPlot->n = yPlot->n = zPlot->n = ref->n;
     987    KapaSetLimits(kapa2, &graphdata);
     988
     989    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
     990    psVectorStats (stats, zPlot, NULL, NULL, 0);
     991    float zero = stats->sampleMedian + 3.0;
     992    float range = 6.0;
     993
     994    for (int i = 0; i < zPlot->n; i++) {
     995        float value = (zero - zPlot->data.F32[i]) / range;
     996        zPlot->data.F32[i] = PS_MAX(0.0, PS_MIN(1.0, value));
     997    }
     998
     999    // the point size will be scaled from the z vector
     1000    graphdata.style = 2;
     1001    graphdata.ptype = 7;
     1002    graphdata.size = -1;
     1003    graphdata.color = KapaColorByName ("black");
     1004
     1005    KapaPrepPlot   (kapa2, xPlot->n, &graphdata);
     1006    KapaPlotVector (kapa2, xPlot->n, xPlot->data.F32, "x");
     1007    KapaPlotVector (kapa2, yPlot->n, yPlot->data.F32, "y");
     1008    KapaPlotVector (kapa2, zPlot->n, zPlot->data.F32, "z");
     1009
     1010    // plot the RAW data (keep previous limits)
     1011    for (int i = 0; i < raw->n; i++) {
     1012        pmAstromObj *obj = raw->data[i];
     1013        xPlot->data.F32[i] = obj->FP->x;
     1014        yPlot->data.F32[i] = obj->FP->y;
     1015        zPlot->data.F32[i] = obj->Mag;
     1016    }
     1017    xPlot->n = yPlot->n = zPlot->n = raw->n;
     1018
     1019    psStatsInit(stats);
     1020    psVectorStats (stats, zPlot, NULL, NULL, 0);
     1021    zero = stats->sampleMedian + 3.0;
     1022    range = 6.0;
     1023
     1024    for (int i = 0; i < zPlot->n; i++) {
     1025        float value = (zero - zPlot->data.F32[i]) / range;
     1026        zPlot->data.F32[i] = PS_MAX(0.0, PS_MIN(1.0, value));
     1027    }
     1028
     1029    // the point size will be scaled from the z vector
     1030    graphdata.style = 2;
     1031    graphdata.ptype = 7;
     1032    graphdata.size = -1;
     1033    graphdata.color = KapaColorByName ("red");
     1034
     1035    KapaPrepPlot   (kapa2, xPlot->n, &graphdata);
     1036    KapaPlotVector (kapa2, xPlot->n, xPlot->data.F32, "x");
     1037    KapaPlotVector (kapa2, yPlot->n, yPlot->data.F32, "y");
     1038    KapaPlotVector (kapa2, zPlot->n, zPlot->data.F32, "z");
     1039
     1040    pmVisualAskUser(&plotGridMatch);
     1041    psFree(xPlot);
     1042    psFree(yPlot);
     1043    psFree(zPlot);
     1044    psFree(stats);
     1045    return true;
     1046}
    9421047
    9431048bool pmAstromVisualPlotTweak (psVector *xHist, // Smoothed Horizontal cut through the histogram

  • trunk/psModules/src/astrom/pmAstrometryVisual.h

    r23487 r26260  
    4545                                  );
    4646
     47
     48bool pmAstromVisualPlotGridMatchOverlay (const psArray *raw,
     49                                         const psArray *ref);
    4750
    4851/**

  • trunk/psModules/src/camera

    • Property svn:mergeinfo deleted

  • trunk/psModules/src/camera/pmReadoutFake.c

    • Property svn:mergeinfo deleted

  • trunk/psModules/src/imcombine

    • Property svn:mergeinfo deleted

  • trunk/psModules/src/imcombine/pmStack.c

    • Property svn:mergeinfo deleted

  • trunk/psModules/src/imcombine/pmStack.h

    • Property svn:mergeinfo deleted

  • trunk/psModules/src/objects

    • Property svn:mergeinfo deleted

  • trunk/psModules/src/objects/pmPSFtryMetric.c

    r25754 r26260  
    7676
    7777    pmSourceVisualPlotPSFMetric (psfTry);
     78    pmSourceVisualPlotPSFMetricSubpix (psfTry);
    7879
    7980    psFree (stats);

  • trunk/psModules/src/objects/pmPSFtryModel.c

    r25819 r26260  
    7272    }
    7373
    74     // XXX set the min number of needed source more carefully (depends on psfTrendMode?)
    75     int orderMax = PS_MAX (options->psfTrendNx, options->psfTrendNy);
    76     if ((sources->n < 15) && (orderMax >= 3)) orderMax = 2;
    77     if ((sources->n < 11) && (orderMax >= 2)) orderMax = 1;
    78     if ((sources->n <  8) && (orderMax >= 1)) orderMax = 0;
     74    // hard limit on minimum number of stars
    7975    if ((sources->n <  3)) {
    8076        psError (PS_ERR_UNKNOWN, true, "failed to determine PSF parameters");
    … …  
    8278    }
    8379
    84     int orderMin;
     80    // this is a bit tricky, because we have two cases (MAP vs POLY), and they have a different
     81    // definition for 'order' (order_MAP = order_POLY + 1).  in addition, we have a
     82    // user-specified MAX order, which we should respect, regardless of the mode
     83
     84    // set the max order (0 = constant) which the number of psf stars can support:
     85    // rule of thumb: require 3 stars per 'cell' (order+1)^2
     86    int MaxOrderForStars = 0;
     87    if (sources->n >= 12) MaxOrderForStars = 1; // 4 cells
     88    if (sources->n >= 27) MaxOrderForStars = 2; // 9 cells
     89    if (sources->n >= 48) MaxOrderForStars = 3; // 16 cells
     90    if (sources->n >  75) MaxOrderForStars = 4; // 25 cells
     91
     92    int orderMax = PS_MAX (options->psfTrendNx, options->psfTrendNy);
     93    int orderMin = 0;
    8594    if (options->psfTrendMode == PM_TREND_MAP) {
    86         orderMin = 1;
    87         orderMax = PS_MAX(orderMax, 1);
    88     } else {
    89         orderMin = 0;
    90     }
     95        MaxOrderForStars ++;
     96        orderMin ++;
     97    }
     98    orderMax = PS_MIN (orderMax, MaxOrderForStars);
    9199
    92100    // save the raw source mask (generated by pmPSFtryFitEXT)

  • trunk/psModules/src/objects/pmSourceIO.c

    r25907 r26260  
    980980        if (!tableHeader) psAbort("cannot read table header");
    981981
     982        char *xtension = psMetadataLookupStr (NULL, tableHeader, "XTENSION");
     983        if (!xtension) psAbort("cannot read table type");
     984        if (strcmp (xtension, "BINTABLE")) {
     985            psWarning ("no binary table in extension %s, skipping\n", dataname);
     986            return false;
     987        }
     988
    982989        char *exttype = psMetadataLookupStr (NULL, tableHeader, "EXTTYPE");
    983990        if (!exttype) psAbort("cannot read table type");

  • trunk/psModules/src/objects/pmSourceIO_MatchedRefs.c

    r24801 r26260  
    120120                        psMetadataAdd (row, PS_LIST_TAIL, "X_CHIP",   PS_DATA_F32, "x coord on chip",              raw->chip->x);
    121121                        psMetadataAdd (row, PS_LIST_TAIL, "Y_CHIP",   PS_DATA_F32, "y coord on chip",              raw->chip->y);
     122                        psMetadataAdd (row, PS_LIST_TAIL, "X_CHIP_FIT",PS_DATA_F32, "x fitted coord on chip",      ref->chip->x);
     123                        psMetadataAdd (row, PS_LIST_TAIL, "Y_CHIP_FIT",PS_DATA_F32, "y fitted coord on chip",      ref->chip->y);
    122124                        psMetadataAdd (row, PS_LIST_TAIL, "X_FPA",    PS_DATA_F32, "x coord on focal plane",       raw->FP->x);
    123125                        psMetadataAdd (row, PS_LIST_TAIL, "Y_FPA",    PS_DATA_F32, "y coord on focal plane",       raw->FP->y);

  • trunk/psModules/src/objects/pmSourceVisual.c

    r25979 r26260  
    7373
    7474    float range;
     75    range = graphdata.xmax - graphdata.xmin;
     76    graphdata.xmax += 0.05*range;
     77    graphdata.xmin -= 0.05*range;
     78    range = graphdata.ymax - graphdata.ymin;
     79    graphdata.ymax += 0.05*range;
     80    graphdata.ymin -= 0.05*range;
     81
     82    // better choice for range?
     83    // graphdata.xmin = -17.0;
     84    // graphdata.xmax =  -9.0;
     85    graphdata.ymin = -0.51;
     86    graphdata.ymax = +0.51;
     87
     88    KapaSetLimits (kapa1, &graphdata);
     89
     90    KapaSetFont (kapa1, "helvetica", 14);
     91    KapaBox (kapa1, &graphdata);
     92    KapaSendLabel (kapa1, "PSF Mag", KAPA_LABEL_XM);
     93    KapaSendLabel (kapa1, "Ap Mag - PSF Mag", KAPA_LABEL_YM);
     94
     95    graphdata.color = KapaColorByName ("black");
     96    graphdata.ptype = 2;
     97    graphdata.size = 0.5;
     98    graphdata.style = 2;
     99    graphdata.etype |= 0x01;
     100
     101    KapaPrepPlot (kapa1, n, &graphdata);
     102    KapaPlotVector (kapa1, n, x->data.F32, "x");
     103    KapaPlotVector (kapa1, n, y->data.F32, "y");
     104    KapaPlotVector (kapa1, n, dy->data.F32, "dym");
     105    KapaPlotVector (kapa1, n, dy->data.F32, "dyp");
     106
     107    psFree (x);
     108    psFree (y);
     109    psFree (dy);
     110
     111    pmVisualAskUser(NULL);
     112    return true;
     113}
     114
     115bool pmSourceVisualPlotPSFMetricSubpix (pmPSFtry *psfTry) {
     116
     117    KapaSection section;  // put the positive profile in one and the residuals in another?
     118    Graphdata graphdata;
     119
     120    if (!pmVisualIsVisual()) return true;
     121
     122    if (kapa1 == -1) {
     123        kapa1 = KapaOpenNamedSocket ("kapa", "pmSource:plots");
     124        if (kapa1 == -1) {
     125            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
     126            pmVisualSetVisual(false);
     127            return false;
     128        }
     129    }
     130
     131    KapaClearSections (kapa1);
     132    KapaInitGraph (&graphdata);
     133
     134    int n;
     135    float range;
     136    psVector *x = psVectorAllocEmpty (psfTry->sources->n, PS_TYPE_F32);
     137    psVector *y = psVectorAllocEmpty (psfTry->sources->n, PS_TYPE_F32);
     138    psVector *dy = psVectorAllocEmpty(psfTry->sources->n, PS_TYPE_F32);
     139
     140    // section a: fractional-x pixel
     141    section.dx = 1.0;
     142    section.dy = 0.5;
     143    section.x = 0.0;
     144    section.y = 0.0;
     145    section.name = NULL;
     146    psStringAppend (&section.name, "a1");
     147    KapaSetSection (kapa1, &section);
     148    psFree (section.name);
     149
     150    graphdata.xmin = +32.0;
     151    graphdata.xmax = -32.0;
     152    graphdata.ymin = +32.0;
     153    graphdata.ymax = -32.0;
     154
     155    // construct the plot vectors
     156    n = 0;
     157    for (int i = 0; i < psfTry->sources->n; i++) {
     158        if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) continue;
     159
     160        pmSource *source = psfTry->sources->data[i];
     161        x->data.F32[n] = source->modelEXT->params->data.F32[PM_PAR_XPOS] - (int)source->modelEXT->params->data.F32[PM_PAR_XPOS];
     162
     163        y->data.F32[n] = psfTry->metric->data.F32[i];
     164        dy->data.F32[n] = psfTry->metricErr->data.F32[i];
     165        graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]);
     166        graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]);
     167        graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[n]);
     168        graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[n]);
     169        n++;
     170    }
     171    x->n = y->n = dy->n = n;
     172
     173    range = graphdata.xmax - graphdata.xmin;
     174    graphdata.xmax += 0.05*range;
     175    graphdata.xmin -= 0.05*range;
     176    range = graphdata.ymax - graphdata.ymin;
     177    graphdata.ymax += 0.05*range;
     178    graphdata.ymin -= 0.05*range;
     179
     180    // better choice for range?
     181    // graphdata.xmin = -17.0;
     182    // graphdata.xmax =  -9.0;
     183    graphdata.ymin = -0.51;
     184    graphdata.ymax = +0.51;
     185
     186    KapaSetLimits (kapa1, &graphdata);
     187
     188    KapaSetFont (kapa1, "helvetica", 14);
     189    KapaBox (kapa1, &graphdata);
     190    KapaSendLabel (kapa1, "PSF Mag", KAPA_LABEL_XM);
     191    KapaSendLabel (kapa1, "Ap Mag - PSF Mag", KAPA_LABEL_YM);
     192
     193    graphdata.color = KapaColorByName ("black");
     194    graphdata.ptype = 2;
     195    graphdata.size = 0.5;
     196    graphdata.style = 2;
     197    graphdata.etype |= 0x01;
     198
     199    KapaPrepPlot (kapa1, n, &graphdata);
     200    KapaPlotVector (kapa1, n, x->data.F32, "x");
     201    KapaPlotVector (kapa1, n, y->data.F32, "y");
     202    KapaPlotVector (kapa1, n, dy->data.F32, "dym");
     203    KapaPlotVector (kapa1, n, dy->data.F32, "dyp");
     204
     205    // *** section b: fractional-x pixel
     206    section.dx = 1.0;
     207    section.dy = 0.5;
     208    section.x = 0.0;
     209    section.y = 0.5;
     210    section.name = NULL;
     211    psStringAppend (&section.name, "a2");
     212    KapaSetSection (kapa1, &section);
     213    psFree (section.name);
     214
     215    graphdata.xmin = +32.0;
     216    graphdata.xmax = -32.0;
     217    graphdata.ymin = +32.0;
     218    graphdata.ymax = -32.0;
     219
     220    // construct the plot vectors
     221    n = 0;
     222    for (int i = 0; i < psfTry->sources->n; i++) {
     223        if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) continue;
     224
     225        pmSource *source = psfTry->sources->data[i];
     226        x->data.F32[n] = source->modelEXT->params->data.F32[PM_PAR_YPOS] - (int)source->modelEXT->params->data.F32[PM_PAR_YPOS];
     227
     228        y->data.F32[n] = psfTry->metric->data.F32[i];
     229        dy->data.F32[n] = psfTry->metricErr->data.F32[i];
     230        graphdata.xmin = PS_MIN(graphdata.xmin, x->data.F32[n]);
     231        graphdata.xmax = PS_MAX(graphdata.xmax, x->data.F32[n]);
     232        graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[n]);
     233        graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[n]);
     234        n++;
     235    }
     236    x->n = y->n = dy->n = n;
     237
    75238    range = graphdata.xmax - graphdata.xmin;
    76239    graphdata.xmax += 0.05*range;

  • trunk/psModules/src/objects/pmSourceVisual.h

    r25754 r26260  
    1919bool pmSourceVisualPSFModelResid (pmTrend2D *trend, psVector *x, psVector *y, psVector *param, psVector *mask);
    2020bool pmSourceVisualPlotPSFMetric (pmPSFtry *try);
     21bool pmSourceVisualPlotPSFMetricSubpix (pmPSFtry *try);
     22
    2123bool pmSourceVisualShowModelFit (pmSource *source);
    2224bool pmSourceVisualShowModelFits (pmPSF *psf, psArray *sources, psImageMaskType maskVal);
Note: See TracChangeset for help on using the changeset viewer.