IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 26259


Ignore:
Timestamp:
Nov 22, 2009, 2:57:41 PM (16 years ago)
Author:
eugene
Message:

various fixes to psastro:

1) added bootstrap resampling to zero point error analysis
2) added iterative clump removal from refstars and rawstars
3) added unique reference match option
4) some improved visualizations
5) improved mosaic iterations

Location:
trunk/psastro
Files:
24 edited
1 copied

Legend:

Unmodified
Added
Removed

  • trunk/psastro/src/psastro.h

    r26173 r26259  
    8888bool              psastroLuminosityFunctionPlot(psVector *lnMag, psVector *Mag, pmLumFunc *lumFunc, pmLumFunc *rawFunc);
    8989psArray          *psastroRemoveClumps (psArray *input, int scale);
     90psArray          *psastroRemoveClumpsIterate (psArray *input, int scale, int nIter);
     91bool              psastroRemoveClumpsRawstars (pmConfig *config);
    9092
    9193
    9294// utility functions:
    93 bool              psastroUpdateChipToFPA (pmFPA *fpa, pmChip *chip, psArray *rawstars, psArray *refstars);
     95bool              psastroUpdateChipToFPA (pmFPA *fpa, pmChip *chip);
    9496bool              psastroWriteStars (char *filename, psArray *sources);
    9597bool              psastroWriteTransform (psPlaneTransform *map);
    … …  
    149151bool              psastroZeroPoint (pmConfig *config);
    150152psVector         *psastroZeroPointReadoutAccum(psVector *dMag, pmReadout *readout, float exptime);
    151 bool              psastroZeroPointAnalysis (psMetadata *header, psVector *dMag, float zeropt, float edgeFraction, bool useMean);
     153bool              psastroZeroPointAnalysis (psMetadata *header, psVector *dMag, float zeropt, psMetadata *recipe);
    152154bool              psastroZeroPointFromRecipe (float *zeropt, float *exptime, float *ghostMaxMag, pmFPA *fpa, psMetadata *recipe);
    153155
    154 psStats          *psastroStatsPercentile (psVector *myVector, float flimit);
     156psStats          *psastroStatsPercentile (psVector *myVector, psMetadata *recipe);
     157float             psastroStatsPercentileValue (psHistogram *histogram, psHistogram *cumulative, psVector *myVector, float flimit);
    155158
    156159// masking functions

  • trunk/psastro/src/psastroAnalysis.c

    r24647 r26259  
    6969    }
    7070
     71    if (!psastroRemoveClumpsRawstars(config)) {
     72        psError (PSASTRO_ERR_UNKNOWN, false, "failed to remove RAWSTAR clumps\n");
     73        return false;
     74    }
     75
    7176    // load the reference stars overlapping the data stars
    7277    psArray *refs = psastroLoadRefstars(config, "PSASTRO.INPUT");
    … …  
    127132    psastroMaskUpdates (config);
    128133
    129     // XXX how do we specify stack astrometry?
    130     // psastroStackAstrom (config, refs);
    131 
    132134    return true;
    133135}

  • trunk/psastro/src/psastroArguments.c

    r23242 r26259  
    6262    // define the reference astrometry file
    6363    status = pmConfigFileSetsMD (config->arguments, &argc, argv, "ASTROM.MODEL", "-astrommodel", "-astrommodellist");
    64     if (status) {
    65         // if supplied, assume -fixchips is desired
    66         psMetadataAddBool (config->arguments, PS_LIST_TAIL, "PSASTRO.FIX.CHIPS", PS_META_REPLACE, "", true);
    67     }
    6864
    6965    // define the reference astrometry file

  • trunk/psastro/src/psastroChipAstrom.c

    r25299 r26259  
    5656
    5757                // select the raw objects for this readout
    58                 psArray *rawstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.RAWSTARS");
     58                psArray *rawstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.RAWSTARS.SUBSET");
    5959                if (rawstars == NULL) { continue; }
    6060
    6161                // select the raw objects for this readout
    62                 psArray *refstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.REFSTARS");
     62                psArray *refstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.REFSTARS.SUBSET");
    6363                if (refstars == NULL) { continue; }
    6464
    … …  
    6666                if ((rawstars->n < 4) || (refstars->n < 4)) {
    6767                    readout->data_exists = false;
    68                     psLogMsg ("psastro", 3, "insufficient rawstars (%ld) or refstars (%ld)",
    69                               rawstars->n, refstars->n);
     68                    psLogMsg ("psastro", 3, "insufficient rawstars (%ld) or refstars (%ld)", rawstars->n, refstars->n);
    7069                    continue;
    7170                }

  • trunk/psastro/src/psastroChooseRefstars.c

    r24806 r26259  
    134134
    135135                psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.REFSTARS", PS_DATA_ARRAY, "astrometry matches", refstars);
     136
     137                // generate a reduced subset excluding the clumps
     138                // XXX do we need both REFSTARS and SUBSET?
     139                psArray *subset = psastroRemoveClumpsIterate(refstars, 150, 3);
     140                psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.REFSTARS.SUBSET", PS_DATA_ARRAY, "astrometry objects", subset);
     141
    136142                psFree (refstars);
     143                psFree (subset);
    137144                psFree (extent);
    138145
    139146                if (matchLumFunc) {
    140                     // in this case, no PSASTRO.REFSTARS is added to readout->analysis
     147                    // limit the total magnitude range of PSASTRO.REFSTARS.SUBSET based on
     148                    // overlapping luminosity functions
    141149                    if (!psastroRefstarSubset (readout)) {
    142150                        psError(PSASTRO_ERR_DATA, false, "Can't determine an appropriate refstar subset\n");

  • trunk/psastro/src/psastroConvert.c

    r23303 r26259  
    127127
    128128    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.RAWSTARS", PS_DATA_ARRAY, "astrometry objects", rawStars);
     129
    129130    psLogMsg ("psastro", 4, "loaded %ld sources, using %ld of %ld good sources (inst mag: %f to %f)\n", sources->n, rawStars->n, inStars->n, mMin, mMax);
    130131    psLogMsg ("psastro", 4, "skip reasons: mode: %d, faint: %d, bright: %d, inf: %d\n", nModeSkip, nFaintSkip, nBrightSkip, nInfSkip);
    … …  
    153154
    154155        psF32 *PAR = model->params->data.F32;
     156        psF32 *dPAR = model->dparams->data.F32;
    155157
    156158        pmAstromObj *obj = pmAstromObjAlloc ();
    157159
    158160        // is the source magnitude calibrated in any sense?
    159         obj->pix->x = PAR[PM_PAR_XPOS];
    160         obj->pix->y = PAR[PM_PAR_YPOS];
     161        obj->pix->x    = PAR[PM_PAR_XPOS];
     162        obj->pix->y    = PAR[PM_PAR_YPOS];
     163        obj->pix->xErr = dPAR[PM_PAR_XPOS];
     164        obj->pix->yErr = dPAR[PM_PAR_YPOS];
    161165        obj->Mag = source->psfMag;
     166        obj->dMag = source->errMag;
    162167
    163168        // XXX do we have the information giving the readout and cell offset?

  • trunk/psastro/src/psastroDemoDump.c

    r21409 r26259  
    160160
    161161                // select the raw objects for this readout
    162                 psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
     162                psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS.SUBSET");
    163163                if (rawstars == NULL) continue;
    164164
    165165                // select the raw objects for this readout
    166                 psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS");
     166                psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS.SUBSET");
    167167                if (refstars == NULL) continue;
    168168                psTrace ("psastro", 4, "Trying %ld refstars\n", refstars->n);

  • trunk/psastro/src/psastroFixChips.c

    r21422 r26259  
    248248        badAstrom |= fabs(obsAngle   - refAngle)   > angleTol;
    249249
    250         fprintf (stderr, "chip %d, angle: %f, pixel: %f,%f\n", view->chip, obsAngle - refAngle, obsCoord.x - refCoord.x, obsCoord.y - refCoord.y);
     250        fprintf (stderr, "chip %d, angle: %f, pixel: %f,%f  : ", view->chip, obsAngle - refAngle, obsCoord.x - refCoord.x, obsCoord.y - refCoord.y);
     251        if (badAstrom) {
     252            fprintf (stderr, "BAD ASTROM\n");
     253        } else {
     254            fprintf (stderr, "GOOD ASTROM\n");
     255        }
    251256
    252257        // XXX for now, just use first readout
    … …  
    265270        psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DY", PS_META_REPLACE, "chip y offset wrt model", obsCoord.y - refCoord.y);
    266271        psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DT", PS_META_REPLACE, "chip rot offset wrt model", obsAngle - refAngle);
     272
     273        continue;
    267274
    268275        // for successful chips, save the measured offsets in the header
    … …  
    319326        }
    320327
    321         psastroUpdateChipToFPA (input->fpa, obsChip, rawstars, refstars);
     328        psastroUpdateChipToFPA (input->fpa, obsChip);
    322329
    323330        // XXX update the header with info to reflect the failure

  • trunk/psastro/src/psastroFixChipsTest.c

    r23989 r26259  
    246246        }
    247247
    248         psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
     248        psastroUpdateChipToFPA (fpa, chip);
    249249
    250250        // XXX skip the re-fit step for a test

  • trunk/psastro/src/psastroLoadGhosts.c

    r24650 r26259  
    119119                if (! readout->data_exists) { continue; }
    120120
    121                 // select the raw objects for this readout (loaded in psastroExtract.c)
     121                // select the raw objects for this readout (loaded in psastroChooseRefstars.c)
     122                // XXX : note that we place limits on the refstar sample in psastroChooseRefstars.c:
     123                // 1) on chip and 2) < PSASTRO.MAX.NREF. magnitude limits and clump exclusion are only
    122124                psArray *refstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.REFSTARS");
    123125                if (refstars == NULL) { continue; }

  • trunk/psastro/src/psastroLoadRefstars.c

    r24805 r26259  
    5050        return false;
    5151    }
     52
     53    // the name in the recipe may be one of:
     54    // (A) the actual directory name
     55    // (B) a reference to the name in the PSASTRO.CATDIRS folder in site.config
     56    // (C) a reference to a folder in the PSASTRO.CATDIRS folder in site.config, containing multiple copy locations
    5257    char *catdir_virtual = psMetadataLookupStr(&status, catdirs, catdir_recipe);
     58
     59    psMetadata *catdir_folder = psMetadataLookupMetadata(&status, catdirs, catdir_recipe);
     60    if (catdir_folder) {
     61        // randomly choose one of the entries
     62        psLogMsg ("psastro", 3, "choosing catdir_folder\n");
     63        int nEntry = catdir_folder->list->n;
     64
     65        psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); 
     66        double frnd = psRandomUniform(rng);
     67        int entry = PS_MIN(nEntry - 1, PS_MAX(0, nEntry * frnd));
     68        psFree(rng);
     69
     70        psMetadataItem *item = psListGet(catdir_folder->list, entry);
     71        if (item->type != PS_DATA_STRING) {
     72            psError(PSASTRO_ERR_CONFIG, true, "Invalid entry in PSASTRO.CATDIR folder: %s\n", item->name);
     73            return false;
     74        }
     75        catdir_virtual = item->data.str;
     76    }
     77
    5378    char *catdir = (catdir_virtual == NULL) ? catdir_recipe : catdir_virtual;
    5479
    … …  
    5681    psString CATDIR = pmConfigConvertFilename(catdir, config, false, false); // Resolved filename
    5782    PS_ASSERT (CATDIR, NULL);
     83
     84    psLogMsg ("psastro", 3, "looking up reference objects in %s\n", CATDIR);
    5885
    5986    char *getstarCommand = psStringCopy(psMetadataLookupStr(NULL, recipe, "DVO.GETSTAR"));

  • trunk/psastro/src/psastroMaskUpdates.c

    r26173 r26259  
    193193
    194194                // select the raw objects for this readout
     195                // XXX : note that we place limits on the refstar sample in psastroChooseRefstars.c:
     196                // 1) on chip and 2) < PSASTRO.MAX.NREF. magnitude limits and clump exclusion are only
     197                // applied to the SUBSETs
    195198                psArray *refstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.REFSTARS");
    196199                if (!refstars) continue;

  • trunk/psastro/src/psastroMosaicAstrom.c

    r21409 r26259  
    2020
    2121    bool status;
    22     char filename[256];
    2322
    2423    // select the current recipe
    … …  
    3837    pmFPA *fpa = input->fpa;
    3938
    40     // before we do object matches, we need to (optionally) fix failed chips.  We compare chips with
    41     // the supplied mosaic model.  Adjust significant outliers to match model.
    42     # if (0)
    43     if (!psastroFixChips (config, recipe)) {
    44         psError(PSASTRO_ERR_UNKNOWN, false, "failed to align problematic chips");
    45         return false;
    46     }
    47     if (!psastroFixChipsTest (config, recipe)) {
    48         psError(PSASTRO_ERR_UNKNOWN, false, "failed to align problematic chips");
    49         return false;
    50     }
    51     # endif
    52 
    5339    char *outroot = psMetadataLookupStr (&status, config->arguments, "OUTPUT");
    5440    if (!status || !outroot) psAbort ("Can't find outroot on config->arguments");
    5541
    56     if (!psastroMosaicFit (fpa, recipe, outroot, 0)) return false;
    57     if (!psastroMosaicFit (fpa, recipe, outroot, 1)) return false;
    58     if (!psastroMosaicFit (fpa, recipe, outroot, 2)) return false;
    59     if (!psastroMosaicFit (fpa, recipe, outroot, 3)) return false;
     42    int nIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NITER");
     43    if (!status) psAbort ("missing config value");
     44
     45    // this should be in a loop with nIter =
     46    for (int iter = 0; iter < nIter; iter++) {
     47        if (!psastroMosaicFit (fpa, recipe, outroot, iter)) return false;
     48    }
    6049
    6150    // now fit the chips under the common distortion with higher-order terms
    6251    // first, re-perform the match with a slightly tighter circle
    63     if (!psastroMosaicSetMatch (fpa, recipe, 4)) {
     52    if (!psastroMosaicSetMatch (fpa, recipe, nIter)) {
    6453        psError(PSASTRO_ERR_UNKNOWN, false, "failed to match raw and ref stars for mosaic (4th pass)");
    6554        return false;
    6655    }
    67     if (!psastroMosaicChipAstrom (fpa, recipe, 4)) {
     56    if (!psastroMosaicChipAstrom (fpa, recipe, nIter)) {
    6857        psError(PSASTRO_ERR_UNKNOWN, false, "failed to measure chip astrometry in mosaic mode (4th pass)");
    6958        return false;
    7059    }
     60
    7161    if (psTraceGetLevel("psastro.dump") > 0) {
    72         snprintf (filename, 256, "%s.10.dat", outroot);
     62        // the last filename (see filenames in psastroMosaicFit)
     63        char filename[256];
     64        snprintf (filename, 256, "%s.%d.dat", outroot, 2*nIter + 2);
    7365        psastroDumpMatches (fpa, filename);
    7466    }

  • trunk/psastro/src/psastroMosaicFPtoTP.c

    r21409 r26259  
    4545
    4646                // select the raw objects for this readout
    47                 psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
     47                psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS.SUBSET");
    4848                if (rawstars == NULL) { continue; }
    4949
    5050                // select the raw objects for this readout
    51                 psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS");
     51                psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS.SUBSET");
    5252                if (refstars == NULL) { continue; }
    5353
    … …  
    126126}
    127127
     128// apply the rotation and scale to all stars in PSASTRO.REFSTARS (also adjusts
     129// PSASTRO.REFSTARS.SUBSET since they are the same pointers)
    128130bool psastroMosaicApplyRotAndScale (pmFPA *fpa, psPlaneTransform *TPtoFP) {
    129131

  • trunk/psastro/src/psastroMosaicGetGrads.c

    r15671 r26259  
    2727
    2828                // select the raw objects for this readout
    29                 psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
     29                psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS.SUBSET");
    3030                if (rawstars == NULL) { continue; }
    3131
    3232                // select the raw objects for this readout
    33                 psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS");
     33                psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS.SUBSET");
    3434                if (refstars == NULL) { continue; }
    3535

  • trunk/psastro/src/psastroMosaicGradients.c

    r21409 r26259  
    4848
    4949                // select the raw objects for this readout
    50                 psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
     50                psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS.SUBSET");
    5151                if (rawstars == NULL) { continue; }
    5252
    5353                // select the raw objects for this readout
    54                 psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS");
     54                psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS.SUBSET");
    5555                if (refstars == NULL) { continue; }
    5656

  • trunk/psastro/src/psastroMosaicOneChip.c

    r24037 r26259  
    3131
    3232    // select the raw objects for this readout
    33     psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
     33    psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS.SUBSET");
    3434    if (rawstars == NULL) return false;
    3535
    36     psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS");
     36    psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS.SUBSET");
    3737    if (refstars == NULL) return false;
    3838
    … …  
    124124    float plateScale = 0.5*(fpa->toSky->Xs + fpa->toSky->Ys)*3600.0*PM_DEG_RAD;
    125125
     126    float rawXstdev = psStatsGetValue (results->xStats, psStatsStdevOption(results->xStats->options));
     127    float rawYstdev = psStatsGetValue (results->yStats, psStatsStdevOption(results->yStats->options));
     128
    126129    // pixError is the average 1D scatter in pixels ('results' are in FPA units = microns)
    127     float pixError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) / pixelScale;
     130    float pixError = 0.5*(rawXstdev + rawYstdev) / pixelScale;
    128131
    129132    // astError is the average 1D scatter in arcsec ('results' are in FPA units = microns)
    130     float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) * plateScale;
     133    float astError = 0.5*(rawXstdev + rawYstdev) * plateScale;
    131134    int astNstar = results->yStats->clippedNvalues;
    132135
    … …  
    161164    psMetadataAddF32 (updates, PS_LIST_TAIL, "EQUINOX",  PS_META_REPLACE, "", 2000.0); // XXX this is bogus: should be defined based on equinox of refstars
    162165
     166    // additional error measurements:
     167    psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_MDX",   PS_META_REPLACE, "mosaic astrometry X stdev (arcsec)", rawXstdev * plateScale);
     168    psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_MSX",   PS_META_REPLACE, "mosaic astrometry X systematic err (arcsec)", results->dXsys * plateScale);
     169    psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_MRX",   PS_META_REPLACE, "mosaic astrometry X 10-90 percentile (arcsec)", results->dXrange * plateScale);
     170
     171    psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_MDY",   PS_META_REPLACE, "mosaic astrometry Y stdev (arcsec)", rawYstdev * plateScale);
     172    psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_MSY",   PS_META_REPLACE, "mosaic astrometry Y systematic err (arcsec)", results->dYsys * plateScale);
     173    psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_MRY",   PS_META_REPLACE, "mosaic astrometry Y 10-90 percentile (arcsec)", results->dYrange * plateScale);
     174
    163175    // determine fromFPA transformation and apply new transformation to raw & ref stars
    164     psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
     176    psastroUpdateChipToFPA (fpa, chip);
    165177
    166178    //plot results

  • trunk/psastro/src/psastroMosaicSetMatch.c

    r21422 r26259  
    2626    double RADIUS = psMetadataLookupF32 (&status, recipe, radiusWord);
    2727    if (!status) {
    28         psError(PS_ERR_IO, false, "Failed to lookup matching radius: %s", radiusWord);
    29         psFree (view);
    30         return false;
     28        psAbort("Failed to lookup matching radius: %s", radiusWord);
     29    }
     30
     31    int uniqIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.MOSAIC.UNIQ.ITER");
     32    if (!status) {
     33        psAbort("Failed to lookup matching PSASTRO.MOSAIC.UNIQ.ITER");
    3134    }
    3235
    … …  
    5861
    5962                // select the raw objects for this readout
    60                 psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
     63                psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS.SUBSET");
    6164                if (rawstars == NULL) { continue; }
    6265
    6366                // select the raw objects for this readout
    64                 psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS");
     67                psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS.SUBSET");
    6568                if (refstars == NULL) { continue; }
    6669                psTrace ("psastro", 4, "Trying %ld refstars\n", refstars->n);
    … …  
    6871                psArray *matches = pmAstromRadiusMatchChip (rawstars, refstars, RADIUS);
    6972                psTrace ("psastro", 4, "Matched %ld refstars\n", matches->n);
     73
     74                if (iteration >= uniqIter) {
     75                    psArray *unique = pmAstromRadiusMatchUniq (rawstars, refstars, matches);
     76                    if (!unique) {
     77                        psLogMsg ("psastro", 3, "failed to generate a uniq set of matched sources\n");
     78                        return false;
     79                    }
     80                    psFree (matches);
     81                    matches = unique;
     82                }
    7083
    7184                pmAstromVisualPlotMosaicMatches(rawstars, refstars, matches, iteration, recipe);

  • trunk/psastro/src/psastroOneChipFit.c

    r24035 r26259  
    2828    REQUIRED_RECIPE_VALUE (int nIter, "PSASTRO.MATCH.FIT.NITER", S32);
    2929
     30    // for iterations >= uniqIter, require a single match per reference
     31    REQUIRED_RECIPE_VALUE (int uniqIter, "PSASTRO.MATCH.UNIQ.ITER", S32);
     32
    3033    // correct radius to FP units (physical pixel scale in microns per pixel)
    3134    REQUIRED_RECIPE_VALUE (double pixelScale, "PSASTRO.PIXEL.SCALE", F32);
    … …  
    6164            return false;
    6265        }
     66
     67        if (iter >= uniqIter) {
     68            psArray *unique = pmAstromRadiusMatchUniq (rawstars, refstars, match);
     69            if (!unique) {
     70                psLogMsg ("psastro", 3, "failed to generate a uniq set of matched sources\n");
     71                return false;
     72            }
     73            psFree (match);
     74            match = unique;
     75        }
    6376
    6477        // modify the order to correspond to the actual number of matched stars:
    … …  
    101114        }
    102115
    103         // XXX allow statitic to be set by the user
     116        // XXX allow statistic to be set by the user
    104117        // fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
    105118        fitStats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    … …  
    117130
    118131        // determine fromFPA transformation and apply new transformation to raw & ref stars
    119         psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
     132        // this applies the transformation to all stars (including subset)
     133        psastroUpdateChipToFPA (fpa, chip); // updates PSASTRO.RAWSTARS and PSASTRO.REFSTARS
    120134
    121135        // toSky converts from FPA & TPA units (microns) to sky units (radians)
    122136        float plateScale = 0.5*(fpa->toSky->Xs + fpa->toSky->Ys)*3600.0*PM_DEG_RAD;
    123         // float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) * plateScale;
    124         float astError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) * plateScale;
     137
     138        float rawXstdev = psStatsGetValue (results->xStats, psStatsStdevOption(results->xStats->options));
     139        float rawYstdev = psStatsGetValue (results->yStats, psStatsStdevOption(results->yStats->options));
     140
     141        float astError = 0.5*(rawXstdev + rawYstdev) * plateScale;
    125142        int astNstar = results->yStats->clippedNvalues;
    126143        psLogMsg ("psastro", PS_LOG_INFO, "pass %d, error: %f arcsec, Nstars: %d", iter, astError, astNstar);
    … …  
    136153    float plateScale = 0.5*(fpa->toSky->Xs + fpa->toSky->Ys)*3600.0*PM_DEG_RAD;
    137154
     155    float rawXstdev = psStatsGetValue (results->xStats, psStatsStdevOption(results->xStats->options));
     156    float rawYstdev = psStatsGetValue (results->yStats, psStatsStdevOption(results->yStats->options));
     157
    138158    // pixError is the average 1D scatter in pixels ('results' are in FPA units = microns)
    139     // float pixError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) / pixelScale;
    140     float pixError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) / pixelScale;
     159    float pixError = 0.5*(rawXstdev + rawYstdev) / pixelScale;
    141160
    142161    // astError is the average 1D scatter in arcsec ('results' are in FPA units = microns)
    143     // float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) * plateScale;
    144     float astError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) * plateScale;
     162    float astError = 0.5*(rawXstdev + rawYstdev) * plateScale;
     163
     164    // x and y are forced to use the same subset of values:
    145165    int astNstar = results->yStats->clippedNvalues;
    146166
    … …  
    170190    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
    171191
    172     // XXX drop from here : determine fromFPA transformation and apply new transformation to raw & ref stars
    173     // psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
     192    // additional error measurements:
     193    psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_CDX",   PS_META_REPLACE, "chip astrometry X stdev (arcsec)", rawXstdev * plateScale);
     194    psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_CSX",   PS_META_REPLACE, "chip astrometry X systematic err (arcsec)", results->dXsys * plateScale);
     195    psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_CRX",   PS_META_REPLACE, "chip astrometry X 10-90 percentile (arcsec)", results->dXrange * plateScale);
     196
     197    psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_CDY",   PS_META_REPLACE, "chip astrometry Y stdev (arcsec)", rawYstdev * plateScale);
     198    psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_CSY",   PS_META_REPLACE, "chip astrometry Y systematic err (arcsec)", results->dYsys * plateScale);
     199    psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_CRY",   PS_META_REPLACE, "chip astrometry Y 10-90 percentile (arcsec)", results->dYrange * plateScale);
    174200
    175201    // XXX check if we correctly applied the new transformation:
    … …  
    189215    psFree (results);
    190216    psFree (fitStats);
     217
    191218    return validSolution;
    192219}

  • trunk/psastro/src/psastroOneChipGrid.c

    r21409 r26259  
    3131
    3232    // generate the bright subset of maxNstar entries (note: rawstars is sorted by S/N)
    33     psArray *subset = psArrayAlloc (PS_MIN (maxNstar, rawstars->n));
     33    psArray *gridStars = psArrayAlloc (PS_MIN (maxNstar, rawstars->n));
    3434    for (int i = 0; (i < maxNstar) && (i < rawstars->n); i++) {
    35         subset->data[i] = psMemIncrRefCounter (rawstars->data[i]);
     35        gridStars->data[i] = psMemIncrRefCounter (rawstars->data[i]);
    3636    }
    3737
    38     // XXX set clump scale from recipe
    39     psArray *gridStars = psastroRemoveClumps (subset, 150);
    40     psFree (subset);
    41 
    42     psArray *refSubset = psastroRemoveClumps (refstars, 150);
    43 
    44     psLogMsg ("psastro", 3, "grid search using %ld raw vs %ld ref stars\n", gridStars->n, refSubset->n);
     38    psLogMsg ("psastro", 3, "grid search using %ld raw vs %ld ref stars\n", gridStars->n, refstars->n);
    4539
    4640    // find initial offset / rotation / scale
    47     pmAstromStats *gridStats = pmAstromGridMatch (gridStars, refSubset, recipe);
     41    pmAstromStats *gridStats = pmAstromGridMatch (gridStars, refstars, recipe);
    4842    if (gridStats == NULL) {
    4943        psLogMsg ("psastro", 3, "failed to find a grid match solution\n");
    5044        psFree (gridStars);
    51         psFree (refSubset);
    5245        return false;
    5346    }
    … …  
    5548
    5649    // tweak the position by finding peak of matches stars
    57     stats = pmAstromGridTweak (gridStars, refSubset, recipe, gridStats);
     50    stats = pmAstromGridTweak (gridStars, refstars, recipe, gridStats);
    5851    if (stats == NULL) {
    5952        psLogMsg ("psastro", 3, "failed to measure tweaked grid solution\n");
    6053        psFree (gridStats);
    6154        psFree (gridStars);
    62         psFree (refSubset);
    6355        return false;
    6456    }
    … …  
    6759    // adjust the chip.toFPA terms only
    6860    pmAstromGridApply (chip->toFPA, stats);
    69     psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
     61    psastroUpdateChipToFPA (fpa, chip); // updates PSASTRO.RAWSTARS and PSASTRO.REFSTARS
    7062    psFree (gridStats);
    7163    psFree (gridStars);
    72     psFree (refSubset);
    7364    psFree (stats);
    7465

  • trunk/psastro/src/psastroRefstarSubset.c

    r21409 r26259  
    1616
    1717  // select the raw objects for this readout
    18   psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
     18  psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS.SUBSET");
    1919  if (rawstars == NULL)  {
    2020    psError(PSASTRO_ERR_DATA, false, "missing rawstars in psastroRefstarSubset\n");
    … …  
    2323
    2424  // select the raw objects for this readout
    25   psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS");
     25  psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS.SUBSET");
    2626  if (refstars == NULL)  {
    2727    psError(PSASTRO_ERR_DATA, false, "missing refstars in psastroRefstarSubset\n");
    … …  
    7474  psLogMsg ("psastro", 4, "keeping %ld of %ld reference stars\n", subset->n, refstars->n);
    7575
    76   psMetadataRemoveKey (readout->analysis, "PSASTRO.REFSTARS");
    77   psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.REFSTARS", PS_DATA_ARRAY, "astrometry matches", subset);
     76  psMetadataRemoveKey (readout->analysis, "PSASTRO.REFSTARS.SUBSET");
     77  psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.REFSTARS.SUBSET", PS_DATA_ARRAY, "astrometry matches", subset);
    7878
    7979  if (psTraceGetLevel("psastro.dump") > 0) {

  • trunk/psastro/src/psastroRemoveClumps.c

    r21422 r26259  
    1212
    1313# include "psastroInternal.h"
     14
     15bool psastroRemoveClumpsRawstars (pmConfig *config) {
     16
     17    bool status;
     18
     19    pmChip *chip = NULL;
     20    pmCell *cell = NULL;
     21    pmReadout *readout = NULL;
     22
     23    // select the current recipe
     24    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSASTRO_RECIPE);
     25    if (!recipe) {
     26        psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe");
     27        return false;
     28    }
     29
     30    // select the input data sources
     31    pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PSASTRO.INPUT");
     32    if (!input) {
     33        psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
     34        return false;
     35    }
     36
     37    pmFPAview *view = pmFPAviewAlloc (0);
     38    pmFPA *fpa = input->fpa;
     39
     40    while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) {
     41        psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
     42        if (!chip->process || !chip->file_exists) { continue; }
     43
     44        while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) {
     45            psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
     46            if (!cell->process || !cell->file_exists) { continue; }
     47            if (!chip->fromFPA) { continue; }
     48
     49            // process each of the readouts
     50            while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) {
     51                if (! readout->data_exists) { continue; }
     52
     53                // select the raw objects for this readout
     54                psArray *rawstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.RAWSTARS");
     55                if (rawstars == NULL) { continue; }
     56
     57                // generate a reduced subset excluding the clumps
     58                // XXX do we need both RAWSTARS and SUBSET?
     59                // XXX put these parameters in the recipe, please
     60                psArray *subset = psastroRemoveClumpsIterate(rawstars, 150, 3);
     61                psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.RAWSTARS.SUBSET", PS_DATA_ARRAY, "astrometry objects", subset);
     62                psFree (subset);
     63            }
     64        }
     65    }
     66    psFree (view);
     67    return true;
     68}
     69
     70psArray *psastroRemoveClumpsIterate (psArray *input, int scale, int nIter) {
     71
     72    psArray *newset = psMemIncrRefCounter (input);
     73    for (int i = 0; i < nIter; i++) {
     74        psArray *subset = psastroRemoveClumps (newset, 150);
     75        psFree (newset);
     76        newset = subset;
     77    }
     78    return newset;
     79}   
    1480
    1581/**

  • trunk/psastro/src/psastroUtils.c

    r21422 r26259  
    106106}
    107107
    108 bool psastroUpdateChipToFPA (pmFPA *fpa, pmChip *chip, psArray *rawstars, psArray *refstars) {
     108bool psastroUpdateChipToFPA (pmFPA *fpa, pmChip *chip) {
    109109
    110110    psRegion *region = pmChipPixels (chip);
    … …  
    114114    psFree (region);
    115115
    116     for (int i = 0; i < rawstars->n; i++) {
    117         pmAstromObj *raw = rawstars->data[i];
    118 
    119         psPlaneTransformApply (raw->FP, chip->toFPA, raw->chip);
    120         psPlaneTransformApply (raw->TP, fpa->toTPA, raw->FP);
    121         psDeproject (raw->sky, raw->TP, fpa->toSky);
    122     }
    123 
    124     for (int i = 0; i < refstars->n; i++) {
    125         pmAstromObj *ref = refstars->data[i];
    126         psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP);
     116    // loop over cells in this chip
     117    for (int nCell = 0; nCell < chip->cells->n; nCell++) {
     118        pmCell *cell = chip->cells->data[nCell];
     119        if (!cell->process || !cell->file_exists) { continue; }
     120
     121        // loop over readouts in this cell
     122        for (int nRead = 0; nRead < cell->readouts->n; nRead++) {
     123            pmReadout *readout = cell->readouts->data[nRead];
     124            if (! readout->data_exists) { continue; }
     125
     126            // select the raw objects for this readout
     127            psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
     128            if (rawstars) {
     129                for (int i = 0; i < rawstars->n; i++) {
     130                    pmAstromObj *raw = rawstars->data[i];
     131                    psPlaneTransformApply (raw->FP, chip->toFPA, raw->chip);
     132                    psPlaneTransformApply (raw->TP, fpa->toTPA, raw->FP);
     133                    psDeproject (raw->sky, raw->TP, fpa->toSky);
     134                }
     135            }
     136
     137            psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS");
     138            if (refstars) {
     139                for (int i = 0; i < refstars->n; i++) {
     140                    pmAstromObj *ref = refstars->data[i];
     141                    psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP);
     142                }
     143            }
     144        }
    127145    }
    128146

  • trunk/psastro/src/psastroZeroPoint.c

    r25906 r26259  
    3030    // recipe options
    3131    bool byExposure = psMetadataLookupBool (&status, recipe, "ZERO.POINT.BY.EXPOSURE");
    32     bool useMean    = psMetadataLookupBool (&status, recipe, "ZERO.POINT.USE.MEAN");
    33     float edgeFraction = psMetadataLookupF32 (&status, recipe, "ZERO.POINT.EDGE.FRACTION");
    3432
    3533    // select the input data sources
    … …  
    8684                    // calculate dMag for the matched stars just for this readout (well, chip)
    8785                    psMetadata *header = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER");
    88                     psastroZeroPointAnalysis (header, dMag, zeropt, edgeFraction, useMean);
     86                    psastroZeroPointAnalysis (header, dMag, zeropt, recipe);
    8987                    psFree (dMag);
    9088                    dMag = NULL;
    … …  
    9694    if (byExposure) {
    9795        psMetadata *header = psMetadataLookupMetadata (&status, fpa->analysis, "PSASTRO.HEADER");
    98         psastroZeroPointAnalysis (header, dMag, zeropt, edgeFraction, useMean);
     96        if (!header) {
     97          header = psMetadataAlloc ();
     98          psMetadataAddMetadata (fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", header);
     99          psFree (header);
     100        }
     101        psastroZeroPointAnalysis (header, dMag, zeropt, recipe);
    99102        psFree (dMag);
    100103        dMag = NULL;
    … …  
    113116
    114117    // select the raw objects for this readout
    115     psArray *rawstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.RAWSTARS");
     118    psArray *rawstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.RAWSTARS.SUBSET");
    116119    if (rawstars == NULL) return dMag;
    117120
    118121    // select the raw objects for this readout
    119     psArray *refstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.REFSTARS");
     122    psArray *refstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.REFSTARS.SUBSET");
    120123    if (refstars == NULL) return dMag;
    121124
    … …  
    141144}
    142145
    143 bool psastroZeroPointAnalysis (psMetadata *header, psVector *dMag, float zeropt, float edgeFraction, bool useMean) {
     146bool psastroZeroPointAnalysis (psMetadata *header, psVector *dMag, float zeropt, psMetadata *recipe) {
    144147
    145148    // XXX make this depend on the mode?
    … …  
    148151      return false;
    149152    }
     153
     154    bool status;
     155    bool useMean = psMetadataLookupBool (&status, recipe, "ZERO.POINT.USE.MEAN");
    150156
    151157    // the zero point analysis depends on the type of desired statistic.  For comparisons
    … …  
    166172        }
    167173    } else {
    168         stats = psastroStatsPercentile (dMag, edgeFraction);
     174        stats = psastroStatsPercentile (dMag, recipe);
    169175        if (!stats) {
    170176            psError(PS_ERR_UNKNOWN, false, "failure to measure zero point by edge");
    … …  
    183189}
    184190
    185 #define MAG_RESOLUTION 0.001
     191#define MAG_RESOLUTION 0.0002
    186192
    187193// set the bin closest to the corresponding value.  if USE_END is +/- 1,
    … …  
    233239// return results in stats->sampleMean, sampleStdev
    234240// XXX this is a misuse of psStats -- make our own structure?
    235 psStats *psastroStatsPercentile (psVector *myVector, float flimit) {
     241psStats *psastroStatsPercentile (psVector *myVector, psMetadata *recipe) {
    236242
    237243    // search for the 'blue' edge of the dMag distribution:
    238244    // the distribution is not a normal population, but instead has a broad range with fairly hard edges.
    239245    // construct a histogram and look for the
     246
     247    bool status;
     248    float edgeFraction = psMetadataLookupF32 (&status, recipe, "ZERO.POINT.EDGE.FRACTION");
     249    int edgeSample = psMetadataLookupS32 (&status, recipe, "ZERO.POINT.EDGE.SAMPLE");
     250    float edgeSampleFraction = psMetadataLookupF32 (&status, recipe, "ZERO.POINT.EDGE.SAMPLE.FRACTION");
    240251
    241252    // stats is first used to find the data range
    … …  
    245256    float min = NAN, max = NAN;         // Mimimum and maximum values
    246257
     258    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); 
     259
    247260    // Get the minimum and maximum values
    248261    if (!psVectorStats(stats, myVector, NULL, NULL, 0)) goto escape;
    … …  
    250263    max = stats->max;
    251264    if (isnan(min) || isnan(max)) goto escape;
    252     psTrace("psastro", 6, "Data min/max is (%.2f, %.2f)\n", min, max);
     265    psTrace("psastro", 5, "Data min/max is (%.2f, %.2f)\n", min, max);
    253266
    254267    // If all data points have the same value, then we set the appropriate members of stats and return.
    … …  
    264277    float binSize = MAG_RESOLUTION;
    265278    long numBins = (max - min) / binSize; // Number of bins
    266     psTrace("psastro", 6, "Numbins is %ld\n", numBins);
    267     psTrace("psastro", 6, "Creating a robust histogram from data range (%.2f - %.2f)\n", min, max);
    268 
    269     // Generate the histogram
     279    psTrace("psastro", 5, "Numbins is %ld\n", numBins);
     280    psTrace("psastro", 5, "Creating a robust histogram from data range (%.2f - %.2f)\n", min, max);
     281
     282    // allocate the histogram containers
    270283    histogram = psHistogramAlloc(min, max, numBins);
    271     if (!psVectorHistogram(histogram, myVector, NULL, NULL, 0)) {
    272       // if psVectorHistogram returns false, we have a programming error
    273       psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram for psastroZeroPointAnalysis\n");
    274       psFree(histogram);
    275       psFree(stats);
    276       return NULL;
    277     }
    278     if (psTraceGetLevel("psastro") >= 8) {
    279       PS_VECTOR_PRINT_F32(histogram->bounds);
    280       PS_VECTOR_PRINT_F32(histogram->nums);
    281     }
    282 
    283     // Convert the specific histogram to a cumulative histogram
    284     // The cumulative histogram data points correspond to the UPPER bound value (N < Bin[i+1])
    285284    cumulative = psHistogramAlloc(min, max, numBins);
    286     cumulative->nums->data.F32[0] = histogram->nums->data.F32[0];
    287     for (long i = 1; i < histogram->nums->n; i++) {
    288       cumulative->nums->data.F32[i] = cumulative->nums->data.F32[i-1] + histogram->nums->data.F32[i];
    289       cumulative->bounds->data.F32[i-1] = histogram->bounds->data.F32[i];
    290     }
    291     if (psTraceGetLevel("psastro") >= 8) {
    292       PS_VECTOR_PRINT_F32(cumulative->bounds);
    293       PS_VECTOR_PRINT_F32(cumulative->nums);
    294     }
    295 
    296     // Find the bin which contains the first data point above the limit
    297     long totalDataPoints = cumulative->nums->data.F32[numBins - 1];
    298     psTrace("psastro", 6, "Total data points is %ld\n", totalDataPoints);
    299 
    300     // find bin which is the lower bound of the limit value (value[bin] < f < value[bin+1]
    301     long bin;
    302     PS_BIN_FOR_VALUE(bin, cumulative->nums, flimit * totalDataPoints, 0);
    303     psTrace("psastro", 6, "The bin is %ld (%.2f to %.2f)\n", bin, cumulative->bounds->data.F32[bin], cumulative->bounds->data.F32[bin+1]);
    304 
    305     // Linear interpolation to the limit value in bin units
    306     float value;
    307     PS_BIN_INTERPOLATE (value, cumulative->nums, cumulative->bounds, bin, totalDataPoints * flimit);
    308     psTrace("psastro", 6, "limit value is %f\n", value);
    309 
    310     stats->clippedMean = value;
    311     stats->clippedStdev = 0.0; // XXX derive correct error value
     285
     286    // find the mean value:
     287    stats->clippedMean = psastroStatsPercentileValue (histogram, cumulative, myVector, edgeFraction);
     288
     289    int nSubset = myVector->n * edgeSampleFraction;
     290    psVector *subset = psVectorAlloc (nSubset, PS_TYPE_F32);
     291
     292    float Sum = 0.0;
     293    float S2 = 0.0;
     294    for (int i = 0; i < edgeSample; i++) {
     295       
     296        // generate the subset vector
     297        for (long i = 0; i < nSubset; i++) {
     298            double frnd = psRandomUniform(rng);
     299            int entry = PS_MIN(myVector->n - 1, PS_MAX(0, myVector->n * frnd));
     300
     301            subset->data.F32[i] = myVector->data.F32[entry];
     302        }
     303
     304        float value = psastroStatsPercentileValue (histogram, cumulative, subset, edgeFraction);
     305
     306        Sum += value;
     307        S2 += value*value;
     308    }
     309    psTrace("psastro", 6, "subset stats: Sum: %f, S2: %f, Npts: %d\n", Sum, S2, edgeSample);
     310
     311    stats->clippedStdev = PS_MAX (sqrt(S2 / edgeSample - PS_SQR(Sum/edgeSample)), MAG_RESOLUTION);
     312    psTrace("psastro", 5, "percentile stats %f +/- %f\n", stats->clippedMean, stats->clippedStdev);
     313
    312314    stats->results |= PS_STAT_CLIPPED_MEAN;
    313315    stats->results |= PS_STAT_CLIPPED_STDEV;
    … …  
    317319    psFree(histogram);
    318320    psFree(cumulative);
     321    psFree(subset);
     322    psFree(rng);
    319323    return stats;
    320324
    … …  
    332336}
    333337
     338
     339// measure the edge of the sample at flimit
     340// return results in stats->sampleMean, sampleStdev
     341// XXX this is a misuse of psStats -- make our own structure?
     342float psastroStatsPercentileValue (psHistogram *histogram, psHistogram *cumulative, psVector *myVector, float flimit) {
     343
     344    // need to initialize the histogram on each pass
     345    psVectorInit (histogram->nums, 0);
     346    if (!psVectorHistogram(histogram, myVector, NULL, NULL, 0)) {
     347        // if psVectorHistogram returns false, we have a programming error
     348        psAbort ("Unable to generate histogram for psastroZeroPointAnalysis");
     349    }
     350    if (psTraceGetLevel("psastro") >= 8) {
     351        PS_VECTOR_PRINT_F32(histogram->bounds);
     352        PS_VECTOR_PRINT_F32(histogram->nums);
     353    }
     354
     355    // Convert the specific histogram to a cumulative histogram
     356    // The cumulative histogram data points correspond to the UPPER bound value (N < Bin[i+1])
     357    cumulative->nums->data.F32[0] = histogram->nums->data.F32[0];
     358    for (long i = 1; i < histogram->nums->n; i++) {
     359        cumulative->nums->data.F32[i] = cumulative->nums->data.F32[i-1] + histogram->nums->data.F32[i];
     360        cumulative->bounds->data.F32[i-1] = histogram->bounds->data.F32[i];
     361    }
     362    if (psTraceGetLevel("psastro") >= 8) {
     363        PS_VECTOR_PRINT_F32(cumulative->bounds);
     364        PS_VECTOR_PRINT_F32(cumulative->nums);
     365    }
     366
     367    // Find the bin which contains the first data point above the limit
     368    long numBins = cumulative->nums->n;
     369    long totalDataPoints = cumulative->nums->data.F32[numBins - 1];
     370    psTrace("psastro", 6, "Total data points is %ld\n", totalDataPoints);
     371
     372    // find bin which is the lower bound of the limit value (value[bin] < f < value[bin+1]
     373    long bin;
     374    PS_BIN_FOR_VALUE(bin, cumulative->nums, flimit * totalDataPoints, 0);
     375    psTrace("psastro", 6, "The bin is %ld (%.4f to %.4f)\n", bin, cumulative->bounds->data.F32[bin], cumulative->bounds->data.F32[bin+1]);
     376
     377    // Linear interpolation to the limit value in bin units
     378    float value;
     379    PS_BIN_INTERPOLATE (value, cumulative->nums, cumulative->bounds, bin, totalDataPoints * flimit);
     380    psTrace("psastro", 6, "limit value is %f\n", value);
     381
     382    return value;
     383}
    334384
    335385# define ESCAPE(MSG) { \
Note: See TracChangeset for help on using the changeset viewer.