Index: /trunk/psModules/src/astrom/Makefile.am =================================================================== --- /trunk/psModules/src/astrom/Makefile.am (revision 36833) +++ /trunk/psModules/src/astrom/Makefile.am (revision 36834) @@ -11,5 +11,6 @@ pmAstrometryRefstars.c \ pmAstrometryWCS.c \ - pmAstrometryVisual.c + pmAstrometryVisual.c \ + pmKHcorrect.c pkginclude_HEADERS = \ @@ -21,5 +22,6 @@ pmAstrometryRefstars.h \ pmAstrometryWCS.h \ - pmAstrometryVisual.h + pmAstrometryVisual.h \ + pmKHcorrect.h CLEANFILES = *~ Index: /trunk/psModules/src/astrom/pmAstrometryModel.c =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryModel.c (revision 36833) +++ /trunk/psModules/src/astrom/pmAstrometryModel.c (revision 36834) @@ -39,4 +39,5 @@ #include "pmFPAExtent.h" #include "pmFPAfileFitsIO.h" +#include "pmConcepts.h" #include "pmAstrometryWCS.h" #include "pmAstrometryUtils.h" @@ -452,28 +453,4 @@ } -int pmConceptsChipNumberFromName (pmFPA *fpa, char *name) { - - for (int i = 0; i < fpa->chips->n; i++) { - pmChip *chip = fpa->chips->data[i]; - if (!chip) continue; - char *thisone = psMetadataLookupStr (NULL, chip->concepts, "CHIP.NAME"); - if (!thisone) continue; - if (!strcmp (name, thisone)) return (i); - } - return -1; -} - -pmChip *pmConceptsChipFromName (pmFPA *fpa, char *name) { - - for (int i = 0; i < fpa->chips->n; i++) { - pmChip *chip = fpa->chips->data[i]; - if (!chip) continue; - char *thisone = psMetadataLookupStr (NULL, chip->concepts, "CHIP.NAME"); - if (!thisone) continue; - if (!strcmp (name, thisone)) return (chip); - } - return NULL; -} - // first layer converts Chip to Focal Plane bool pmAstromModelReadChips (pmFPAfile *file) { Index: /trunk/psModules/src/astrom/pmAstrometryModel.h =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryModel.h (revision 36833) +++ /trunk/psModules/src/astrom/pmAstrometryModel.h (revision 36834) @@ -27,7 +27,4 @@ bool pmAstromModelWriteChips (pmFPAfile *file); -int pmConceptsChipNumberFromName (pmFPA *fpa, char *name); -pmChip *pmConceptsChipFromName (pmFPA *fpa, char *name); - bool pmAstromModelReadForView (const pmFPAview *view, pmFPAfile *file, const pmConfig *config); bool pmAstromModelReadFPA (pmFPAfile *file); Index: /trunk/psModules/src/astrom/pmAstrometryObjects.c =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryObjects.c (revision 36833) +++ /trunk/psModules/src/astrom/pmAstrometryObjects.c (revision 36834) @@ -571,13 +571,21 @@ psMemSetDeallocator (stats, (psFreeFunc)pmAstromStatsFree); - // stats->center = {0, 0, 0, 0}; - // stats->offset = {0, 0, 0, 0}; - stats->angle = 0.0; - stats->scale = 1.0; - stats->minMetric = 0.0; - stats->minVar = 0.0; - stats->nMatch = 0; - stats->nTest = 0; - stats->nSigma = 0; + stats->center.x = 0; + stats->center.y = 0; + stats->center.xErr = 0; + stats->center.yErr = 0; + + stats->offset.x = 0; + stats->offset.y = 0; + stats->offset.xErr = 0; + stats->offset.yErr = 0; + + stats->angle = 0.0; + stats->scale = 1.0; + stats->minMetric = 0.0; + stats->minVar = 0.0; + stats->nMatch = 0; + stats->nTest = 0; + stats->nSigma = 0; return (stats); @@ -914,6 +922,5 @@ // fprintf (stderr, "sigma: nMatch: %d, nTest: %d, nTen: %d\n", stats->nMatch, stats->nTest, sort->data.U32[sort->n - 10]); - - psFree (sort); + psFree (sort); psFree (listNP); psFree (gridNP); Index: /trunk/psModules/src/astrom/pmAstrometryObjects.h =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryObjects.h (revision 36833) +++ /trunk/psModules/src/astrom/pmAstrometryObjects.h (revision 36834) @@ -42,4 +42,5 @@ float Color; ///< object color float dMag; ///< error on object magnitude + float SBinst; ///< surface brightness, used for Koppenhoefer correction } pmAstromObj; Index: /trunk/psModules/src/astrom/pmAstrometryVisual.c =================================================================== --- /trunk/psModules/src/astrom/pmAstrometryVisual.c (revision 36833) +++ /trunk/psModules/src/astrom/pmAstrometryVisual.c (revision 36834) @@ -946,5 +946,5 @@ KapaSendLabel (kapa2, "X (FP)", KAPA_LABEL_XM); KapaSendLabel (kapa2, "Y (FP)", KAPA_LABEL_YM); - KapaSendLabel (kapa2, "pmAstromGridAngle residuals. Box: Correlation Peak.", KAPA_LABEL_XP); + KapaSendLabel (kapa2, "pmAstromGridAngle red: raw, black: ref.", KAPA_LABEL_XP); // plot the REF data. (also calculate the plot ranges, accumulate the plot vectors) @@ -966,13 +966,18 @@ KapaSetLimits(kapa2, &graphdata); - psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); + psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN | PS_STAT_MIN | PS_STAT_MAX ); psVectorStats (stats, zPlot, NULL, NULL, 0); - float zero = stats->sampleMedian + 3.0; - float range = 6.0; - + float range = stats->max - stats->min; + range = PS_MAX (0.5, PS_MIN (6.0, range)); + float zero = stats->sampleMedian + 0.25*range; + + float maxZ = zPlot->data.F32[0], minZ = zPlot->data.F32[0]; for (int i = 0; i < zPlot->n; i++) { + maxZ = PS_MAX (maxZ, zPlot->data.F32[i]); + minZ = PS_MIN (minZ, zPlot->data.F32[i]); float value = (zero - zPlot->data.F32[i]) / range; zPlot->data.F32[i] = PS_MAX(0.0, PS_MIN(1.0, value)); } + fprintf (stderr, "ref mags: %f to %f (%f median)\n", minZ, maxZ, stats->sampleMedian); // the point size will be scaled from the z vector @@ -998,11 +1003,18 @@ psStatsInit(stats); psVectorStats (stats, zPlot, NULL, NULL, 0); - zero = stats->sampleMedian + 3.0; - range = 6.0; - + range = stats->max - stats->min; + range = PS_MAX (0.5, PS_MIN (6.0, range)); + zero = stats->sampleMedian + 0.25*range; + // zero = stats->sampleMedian + 1.0; + // range = 6.0; + + maxZ = zPlot->data.F32[0], minZ = zPlot->data.F32[0]; for (int i = 0; i < zPlot->n; i++) { + maxZ = PS_MAX (maxZ, zPlot->data.F32[i]); + minZ = PS_MIN (minZ, zPlot->data.F32[i]); float value = (zero - zPlot->data.F32[i]) / range; zPlot->data.F32[i] = PS_MAX(0.0, PS_MIN(1.0, value)); } + fprintf (stderr, "raw mags: %f to %f (%f median)\n", minZ, maxZ, stats->sampleMedian); // the point size will be scaled from the z vector Index: /trunk/psModules/src/astrom/pmKHcorrect.c =================================================================== --- /trunk/psModules/src/astrom/pmKHcorrect.c (revision 36834) +++ /trunk/psModules/src/astrom/pmKHcorrect.c (revision 36834) @@ -0,0 +1,252 @@ +/** @file pmKHcorrect.c + * @brief Functions to read (and write?) Koppenhoefer correction file + * + * The Koppenhoefer correction is needed for some chips of gpc1 before the camera voltages were adjusted 2011/05/11. + * The correction is a modification of the X (and possibly Y) coordinate of a star which depends on the instrumental + * surface brightness, defined as -2.5 log_10 (DN) + 5.0 log_10 fwhm_maj [XXX be careful about the definition of fwhm_maj] + * + * @ingroup AstroImage + * @author EAM, IfA + * + * Copyright 2014 Institute for Astronomy, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +/******************************************************************************/ +/* INCLUDE FILES */ +/******************************************************************************/ +#include +#include +#include +#include +#include +#include // for unlink +#include + +#include "pmConfig.h" +#include "pmDetrendDB.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmFPALevel.h" +#include "pmFPAview.h" +#include "pmFPAfile.h" +#include "pmFPAExtent.h" +#include "pmFPAfileFitsIO.h" +#include "pmConcepts.h" +#include "pmKHcorrect.h" + +static void KHcorrectDataFree (KHcorrectData *spline) { + + if (!spline) return; + + psFree (spline->xk); + psFree (spline->yk); + psFree (spline->y2); + return; +} + +KHcorrectData *KHcorrectDataAlloc (int Nrow) { + + // allocate xk[Nrow], etc + + KHcorrectData *spline = (KHcorrectData *) psAlloc(sizeof(KHcorrectData)); + psMemSetDeallocator(spline, (psFreeFunc) KHcorrectDataFree); + + spline->N = Nrow; + spline->xk = (float *) psAlloc(Nrow*sizeof(float)); + spline->yk = (float *) psAlloc(Nrow*sizeof(float)); + spline->y2 = (float *) psAlloc(Nrow*sizeof(float)); + + return spline; +} + +// the KH correction is a function of the instrumental surface brightness, SBinst +float KHcorrectApply (KHcorrectData *spline, float X) { + + int N = spline->N; + + // saturate correction at high and low ends + if (X < spline->xk[ 0]) return spline->yk[ 0]; + if (X > spline->xk[N-1]) return spline->yk[N-1]; + + float *xk = spline->xk; + float *yk = spline->yk; + float *y2 = spline->y2; + + /* find correct element in array (x must be sorted) */ + int lo = 0; + int hi = N-1; + while (hi - lo > 1) { + int i = 0.5*(hi+lo); + if (xk[i] > X) { + hi = i; + } else { + lo = i; + } + } + + /* error condition: duplicate abssisca */ + float dx = xk[hi] - xk[lo]; + if (dx == 0.0) { + return (0.0); + } + + /* evaluate spline */ + float a = (xk[hi] - X) / dx; + float b = (X - xk[lo]) / dx; + + float value = a*yk[lo] + b*yk[hi] + ((a*a - 1.0)*a*y2[lo] + (b*b - 1.0)*b*y2[hi])*(dx*dx) / 6.0; + return (value); +} + +/********************* CheckDataStatus functions *****************************/ + +bool pmKHcorrectCheckDataStatusForView (const pmFPAview *view, pmFPAfile *file) { + psError(PS_ERR_IO, false, "Check Data Status not defined"); + return false; +} + +bool pmKHcorrectCheckDataStatusForFPA (const pmFPA *fpa) { + psError(PS_ERR_IO, false, "Check Data Status not defined"); + return false; +} + +bool pmKHcorrectCheckDataStatusForChip (const pmChip *chip) { + psError(PS_ERR_IO, false, "Check Data Status not defined"); + return false; +} + +/********************* Write Data functions *****************************/ + +// NOTE : these are not exposed because I don't think we need them (we do not create KHcorrect +// in psModules based tool, but in DVO. + +bool pmKHcorrectWriteFPA (pmFPAfile *file, const pmFPA *fpa); +bool pmKHcorrectWriteForView (const pmFPAview *view, pmFPAfile *file, pmConfig *config) +{ + // write the full model in one pass: require the level to be FPA + if (view->chip != -1) { + psError(PS_ERR_IO, false, "Koppenhoefer Correction must be written at the FPA level"); + return false; + } + + pmFPA *fpa = pmFPAfileSuitableFPA(file, view, config, false); // Suitable FPA for writing + + if (!pmKHcorrectWriteFPA(file, fpa)) { + psError(PS_ERR_IO, false, "Failed to write KH Correction for fpa"); + psFree(fpa); + return false; + } + + psFree(fpa); + + return true; +} + +// write out all chip-level KH Correction data for this FPA +bool pmKHcorrectWriteFPA (pmFPAfile *file, const pmFPA *fpa) +{ + psError(PS_ERR_IO, false, "output for KH Correction is not defined"); + return false; +} + +/********************* Read Data functions *****************************/ + +bool pmKHcorrectReadForView (const pmFPAview *view, pmFPAfile *file, const pmConfig *config) + { + // read the full model in one pass: require the level to be FPA + if (view->chip != -1) { + psError(PS_ERR_IO, false, "KH Correction must be read at the FPA level"); + return false; + } + + if (!pmKHcorrectReadFPA (file)) { + psError(PS_ERR_IO, false, "Failed to read KH Correction for fpa"); + return false; + } + return true; + } + +// read in all chip-level KH Correction data for this FPA +bool pmKHcorrectReadFPA (pmFPAfile *file) { + + if (!pmKHcorrectReadChips (file)) { + psError(PS_ERR_IO, false, "Failed to read KH Correction for chips"); + return false; + } + + return true; +} + +// read the set of tables, one for each chip +bool pmKHcorrectReadChips (pmFPAfile *file) { + + bool haveData, status; + + // loop over the extensions + // for each extension, use the extname (eg, XY01.DX.T0) to assign to a chip + + // move to the start of the file + haveData = psFitsMoveExtNum (file->fits, 1, false); + if (!haveData) { + psError(PS_ERR_IO, false, "Failed to read even the first extension?"); + return false; + } + + while (haveData) { + + // load the header + psMetadata *header = psFitsReadHeader(NULL, file->fits); // The FITS header + if (!header) psAbort("cannot read model header"); + + // load the full model in one shot + psArray *model = psFitsReadTable (file->fits); + if (!model) psAbort("cannot read model"); + + // determine the chip: + char *extname = psMetadataLookupStr (&status, header, "EXTNAME"); + psLogMsg ("psModules.astrom", 4, "read %ld rows from Koppenhoefer correction file, extname %s\n", model->n, extname); + + // I expect to find a name of the form: chipName.dir.tset (eg, XY01.DX.T0) + // where chipName like 'XY01' + // dir = 'DX' + // tset = 'T0' + psAssert (strlen(extname) == 10, "invalid extension %s", extname); + psAssert (extname[5] == 'D', "invalid extension %s", extname); + psAssert (extname[6] == 'X', "invalid extension %s", extname); + psAssert (extname[8] == 'T', "invalid extension %s", extname); + psAssert (extname[9] == '0', "invalid extension %s", extname); + + char chipName[5]; + strncpy (chipName, extname, 4); + chipName[4] = 0; + + pmChip *chip = pmConceptsChipFromName (file->fpa, chipName); + if (!chip) psAbort ("invalid chip?"); + + KHcorrectData *spline = KHcorrectDataAlloc (model->n); + + // parse the model entries + for (int i = 0; i < model->n; i++) { + psMetadata *row = model->data[i]; + + spline->xk[i] = psMetadataLookupF32(&status, row, "X_KNOT"); + spline->yk[i] = psMetadataLookupF32(&status, row, "Y_KNOT"); + spline->y2[i] = psMetadataLookupF32(&status, row, "DY2_DX"); + } + psMetadataAddUnknown (chip->analysis, PS_LIST_TAIL, "KH.CORRECT", PS_META_REPLACE, "", spline); + psFree (spline); + + psFree (model); + psFree (header); + + // move to the next extension + haveData = psFitsMoveExtNum (file->fits, 1, true); + } + + return true; +} + Index: /trunk/psModules/src/astrom/pmKHcorrect.h =================================================================== --- /trunk/psModules/src/astrom/pmKHcorrect.h (revision 36834) +++ /trunk/psModules/src/astrom/pmKHcorrect.h (revision 36834) @@ -0,0 +1,35 @@ +/** @file pmKHcorrect.h + * @brief Functions to read (and write?) Koppenhoefer correction file + * + * The Koppenhoefer correction is needed for some chips of gpc1 before the camera voltages were adjusted 2011/05/11. + * The correction is a modification of the X (and possibly Y) coordinate of a star which depends on the instrumental + * surface brightness, defined as -2.5 log_10 (DN) + 5.0 log_10 fwhm_maj [XXX be careful about the definition of fwhm_maj] + * + * @ingroup AstroImage + * @author EAM, IfA + * + * Copyright 2014 Institute for Astronomy, University of Hawaii + */ + +#ifndef PM_KH_CORRECT_H +#define PM_KH_CORRECT_H + +/// @addtogroup Astrometry +/// @{ + +typedef struct { + int N; + float *xk; + float *yk; + float *y2; +} KHcorrectData; + +KHcorrectData *KHcorrectDataAlloc (int Nrow); +float KHcorrectApply (KHcorrectData *spline, float X); + +bool pmKHcorrectReadForView (const pmFPAview *view, pmFPAfile *file, const pmConfig *config); +bool pmKHcorrectReadFPA (pmFPAfile *file); +bool pmKHcorrectReadChips (pmFPAfile *file); + +/// @} +#endif Index: /trunk/psModules/src/camera/pmFPAfile.c =================================================================== --- /trunk/psModules/src/camera/pmFPAfile.c (revision 36833) +++ /trunk/psModules/src/camera/pmFPAfile.c (revision 36834) @@ -531,7 +531,6 @@ return PM_FPA_FILE_LINEARITY; } - // deprecate this? if (!strcasecmp(type, "ASTROM")) { - return PM_FPA_FILE_ASTROM_MODEL; + return PM_FPA_FILE_ASTROM_MODEL; } if (!strcasecmp(type, "ASTROM.MODEL")) { @@ -540,4 +539,7 @@ if (!strcasecmp(type, "ASTROM.REFSTARS")) { return PM_FPA_FILE_ASTROM_REFSTARS; + } + if (!strcasecmp(type, "KH.CORRECT")) { + return PM_FPA_FILE_KH_CORRECT; } if (!strcasecmp(type, "SUBKERNEL")) { @@ -593,4 +595,6 @@ case PM_FPA_FILE_ASTROM_REFSTARS: return ("ASTROM.REFSTARS"); + case PM_FPA_FILE_KH_CORRECT: + return ("KH.CORRECT"); case PM_FPA_FILE_SUBKERNEL: return ("SUBKERNEL"); Index: /trunk/psModules/src/camera/pmFPAfile.h =================================================================== --- /trunk/psModules/src/camera/pmFPAfile.h (revision 36833) +++ /trunk/psModules/src/camera/pmFPAfile.h (revision 36834) @@ -48,4 +48,5 @@ PM_FPA_FILE_ASTROM_MODEL, PM_FPA_FILE_ASTROM_REFSTARS, + PM_FPA_FILE_KH_CORRECT, PM_FPA_FILE_SUBKERNEL, PM_FPA_FILE_SRCTEXT, Index: /trunk/psModules/src/camera/pmFPAfileFitsIO.c =================================================================== --- /trunk/psModules/src/camera/pmFPAfileFitsIO.c (revision 36833) +++ /trunk/psModules/src/camera/pmFPAfileFitsIO.c (revision 36834) @@ -150,5 +150,7 @@ case PM_FPA_FILE_PSF: case PM_FPA_FILE_ASTROM_MODEL: - case PM_FPA_FILE_ASTROM_REFSTARS: { + case PM_FPA_FILE_ASTROM_REFSTARS: + case PM_FPA_FILE_KH_CORRECT: + { pmHDU *hdu = pmFPAviewThisHDU(view, fpa); if (hdu) { Index: /trunk/psModules/src/camera/pmFPAfileIO.c =================================================================== --- /trunk/psModules/src/camera/pmFPAfileIO.c (revision 36833) +++ /trunk/psModules/src/camera/pmFPAfileIO.c (revision 36834) @@ -48,4 +48,5 @@ #include "pmPSF_IO.h" +#include "pmKHcorrect.h" #include "pmAstrometryModel.h" #include "pmAstrometryRefstars.h" @@ -234,4 +235,7 @@ status = pmAstromModelReadForView (view, file, config); break; + case PM_FPA_FILE_KH_CORRECT: + status = pmKHcorrectReadForView (view, file, config); + break; case PM_FPA_FILE_EXPNUM: status = pmExpNumRead(view, file, config); @@ -324,4 +328,5 @@ case PM_FPA_FILE_ASTROM_MODEL: case PM_FPA_FILE_ASTROM_REFSTARS: + case PM_FPA_FILE_KH_CORRECT: case PM_FPA_FILE_JPEG: case PM_FPA_FILE_KAPA: @@ -408,4 +413,8 @@ } } + if (file->type == PM_FPA_FILE_KH_CORRECT) { + psTrace("psModules.camera", 6, "skip write for %s, no write function defined", file->name); + return true; + } // open the file if not yet opened @@ -500,4 +509,8 @@ case PM_FPA_FILE_ASTROM_REFSTARS: status = pmAstromRefstarsWriteForView (view, file, config); + break; + + case PM_FPA_FILE_KH_CORRECT: + psError(PS_ERR_IO, true, "cannot write type KH.CORRECT (%s)", file->name); break; @@ -576,4 +589,5 @@ case PM_FPA_FILE_ASTROM_MODEL: case PM_FPA_FILE_ASTROM_REFSTARS: + case PM_FPA_FILE_KH_CORRECT: case PM_FPA_FILE_LINEARITY: case PM_FPA_FILE_EXPNUM: @@ -653,4 +667,5 @@ case PM_FPA_FILE_ASTROM_MODEL: case PM_FPA_FILE_ASTROM_REFSTARS: + case PM_FPA_FILE_KH_CORRECT: case PM_FPA_FILE_EXPNUM: psTrace ("psModules.camera", 6, "NOT freeing %s (%s) : save for further analysis\n", file->filename, file->name); @@ -815,4 +830,5 @@ case PM_FPA_FILE_ASTROM_MODEL: case PM_FPA_FILE_ASTROM_REFSTARS: + case PM_FPA_FILE_KH_CORRECT: case PM_FPA_FILE_LINEARITY: case PM_FPA_FILE_EXPNUM: @@ -1018,4 +1034,5 @@ case PM_FPA_FILE_EXPNUM: case PM_FPA_FILE_ASTROM_MODEL: + case PM_FPA_FILE_KH_CORRECT: case PM_FPA_FILE_SX: case PM_FPA_FILE_RAW: Index: /trunk/psModules/src/concepts/pmConcepts.c =================================================================== --- /trunk/psModules/src/concepts/pmConcepts.c (revision 36833) +++ /trunk/psModules/src/concepts/pmConcepts.c (revision 36834) @@ -578,2 +578,26 @@ } +int pmConceptsChipNumberFromName (pmFPA *fpa, char *name) { + + for (int i = 0; i < fpa->chips->n; i++) { + pmChip *chip = fpa->chips->data[i]; + if (!chip) continue; + char *thisone = psMetadataLookupStr (NULL, chip->concepts, "CHIP.NAME"); + if (!thisone) continue; + if (!strcmp (name, thisone)) return (i); + } + return -1; +} + +pmChip *pmConceptsChipFromName (pmFPA *fpa, char *name) { + + for (int i = 0; i < fpa->chips->n; i++) { + pmChip *chip = fpa->chips->data[i]; + if (!chip) continue; + char *thisone = psMetadataLookupStr (NULL, chip->concepts, "CHIP.NAME"); + if (!thisone) continue; + if (!strcmp (name, thisone)) return (chip); + } + return NULL; +} + Index: /trunk/psModules/src/concepts/pmConcepts.h =================================================================== --- /trunk/psModules/src/concepts/pmConcepts.h (revision 36833) +++ /trunk/psModules/src/concepts/pmConcepts.h (revision 36834) @@ -161,4 +161,8 @@ ); +// some utility functions: +int pmConceptsChipNumberFromName (pmFPA *fpa, char *name); +pmChip *pmConceptsChipFromName (pmFPA *fpa, char *name); + /// @} #endif Index: /trunk/psModules/src/detrend/pmDetrendDB.c =================================================================== --- /trunk/psModules/src/detrend/pmDetrendDB.c (revision 36833) +++ /trunk/psModules/src/detrend/pmDetrendDB.c (revision 36834) @@ -111,4 +111,5 @@ DETREND_STRING_CASE(LINEARITY); DETREND_STRING_CASE(AUXMASK); + DETREND_STRING_CASE(KH_CORRECT); default: return NULL; Index: /trunk/psModules/src/detrend/pmDetrendDB.h =================================================================== --- /trunk/psModules/src/detrend/pmDetrendDB.h (revision 36833) +++ /trunk/psModules/src/detrend/pmDetrendDB.h (revision 36834) @@ -40,4 +40,5 @@ PM_DETREND_TYPE_LINEARITY, PM_DETREND_TYPE_AUXMASK, + PM_DETREND_TYPE_KH_CORRECT, } pmDetrendType; Index: /trunk/psModules/src/extras/pmVisual.c =================================================================== --- /trunk/psModules/src/extras/pmVisual.c (revision 36833) +++ /trunk/psModules/src/extras/pmVisual.c (revision 36834) @@ -162,4 +162,77 @@ +// ask the user to continue or not. give up after 2 seconds. +bool pmVisualAskUserOrDump(bool *plotFlag, bool *dumpData) +{ + struct timeval timeout; + fd_set fdSet; + int status; + + if (dumpData) *dumpData = false; + + char key[10]; + if (plotFlag && dumpData) { + fprintf (stderr, "[p]ause? [c]ontinue? [s]kip the rest of these plots? [d]ump the data? [a]bort all visual plots? (c) "); + } + if (plotFlag && !dumpData) { + fprintf (stderr, "[p]ause? [c]ontinue? [s]kip the rest of these plots? [a]bort all visual plots? (c) "); + } + if (!plotFlag && dumpData) { + fprintf (stderr, "[p]ause? [c]ontinue? [d]ump the data? [a]bort all visual plots? (c) "); + } + if (!plotFlag && !dumpData) { + fprintf (stderr, "[p]ause? [c]ontinue? [a]bort all visual plots? (c) "); + } + + /* Wait up to 1.0 second for a response, then continue */ + timeout.tv_sec = 10; + timeout.tv_usec = 0; + + FD_ZERO (&fdSet); + FD_SET (STDIN_FILENO, &fdSet); + + status = select (1, &fdSet, NULL, NULL, &timeout); + if (status <= 0) { + fprintf (stderr, "\n"); + return true; // if no data, give up + } + + while (true) { + if (!fgets(key, 8, stdin)) { + psWarning("Unable to read option"); + } + switch (key[0]) { + case 's': + if (plotFlag) *plotFlag = false; + return true; + case 'd': + if (dumpData) *dumpData = true; + return true; + case 'a': + isVisual = false; + return true; + case 'c': + case '\n': + return true; + default: + break; + } + + if (plotFlag && dumpData) { + fprintf (stderr, "[p]ause? [c]ontinue? [s]kip the rest of these plots? [d]ump the data? [a]bort all visual plots? (c) "); + } + if (plotFlag && !dumpData) { + fprintf (stderr, "[p]ause? [c]ontinue? [s]kip the rest of these plots? [a]bort all visual plots? (c) "); + } + if (!plotFlag && dumpData) { + fprintf (stderr, "[p]ause? [c]ontinue? [d]ump the data? [a]bort all visual plots? (c) "); + } + if (!plotFlag && !dumpData) { + fprintf (stderr, "[p]ause? [c]ontinue? [a]bort all visual plots? (c) "); + } + } + return true; +} + bool pmVisualImStats(psImage *image, double *mean, double *stdev, double *min, double *max) { Index: /trunk/psModules/src/extras/pmVisual.h =================================================================== --- /trunk/psModules/src/extras/pmVisual.h (revision 36833) +++ /trunk/psModules/src/extras/pmVisual.h (revision 36834) @@ -52,4 +52,12 @@ */ bool pmVisualAskUser(bool *plotFlag); + + +/** Ask the user how to proceed. + * At the user's request, this will disable diagnostic plotting. + * @param plotFlag, set to false if this plot should be disabled in the future + * @param dumpData, set to true if user requests a data dump + */ +bool pmVisualAskUserOrDump(bool *plotFlag, bool *dumpData); @@ -138,4 +146,5 @@ bool pmVisualInitGraph (int kapa, void *section, void *graphdata); bool pmVisualAskUser(bool *plotFlag); +bool pmVisualAskUserOrDump(bool *plotFlag, bool *dumpData); bool pmVisualScaleImage(int kapaFD, psImage *inImage, const char *name, int channel, bool clip); Index: /trunk/psModules/src/objects/pmModel.c =================================================================== --- /trunk/psModules/src/objects/pmModel.c (revision 36833) +++ /trunk/psModules/src/objects/pmModel.c (revision 36834) @@ -39,4 +39,6 @@ { psTrace("psModules.objects", 10, "---- %s() begin ----\n", __func__); + if (!tmp) return; + psFree(tmp->params); psFree(tmp->dparams); Index: /trunk/psModules/src/objects/pmPSFtryMakePSF.c =================================================================== --- /trunk/psModules/src/objects/pmPSFtryMakePSF.c (revision 36833) +++ /trunk/psModules/src/objects/pmPSFtryMakePSF.c (revision 36834) @@ -220,4 +220,26 @@ } + // weed out extreme e0 outliers here: find the median and exclude points not in the + // range MEDIAN / 5 < e0 < 5 * MEDIAN + { + psStats *e0stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); + if (psVectorStats (e0stats, e0, NULL, srcMask, 0xff)) { + float e0med = e0stats->sampleMedian; + + for (int i = 0; i < sources->n; i++) { + // skip any masked sources (failed to fit one of the model steps or get a magnitude) + if (srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue; + + if (e0->data.F32[i] < 0.2*e0med) { + srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_OUTLIER; + } + if (e0->data.F32[i] > 5.0*e0med) { + srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_OUTLIER; + } + } + } + psFree (e0stats); + } + // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each. // This way, the parameters masked by one of the fits will be applied to the others Index: /trunk/psModules/src/objects/pmSourceIO_MatchedRefs.c =================================================================== --- /trunk/psModules/src/objects/pmSourceIO_MatchedRefs.c (revision 36833) +++ /trunk/psModules/src/objects/pmSourceIO_MatchedRefs.c (revision 36834) @@ -99,22 +99,11 @@ // 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); - -# if (0) - // XXX test - FILE *outfile = fopen ("refstars.dat", "w"); - assert (outfile); - for (int nn = 0; nn < refstars->n; nn++) { - pmAstromObj *ref = refstars->data[nn]; - fprintf (outfile, "%lf %lf\n", ref->sky->r*PS_DEG_RAD, ref->sky->d*PS_DEG_RAD); - } - fclose (outfile); -# endif psArray *matches = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.MATCH"); @@ -128,16 +117,18 @@ psMetadata *row = psMetadataAlloc (); - psMetadataAdd (row, PS_LIST_TAIL, "RA", PS_DATA_F64, "right ascension (deg, J2000)", PM_DEG_RAD*ref->sky->r); - psMetadataAdd (row, PS_LIST_TAIL, "DEC", PS_DATA_F64, "declination (deg, J2000)", PM_DEG_RAD*ref->sky->d); - psMetadataAdd (row, PS_LIST_TAIL, "X_CHIP", PS_DATA_F32, "x coord on chip", raw->chip->x); - psMetadataAdd (row, PS_LIST_TAIL, "Y_CHIP", PS_DATA_F32, "y coord on chip", raw->chip->y); - psMetadataAdd (row, PS_LIST_TAIL, "X_CHIP_FIT",PS_DATA_F32, "x fitted coord on chip", ref->chip->x); - psMetadataAdd (row, PS_LIST_TAIL, "Y_CHIP_FIT",PS_DATA_F32, "y fitted coord on chip", ref->chip->y); - psMetadataAdd (row, PS_LIST_TAIL, "X_FPA", PS_DATA_F32, "x coord on focal plane", raw->FP->x); - psMetadataAdd (row, PS_LIST_TAIL, "Y_FPA", PS_DATA_F32, "y coord on focal plane", raw->FP->y); - psMetadataAdd (row, PS_LIST_TAIL, "MAG_INST", PS_DATA_F32, "instrumental magnitude", raw->Mag); - psMetadataAdd (row, PS_LIST_TAIL, "MAG_REF", PS_DATA_F32, "reference star magnitude", ref->Mag); - psMetadataAdd (row, PS_LIST_TAIL, "COLOR_REF",PS_DATA_F32, "reference star color", ref->Color); - psMetadataAdd (row, PS_LIST_TAIL, "CHIP_ID", PS_DATA_STRING, "chip identifier", chipName); + psMetadataAdd (row, PS_LIST_TAIL, "RA_REF", PS_DATA_F64, "right ascension (deg, J2000)", PM_DEG_RAD*ref->sky->r); + psMetadataAdd (row, PS_LIST_TAIL, "DEC_REF", PS_DATA_F64, "declination (deg, J2000)", PM_DEG_RAD*ref->sky->d); + psMetadataAdd (row, PS_LIST_TAIL, "X_CHIP_REF", PS_DATA_F32, "x fitted coord on chip", ref->chip->x); + psMetadataAdd (row, PS_LIST_TAIL, "Y_CHIP_REF", PS_DATA_F32, "y fitted coord on chip", ref->chip->y); + psMetadataAdd (row, PS_LIST_TAIL, "X_CHIP_RAW", PS_DATA_F32, "x coord on chip", raw->chip->x); + psMetadataAdd (row, PS_LIST_TAIL, "Y_CHIP_RAW", PS_DATA_F32, "y coord on chip", raw->chip->y); + psMetadataAdd (row, PS_LIST_TAIL, "X_FPA_RAW", PS_DATA_F32, "x coord on focal plane", raw->FP->x); + psMetadataAdd (row, PS_LIST_TAIL, "Y_FPA_RAW", PS_DATA_F32, "y coord on focal plane", raw->FP->y); + psMetadataAdd (row, PS_LIST_TAIL, "X_TPA_RAW", PS_DATA_F32, "x coord on focal plane", raw->TP->x); + psMetadataAdd (row, PS_LIST_TAIL, "Y_TPA_RAW", PS_DATA_F32, "y coord on focal plane", raw->TP->y); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_INST", PS_DATA_F32, "instrumental magnitude", raw->Mag); + psMetadataAdd (row, PS_LIST_TAIL, "MAG_REF", PS_DATA_F32, "reference star magnitude", ref->Mag); + psMetadataAdd (row, PS_LIST_TAIL, "COLOR_REF", PS_DATA_F32, "reference star color", ref->Color); + psMetadataAdd (row, PS_LIST_TAIL, "CHIP_ID", PS_DATA_STRING, "chip identifier", chipName); // XXX need to add the reference color, but this needs getstar / dvo.photcodes for the reference to be refined. Index: /trunk/psModules/src/objects/pmSourceVisual.c =================================================================== --- /trunk/psModules/src/objects/pmSourceVisual.c (revision 36833) +++ /trunk/psModules/src/objects/pmSourceVisual.c (revision 36834) @@ -545,7 +545,31 @@ psFree (model); + bool dumpData = false; + // pause and wait for user input: // continue, save (provide name), ?? - pmVisualAskUser(&plotPSF); +retry: + pmVisualAskUserOrDump(&plotPSF, &dumpData); + if (dumpData) { + char name[128]; + fprintf (stderr, "filename: "); + int status = fscanf (stdin, "%127s", name); + if (status != 1) { + fprintf (stderr, "odd response\n"); + goto retry; + } + + FILE *f = fopen (name, "w"); + if (!f) { + fprintf (stderr, "cannot open %s for output\n", name); + goto retry; + } + for (int i = 0; i < x->n; i++) { + float vModel = pmTrend2DEval (trend, x->data.F32[i], y->data.F32[i]); + fprintf (f, "%f %f %f %f %d\n", x->data.F32[i], y->data.F32[i], param->data.F32[i], vModel, mask->data.PS_TYPE_VECTOR_MASK_DATA[i]); + } + fclose (f); + goto retry; + } return true; Index: /trunk/psModules/src/psmodules.h =================================================================== --- /trunk/psModules/src/psmodules.h (revision 36833) +++ /trunk/psModules/src/psmodules.h (revision 36834) @@ -95,4 +95,5 @@ #include #include +#include // the following headers are from psModule:imcombine