Index: trunk/psastro/doc/outline.txt =================================================================== --- trunk/psastro/doc/outline.txt (revision 26259) +++ trunk/psastro/doc/outline.txt (revision 26259) @@ -0,0 +1,45 @@ + +psastro outline: + + psastroParseCamera + psastroDataLoad + psastroConvertReadout (save PSASTRO.RAWSTARS & PSASTRO.RAWSTARS.SUBSET on readout->analysis) + psastroAnalysis + psastroUseModel + psastroAstromGuess (plots PSASTRO.RAWSTARS) + psastroLoadRefstars (loads full reference star set) + psastroChooseRefstars (save PSASTRO.REFSTARS & PSASTRO.REFSTARS.SUBSET on readout->analysis) + psastroRefstarSubset (adjust PSASTRO.REFSTARS.SUBSET on readout->analysis) + psastroChooseGlintStars (uses full reference star set) + psastroChipAstrom (uses PSASTRO.REFSTARS.SUBSET & PSASTRO.RAWSTARS.SUBSET) + psastroOneChipGrid + psastroOneChipFit + psastroMosaicAstrom + psastroFixChips + psastroFixChipsTest + psastroMosaicFit + psastroMosaicSetMatch (uses PSASTRO.REFSTARS.SUBSET & PSASTRO.RAWSTARS.SUBSET) + psastroMosaicCommonScale + psastroMosaicDistortion + psastroMosaicChipAstrom + psastroZeroPoint + psastroAstromGuessCheck + psastroMaskUpdates + psastroDataSave + +psastroExtract outline: + psastroExtractParseCamera + psastroExtractDataLoad + psastroExtractAnalysis + psastroLoadRefstars (loads full reference star set) + psastroChooseRefstars (save PSASTRO.REFSTARS & PSASTRO.REFSTARS.SUBSET on readout->analysis) + psastroRefstarSubset (adjust PSASTRO.REFSTARS.SUBSET on readout->analysis) + psastroExtractGhosts + psastroExtractStars + psastroDataSave + + +psastroModel outline: + +psastroModelFit outline: + Index: trunk/psastro/src/psastro.h =================================================================== --- trunk/psastro/src/psastro.h (revision 26258) +++ trunk/psastro/src/psastro.h (revision 26259) @@ -88,8 +88,10 @@ bool psastroLuminosityFunctionPlot(psVector *lnMag, psVector *Mag, pmLumFunc *lumFunc, pmLumFunc *rawFunc); psArray *psastroRemoveClumps (psArray *input, int scale); +psArray *psastroRemoveClumpsIterate (psArray *input, int scale, int nIter); +bool psastroRemoveClumpsRawstars (pmConfig *config); // utility functions: -bool psastroUpdateChipToFPA (pmFPA *fpa, pmChip *chip, psArray *rawstars, psArray *refstars); +bool psastroUpdateChipToFPA (pmFPA *fpa, pmChip *chip); bool psastroWriteStars (char *filename, psArray *sources); bool psastroWriteTransform (psPlaneTransform *map); @@ -149,8 +151,9 @@ bool psastroZeroPoint (pmConfig *config); psVector *psastroZeroPointReadoutAccum(psVector *dMag, pmReadout *readout, float exptime); -bool psastroZeroPointAnalysis (psMetadata *header, psVector *dMag, float zeropt, float edgeFraction, bool useMean); +bool psastroZeroPointAnalysis (psMetadata *header, psVector *dMag, float zeropt, psMetadata *recipe); bool psastroZeroPointFromRecipe (float *zeropt, float *exptime, float *ghostMaxMag, pmFPA *fpa, psMetadata *recipe); -psStats *psastroStatsPercentile (psVector *myVector, float flimit); +psStats *psastroStatsPercentile (psVector *myVector, psMetadata *recipe); +float psastroStatsPercentileValue (psHistogram *histogram, psHistogram *cumulative, psVector *myVector, float flimit); // masking functions Index: trunk/psastro/src/psastroAnalysis.c =================================================================== --- trunk/psastro/src/psastroAnalysis.c (revision 26258) +++ trunk/psastro/src/psastroAnalysis.c (revision 26259) @@ -69,4 +69,9 @@ } + if (!psastroRemoveClumpsRawstars(config)) { + psError (PSASTRO_ERR_UNKNOWN, false, "failed to remove RAWSTAR clumps\n"); + return false; + } + // load the reference stars overlapping the data stars psArray *refs = psastroLoadRefstars(config, "PSASTRO.INPUT"); @@ -127,7 +132,4 @@ psastroMaskUpdates (config); - // XXX how do we specify stack astrometry? - // psastroStackAstrom (config, refs); - return true; } Index: trunk/psastro/src/psastroArguments.c =================================================================== --- trunk/psastro/src/psastroArguments.c (revision 26258) +++ trunk/psastro/src/psastroArguments.c (revision 26259) @@ -62,8 +62,4 @@ // define the reference astrometry file status = pmConfigFileSetsMD (config->arguments, &argc, argv, "ASTROM.MODEL", "-astrommodel", "-astrommodellist"); - if (status) { - // if supplied, assume -fixchips is desired - psMetadataAddBool (config->arguments, PS_LIST_TAIL, "PSASTRO.FIX.CHIPS", PS_META_REPLACE, "", true); - } // define the reference astrometry file Index: trunk/psastro/src/psastroChipAstrom.c =================================================================== --- trunk/psastro/src/psastroChipAstrom.c (revision 26258) +++ trunk/psastro/src/psastroChipAstrom.c (revision 26259) @@ -56,9 +56,9 @@ // select the raw objects for this readout - psArray *rawstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.RAWSTARS"); + psArray *rawstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.RAWSTARS.SUBSET"); if (rawstars == NULL) { continue; } // select the raw objects for this readout - psArray *refstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.REFSTARS"); + psArray *refstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.REFSTARS.SUBSET"); if (refstars == NULL) { continue; } @@ -66,6 +66,5 @@ if ((rawstars->n < 4) || (refstars->n < 4)) { readout->data_exists = false; - psLogMsg ("psastro", 3, "insufficient rawstars (%ld) or refstars (%ld)", - rawstars->n, refstars->n); + psLogMsg ("psastro", 3, "insufficient rawstars (%ld) or refstars (%ld)", rawstars->n, refstars->n); continue; } Index: trunk/psastro/src/psastroChooseRefstars.c =================================================================== --- trunk/psastro/src/psastroChooseRefstars.c (revision 26258) +++ trunk/psastro/src/psastroChooseRefstars.c (revision 26259) @@ -134,9 +134,17 @@ psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.REFSTARS", PS_DATA_ARRAY, "astrometry matches", refstars); + + // generate a reduced subset excluding the clumps + // XXX do we need both REFSTARS and SUBSET? + psArray *subset = psastroRemoveClumpsIterate(refstars, 150, 3); + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.REFSTARS.SUBSET", PS_DATA_ARRAY, "astrometry objects", subset); + psFree (refstars); + psFree (subset); psFree (extent); if (matchLumFunc) { - // in this case, no PSASTRO.REFSTARS is added to readout->analysis + // limit the total magnitude range of PSASTRO.REFSTARS.SUBSET based on + // overlapping luminosity functions if (!psastroRefstarSubset (readout)) { psError(PSASTRO_ERR_DATA, false, "Can't determine an appropriate refstar subset\n"); Index: trunk/psastro/src/psastroConvert.c =================================================================== --- trunk/psastro/src/psastroConvert.c (revision 26258) +++ trunk/psastro/src/psastroConvert.c (revision 26259) @@ -127,4 +127,5 @@ psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.RAWSTARS", PS_DATA_ARRAY, "astrometry objects", rawStars); + 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); psLogMsg ("psastro", 4, "skip reasons: mode: %d, faint: %d, bright: %d, inf: %d\n", nModeSkip, nFaintSkip, nBrightSkip, nInfSkip); @@ -153,11 +154,15 @@ psF32 *PAR = model->params->data.F32; + psF32 *dPAR = model->dparams->data.F32; pmAstromObj *obj = pmAstromObjAlloc (); // is the source magnitude calibrated in any sense? - obj->pix->x = PAR[PM_PAR_XPOS]; - obj->pix->y = PAR[PM_PAR_YPOS]; + obj->pix->x = PAR[PM_PAR_XPOS]; + obj->pix->y = PAR[PM_PAR_YPOS]; + obj->pix->xErr = dPAR[PM_PAR_XPOS]; + obj->pix->yErr = dPAR[PM_PAR_YPOS]; obj->Mag = source->psfMag; + obj->dMag = source->errMag; // XXX do we have the information giving the readout and cell offset? Index: trunk/psastro/src/psastroDemoDump.c =================================================================== --- trunk/psastro/src/psastroDemoDump.c (revision 26258) +++ trunk/psastro/src/psastroDemoDump.c (revision 26259) @@ -160,9 +160,9 @@ // select the raw objects for this readout - psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS"); + psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS.SUBSET"); if (rawstars == NULL) continue; // select the raw objects for this readout - psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS"); + psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS.SUBSET"); if (refstars == NULL) continue; psTrace ("psastro", 4, "Trying %ld refstars\n", refstars->n); Index: trunk/psastro/src/psastroFixChips.c =================================================================== --- trunk/psastro/src/psastroFixChips.c (revision 26258) +++ trunk/psastro/src/psastroFixChips.c (revision 26259) @@ -248,5 +248,10 @@ badAstrom |= fabs(obsAngle - refAngle) > angleTol; - fprintf (stderr, "chip %d, angle: %f, pixel: %f,%f\n", view->chip, obsAngle - refAngle, obsCoord.x - refCoord.x, obsCoord.y - refCoord.y); + fprintf (stderr, "chip %d, angle: %f, pixel: %f,%f : ", view->chip, obsAngle - refAngle, obsCoord.x - refCoord.x, obsCoord.y - refCoord.y); + if (badAstrom) { + fprintf (stderr, "BAD ASTROM\n"); + } else { + fprintf (stderr, "GOOD ASTROM\n"); + } // XXX for now, just use first readout @@ -265,4 +270,6 @@ psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DY", PS_META_REPLACE, "chip y offset wrt model", obsCoord.y - refCoord.y); psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_DT", PS_META_REPLACE, "chip rot offset wrt model", obsAngle - refAngle); + + continue; // for successful chips, save the measured offsets in the header @@ -319,5 +326,5 @@ } - psastroUpdateChipToFPA (input->fpa, obsChip, rawstars, refstars); + psastroUpdateChipToFPA (input->fpa, obsChip); // XXX update the header with info to reflect the failure Index: trunk/psastro/src/psastroFixChipsTest.c =================================================================== --- trunk/psastro/src/psastroFixChipsTest.c (revision 26258) +++ trunk/psastro/src/psastroFixChipsTest.c (revision 26259) @@ -246,5 +246,5 @@ } - psastroUpdateChipToFPA (fpa, chip, rawstars, refstars); + psastroUpdateChipToFPA (fpa, chip); // XXX skip the re-fit step for a test Index: trunk/psastro/src/psastroLoadGhosts.c =================================================================== --- trunk/psastro/src/psastroLoadGhosts.c (revision 26258) +++ trunk/psastro/src/psastroLoadGhosts.c (revision 26259) @@ -119,5 +119,7 @@ if (! readout->data_exists) { continue; } - // select the raw objects for this readout (loaded in psastroExtract.c) + // select the raw objects for this readout (loaded in psastroChooseRefstars.c) + // XXX : note that we place limits on the refstar sample in psastroChooseRefstars.c: + // 1) on chip and 2) < PSASTRO.MAX.NREF. magnitude limits and clump exclusion are only psArray *refstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.REFSTARS"); if (refstars == NULL) { continue; } Index: trunk/psastro/src/psastroLoadRefstars.c =================================================================== --- trunk/psastro/src/psastroLoadRefstars.c (revision 26258) +++ trunk/psastro/src/psastroLoadRefstars.c (revision 26259) @@ -50,5 +50,30 @@ return false; } + + // the name in the recipe may be one of: + // (A) the actual directory name + // (B) a reference to the name in the PSASTRO.CATDIRS folder in site.config + // (C) a reference to a folder in the PSASTRO.CATDIRS folder in site.config, containing multiple copy locations char *catdir_virtual = psMetadataLookupStr(&status, catdirs, catdir_recipe); + + psMetadata *catdir_folder = psMetadataLookupMetadata(&status, catdirs, catdir_recipe); + if (catdir_folder) { + // randomly choose one of the entries + psLogMsg ("psastro", 3, "choosing catdir_folder\n"); + int nEntry = catdir_folder->list->n; + + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); + double frnd = psRandomUniform(rng); + int entry = PS_MIN(nEntry - 1, PS_MAX(0, nEntry * frnd)); + psFree(rng); + + psMetadataItem *item = psListGet(catdir_folder->list, entry); + if (item->type != PS_DATA_STRING) { + psError(PSASTRO_ERR_CONFIG, true, "Invalid entry in PSASTRO.CATDIR folder: %s\n", item->name); + return false; + } + catdir_virtual = item->data.str; + } + char *catdir = (catdir_virtual == NULL) ? catdir_recipe : catdir_virtual; @@ -56,4 +81,6 @@ psString CATDIR = pmConfigConvertFilename(catdir, config, false, false); // Resolved filename PS_ASSERT (CATDIR, NULL); + + psLogMsg ("psastro", 3, "looking up reference objects in %s\n", CATDIR); char *getstarCommand = psStringCopy(psMetadataLookupStr(NULL, recipe, "DVO.GETSTAR")); Index: trunk/psastro/src/psastroMaskUpdates.c =================================================================== --- trunk/psastro/src/psastroMaskUpdates.c (revision 26258) +++ trunk/psastro/src/psastroMaskUpdates.c (revision 26259) @@ -193,4 +193,7 @@ // select the raw objects for this readout + // XXX : note that we place limits on the refstar sample in psastroChooseRefstars.c: + // 1) on chip and 2) < PSASTRO.MAX.NREF. magnitude limits and clump exclusion are only + // applied to the SUBSETs psArray *refstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.REFSTARS"); if (!refstars) continue; Index: trunk/psastro/src/psastroMosaicAstrom.c =================================================================== --- trunk/psastro/src/psastroMosaicAstrom.c (revision 26258) +++ trunk/psastro/src/psastroMosaicAstrom.c (revision 26259) @@ -20,5 +20,4 @@ bool status; - char filename[256]; // select the current recipe @@ -38,37 +37,30 @@ pmFPA *fpa = input->fpa; - // before we do object matches, we need to (optionally) fix failed chips. We compare chips with - // the supplied mosaic model. Adjust significant outliers to match model. - # if (0) - if (!psastroFixChips (config, recipe)) { - psError(PSASTRO_ERR_UNKNOWN, false, "failed to align problematic chips"); - return false; - } - if (!psastroFixChipsTest (config, recipe)) { - psError(PSASTRO_ERR_UNKNOWN, false, "failed to align problematic chips"); - return false; - } - # endif - char *outroot = psMetadataLookupStr (&status, config->arguments, "OUTPUT"); if (!status || !outroot) psAbort ("Can't find outroot on config->arguments"); - if (!psastroMosaicFit (fpa, recipe, outroot, 0)) return false; - if (!psastroMosaicFit (fpa, recipe, outroot, 1)) return false; - if (!psastroMosaicFit (fpa, recipe, outroot, 2)) return false; - if (!psastroMosaicFit (fpa, recipe, outroot, 3)) return false; + int nIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.NITER"); + if (!status) psAbort ("missing config value"); + + // this should be in a loop with nIter = + for (int iter = 0; iter < nIter; iter++) { + if (!psastroMosaicFit (fpa, recipe, outroot, iter)) return false; + } // now fit the chips under the common distortion with higher-order terms // first, re-perform the match with a slightly tighter circle - if (!psastroMosaicSetMatch (fpa, recipe, 4)) { + if (!psastroMosaicSetMatch (fpa, recipe, nIter)) { psError(PSASTRO_ERR_UNKNOWN, false, "failed to match raw and ref stars for mosaic (4th pass)"); return false; } - if (!psastroMosaicChipAstrom (fpa, recipe, 4)) { + if (!psastroMosaicChipAstrom (fpa, recipe, nIter)) { psError(PSASTRO_ERR_UNKNOWN, false, "failed to measure chip astrometry in mosaic mode (4th pass)"); return false; } + if (psTraceGetLevel("psastro.dump") > 0) { - snprintf (filename, 256, "%s.10.dat", outroot); + // the last filename (see filenames in psastroMosaicFit) + char filename[256]; + snprintf (filename, 256, "%s.%d.dat", outroot, 2*nIter + 2); psastroDumpMatches (fpa, filename); } Index: trunk/psastro/src/psastroMosaicFPtoTP.c =================================================================== --- trunk/psastro/src/psastroMosaicFPtoTP.c (revision 26258) +++ trunk/psastro/src/psastroMosaicFPtoTP.c (revision 26259) @@ -45,9 +45,9 @@ // select the raw objects for this readout - psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS"); + psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS.SUBSET"); if (rawstars == NULL) { continue; } // select the raw objects for this readout - psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS"); + psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS.SUBSET"); if (refstars == NULL) { continue; } @@ -126,4 +126,6 @@ } +// apply the rotation and scale to all stars in PSASTRO.REFSTARS (also adjusts +// PSASTRO.REFSTARS.SUBSET since they are the same pointers) bool psastroMosaicApplyRotAndScale (pmFPA *fpa, psPlaneTransform *TPtoFP) { Index: trunk/psastro/src/psastroMosaicGetGrads.c =================================================================== --- trunk/psastro/src/psastroMosaicGetGrads.c (revision 26258) +++ trunk/psastro/src/psastroMosaicGetGrads.c (revision 26259) @@ -27,9 +27,9 @@ // select the raw objects for this readout - psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS"); + psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS.SUBSET"); if (rawstars == NULL) { continue; } // select the raw objects for this readout - psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS"); + psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS.SUBSET"); if (refstars == NULL) { continue; } Index: trunk/psastro/src/psastroMosaicGradients.c =================================================================== --- trunk/psastro/src/psastroMosaicGradients.c (revision 26258) +++ trunk/psastro/src/psastroMosaicGradients.c (revision 26259) @@ -48,9 +48,9 @@ // select the raw objects for this readout - psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS"); + psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS.SUBSET"); if (rawstars == NULL) { continue; } // select the raw objects for this readout - psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS"); + psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS.SUBSET"); if (refstars == NULL) { continue; } Index: trunk/psastro/src/psastroMosaicOneChip.c =================================================================== --- trunk/psastro/src/psastroMosaicOneChip.c (revision 26258) +++ trunk/psastro/src/psastroMosaicOneChip.c (revision 26259) @@ -31,8 +31,8 @@ // select the raw objects for this readout - psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS"); + psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS.SUBSET"); if (rawstars == NULL) return false; - psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS"); + psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS.SUBSET"); if (refstars == NULL) return false; @@ -124,9 +124,12 @@ float plateScale = 0.5*(fpa->toSky->Xs + fpa->toSky->Ys)*3600.0*PM_DEG_RAD; + float rawXstdev = psStatsGetValue (results->xStats, psStatsStdevOption(results->xStats->options)); + float rawYstdev = psStatsGetValue (results->yStats, psStatsStdevOption(results->yStats->options)); + // pixError is the average 1D scatter in pixels ('results' are in FPA units = microns) - float pixError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) / pixelScale; + float pixError = 0.5*(rawXstdev + rawYstdev) / pixelScale; // astError is the average 1D scatter in arcsec ('results' are in FPA units = microns) - float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) * plateScale; + float astError = 0.5*(rawXstdev + rawYstdev) * plateScale; int astNstar = results->yStats->clippedNvalues; @@ -161,6 +164,15 @@ psMetadataAddF32 (updates, PS_LIST_TAIL, "EQUINOX", PS_META_REPLACE, "", 2000.0); // XXX this is bogus: should be defined based on equinox of refstars + // additional error measurements: + psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_MDX", PS_META_REPLACE, "mosaic astrometry X stdev (arcsec)", rawXstdev * plateScale); + psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_MSX", PS_META_REPLACE, "mosaic astrometry X systematic err (arcsec)", results->dXsys * plateScale); + psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_MRX", PS_META_REPLACE, "mosaic astrometry X 10-90 percentile (arcsec)", results->dXrange * plateScale); + + psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_MDY", PS_META_REPLACE, "mosaic astrometry Y stdev (arcsec)", rawYstdev * plateScale); + psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_MSY", PS_META_REPLACE, "mosaic astrometry Y systematic err (arcsec)", results->dYsys * plateScale); + psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_MRY", PS_META_REPLACE, "mosaic astrometry Y 10-90 percentile (arcsec)", results->dYrange * plateScale); + // determine fromFPA transformation and apply new transformation to raw & ref stars - psastroUpdateChipToFPA (fpa, chip, rawstars, refstars); + psastroUpdateChipToFPA (fpa, chip); //plot results Index: trunk/psastro/src/psastroMosaicSetMatch.c =================================================================== --- trunk/psastro/src/psastroMosaicSetMatch.c (revision 26258) +++ trunk/psastro/src/psastroMosaicSetMatch.c (revision 26259) @@ -26,7 +26,10 @@ double RADIUS = psMetadataLookupF32 (&status, recipe, radiusWord); if (!status) { - psError(PS_ERR_IO, false, "Failed to lookup matching radius: %s", radiusWord); - psFree (view); - return false; + psAbort("Failed to lookup matching radius: %s", radiusWord); + } + + int uniqIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.MOSAIC.UNIQ.ITER"); + if (!status) { + psAbort("Failed to lookup matching PSASTRO.MOSAIC.UNIQ.ITER"); } @@ -58,9 +61,9 @@ // select the raw objects for this readout - psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS"); + psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS.SUBSET"); if (rawstars == NULL) { continue; } // select the raw objects for this readout - psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS"); + psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS.SUBSET"); if (refstars == NULL) { continue; } psTrace ("psastro", 4, "Trying %ld refstars\n", refstars->n); @@ -68,4 +71,14 @@ psArray *matches = pmAstromRadiusMatchChip (rawstars, refstars, RADIUS); psTrace ("psastro", 4, "Matched %ld refstars\n", matches->n); + + if (iteration >= uniqIter) { + psArray *unique = pmAstromRadiusMatchUniq (rawstars, refstars, matches); + if (!unique) { + psLogMsg ("psastro", 3, "failed to generate a uniq set of matched sources\n"); + return false; + } + psFree (matches); + matches = unique; + } pmAstromVisualPlotMosaicMatches(rawstars, refstars, matches, iteration, recipe); Index: trunk/psastro/src/psastroOneChipFit.c =================================================================== --- trunk/psastro/src/psastroOneChipFit.c (revision 26258) +++ trunk/psastro/src/psastroOneChipFit.c (revision 26259) @@ -28,4 +28,7 @@ REQUIRED_RECIPE_VALUE (int nIter, "PSASTRO.MATCH.FIT.NITER", S32); + // for iterations >= uniqIter, require a single match per reference + REQUIRED_RECIPE_VALUE (int uniqIter, "PSASTRO.MATCH.UNIQ.ITER", S32); + // correct radius to FP units (physical pixel scale in microns per pixel) REQUIRED_RECIPE_VALUE (double pixelScale, "PSASTRO.PIXEL.SCALE", F32); @@ -61,4 +64,14 @@ return false; } + + if (iter >= uniqIter) { + psArray *unique = pmAstromRadiusMatchUniq (rawstars, refstars, match); + if (!unique) { + psLogMsg ("psastro", 3, "failed to generate a uniq set of matched sources\n"); + return false; + } + psFree (match); + match = unique; + } // modify the order to correspond to the actual number of matched stars: @@ -101,5 +114,5 @@ } - // XXX allow statitic to be set by the user + // XXX allow statistic to be set by the user // fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); fitStats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); @@ -117,10 +130,14 @@ // determine fromFPA transformation and apply new transformation to raw & ref stars - psastroUpdateChipToFPA (fpa, chip, rawstars, refstars); + // this applies the transformation to all stars (including subset) + psastroUpdateChipToFPA (fpa, chip); // updates PSASTRO.RAWSTARS and PSASTRO.REFSTARS // toSky converts from FPA & TPA units (microns) to sky units (radians) float plateScale = 0.5*(fpa->toSky->Xs + fpa->toSky->Ys)*3600.0*PM_DEG_RAD; - // float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) * plateScale; - float astError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) * plateScale; + + float rawXstdev = psStatsGetValue (results->xStats, psStatsStdevOption(results->xStats->options)); + float rawYstdev = psStatsGetValue (results->yStats, psStatsStdevOption(results->yStats->options)); + + float astError = 0.5*(rawXstdev + rawYstdev) * plateScale; int astNstar = results->yStats->clippedNvalues; psLogMsg ("psastro", PS_LOG_INFO, "pass %d, error: %f arcsec, Nstars: %d", iter, astError, astNstar); @@ -136,11 +153,14 @@ float plateScale = 0.5*(fpa->toSky->Xs + fpa->toSky->Ys)*3600.0*PM_DEG_RAD; + float rawXstdev = psStatsGetValue (results->xStats, psStatsStdevOption(results->xStats->options)); + float rawYstdev = psStatsGetValue (results->yStats, psStatsStdevOption(results->yStats->options)); + // pixError is the average 1D scatter in pixels ('results' are in FPA units = microns) - // float pixError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) / pixelScale; - float pixError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) / pixelScale; + float pixError = 0.5*(rawXstdev + rawYstdev) / pixelScale; // astError is the average 1D scatter in arcsec ('results' are in FPA units = microns) - // float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) * plateScale; - float astError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) * plateScale; + float astError = 0.5*(rawXstdev + rawYstdev) * plateScale; + + // x and y are forced to use the same subset of values: int astNstar = results->yStats->clippedNvalues; @@ -170,6 +190,12 @@ 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 - // XXX drop from here : determine fromFPA transformation and apply new transformation to raw & ref stars - // psastroUpdateChipToFPA (fpa, chip, rawstars, refstars); + // additional error measurements: + psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_CDX", PS_META_REPLACE, "chip astrometry X stdev (arcsec)", rawXstdev * plateScale); + psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_CSX", PS_META_REPLACE, "chip astrometry X systematic err (arcsec)", results->dXsys * plateScale); + psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_CRX", PS_META_REPLACE, "chip astrometry X 10-90 percentile (arcsec)", results->dXrange * plateScale); + + psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_CDY", PS_META_REPLACE, "chip astrometry Y stdev (arcsec)", rawYstdev * plateScale); + psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_CSY", PS_META_REPLACE, "chip astrometry Y systematic err (arcsec)", results->dYsys * plateScale); + psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_CRY", PS_META_REPLACE, "chip astrometry Y 10-90 percentile (arcsec)", results->dYrange * plateScale); // XXX check if we correctly applied the new transformation: @@ -189,4 +215,5 @@ psFree (results); psFree (fitStats); + return validSolution; } Index: trunk/psastro/src/psastroOneChipGrid.c =================================================================== --- trunk/psastro/src/psastroOneChipGrid.c (revision 26258) +++ trunk/psastro/src/psastroOneChipGrid.c (revision 26259) @@ -31,23 +31,16 @@ // generate the bright subset of maxNstar entries (note: rawstars is sorted by S/N) - psArray *subset = psArrayAlloc (PS_MIN (maxNstar, rawstars->n)); + psArray *gridStars = psArrayAlloc (PS_MIN (maxNstar, rawstars->n)); for (int i = 0; (i < maxNstar) && (i < rawstars->n); i++) { - subset->data[i] = psMemIncrRefCounter (rawstars->data[i]); + gridStars->data[i] = psMemIncrRefCounter (rawstars->data[i]); } - // XXX set clump scale from recipe - psArray *gridStars = psastroRemoveClumps (subset, 150); - psFree (subset); - - psArray *refSubset = psastroRemoveClumps (refstars, 150); - - psLogMsg ("psastro", 3, "grid search using %ld raw vs %ld ref stars\n", gridStars->n, refSubset->n); + psLogMsg ("psastro", 3, "grid search using %ld raw vs %ld ref stars\n", gridStars->n, refstars->n); // find initial offset / rotation / scale - pmAstromStats *gridStats = pmAstromGridMatch (gridStars, refSubset, recipe); + pmAstromStats *gridStats = pmAstromGridMatch (gridStars, refstars, recipe); if (gridStats == NULL) { psLogMsg ("psastro", 3, "failed to find a grid match solution\n"); psFree (gridStars); - psFree (refSubset); return false; } @@ -55,10 +48,9 @@ // tweak the position by finding peak of matches stars - stats = pmAstromGridTweak (gridStars, refSubset, recipe, gridStats); + stats = pmAstromGridTweak (gridStars, refstars, recipe, gridStats); if (stats == NULL) { psLogMsg ("psastro", 3, "failed to measure tweaked grid solution\n"); psFree (gridStats); psFree (gridStars); - psFree (refSubset); return false; } @@ -67,8 +59,7 @@ // adjust the chip.toFPA terms only pmAstromGridApply (chip->toFPA, stats); - psastroUpdateChipToFPA (fpa, chip, rawstars, refstars); + psastroUpdateChipToFPA (fpa, chip); // updates PSASTRO.RAWSTARS and PSASTRO.REFSTARS psFree (gridStats); psFree (gridStars); - psFree (refSubset); psFree (stats); Index: trunk/psastro/src/psastroRefstarSubset.c =================================================================== --- trunk/psastro/src/psastroRefstarSubset.c (revision 26258) +++ trunk/psastro/src/psastroRefstarSubset.c (revision 26259) @@ -16,5 +16,5 @@ // select the raw objects for this readout - psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS"); + psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS.SUBSET"); if (rawstars == NULL) { psError(PSASTRO_ERR_DATA, false, "missing rawstars in psastroRefstarSubset\n"); @@ -23,5 +23,5 @@ // select the raw objects for this readout - psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS"); + psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS.SUBSET"); if (refstars == NULL) { psError(PSASTRO_ERR_DATA, false, "missing refstars in psastroRefstarSubset\n"); @@ -74,6 +74,6 @@ psLogMsg ("psastro", 4, "keeping %ld of %ld reference stars\n", subset->n, refstars->n); - psMetadataRemoveKey (readout->analysis, "PSASTRO.REFSTARS"); - psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.REFSTARS", PS_DATA_ARRAY, "astrometry matches", subset); + psMetadataRemoveKey (readout->analysis, "PSASTRO.REFSTARS.SUBSET"); + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.REFSTARS.SUBSET", PS_DATA_ARRAY, "astrometry matches", subset); if (psTraceGetLevel("psastro.dump") > 0) { Index: trunk/psastro/src/psastroRemoveClumps.c =================================================================== --- trunk/psastro/src/psastroRemoveClumps.c (revision 26258) +++ trunk/psastro/src/psastroRemoveClumps.c (revision 26259) @@ -12,4 +12,70 @@ # include "psastroInternal.h" + +bool psastroRemoveClumpsRawstars (pmConfig *config) { + + bool status; + + pmChip *chip = NULL; + pmCell *cell = NULL; + pmReadout *readout = NULL; + + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSASTRO_RECIPE); + if (!recipe) { + psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe"); + return false; + } + + // select the input data sources + pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PSASTRO.INPUT"); + if (!input) { + psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); + return false; + } + + pmFPAview *view = pmFPAviewAlloc (0); + pmFPA *fpa = input->fpa; + + while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) { + psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); + if (!chip->process || !chip->file_exists) { continue; } + + while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) { + psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + if (!cell->process || !cell->file_exists) { continue; } + if (!chip->fromFPA) { continue; } + + // process each of the readouts + while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) { + if (! readout->data_exists) { continue; } + + // select the raw objects for this readout + psArray *rawstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.RAWSTARS"); + if (rawstars == NULL) { continue; } + + // generate a reduced subset excluding the clumps + // XXX do we need both RAWSTARS and SUBSET? + // XXX put these parameters in the recipe, please + psArray *subset = psastroRemoveClumpsIterate(rawstars, 150, 3); + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.RAWSTARS.SUBSET", PS_DATA_ARRAY, "astrometry objects", subset); + psFree (subset); + } + } + } + psFree (view); + return true; +} + +psArray *psastroRemoveClumpsIterate (psArray *input, int scale, int nIter) { + + psArray *newset = psMemIncrRefCounter (input); + for (int i = 0; i < nIter; i++) { + psArray *subset = psastroRemoveClumps (newset, 150); + psFree (newset); + newset = subset; + } + return newset; +} /** Index: trunk/psastro/src/psastroUtils.c =================================================================== --- trunk/psastro/src/psastroUtils.c (revision 26258) +++ trunk/psastro/src/psastroUtils.c (revision 26259) @@ -106,5 +106,5 @@ } -bool psastroUpdateChipToFPA (pmFPA *fpa, pmChip *chip, psArray *rawstars, psArray *refstars) { +bool psastroUpdateChipToFPA (pmFPA *fpa, pmChip *chip) { psRegion *region = pmChipPixels (chip); @@ -114,15 +114,33 @@ psFree (region); - for (int i = 0; i < rawstars->n; i++) { - pmAstromObj *raw = rawstars->data[i]; - - psPlaneTransformApply (raw->FP, chip->toFPA, raw->chip); - psPlaneTransformApply (raw->TP, fpa->toTPA, raw->FP); - psDeproject (raw->sky, raw->TP, fpa->toSky); - } - - for (int i = 0; i < refstars->n; i++) { - pmAstromObj *ref = refstars->data[i]; - psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP); + // loop over cells in this chip + for (int nCell = 0; nCell < chip->cells->n; nCell++) { + pmCell *cell = chip->cells->data[nCell]; + if (!cell->process || !cell->file_exists) { continue; } + + // loop over readouts in this cell + for (int nRead = 0; nRead < cell->readouts->n; nRead++) { + pmReadout *readout = cell->readouts->data[nRead]; + if (! readout->data_exists) { continue; } + + // select the raw objects for this readout + psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS"); + if (rawstars) { + for (int i = 0; i < rawstars->n; i++) { + pmAstromObj *raw = rawstars->data[i]; + psPlaneTransformApply (raw->FP, chip->toFPA, raw->chip); + psPlaneTransformApply (raw->TP, fpa->toTPA, raw->FP); + psDeproject (raw->sky, raw->TP, fpa->toSky); + } + } + + psArray *refstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.REFSTARS"); + if (refstars) { + for (int i = 0; i < refstars->n; i++) { + pmAstromObj *ref = refstars->data[i]; + psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP); + } + } + } } Index: trunk/psastro/src/psastroZeroPoint.c =================================================================== --- trunk/psastro/src/psastroZeroPoint.c (revision 26258) +++ trunk/psastro/src/psastroZeroPoint.c (revision 26259) @@ -30,6 +30,4 @@ // recipe options bool byExposure = psMetadataLookupBool (&status, recipe, "ZERO.POINT.BY.EXPOSURE"); - bool useMean = psMetadataLookupBool (&status, recipe, "ZERO.POINT.USE.MEAN"); - float edgeFraction = psMetadataLookupF32 (&status, recipe, "ZERO.POINT.EDGE.FRACTION"); // select the input data sources @@ -86,5 +84,5 @@ // calculate dMag for the matched stars just for this readout (well, chip) psMetadata *header = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER"); - psastroZeroPointAnalysis (header, dMag, zeropt, edgeFraction, useMean); + psastroZeroPointAnalysis (header, dMag, zeropt, recipe); psFree (dMag); dMag = NULL; @@ -96,5 +94,10 @@ if (byExposure) { psMetadata *header = psMetadataLookupMetadata (&status, fpa->analysis, "PSASTRO.HEADER"); - psastroZeroPointAnalysis (header, dMag, zeropt, edgeFraction, useMean); + if (!header) { + header = psMetadataAlloc (); + psMetadataAddMetadata (fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER", PS_META_REPLACE, "psastro header stats", header); + psFree (header); + } + psastroZeroPointAnalysis (header, dMag, zeropt, recipe); psFree (dMag); dMag = NULL; @@ -113,9 +116,9 @@ // select the raw objects for this readout - psArray *rawstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.RAWSTARS"); + psArray *rawstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.RAWSTARS.SUBSET"); if (rawstars == NULL) return dMag; // select the raw objects for this readout - psArray *refstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.REFSTARS"); + psArray *refstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.REFSTARS.SUBSET"); if (refstars == NULL) return dMag; @@ -141,5 +144,5 @@ } -bool psastroZeroPointAnalysis (psMetadata *header, psVector *dMag, float zeropt, float edgeFraction, bool useMean) { +bool psastroZeroPointAnalysis (psMetadata *header, psVector *dMag, float zeropt, psMetadata *recipe) { // XXX make this depend on the mode? @@ -148,4 +151,7 @@ return false; } + + bool status; + bool useMean = psMetadataLookupBool (&status, recipe, "ZERO.POINT.USE.MEAN"); // the zero point analysis depends on the type of desired statistic. For comparisons @@ -166,5 +172,5 @@ } } else { - stats = psastroStatsPercentile (dMag, edgeFraction); + stats = psastroStatsPercentile (dMag, recipe); if (!stats) { psError(PS_ERR_UNKNOWN, false, "failure to measure zero point by edge"); @@ -183,5 +189,5 @@ } -#define MAG_RESOLUTION 0.001 +#define MAG_RESOLUTION 0.0002 // set the bin closest to the corresponding value. if USE_END is +/- 1, @@ -233,9 +239,14 @@ // return results in stats->sampleMean, sampleStdev // XXX this is a misuse of psStats -- make our own structure? -psStats *psastroStatsPercentile (psVector *myVector, float flimit) { +psStats *psastroStatsPercentile (psVector *myVector, psMetadata *recipe) { // search for the 'blue' edge of the dMag distribution: // the distribution is not a normal population, but instead has a broad range with fairly hard edges. // construct a histogram and look for the + + bool status; + float edgeFraction = psMetadataLookupF32 (&status, recipe, "ZERO.POINT.EDGE.FRACTION"); + int edgeSample = psMetadataLookupS32 (&status, recipe, "ZERO.POINT.EDGE.SAMPLE"); + float edgeSampleFraction = psMetadataLookupF32 (&status, recipe, "ZERO.POINT.EDGE.SAMPLE.FRACTION"); // stats is first used to find the data range @@ -245,4 +256,6 @@ float min = NAN, max = NAN; // Mimimum and maximum values + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); + // Get the minimum and maximum values if (!psVectorStats(stats, myVector, NULL, NULL, 0)) goto escape; @@ -250,5 +263,5 @@ max = stats->max; if (isnan(min) || isnan(max)) goto escape; - psTrace("psastro", 6, "Data min/max is (%.2f, %.2f)\n", min, max); + psTrace("psastro", 5, "Data min/max is (%.2f, %.2f)\n", min, max); // If all data points have the same value, then we set the appropriate members of stats and return. @@ -264,50 +277,39 @@ float binSize = MAG_RESOLUTION; long numBins = (max - min) / binSize; // Number of bins - psTrace("psastro", 6, "Numbins is %ld\n", numBins); - psTrace("psastro", 6, "Creating a robust histogram from data range (%.2f - %.2f)\n", min, max); - - // Generate the histogram + psTrace("psastro", 5, "Numbins is %ld\n", numBins); + psTrace("psastro", 5, "Creating a robust histogram from data range (%.2f - %.2f)\n", min, max); + + // allocate the histogram containers histogram = psHistogramAlloc(min, max, numBins); - if (!psVectorHistogram(histogram, myVector, NULL, NULL, 0)) { - // if psVectorHistogram returns false, we have a programming error - psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram for psastroZeroPointAnalysis\n"); - psFree(histogram); - psFree(stats); - return NULL; - } - if (psTraceGetLevel("psastro") >= 8) { - PS_VECTOR_PRINT_F32(histogram->bounds); - PS_VECTOR_PRINT_F32(histogram->nums); - } - - // Convert the specific histogram to a cumulative histogram - // The cumulative histogram data points correspond to the UPPER bound value (N < Bin[i+1]) cumulative = psHistogramAlloc(min, max, numBins); - cumulative->nums->data.F32[0] = histogram->nums->data.F32[0]; - for (long i = 1; i < histogram->nums->n; i++) { - cumulative->nums->data.F32[i] = cumulative->nums->data.F32[i-1] + histogram->nums->data.F32[i]; - cumulative->bounds->data.F32[i-1] = histogram->bounds->data.F32[i]; - } - if (psTraceGetLevel("psastro") >= 8) { - PS_VECTOR_PRINT_F32(cumulative->bounds); - PS_VECTOR_PRINT_F32(cumulative->nums); - } - - // Find the bin which contains the first data point above the limit - long totalDataPoints = cumulative->nums->data.F32[numBins - 1]; - psTrace("psastro", 6, "Total data points is %ld\n", totalDataPoints); - - // find bin which is the lower bound of the limit value (value[bin] < f < value[bin+1] - long bin; - PS_BIN_FOR_VALUE(bin, cumulative->nums, flimit * totalDataPoints, 0); - psTrace("psastro", 6, "The bin is %ld (%.2f to %.2f)\n", bin, cumulative->bounds->data.F32[bin], cumulative->bounds->data.F32[bin+1]); - - // Linear interpolation to the limit value in bin units - float value; - PS_BIN_INTERPOLATE (value, cumulative->nums, cumulative->bounds, bin, totalDataPoints * flimit); - psTrace("psastro", 6, "limit value is %f\n", value); - - stats->clippedMean = value; - stats->clippedStdev = 0.0; // XXX derive correct error value + + // find the mean value: + stats->clippedMean = psastroStatsPercentileValue (histogram, cumulative, myVector, edgeFraction); + + int nSubset = myVector->n * edgeSampleFraction; + psVector *subset = psVectorAlloc (nSubset, PS_TYPE_F32); + + float Sum = 0.0; + float S2 = 0.0; + for (int i = 0; i < edgeSample; i++) { + + // generate the subset vector + for (long i = 0; i < nSubset; i++) { + double frnd = psRandomUniform(rng); + int entry = PS_MIN(myVector->n - 1, PS_MAX(0, myVector->n * frnd)); + + subset->data.F32[i] = myVector->data.F32[entry]; + } + + float value = psastroStatsPercentileValue (histogram, cumulative, subset, edgeFraction); + + Sum += value; + S2 += value*value; + } + psTrace("psastro", 6, "subset stats: Sum: %f, S2: %f, Npts: %d\n", Sum, S2, edgeSample); + + stats->clippedStdev = PS_MAX (sqrt(S2 / edgeSample - PS_SQR(Sum/edgeSample)), MAG_RESOLUTION); + psTrace("psastro", 5, "percentile stats %f +/- %f\n", stats->clippedMean, stats->clippedStdev); + stats->results |= PS_STAT_CLIPPED_MEAN; stats->results |= PS_STAT_CLIPPED_STDEV; @@ -317,4 +319,6 @@ psFree(histogram); psFree(cumulative); + psFree(subset); + psFree(rng); return stats; @@ -332,4 +336,50 @@ } + +// measure the edge of the sample at flimit +// return results in stats->sampleMean, sampleStdev +// XXX this is a misuse of psStats -- make our own structure? +float psastroStatsPercentileValue (psHistogram *histogram, psHistogram *cumulative, psVector *myVector, float flimit) { + + // need to initialize the histogram on each pass + psVectorInit (histogram->nums, 0); + if (!psVectorHistogram(histogram, myVector, NULL, NULL, 0)) { + // if psVectorHistogram returns false, we have a programming error + psAbort ("Unable to generate histogram for psastroZeroPointAnalysis"); + } + if (psTraceGetLevel("psastro") >= 8) { + PS_VECTOR_PRINT_F32(histogram->bounds); + PS_VECTOR_PRINT_F32(histogram->nums); + } + + // Convert the specific histogram to a cumulative histogram + // The cumulative histogram data points correspond to the UPPER bound value (N < Bin[i+1]) + cumulative->nums->data.F32[0] = histogram->nums->data.F32[0]; + for (long i = 1; i < histogram->nums->n; i++) { + cumulative->nums->data.F32[i] = cumulative->nums->data.F32[i-1] + histogram->nums->data.F32[i]; + cumulative->bounds->data.F32[i-1] = histogram->bounds->data.F32[i]; + } + if (psTraceGetLevel("psastro") >= 8) { + PS_VECTOR_PRINT_F32(cumulative->bounds); + PS_VECTOR_PRINT_F32(cumulative->nums); + } + + // Find the bin which contains the first data point above the limit + long numBins = cumulative->nums->n; + long totalDataPoints = cumulative->nums->data.F32[numBins - 1]; + psTrace("psastro", 6, "Total data points is %ld\n", totalDataPoints); + + // find bin which is the lower bound of the limit value (value[bin] < f < value[bin+1] + long bin; + PS_BIN_FOR_VALUE(bin, cumulative->nums, flimit * totalDataPoints, 0); + psTrace("psastro", 6, "The bin is %ld (%.4f to %.4f)\n", bin, cumulative->bounds->data.F32[bin], cumulative->bounds->data.F32[bin+1]); + + // Linear interpolation to the limit value in bin units + float value; + PS_BIN_INTERPOLATE (value, cumulative->nums, cumulative->bounds, bin, totalDataPoints * flimit); + psTrace("psastro", 6, "limit value is %f\n", value); + + return value; +} # define ESCAPE(MSG) { \