Index: /trunk/psphot/doc/notes.20101222.txt =================================================================== --- /trunk/psphot/doc/notes.20101222.txt (revision 30624) +++ /trunk/psphot/doc/notes.20101222.txt (revision 30624) @@ -0,0 +1,18 @@ + +20101222 : EAM + + continuing work on psphot stack. Some notes: + + * visualization now has duplicate image display calls + * visual needs to do a better jobs scaling the smoothed, subtracted images. + * the fit radius for the convolved images should be a bit bigger than the unconvolved one + * *** fitRadius vs apRadius *** + I have been using 'markVal' to define the valid pixels for different operations + I think this has been inconsistent. + + ** Have pmSourceMagnitudes apply the defined apRadius on the fly + when measuring aperture photometry + + ** Set the fitRadius and keep it until redefined? + + ** RadialApertures does not use the marked region? Index: /trunk/psphot/doc/stack.txt =================================================================== --- /trunk/psphot/doc/stack.txt (revision 30623) +++ /trunk/psphot/doc/stack.txt (revision 30624) @@ -1,2 +1,92 @@ + +20101221 + + psphotStackReadout is now correctly subtracting the PSF models from + the images as it measures the radial aperture fluxes. + + Some issues: + + * the source properties get buggered up by the radial aperture flux + analysis (we need to re-determine the psf, re-create the source + pixels, and re-fit the models (linearly) to subtract them + correctly). + + * the standard analysis sequence is not doing a second pass + + * the choice of the best model is ambiguous + + * the radial aperture (and source addition / subtraction) is only + using the PSF model + + * make sure psphotEfficiency actually subtracts the inserted fake + sources + + * XXX how does the PSF-convolved model work with changing psf models? + +20101207 + + header fields for PSPS: + + * stack needs to count # of inputs + * propagate the input header to the output image (DONE) + * save PSF parameters in header + * stack type is only known to the launcher + + * what are we doing for the stack calibrations?? + +skycellID : code +surveyID : code +filterID : code +stackMetaID : stack_id +photoCalID : photcode -> number +magSat : FSATUR saturation magnitude level ** not correctly set +completMag : FLIMIT 95% completion level in mag ** not correctly set +stackTypeID : STK_TYPE stack type identifier ** deep stack, nightly stack, best IQ stack? +refImageID : identifier of image used as reference for analysis +subtrImageID : N/A (stack) identifier of image subtracted to generate difference image +analVer : analysis version index ** index for tess_id + skycell_id + filter? +nP2Images : NINPUTS number of P2 images contributing to this cell ** missing from input stack +astroScat : astrometric scatter for chip ** measure scatter on stack creation? +photoScat : photometric scatter for chip ** internal scatter? +nAstroRef : number of astrometric reference sources ** (connected to above) +nPhoRef : number of photometric reference sources ** same +psfFwhm : PSF full width at half maximum ** 0.5(FWHM_MAJ + FWHM_MIN) +psfmodelID : * PSFMODEL PSF model identifier ** save the PSF model name in the header +psfWidMajor : * FWHM_MAJ PSF parameters +psfWidMinor : * FWHM_MIN PSF parameters +psfTheta : * ANGLE PSF parameters +psfExtra1 : * PSF_EXT1 PSF parameters ** (at field center?) +psfExtra2 : * PSF_EXT2 PSF parameters ** (not set for all models) +photoZero : local derived photometric zero point +photoColor : local derived photometric color term +ctype1 : * CTYPE1 name of astrometric projection in RA ** propagate from input stacks +ctype2 : * CTYPE2 name of astrometric projection in DEC +crval1 : * CRVAL1 RA corresponding to reference pixel +crval2 : * CRVAL2 DEC corresponding to reference pixel +crpix1 : * CRPIX1 reference pixel value for RA +crpix2 : * CRPIX2 reference pixel value for DEC +cdelt1 : * CDELT1 scale factor for RA +cdelt2 : * CDELT2 scale factor for DEC +pc001001 : * PC001001 elements of rotation/Dcale matrix +pc001002 : * PC001002 elements of rotation/Dcale matrix +pc002001 : * PC002001 elements of rotation/Dcale matrix +pc002002 : * PC002002 elements of rotation/Dcale matrix +polyOrder : * NPLYTERM polynomial order of astrometry fit ** default to 1 +pca1x3y0 : * PCA1X3Y0 polynomial coefficients for the astrometric fit +pca1x2y1 : * PCA1X2Y1 polynomial coefficients for the astrometric fit +pca1x1y2 : * PCA1X1Y2 polynomial coefficients for the astrometric fit +pca1x0y3 : * PCA1X0Y3 polynomial coefficients for the astrometric fit +pca1x2y0 : * PCA1X2Y0 polynomial coefficients for the astrometric fit +pca1x1y1 : * PCA1X1Y1 polynomial coefficients for the astrometric fit +pca1x0y2 : * PCA1X0Y2 polynomial coefficients for the astrometric fit +pca2x3y0 : * PCA2X3Y0 polynomial coefficients for the astrometric fit +pca2x2y1 : * PCA2X2Y1 polynomial coefficients for the astrometric fit +pca2x1y2 : * PCA2X1Y2 polynomial coefficients for the astrometric fit +pca2x0y3 : * PCA2X0Y3 polynomial coefficients for the astrometric fit +pca2x2y0 : * PCA2X2Y0 polynomial coefficients for the astrometric fit +pca2x1y1 : * PCA2X1Y1 polynomial coefficients for the astrometric fit +pca2x0y2 : * PCA2X0Y2 polynomial coefficients for the astrometric fit +calibModNum : calibration modification number +dataRelease : Data release 20100506: Index: /trunk/psphot/src/Makefile.am =================================================================== --- /trunk/psphot/src/Makefile.am (revision 30623) +++ /trunk/psphot/src/Makefile.am (revision 30624) @@ -98,4 +98,5 @@ psphotStackMatchPSFsUtils.c \ psphotStackMatchPSFsPrepare.c \ + psphotStackMatchPSFsNext.c \ psphotStackOptions.c \ psphotStackObjects.c \ Index: /trunk/psphot/src/psphot.h =================================================================== --- /trunk/psphot/src/psphot.h (revision 30623) +++ /trunk/psphot/src/psphot.h (revision 30624) @@ -75,6 +75,6 @@ bool psphotImageQualityReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); -bool psphotChoosePSF (pmConfig *config, const pmFPAview *view, const char *filerule); -bool psphotChoosePSFReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); +bool psphotChoosePSF (pmConfig *config, const pmFPAview *view, const char *filerule, bool newSources); +bool psphotChoosePSFReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool newSources); bool psphotGuessModels (pmConfig *config, const pmFPAview *view, const char *filerule); @@ -260,5 +260,6 @@ bool psphotVisualShowSourceSize (pmReadout *readout, psArray *sources); bool psphotVisualPlotSourceSizeAlt (psMetadata *recipe, psMetadata *analysis, psArray *sources); -bool psphotVisualShowResidualImage (pmReadout *readout); +bool psphotVisualShowResidualImage (pmReadout *readout, bool reshow); +bool psphotVisualShowObjectRegions (pmReadout *readout, psMetadata *recipe, psArray *sources); bool psphotVisualPlotApResid (psArray *sources, float mean, float error, bool useApMag); bool psphotVisualPlotChisq (psArray *sources); @@ -312,4 +313,5 @@ int psphotFileruleCount(const pmConfig *config, const char *filerule); +bool psphotFileruleCountSet(const pmConfig *config, const char *filerule, int num); bool psphotAddKnownSources (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *inSources); @@ -403,4 +405,6 @@ bool psphotStackMatchPSFsReadout (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index); bool psphotStackMatchPSFsPrepare (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index); +bool psphotStackMatchPSFsNext (bool *smoothAgain, pmConfig *config, const pmFPAview *view, const char *filerule, int lastSize); +bool psphotStackMatchPSFsNextReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, psVector *fwhmValues, int lastSize); // psphotStackMatchPSFsUtils @@ -426,8 +430,8 @@ bool psphotRadialApertures (pmConfig *config, const pmFPAview *view, const char *filerule); bool psphotRadialAperturesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); -bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *radMax); +bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *radMax, int entry); bool psphotExtendedSourceAnalysisByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule); -bool psphotRadialAperturesByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule); +bool psphotRadialAperturesByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule, int nMatchedPSF); bool psphotStackObjectsUnifyPosition (psArray *objects); @@ -440,3 +444,13 @@ bool psphotCleanInputs (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotResetModels (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotResetModelsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); + +bool psphotRedefinePixels (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotRedefinePixelsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); + +bool psphotSourceChildren (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc); +bool psphotSourceChildrenReadout (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc, int index); +psArray *psphotSourceChildrenByObject (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objectsSrc); + #endif Index: /trunk/psphot/src/psphotBlendFit.c =================================================================== --- /trunk/psphot/src/psphotBlendFit.c (revision 30623) +++ /trunk/psphot/src/psphotBlendFit.c (revision 30624) @@ -187,5 +187,6 @@ psLogMsg ("psphot.psphotBlendFit", PS_LOG_INFO, "fit models: %f sec for %d objects (%d psf, %d ext, %d failed, %ld skipped)\n", psTimerMark ("psphot.fit.nonlinear"), Nfit, Npsf, Next, Nfail, sources->n - Nfit); - psphotVisualShowResidualImage (readout); + psphotVisualShowResidualImage (readout, false); + psphotVisualShowObjectRegions (readout, recipe, sources); psphotVisualShowFlags (sources); Index: /trunk/psphot/src/psphotCheckStarDistribution.c =================================================================== --- /trunk/psphot/src/psphotCheckStarDistribution.c (revision 30623) +++ /trunk/psphot/src/psphotCheckStarDistribution.c (revision 30624) @@ -34,5 +34,5 @@ pmSource *source = stars->data[i]; if (source->peak == NULL) continue; - if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue; + if (!(source->tmpFlags & PM_SOURCE_TMPF_CANDIDATE_PSFSTAR)) continue; int binX = psImageBinningGetRuffX (binning, source->peak->xf); @@ -85,5 +85,5 @@ if (source->peak == NULL) continue; if (source->moments == NULL) continue; - if (source->mode & PM_SOURCE_MODE_PSFSTAR) continue; + if (source->tmpFlags & PM_SOURCE_TMPF_CANDIDATE_PSFSTAR) continue; if (source->mode & PM_SOURCE_MODE_SATSTAR) continue; if (source->type != PM_SOURCE_TYPE_STAR) continue; @@ -97,5 +97,5 @@ if (y > Ye) continue; - source->mode |= PM_SOURCE_MODE_PSFSTAR; + source->tmpFlags |= PM_SOURCE_TMPF_CANDIDATE_PSFSTAR; psArrayAdd (stars, 200, source); Index: /trunk/psphot/src/psphotChoosePSF.c =================================================================== --- /trunk/psphot/src/psphotChoosePSF.c (revision 30623) +++ /trunk/psphot/src/psphotChoosePSF.c (revision 30624) @@ -2,5 +2,5 @@ // generate a PSF model for inputs without PSF models already loaded -bool psphotChoosePSF (pmConfig *config, const pmFPAview *view, const char *filerule) +bool psphotChoosePSF (pmConfig *config, const pmFPAview *view, const char *filerule, bool newSources) { bool status = true; @@ -19,5 +19,5 @@ for (int i = 0; i < num; i++) { if (i == chisqNum) continue; // skip chisq image - if (!psphotChoosePSFReadout (config, view, filerule, i, recipe)) { + if (!psphotChoosePSFReadout (config, view, filerule, i, recipe, newSources)) { psError (PSPHOT_ERR_CONFIG, false, "failed to choose a psf model for %s entry %d", filerule, i); return false; @@ -28,5 +28,5 @@ // try PSF models and select best option -bool psphotChoosePSFReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) { +bool psphotChoosePSFReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool newSources) { bool status; @@ -50,5 +50,5 @@ psAssert (detections, "missing detections?"); - psArray *sources = detections->newSources; + psArray *sources = newSources ? detections->newSources : detections->allSources; psAssert (sources, "missing sources?"); @@ -166,10 +166,8 @@ for (int i = 0; i < sources->n; i++) { pmSource *source = sources->data[i]; - if (source->mode & PM_SOURCE_MODE_PSFSTAR) { - // keep NSTARS PSF stars, unmark the rest + if (source->tmpFlags & PM_SOURCE_TMPF_CANDIDATE_PSFSTAR) { + // keep NSTARS PSF stars if (stars->n < NSTARS) { psArrayAdd (stars, 200, source); - } else { - source->mode &= ~PM_SOURCE_MODE_PSFSTAR; } } @@ -214,9 +212,5 @@ psFree(options); - // unset the PSFSTAR flags (none are used): - for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - source->mode &= ~PM_SOURCE_MODE_PSFSTAR; - } + // no sources are used as PSF stars // XXX set sxx, etc from FWHM in recipe @@ -292,9 +286,5 @@ psLogMsg ("psphot.pspsf", PS_LOG_INFO, "Using guess PSF model"); - // unset the PSFSTAR flags (none are used): - for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - source->mode &= ~PM_SOURCE_MODE_PSFSTAR; - } + // no sources are used as PSF stars // XXX set sxx, etc from FWHM in recipe @@ -322,45 +312,21 @@ pmPSFtry *try = models->data[bestN]; - // unset the PSFSTAR flag for stars not used for PSF model - // XXX a more efficient way of achieving this would be to record a pair of arrays - // of the source index and the source id for the psf stars. but that would require we do - // not re-sort the source list in the meanwhile - int nDrop = 0; + // set the PSFSTAR flag for stars used for the PSF model + int nKeep = 0; for (int i = 0; i < try->sources->n; i++) { pmSource *source = try->sources->data[i]; - if (try->mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) { - // need to find this source in the original list (these are copies, not pointers) - for (int j = 0; j < sources->n; j++) { - pmSource *realSource = sources->data[j]; - if (realSource->id != source->id) continue; - realSource->mode &= ~PM_SOURCE_MODE_PSFSTAR; - source->mode &= ~PM_SOURCE_MODE_PSFSTAR; - nDrop ++; - break; - } - } - } - // fprintf (stderr, "drop %d stars as PSF stars\n", nDrop); - - // XXX is this working? - // int N1 = 0; - // for (int i = 0; i < try->sources->n; i++) { - // pmSource *source = try->sources->data[i]; - // fprintf (stderr, "%llx : %d\n", (long long int) source, (source->mode & PM_SOURCE_MODE_PSFSTAR)); - // if (source->mode & PM_SOURCE_MODE_PSFSTAR) { - // N1 ++; - // } - // } - // int N2 = 0; - // for (int i = 0; i < sources->n; i++) { - // pmSource *source = sources->data[i]; - // fprintf (stderr, "%llx : %d\n", (long long int) source, (source->mode & PM_SOURCE_MODE_PSFSTAR)); - // if (source->mode & PM_SOURCE_MODE_PSFSTAR) { - // N2 ++; - // } - // } - // fprintf (stderr, "N1: %d, N2: %d\n", N1, N2); + if (try->mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue; + source->mode |= PM_SOURCE_MODE_PSFSTAR; + + // this source was used: find the parent and set its PSFSTAR flag + pmSource *realSource = source->parent; + psAssert (realSource, "pmPSFtryAlloc should have set the parent pointers"); + realSource->mode |= PM_SOURCE_MODE_PSFSTAR; + nKeep ++; + } + psLogMsg ("psphot.pspsf", PS_LOG_DETAIL, "used %d of %ld candidate PSF objects\n", nKeep, try->sources->n); // build a PSF residual image + // we need to use the 'try' set because they have the associated models defined if (!psphotMakeResiduals (try->sources, recipe, try->psf, maskVal)) { psError(PSPHOT_ERR_PSF, false, "Unable to construct residual table for PSF"); @@ -386,5 +352,5 @@ } - // XXX test dump of psf star data and psf-subtracted image + // test dump of psf star data and psf-subtracted image if (psTraceGetLevel("psphot.psfstars") > 5) { psphotDumpPSFStars (readout, try, options->fitRadius, maskVal, markVal); @@ -413,7 +379,6 @@ return false; } - psFree (psf); // XXX double-check - - // move into psphotChoosePSF + psFree (psf); + psphotVisualShowPSFModel (readout, psf); @@ -436,4 +401,6 @@ psVector *fwhmMajor = psVectorAllocEmpty (100, PS_DATA_F32); psVector *fwhmMinor = psVectorAllocEmpty (100, PS_DATA_F32); + psVector *psfExtra1 = psVectorAllocEmpty (100, PS_DATA_F32); + psVector *psfExtra2 = psVectorAllocEmpty (100, PS_DATA_F32); for (float ix = -0.4; ix <= +0.4; ix += 0.1) { @@ -459,13 +426,21 @@ shape.sxy = modelPSF->params->data.F32[PM_PAR_SXY]; axes = psEllipseShapeToAxes (shape, 20.0); - psFree (modelPSF); float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major); if (!isfinite(FWHM_MAJOR) || !isfinite(FWHM_MINOR)) { fprintf (stderr, "!"); + psFree (modelPSF); continue; } psVectorAppend (fwhmMajor, FWHM_MAJOR); psVectorAppend (fwhmMinor, FWHM_MINOR); + + if (modelPSF->params->n >= 7) { + psVectorAppend (psfExtra1, modelPSF->params->data.F32[7]); + } + if (modelPSF->params->n >= 8) { + psVectorAppend (psfExtra2, modelPSF->params->data.F32[8]); + } + psFree (modelPSF); } } @@ -474,5 +449,5 @@ if (!psVectorStats (stats, fwhmMajor, NULL, NULL, 0)) { psError(PS_ERR_UNKNOWN, false, "failure to measure stats for FWHM MAJOR"); - return false; + goto escape; } @@ -486,5 +461,5 @@ if (!psVectorStats (stats, fwhmMinor, NULL, NULL, 0)) { psError(PS_ERR_UNKNOWN, false, "failure to measure stats for FWHM MINOR"); - return false; + goto escape; } psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FWHM_MIN", PS_META_REPLACE, "PSF FWHM Minor axis (mean)", stats->sampleMean); @@ -494,19 +469,64 @@ float fwhmMin = stats->sampleMean; // FWHM on minor axis + if (readout->parent) { + + // we now have 2 possible measurements of the seeing : the PSF model based version and + // the Moments based value we need to define a definitive "CHIP.SEEING" value. This + // value is used by the PSF-matching programs (ppSub and ppStack). The moments-based + // value is probably more representative of the value needed by those tools (it is also + // more stable?) + bool status = false; + float fwhmMajorMoments = psMetadataLookupF32(&status, readout->analysis, "IQ_FW1"); + float fwhmMinorMoments = psMetadataLookupF32(&status, readout->analysis, "IQ_FW1"); + pmChip *chip = readout->parent->parent; // Parent chip psAssert(chip, "Cell should be attached to a chip."); psMetadataItem *item = psMetadataLookup(chip->concepts, "CHIP.SEEING"); // Item with chip - item->data.F32 = 0.5 * (fwhmMaj + fwhmMin); + item->data.F32 = 0.5 * (fwhmMajorMoments + fwhmMinorMoments); + + psLogMsg ("psphot", PS_LOG_DETAIL, "fwhm (psf): %f,%f (moments): %f,%f", fwhmMaj, fwhmMin, fwhmMajorMoments, fwhmMinorMoments); + } + + if (psfExtra1->n) { + if (!psVectorStats (stats, psfExtra1, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats for PSF EXTRA 1"); + goto escape; + } + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PSF_EXT1", PS_META_REPLACE, "PSF extra param 1", stats->sampleMean); + } + + if (psfExtra2->n) { + if (!psVectorStats (stats, psfExtra2, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats for PSF EXTRA 2"); + goto escape; + } + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PSF_EXT2", PS_META_REPLACE, "PSF extra param 2", stats->sampleMean); } psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "ANGLE", PS_META_REPLACE, "PSF angle", axes.theta); psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "NPSFSTAR", PS_META_REPLACE, "Number of stars used to make PSF", psf->nPSFstars); - psMetadataAddBool(readout->analysis, PS_LIST_TAIL, "PSFMODEL", PS_META_REPLACE, "Valid PSF Model?", true); + + char *psfModelName = pmModelClassGetName(psf->type); + psMetadataAddStr(readout->analysis, PS_LIST_TAIL, "PSFMODEL", PS_META_REPLACE, "PSF Model Name", psfModelName); + psMetadataAddBool(readout->analysis, PS_LIST_TAIL, "PSF_OK", PS_META_REPLACE, "Valid PSF Model?", true); + + int nParams = pmModelClassParameterCount(psf->type); + psMetadataAddS32(readout->analysis, PS_LIST_TAIL, "PSF_NPAR", PS_META_REPLACE, "Number of PSF parameters", nParams); psFree (fwhmMajor); psFree (fwhmMinor); + psFree (psfExtra1); + psFree (psfExtra2); psFree (stats); return true; + +escape: + psFree (fwhmMajor); + psFree (fwhmMinor); + psFree (psfExtra1); + psFree (psfExtra2); + psFree (stats); + return false; } @@ -568,7 +588,8 @@ psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MN_LQ", PS_META_REPLACE, "PSF FWHM Minor axis (lower quartile)", 0); psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MN_UQ", PS_META_REPLACE, "PSF FWHM Minor axis (upper quartile)", 0); - psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "ANGLE", PS_META_REPLACE, "PSF angle", FWHM_T); - psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "NPSFSTAR", PS_META_REPLACE, "Number of stars used to make PSF", 0); - psMetadataAddBool(readout->analysis, PS_LIST_TAIL, "PSFMODEL", PS_META_REPLACE, "Valid PSF Model?", false); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "ANGLE", PS_META_REPLACE, "PSF angle", FWHM_T); + psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "NPSFSTAR", PS_META_REPLACE, "Number of stars used to make PSF", 0); + psMetadataAddStr(readout->analysis, PS_LIST_TAIL, "PSFMODEL", PS_META_REPLACE, "PSF Model Name", "NONE"); + psMetadataAddBool(readout->analysis, PS_LIST_TAIL, "PSF_OK", PS_META_REPLACE, "Valid PSF Model?", false); return true; Index: /trunk/psphot/src/psphotExtendedSourceAnalysis.c =================================================================== --- /trunk/psphot/src/psphotExtendedSourceAnalysis.c (revision 30623) +++ /trunk/psphot/src/psphotExtendedSourceAnalysis.c (revision 30624) @@ -173,5 +173,5 @@ psLogMsg ("psphot", PS_LOG_INFO, " %d kron\n", Nkron); - psphotVisualShowResidualImage (readout); + psphotVisualShowResidualImage (readout, false); if (doPetrosian) { Index: /trunk/psphot/src/psphotExtendedSourceAnalysisByObject.c =================================================================== --- /trunk/psphot/src/psphotExtendedSourceAnalysisByObject.c (revision 30623) +++ /trunk/psphot/src/psphotExtendedSourceAnalysisByObject.c (revision 30624) @@ -53,4 +53,7 @@ pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); psAssert (readout, "missing readout?"); + + psLogMsg("psphot", PS_LOG_INFO, "petrosians for image %d", i); + psphotVisualShowImage(readout); readouts->data[i] = psMemIncrRefCounter(readout); Index: /trunk/psphot/src/psphotExtendedSourceFits.c =================================================================== --- /trunk/psphot/src/psphotExtendedSourceFits.c (revision 30623) +++ /trunk/psphot/src/psphotExtendedSourceFits.c (revision 30624) @@ -18,6 +18,11 @@ int num = psphotFileruleCount(config, filerule); + // skip the chisq image (optionally?) + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + // loop over the available readouts for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; // skip chisq image if (!psphotExtendedSourceFitsReadout (config, view, filerule, i, recipe)) { psError (PSPHOT_ERR_CONFIG, false, "failed on to fit extended sources for %s entry %d", filerule, i); @@ -37,4 +42,5 @@ int Nplain = 0; int NplainPass = 0; + int Nfaint = 0; psTimerStart ("psphot.extended"); @@ -46,4 +52,7 @@ pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); psAssert (readout, "missing readout?"); + + psLogMsg("psphot", PS_LOG_INFO, "extended source fits for image %d", index); + psphotVisualShowImage(readout); pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); @@ -167,4 +176,5 @@ PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nplain PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for NplainPass + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfain if (false && !psThreadJobAddPending(job)) { @@ -189,4 +199,6 @@ scalar = job->args->data[11]; NplainPass += scalar->data.S32; + scalar = job->args->data[12]; + Nfaint += scalar->data.S32; psFree(job); } @@ -217,4 +229,6 @@ scalar = job->args->data[11]; NplainPass += scalar->data.S32; + scalar = job->args->data[12]; + Nfaint += scalar->data.S32; } psFree(job); @@ -227,4 +241,5 @@ psLogMsg ("psphot", PS_LOG_INFO, " %d convolved models (%d passed)\n", Nconvolve, NconvolvePass); psLogMsg ("psphot", PS_LOG_INFO, " %d plain models (%d passed)\n", Nplain, NplainPass); + psLogMsg ("psphot", PS_LOG_INFO, " %d too faint to fit\n", Nfaint); return true; } @@ -238,4 +253,5 @@ int NconvolvePass = 0; int Nplain = 0; + int Nfaint = 0; int NplainPass = 0; bool savePics = false; @@ -326,9 +342,9 @@ // array to store the pointers to the model flux images while the models are being fitted - psArray *modelFluxes = psArrayAllocEmpty (4); + psArray *modelFluxes = psArrayAllocEmpty (models->list->n); // allocate the array to store the model fits if (source->modelFits == NULL) { - source->modelFits = psArrayAllocEmpty (4); + source->modelFits = psArrayAllocEmpty (models->list->n); } @@ -350,5 +366,8 @@ // limit selection to some SN limit assert (source->peak); // how can a source not have a peak? - if (source->peak->SN < SNlim) continue; + if (source->peak->SN < SNlim) { + Nfaint ++; + continue; + } // check on the model type @@ -497,4 +516,7 @@ scalar->data.S32 = NplainPass; + scalar = job->args->data[12]; + scalar->data.S32 = Nfaint; + return true; } Index: /trunk/psphot/src/psphotFindDetections.c =================================================================== --- /trunk/psphot/src/psphotFindDetections.c (revision 30623) +++ /trunk/psphot/src/psphotFindDetections.c (revision 30624) @@ -86,9 +86,9 @@ psImage *significance = psphotSignificanceImage (readout, recipe, maskVal); - // display the significance image - psphotVisualShowSignificance (significance, -1.0, PS_SQR(3.0*NSIGMA_PEAK)); - // display the log significance image psphotVisualShowLogSignificance (significance, 0.0, 4.5); + + // display the significance image + psphotVisualShowSignificance (significance, 0.98*threshold, 1.02*threshold); // detect the peaks in the significance image @@ -108,4 +108,6 @@ } + // XXX do a second (or third?) pass with rebinning (to detected more extended sources) + psFree (significance); Index: /trunk/psphot/src/psphotFindFootprints.c =================================================================== --- /trunk/psphot/src/psphotFindFootprints.c (revision 30623) +++ /trunk/psphot/src/psphotFindFootprints.c (revision 30624) @@ -21,6 +21,7 @@ psArray *footprints = pmFootprintsFind (significance, threshold, npixMin); - pmFootprintsAssignPeaks(footprints, detections->peaks); - // XXX handle the error conditions here + if (pmFootprintsAssignPeaks(footprints, detections->peaks) != PS_ERR_NONE) { + psAbort ("inconsistent peaks and footprints"); + } // footprints now owns the peaks; after culling (below), we will rebuild the peaks array Index: /trunk/psphot/src/psphotFitSourcesLinear.c =================================================================== --- /trunk/psphot/src/psphotFitSourcesLinear.c (revision 30623) +++ /trunk/psphot/src/psphotFitSourcesLinear.c (revision 30624) @@ -46,4 +46,8 @@ return false; } + + psphotVisualShowResidualImage (readout, (num > 0)); + psphotVisualShowPeaks (detections); + psphotVisualShowObjectRegions (readout, recipe, sources); } return true; @@ -151,4 +155,23 @@ if (x > AnalysisRegion.x1) continue; if (y > AnalysisRegion.y1) continue; + + // check the integral of the model : is it large enough? + float modelSum = 0.0; + for (int iy = 0; iy < source->modelFlux->numRows; iy++) { + for (int ix = 0; ix < source->modelFlux->numCols; ix++) { + modelSum += source->modelFlux->data.F32[iy][ix]; + } + } + if (modelSum < 0.5) continue; // skip sources with no model constraint (somewhat arbitrary limit) + // if (modelSum < 0.01) continue; // skip sources with no model constraint (somewhat arbitrary limit) + if (modelSum < 0.8) { + fprintf (stderr, "low-sig model @ %f, %f (%f sum, %f peak)\n", + source->peak->xf, source->peak->yf, modelSum, source->peak->flux); + } + + pmModel *model = pmSourceGetModel (NULL, source); + + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); + psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, model->fitRadius, "OR", markVal); source->mode |= PM_SOURCE_MODE_LINEAR_FIT; @@ -270,5 +293,4 @@ model->params->data.F32[PM_PAR_I0] = norm->data.F32[i]; model->dparams->data.F32[PM_PAR_I0] = errors->data.F32[i]; - // XXX is the value of 'errors' modified by the sky fit? // subtract object @@ -297,5 +319,4 @@ psLogMsg ("psphot.ensemble", PS_LOG_INFO, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot.linear")); - psphotVisualShowResidualImage (readout); psphotVisualPlotChisq (sources); // psphotVisualShowFlags (sources); Index: /trunk/psphot/src/psphotFitSourcesLinearStack.c =================================================================== --- /trunk/psphot/src/psphotFitSourcesLinearStack.c (revision 30623) +++ /trunk/psphot/src/psphotFitSourcesLinearStack.c (revision 30624) @@ -176,5 +176,4 @@ psLogMsg ("psphot.ensemble", PS_LOG_INFO, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot.linear")); - return true; } Index: /trunk/psphot/src/psphotImageQuality.c =================================================================== --- /trunk/psphot/src/psphotImageQuality.c (revision 30623) +++ /trunk/psphot/src/psphotImageQuality.c (revision 30624) @@ -79,5 +79,5 @@ // select by PSFSTAR mode? // ?? - if (source->type != PM_SOURCE_TYPE_STAR || !(source->mode & PM_SOURCE_MODE_PSFSTAR)) { + if (source->type != PM_SOURCE_TYPE_STAR || !(source->tmpFlags & PM_SOURCE_TMPF_CANDIDATE_PSFSTAR)) { psTrace("psphot", 10, "Ignoring source for image quality because not a good star"); continue; @@ -133,4 +133,12 @@ if (num == 0) { psLogMsg ("psphot", PS_LOG_INFO, "no valid sources for image quality, skipping"); + psFree(FWHM_MAJOR); + psFree(FWHM_MINOR); + psFree(M2); + psFree(M2c); + psFree(M2s); + psFree(M3); + psFree(M4); + return true; } Index: /trunk/psphot/src/psphotMakePSFReadout.c =================================================================== --- /trunk/psphot/src/psphotMakePSFReadout.c (revision 30623) +++ /trunk/psphot/src/psphotMakePSFReadout.c (revision 30624) @@ -60,5 +60,5 @@ // Use bright stellar objects to measure PSF. If we do not have enough stars to generate // the PSF, build one from the SEEING guess and model class - if (!psphotChoosePSF (config, view, filerule)) { + if (!psphotChoosePSF (config, view, filerule, true)) { psLogMsg ("psphot", 3, "failure to construct a psf model"); return psphotReadoutCleanup (config, view, filerule); Index: /trunk/psphot/src/psphotMakeResiduals.c =================================================================== --- /trunk/psphot/src/psphotMakeResiduals.c (revision 30623) +++ /trunk/psphot/src/psphotMakeResiduals.c (revision 30624) @@ -136,4 +136,7 @@ psFree (interp); } + xSize = PS_MIN(xSize, 2*radiusMax+3); + ySize = PS_MIN(ySize, 2*radiusMax+3); + pmResiduals *resid = pmResidualsAlloc (xSize, ySize, xBin, yBin); psImageInit (resid->mask, 0); @@ -316,5 +319,5 @@ psFree (B); - psLogMsg ("psphot.pspsf", PS_LOG_MINUTIA, "generate residuals for %ld objects: %f sec\n", input->n, psTimerMark ("psphot.residuals")); + psLogMsg ("psphot.pspsf", PS_LOG_DETAIL, "generate residuals for %ld objects: %f sec\n", input->n, psTimerMark ("psphot.residuals")); psFree (xC); Index: /trunk/psphot/src/psphotMergeSources.c =================================================================== --- /trunk/psphot/src/psphotMergeSources.c (revision 30623) +++ /trunk/psphot/src/psphotMergeSources.c (revision 30624) @@ -513,5 +513,236 @@ } - return true; -} - + // loop over the sources, redefine their pixels to point at the new filerule image, + // copy the source data, and add a reference back to the original source + + + return true; +} + +// create source children from ruleSrc for ruleOut +bool psphotSourceChildren (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc) +{ + bool status = true; + + int num = psphotFileruleCount(config, ruleSrc); + + // skip the chisq image because it is a duplicate of the detection version + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; // skip chisq image + if (!psphotSourceChildrenReadout (config, view, ruleOut, ruleSrc, i)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to copy sources from %s to %s entry %d", ruleSrc, ruleOut, i); + return false; + } + } + return true; +} + +// create source children from ruleSrc for ruleOut for this entry. XXX currently, this is only +// used by psphotStackReadout (sources go on allSources so that psphotChoosePSF can be called +// repeatedly) +bool psphotSourceChildrenReadout (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc, int index) { + + bool status; + + // find the currently selected readout + pmFPAfile *fileSrc = pmFPAfileSelectSingle(config->files, ruleSrc, index); // File of interest + psAssert (fileSrc, "missing file?"); + + pmReadout *readoutSrc = pmFPAviewThisReadout(view, fileSrc->fpa); + psAssert (readoutSrc, "missing readout?"); + + pmDetections *detectionsSrc = psMetadataLookupPtr (&status, readoutSrc->analysis, "PSPHOT.DETECTIONS"); + psAssert (detectionsSrc, "missing detections?"); + + psArray *sourcesSrc = detectionsSrc->allSources; + psAssert (sourcesSrc, "missing sources?"); + + // find the currently selected readout + pmFPAfile *fileOut = pmFPAfileSelectSingle(config->files, ruleOut, index); // File of interest + psAssert (fileOut, "missing file?"); + + pmReadout *readoutOut = pmFPAviewThisReadout(view, fileOut->fpa); + psAssert (readoutOut, "missing readout?"); + + // generate a new detection structure for the output filerule + pmDetections *detectionsOut = psMetadataLookupPtr (&status, readoutOut->analysis, "PSPHOT.DETECTIONS"); + if (!detectionsOut) { + detectionsOut = pmDetectionsAlloc(); + detectionsOut->allSources = psArrayAllocEmpty (100); + // save detections on the readout->analysis + if (!psMetadataAddPtr (readoutOut->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot detections", detectionsOut)) { + psError (PSPHOT_ERR_CONFIG, false, "problem saving detections on readout"); + return false; + } + psFree(detectionsOut); // a copy remains on the analysis metadata + } + + // loop over the sources, redefine their pixels to point at the new filerule image, + // copy the source data, and add a reference back to the original source + + // copy the sources from sourceSrcs to the new detection structure + for (int i = 0; i < sourcesSrc->n; i++) { + pmSource *sourceSrc = sourcesSrc->data[i]; + + pmSource *sourceOut = pmSourceCopy(sourceSrc); + sourceOut->parent = sourceSrc; + + // keep the original source flags + sourceOut->type = sourceSrc->type; + sourceOut->mode = sourceSrc->mode; + sourceOut->mode2 = sourceSrc->mode2; + sourceOut->tmpFlags = sourceSrc->tmpFlags; + + // does this copy all model data? (NO) + sourceOut->modelPSF = pmModelCopy(sourceSrc->modelPSF); + sourceOut->modelEXT = pmModelCopy(sourceSrc->modelEXT); + + if (sourceSrc->modelFits) { + sourceOut->modelFits = psArrayAlloc(sourceSrc->modelFits->n); + for (int j = 0; j < sourceSrc->modelFits->n; j++) { + sourceOut->modelFits->data[j] = pmModelCopy(sourceSrc->modelFits->data[j]); + } + } + + // drop the references to the original image pixels: + pmSourceFreePixels (sourceOut); + + // allocate image, weight, mask for the new image for each peak + pmSourceRedefinePixels (sourceOut, readoutOut, sourceOut->peak->x, sourceOut->peak->y, sourceOut->modelPSF->fitRadius); + + // child sources have not been subtracted in this image, but this flag may be raised if + // they were subtracted in the parent's image + sourceOut->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED; + + psArrayAdd (detectionsOut->allSources, 100, sourceOut); + psFree (sourceOut); + } + psLogMsg ("psphot", 3, "%ld known sources supplied", detectionsOut->allSources->n); + + + return true; +} + +// create source children associated with 'filerule' from the objectsSrc. XXX currently, this +// is only used by psphotStackReadout (sources go on allSources so that psphotChoosePSF can be +// called repeatedly) +psArray *psphotSourceChildrenByObject (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objectsSrc) { + + bool status; + + int nImages = psphotFileruleCount(config, filerule); + + // generate look-up arrays for detections and readouts + psArray *detArrays = psArrayAlloc(nImages); + psArray *readouts = psArrayAlloc(nImages); + + for (int i = 0; i < nImages; i++) { + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + if (!detections) { + detections = pmDetectionsAlloc(); + detections->allSources = psArrayAllocEmpty (100); + detections->peaks = psArrayAllocEmpty (100); + // save detections on the readout->analysis + if (!psMetadataAddPtr (readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot detections", detections)) { + psError (PSPHOT_ERR_CONFIG, false, "problem saving detections on readout"); + return NULL; + } + psFree(detections); // a copy remains on the analysis metadata + psAssert (detections, "missing detections?"); + } + detArrays->data[i] = psMemIncrRefCounter(detections); + readouts->data[i] = psMemIncrRefCounter(readout); + } + + psArray *objectsOut = psArrayAlloc(objectsSrc->n); + + // copy all sources for each object + for (int k = 0; k < objectsSrc->n; k++) { + + pmPhotObj *objectSrc = objectsSrc->data[k]; + if (!objectSrc) continue; + if (!objectSrc->sources) continue; + + pmPhotObj *objectOut = pmPhotObjAlloc(); + objectsOut->data[k] = objectOut; + + objectOut->SN = objectSrc->SN; + objectOut->x = objectSrc->x; + objectOut->y = objectSrc->y; + + objectOut->sources = psArrayAlloc(objectSrc->sources->n); + + // copy the sources from sourceSrcs to the new detection structure + // loop over the sources, redefine their pixels to point at the new filerule image, + // copy the source data, and add a reference back to the original source + for (int i = 0; i < objectSrc->sources->n; i++) { + + pmSource *sourceSrc = objectSrc->sources->data[i]; + + pmSource *sourceOut = pmSourceCopy(sourceSrc); + objectOut->sources->data[i] = sourceOut; + + sourceOut->parent = sourceSrc; + + // keep the original source flags + sourceOut->type = sourceSrc->type; + sourceOut->mode = sourceSrc->mode; + sourceOut->mode2 = sourceSrc->mode2; + sourceOut->tmpFlags = sourceSrc->tmpFlags; + + // does this copy all model data? (NO) + sourceOut->modelPSF = pmModelCopy(sourceSrc->modelPSF); + sourceOut->modelEXT = pmModelCopy(sourceSrc->modelEXT); + + if (sourceSrc->modelFits) { + sourceOut->modelFits = psArrayAlloc(sourceSrc->modelFits->n); + for (int j = 0; j < sourceSrc->modelFits->n; j++) { + sourceOut->modelFits->data[j] = pmModelCopy(sourceSrc->modelFits->data[j]); + } + } + + // drop the references to the original image pixels: + pmSourceFreePixels (sourceOut); + + // set the output readotu + int index = sourceOut->imageID; + if (index >= readouts->n) continue; // skip the sources generated by the chisq image + pmReadout *readout = readouts->data[index]; + + // allocate image, weight, mask for the new image for each peak + pmSourceRedefinePixels (sourceOut, readout, sourceOut->peak->x, sourceOut->peak->y, sourceOut->modelPSF->fitRadius); + + // child sources have not been subtracted in this image, but this flag may be raised if + // they were subtracted in the parent's image + sourceOut->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED; + + // set the output detections: + pmDetections *detectionsOut = detArrays->data[index]; + psArrayAdd (detectionsOut->allSources, 100, sourceOut); + psArrayAdd (detectionsOut->peaks, 100, sourceOut->peak); + } + } + + for (int i = 0; i < nImages; i++) { + pmDetections *detections = detArrays->data[i]; + psLogMsg ("psphot", 3, "%ld source children for image %d", detections->allSources->n, i); + } + + psFree (detArrays); + psFree (readouts); + + return objectsOut; +} + Index: /trunk/psphot/src/psphotOutput.c =================================================================== --- /trunk/psphot/src/psphotOutput.c (revision 30623) +++ /trunk/psphot/src/psphotOutput.c (revision 30624) @@ -15,4 +15,16 @@ psFree (name); return num; +} + +// convert filerule to filerule.NUM and look up in the config->arguments metadata +bool psphotFileruleCountSet(const pmConfig *config, const char *filerule, int num) { + + psString name = NULL; + psStringAppend(&name, "%s.NUM", filerule); + + bool status = psMetadataAddS32(config->arguments, PS_LIST_TAIL, name, PS_META_REPLACE, "", num); + psFree (name); + + return status; } @@ -261,4 +273,7 @@ psMetadataItemSupplement (&status, header, analysis, "ANGLE"); + psMetadataItemSupplement (&status, header, analysis, "PSFMODEL"); + psMetadataItemSupplement (&status, header, analysis, "PSF_OK"); + // Image Quality measurements psMetadataItemSupplement (&status, header, analysis, "IQ_NSTAR"); Index: /trunk/psphot/src/psphotPetrosianAnalysis.c =================================================================== --- /trunk/psphot/src/psphotPetrosianAnalysis.c (revision 30623) +++ /trunk/psphot/src/psphotPetrosianAnalysis.c (revision 30624) @@ -64,5 +64,5 @@ } - psphotVisualShowResidualImage (readout); + psphotVisualShowResidualImage (readout, false); return true; } Index: /trunk/psphot/src/psphotPetrosianProfile.c =================================================================== --- /trunk/psphot/src/psphotPetrosianProfile.c (revision 30623) +++ /trunk/psphot/src/psphotPetrosianProfile.c (revision 30624) @@ -68,5 +68,4 @@ // XXX this will only work in the psphot context, not the psphotPetrosianStudy... // XXX add the petrosian to the pmSource structure... - // psphotVisualShowResidualImage (readout); psphotVisualShowPetrosian (source, petrosian); Index: /trunk/psphot/src/psphotRadialApertures.c =================================================================== --- /trunk/psphot/src/psphotRadialApertures.c (revision 30623) +++ /trunk/psphot/src/psphotRadialApertures.c (revision 30624) @@ -102,4 +102,9 @@ if (source->peak->x > AnalysisRegion.x1) continue; if (source->peak->y > AnalysisRegion.y1) continue; + + // allocate pmSourceExtendedParameters, if not already defined + if (!source->radialAper) { + source->radialAper = psArrayAlloc(1); + } // replace object in image @@ -116,5 +121,5 @@ pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, 1.5*radius); - if (!psphotRadialApertureSource (source, recipe, skynoise, maskVal, radMax)) { + if (!psphotRadialApertureSource (source, recipe, skynoise, maskVal, radMax, 0)) { psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); } else { @@ -130,12 +135,122 @@ } -bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *radMax) { - - // allocate pmSourceExtendedParameters, if not already defined - if (!source->radial) { - source->radial = pmSourceRadialAperturesAlloc (); +bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *aperRadii, int entry) { + + // if we are a child source, save the results to the parent source radial aperture array + psArray *radialAperSet = source->radialAper; + if (source->parent) { + radialAperSet = source->parent->radialAper; + } + psAssert(radialAperSet, "this should be defined before calling"); + psAssert(radialAperSet->data[entry] == NULL, "why is this already defined?"); + + pmSourceRadialApertures *radialAper = pmSourceRadialAperturesAlloc (); + radialAperSet->data[entry] = radialAper; + + // storage for the derived pixel values + psVector *pixRadius2 = psVectorAllocEmpty(100, PS_TYPE_F32); + psVector *pixFlux = psVectorAllocEmpty(100, PS_TYPE_F32); + psVector *pixVar = psVectorAllocEmpty(100, PS_TYPE_F32); + + // outer-most radius for initial truncation + float Rmax = aperRadii->data.F32[aperRadii->n - 1]; + float Rmax2 = PS_SQR(Rmax); + + // store the R^2 values for the apertures + psVector *aperRadii2 = psVectorAlloc(aperRadii->n, PS_TYPE_F32); + for (int i = 0; i < aperRadii->n; i++) { + aperRadii2->data.F32[i] = PS_SQR(aperRadii->data.F32[i]); + } + + // center of the apertures + float xCM = source->moments->Mx - 0.5 - source->pixels->col0; // coord of peak in subimage + float yCM = source->moments->My - 0.5 - source->pixels->row0; // coord of peak in subimage + + // one pass through the pixels to select the valid pixels and calculate R^2 + for (int iy = 0; iy < source->pixels->numRows; iy++) { + + float yDiff = iy - yCM; + if (fabs(yDiff) > Rmax) continue; + + float *vPix = source->pixels->data.F32[iy]; + float *vWgt = source->variance->data.F32[iy]; + psImageMaskType *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy]; + + for (int ix = 0; ix < source->pixels->numCols; ix++, vPix++, vWgt++) { + + if (vMsk) { + if (*vMsk & maskVal) { + vMsk++; + continue; + } + vMsk++; + } + if (isnan(*vPix)) continue; + + float xDiff = ix - xCM; + if (fabs(xDiff) > Rmax) continue; + + // radius is just a function of (xDiff, yDiff) + float r2 = PS_SQR(xDiff) + PS_SQR(yDiff); + if (r2 > Rmax2) continue; + + psVectorAppend(pixRadius2, r2); + psVectorAppend(pixFlux, *vPix); + psVectorAppend(pixVar, *vWgt); + } + } + + psVector *flux = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin + psVector *fluxErr = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin + psVector *fill = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin + + psVectorInit (flux, 0.0); + psVectorInit (fluxErr, 0.0); + psVectorInit (fill, 0.0); + + float *rPix2 = pixRadius2->data.F32; + for (int i = 0; i < pixRadius2->n; i++, rPix2++) { + + float *aRad2 = aperRadii2->data.F32; + for (int j = 0; (*rPix2 < *aRad2) && (j < aperRadii2->n); j++, aRad2++) { + flux->data.F32[j] += pixFlux->data.F32[i]; + fluxErr->data.F32[j] += pixVar->data.F32[i]; + fill->data.F32[j] += 1.0; + } + } + + for (int i = 0; i < flux->n; i++) { + // calculate the total flux for bin 'nOut' + float Area = M_PI*aperRadii2->data.F32[i]; + fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]); + fill->data.F32[i] /= Area; + psTrace ("psphot", 5, "radial bins: %3d %5.1f : %8.1f +/- %7.2f : %4.2f %6.1f\n", + i, aperRadii->data.F32[i], flux->data.F32[i], fluxErr->data.F32[i], fill->data.F32[i], Area); } - psVector *radius = psVectorAllocEmpty(100, PS_TYPE_F32); + radialAper->flux = flux; + radialAper->fluxErr = fluxErr; + radialAper->fill = fill; + + psFree (aperRadii2); + psFree (pixRadius2); + psFree (pixFlux); + psFree (pixVar); + + return true; +} + +static int nCalls = 0; +static int nPass = 0; +static int nPix = 0; + +bool psphotRadialApertureSource_With_Sort (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *radMax, int entry) { + + psAssert(source->radialAper->data[entry] == NULL, "why is this already defined?"); + + pmSourceRadialApertures *radialAper = pmSourceRadialAperturesAlloc (); + source->radialAper->data[entry] = radialAper; + + psVector *pixRadius = psVectorAllocEmpty(100, PS_TYPE_F32); psVector *pixFlux = psVectorAllocEmpty(100, PS_TYPE_F32); psVector *pixVar = psVectorAllocEmpty(100, PS_TYPE_F32); @@ -144,5 +259,5 @@ for (int ix = 0; ix < source->pixels->numCols; ix++) { - // 0.5 PIX: get radius as a function of pixel coord + // 0.5 PIX: get pixRadius as a function of pixel coord float x = ix + 0.5 - source->peak->xf + source->pixels->col0; float y = iy + 0.5 - source->peak->yf + source->pixels->row0; @@ -150,10 +265,12 @@ float r = hypot(x, y); - psVectorAppend(radius, r); + psVectorAppend(pixRadius, r); psVectorAppend(pixFlux, source->pixels->data.F32[iy][ix]); psVectorAppend(pixVar, source->variance->data.F32[iy][ix]); - } - } - psphotRadialAperturesSortFlux(radius, pixFlux, pixVar); + nPix ++; + // if (nPix % 10000 == 0) {fprintf (stderr, "?");} + } + } + psphotRadialAperturesSortFlux(pixRadius, pixFlux, pixVar); psVector *flux = psVectorAllocEmpty(radMax->n, PS_TYPE_F32); // surface brightness of radial bin @@ -174,6 +291,6 @@ // XXX assume (or enforce) that the bins are contiguous and non-overlapping (Rmax[i] = Rmin[i+1]) - for (int i = 0; !done && (i < radius->n); i++) { - if (radius->data.F32[i] > Rmax) { + for (int i = 0; !done && (i < pixRadius->n); i++) { + if (pixRadius->data.F32[i] > Rmax) { // calculate the total flux for bin 'nOut' float Area = M_PI*PS_SQR(Rmax); @@ -185,4 +302,7 @@ nOut, radMax->data.F32[nOut], flux->data.F32[nOut], fluxErr->data.F32[nOut], fill->data.F32[nOut], Area); + nPass ++; + // if (nPass % 1000 == 0) {fprintf (stderr, "!");} + nOut ++; if (nOut >= radMax->n) break; @@ -195,16 +315,14 @@ flux->n = fluxErr->n = fill->n = nOut; - psFree(source->radial->flux); - psFree(source->radial->fluxErr); - psFree(source->radial->fill); - - source->radial->flux = flux; - source->radial->fluxErr = fluxErr; - source->radial->fill = fill; - - psFree (radius); + radialAper->flux = flux; + radialAper->fluxErr = fluxErr; + radialAper->fill = fill; + + psFree (pixRadius); psFree (pixFlux); psFree (pixVar); + nCalls ++; + // if (nCalls % 100 == 0) {fprintf (stderr, "*");} return true; } Index: /trunk/psphot/src/psphotRadialAperturesByObject.c =================================================================== --- /trunk/psphot/src/psphotRadialAperturesByObject.c (revision 30623) +++ /trunk/psphot/src/psphotRadialAperturesByObject.c (revision 30624) @@ -3,5 +3,14 @@ // aperture-like measurements for extended sources // flux in simple, circular apertures -bool psphotRadialAperturesByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule) { + +// **** it looks like this function will re-point the source pixels at the specified FILERULE +// **** I need to distinguish PSF-matched images from raw +// **** save (somewhere) the PSF-matched PSF values + +// this function measures the radial aperture fluxes for the set of readouts. this function +// may be called multiple times (presumably with different matched PSF sizes). we must have +// already added an entry to the readout->analysis identifying the FWHM of this version. + +bool psphotRadialAperturesByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule, int nMatchedPSF) { bool status; @@ -28,4 +37,5 @@ psAssert (radMax, "annular bins (RADIAL.ANNULAR.BINS.UPPER) are not defined in the recipe"); psAssert (radMax->n, "no valid annular bins (RADIAL.ANNULAR.BINS.UPPER) are define"); + float outerRadius = radMax->data.F32[radMax->n - 1]; // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) @@ -39,4 +49,17 @@ float SN_LIM = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM"); + // how many target PSFs do we want? + int nPSFsizes = 0; + { + psMetadataLookupF32 (&status, recipe, "PSPHOT.STACK.TARGET.PSF.FWHM"); + if (status) { + nPSFsizes = 1; + } else { + psVector *fwhmValues = psMetadataLookupVector(&status, recipe, "PSPHOT.STACK.TARGET.PSF.FWHM"); // Magnitude offsets + psAssert (status, "missing psphot recipe value PSPHOT.STACK.TARGET.PSF.FWHM"); + nPSFsizes = fwhmValues->n; + } + } + // source analysis is done in S/N order (brightest first) objects = psArraySort (objects, pmPhotObjSortBySN); @@ -52,4 +75,15 @@ pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); psAssert (readout, "missing readout?"); + + psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES"); + if (!fwhmValues) { + psError (PSPHOT_ERR_CONFIG, true, "convolved or measured FWHM is not defined for this readout"); + return false; + } + if (fwhmValues->n != nMatchedPSF + 1) { + psError (PSPHOT_ERR_CONFIG, true, "convolved or measured FWHM sequence is inconsistent this readout"); + return false; + } + psLogMsg ("psphot", PS_LOG_DETAIL, "PSF FWHM of %s : %f pixels\n", file->name, fwhmValues->data.F32[nMatchedPSF]); readouts->data[i] = psMemIncrRefCounter(readout); @@ -77,21 +111,47 @@ if (source->peak->SN < SN_LIM) continue; + int index = source->imageID; + if (index >= readouts->n) continue; // skip the sources generated by the chisq image + pmReadout *readout = readouts->data[index]; + + // psLogMsg("psphot", PS_LOG_INFO, "radial apertures for %d", index); + // psphotVisualShowImage(readout); + + // allocate pmSourceExtendedParameters, if not already defined + if (source->parent) { + if (!source->parent->radialAper) { + source->parent->radialAper = psArrayAlloc(nPSFsizes); + } + } else { + if (!source->radialAper) { + source->radialAper = psArrayAlloc(nPSFsizes); + } + } + // replace object in image if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); } + + // we need to change the view for the radial aperture analysis, but we want to recover exactly + // the original view; the following elements get destroyed by pmSourceRedefinePixels so save them: + psImage *oldMaskObj = psMemIncrRefCounter(source->maskObj); + psImage *oldModelFlux = psMemIncrRefCounter(source->modelFlux); + psImage *oldPSFimage = psMemIncrRefCounter(source->psfImage); + psRegion oldRegion = source->region; + Nradial ++; - int index = source->imageID; - pmReadout *readout = readouts->data[index]; + // psLogMsg("psphot", PS_LOG_INFO, "radial apertures for %d", index); + // psphotVisualShowImage(readout); // force source image to be a bit larger... - float radius = source->peak->xf - source->pixels->col0; - radius = PS_MAX (radius, source->peak->yf - source->pixels->row0); - radius = PS_MAX (radius, source->pixels->numRows - source->peak->yf + source->pixels->row0); - radius = PS_MAX (radius, source->pixels->numCols - source->peak->xf + source->pixels->col0); - pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, 1.5*radius); + // float radius = source->peak->xf - source->pixels->col0; + // radius = PS_MAX (radius, source->peak->yf - source->pixels->row0); + // radius = PS_MAX (radius, source->pixels->numRows - source->peak->yf + source->pixels->row0); + // radius = PS_MAX (radius, source->pixels->numCols - source->peak->xf + source->pixels->col0); + pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, outerRadius + 2); - if (!psphotRadialApertureSource (source, recipe, skynoise, maskVal, radMax)) { + if (!psphotRadialApertureSource (source, recipe, skynoise, maskVal, radMax, nMatchedPSF)) { psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); } else { @@ -99,6 +159,14 @@ } + pmSourceRedefinePixelsByRegion (source, readout, oldRegion); + psFree(source->maskObj); source->maskObj = oldMaskObj; + psFree(source->modelFlux); source->modelFlux = oldModelFlux; + psFree(source->psfImage); source->psfImage = oldPSFimage; + // re-subtract the object, leave local sky pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + + // psLogMsg("psphot", PS_LOG_INFO, "radial apertures for %d", index); + // psphotVisualShowImage(readout); } } Index: /trunk/psphot/src/psphotReadout.c =================================================================== --- /trunk/psphot/src/psphotReadout.c (revision 30623) +++ /trunk/psphot/src/psphotReadout.c (revision 30624) @@ -112,5 +112,5 @@ // use bright stellar objects to measure PSF if we were supplied a PSF for any input file, // this step is skipped - if (!psphotChoosePSF (config, view, filerule)) { // pass 1 + if (!psphotChoosePSF (config, view, filerule, true)) { // pass 1 psLogMsg ("psphot", 3, "failure to construct a psf model"); return psphotReadoutCleanup (config, view, filerule); Index: /trunk/psphot/src/psphotReadoutFindPSF.c =================================================================== --- /trunk/psphot/src/psphotReadoutFindPSF.c (revision 30623) +++ /trunk/psphot/src/psphotReadoutFindPSF.c (revision 30624) @@ -49,5 +49,5 @@ } - if (!psphotChoosePSF(config, view, filerule)) { + if (!psphotChoosePSF(config, view, filerule, true)) { psError(PSPHOT_ERR_PSF, false, "Failed to construct a psf model"); return psphotReadoutCleanup (config, view, filerule); Index: /trunk/psphot/src/psphotReadoutKnownSources.c =================================================================== --- /trunk/psphot/src/psphotReadoutKnownSources.c (revision 30623) +++ /trunk/psphot/src/psphotReadoutKnownSources.c (revision 30624) @@ -43,5 +43,5 @@ } - if (!psphotChoosePSF (config, view, filerule)) { + if (!psphotChoosePSF (config, view, filerule, true)) { psError(PSPHOT_ERR_PSF, false, "Failed to construct a psf model"); return psphotReadoutCleanup (config, view, filerule); Index: /trunk/psphot/src/psphotReplaceUnfit.c =================================================================== --- /trunk/psphot/src/psphotReplaceUnfit.c (revision 30623) +++ /trunk/psphot/src/psphotReplaceUnfit.c (revision 30624) @@ -75,4 +75,6 @@ pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); } + + psphotVisualShowImage(readout); psLogMsg ("psphot.replace", PS_LOG_INFO, "replaced models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.replace")); return true; @@ -101,2 +103,182 @@ return true; } + +// modify the sources to point at the corresponding pixels for the given filerule +bool psphotRedefinePixels (pmConfig *config, const pmFPAview *view, const char *filerule) +{ + bool status = true; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psphotFileruleCount(config, filerule); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotRedefinePixelsReadout (config, view, filerule, i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to replace all sources for %s entry %d", filerule, i); + return false; + } + } + return true; +} + +bool psphotRedefinePixelsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) { + + bool status; + pmSource *source; + + psTimerStart ("psphot.replace"); + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + // XXX the sources have already been copied (merge into here?) + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + psAssert (maskVal, "missing mask value?"); + + for (int i = 0; i < sources->n; i++) { + source = sources->data[i]; + + // sources have not yet been subtracted in this image (but this flag may be raised) + source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED; + + float Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS]; + float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS]; + float radius = source->modelPSF->fitRadius; + + // force a redefine to this image + pmSourceFreePixels(source); + pmSourceRedefinePixels (source, readout, Xo, Yo, radius); + } + return true; +} + +// for now, let's store the detections on the readout->analysis for each readout +bool psphotResetModels (pmConfig *config, const pmFPAview *view, const char *filerule) +{ + bool status = true; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psphotFileruleCount(config, filerule); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotResetModelsReadout (config, view, filerule, i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to replace all sources for %s entry %d", filerule, i); + return false; + } + } + return true; +} + +bool psphotResetModelsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) { + + bool status; + pmSource *source; + + psTimerStart ("psphot.replace"); + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + // XXX the sources have already been copied (merge into here?) + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + psAssert (maskVal, "missing mask value?"); + + // what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box) + int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE"); + assert (status); + + pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); + psAssert (psf, "missing psf?"); + + for (int i = 0; i < sources->n; i++) { + source = sources->data[i]; + + // *** we need to cache the 'best' model, and we have 3 cases: + // 1) model is the psf model --> generate from the new psf + // 2) model is an unconvolved extended model --> just cache the copy (not perfect) + // 3) model is a convolved extended model --> re-generate + + // use the 'best' model to cache the model (PSF or EXT : EXT may point at one of modelFits + bool isPSF = false; + pmModel *model = pmSourceGetModel(&isPSF, source); + float radius = model->fitRadius; // save for future use below + + // regenerate the PSF if the model is a PSF, or if we need the PSF for a PCM + if (isPSF || model->isPCM) { + // the guess central intensity comes from the peak: + float Io = source->peak->flux; + float Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS]; + float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS]; + + // generate a model for this object with Io = 1.0 + pmModel *modelPSF = pmModelFromPSFforXY(psf, Xo, Yo, Io); + if (modelPSF == NULL) { + psWarning ("Failed to determine PSF model for source at (%f,%f), skipping", Xo, Yo); + continue; + } + + // set the source PSF model + psFree (source->modelPSF); + source->modelPSF = modelPSF; + source->modelPSF->fitRadius = radius; + } + + if (model->isPCM) { + psAssert(false, "this section is not complete"); + + pmSourceCachePSF (source, maskVal); + + psKernel *psfKernel = pmPCMkernelFromPSF(source, psfSize); + if (!psfKernel) { + psWarning ("no psf kernel"); + } + + // generate an image of the right size + psImage *rawModelFlux = psImageCopy (NULL, source->pixels, PS_TYPE_F32); + psImageInit (rawModelFlux, 0.0); + + // insert the model image normalized to 1.0 + pmModelAdd (rawModelFlux, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM, maskVal); + + psImageConvolveFFT (source->modelFlux, rawModelFlux, NULL, 0, psfKernel); + + psFree (psfKernel); + psFree (rawModelFlux); + } + + pmSourceCacheModel (source, maskVal); // ALLOC x14 (!) + } + + psLogMsg ("psphot.replace", PS_LOG_INFO, "subtracted models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.replace")); + return true; +} + Index: /trunk/psphot/src/psphotSetThreads.c =================================================================== --- /trunk/psphot/src/psphotSetThreads.c (revision 30623) +++ /trunk/psphot/src/psphotSetThreads.c (revision 30624) @@ -35,5 +35,5 @@ psFree(task); - task = psThreadTaskAlloc("PSPHOT_EXTENDED_FIT", 12); + task = psThreadTaskAlloc("PSPHOT_EXTENDED_FIT", 13); task->function = &psphotExtendedSourceFits_Threaded; psThreadTaskAdd(task); Index: /trunk/psphot/src/psphotSourceFits.c =================================================================== --- /trunk/psphot/src/psphotSourceFits.c (revision 30623) +++ /trunk/psphot/src/psphotSourceFits.c (revision 30624) @@ -365,5 +365,5 @@ // copy most data from the primary source (modelEXT, blends stay NULL) - pmSource *newSrc = pmSourceCopyData (source); + pmSource *newSrc = pmSourceCopy (source); newSrc->modelPSF = psMemIncrRefCounter (DBL->data[1]); @@ -385,5 +385,5 @@ psLogMsg ("psphot", 1, "PAR %d : %f +/- %f\n", i, source->modelPSF->params->data.F32[i], source->modelPSF->dparams->data.F32[i]); } - psphotVisualShowResidualImage (readout); + psphotVisualShowResidualImage (readout, false); } # endif Index: /trunk/psphot/src/psphotSourceMatch.c =================================================================== --- /trunk/psphot/src/psphotSourceMatch.c (revision 30623) +++ /trunk/psphot/src/psphotSourceMatch.c (revision 30624) @@ -187,4 +187,8 @@ pmPhotObj *obj = objects->data[i]; + // we will find the input source with the max number of spans and reproduce that footprint + int nSpansMax = 0; + int iSpansMax = -1; + // mark the images for which sources have been found psVectorInit (found, 0); @@ -196,5 +200,22 @@ psAssert (index < found->n, "invalid index"); + if (src->peak && src->peak->footprint && src->peak->footprint->nspans > nSpansMax) { + nSpansMax = src->peak->footprint->nspans; + iSpansMax = j; + } + found->data.U8[index] = 1; + } + + // we make a copy of the largest footprint; this will be used for all new sources associated with this object + pmFootprint *footprint = NULL; + if (iSpansMax != -1) { // copy the footprint info + pmSource *src = obj->sources->data[iSpansMax]; + psAssert(src->peak, "source does not exist?"); + psAssert(src->peak->footprint, "footprint does not exist"); + psAssert(src->peak->footprint->nspans == nSpansMax, "wrong footprint?"); + + // we only care about the spans, do not worry about the image of this footprint + footprint = pmFootprintCopyData(src->peak->footprint, NULL); } @@ -219,6 +240,12 @@ peak->dy = NAN; - // XXX assign to a footprint? - + // assign to a footprint on this readout->image + peak->footprint = pmFootprintCopyData(footprint, readout->image); + + // the peak does not claim ownership of the footprint (it does not free it). save a copy of this + // footprint on detections->footprints so we can free it later + psArrayAdd(detections->footprints, 100, peak->footprint); + psFree (peak->footprint); + // create a new source pmSource *source = pmSourceAlloc(); @@ -227,5 +254,5 @@ // add the peak - source->peak = psMemIncrRefCounter(peak); + source->peak = peak; // allocate space for moments @@ -239,6 +266,6 @@ psArrayAdd (detections->newSources, 100, source); psFree (source); - psFree (peak); - } + } + psFree (footprint); } Index: /trunk/psphot/src/psphotSourceSize.c =================================================================== --- /trunk/psphot/src/psphotSourceSize.c (revision 30623) +++ /trunk/psphot/src/psphotSourceSize.c (revision 30624) @@ -499,11 +499,10 @@ // clear the mask bit and set the circular mask pixels - psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal)); - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", options->markVal); - + // psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal)); + // psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", options->markVal); pmSourceMagnitudes (source, psf, photMode, maskVal, markVal); // clear the mask bit - psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal)); + // psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(options->markVal)); // re-subtract the object, leave local sky Index: /trunk/psphot/src/psphotStackArguments.c =================================================================== --- /trunk/psphot/src/psphotStackArguments.c (revision 30623) +++ /trunk/psphot/src/psphotStackArguments.c (revision 30624) @@ -41,5 +41,5 @@ if ((N = psArgumentGet (argc, argv, "-break"))) { if (argc <= N+1) { - psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1); + psErrorStackPrint(stderr, "Expected to see an argument for -break"); exit(PS_EXIT_CONFIG_ERROR); } Index: /trunk/psphot/src/psphotStackChisqImage.c =================================================================== --- /trunk/psphot/src/psphotStackChisqImage.c (revision 30623) +++ /trunk/psphot/src/psphotStackChisqImage.c (revision 30624) @@ -6,5 +6,5 @@ // XXX supply filename or keep PSPHOT.INPUT fixed? -bool psphotStackChisqImage (pmConfig *config, const pmFPAview *view, const char *ruleDet, const char *ruleCnv) +bool psphotStackChisqImage (pmConfig *config, const pmFPAview *view, const char *ruleDet, const char *ruleSrc) { psTimerStart ("psphot.chisq.image"); @@ -27,6 +27,7 @@ psMetadataAddS32(config->arguments, PS_LIST_TAIL, "PSPHOT.CHISQ.NUM", PS_META_REPLACE, "", num); + + // we need to increment the counter for ruleDet and ruleSrc: num++; - psMetadataAddS32(config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.NUM", PS_META_REPLACE, "", num); // save the resulting image in the 'detection' set @@ -35,9 +36,13 @@ return false; } + psphotFileruleCountSet(config, ruleDet, num); - // save the resulting image in the 'convolved' set - if (!psMetadataAddPtr(config->files, PS_LIST_TAIL, ruleCnv, PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "", chisqFile)) { - psError(PM_ERR_CONFIG, false, "could not add chisqFPA to config files"); - return false; + // also save the resulting image in the 'source' set (analysis set) + if (strcmp(ruleDet, ruleSrc)) { + if (!psMetadataAddPtr(config->files, PS_LIST_TAIL, ruleSrc, PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "", chisqFile)) { + psError(PM_ERR_CONFIG, false, "could not add chisqFPA to config files"); + return false; + } + psphotFileruleCountSet(config, ruleSrc, num); } @@ -119,10 +124,10 @@ psAssert (status, "programming error: must define PSPHOT.CHISQ.NUM"); - int inputNum = psphotFileruleCount(config, "PSPHOT.INPUT"); + int inputNum = psphotFileruleCount(config, filerule); pmFPAfileRemoveSingle (config->files, filerule, chisqNum); inputNum --; - psMetadataAddS32(config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.NUM", PS_META_REPLACE, "", inputNum); + psphotFileruleCountSet(config, filerule, inputNum); return true; Index: /trunk/psphot/src/psphotStackImageLoop.c =================================================================== --- /trunk/psphot/src/psphotStackImageLoop.c (revision 30623) +++ /trunk/psphot/src/psphotStackImageLoop.c (revision 30624) @@ -101,5 +101,5 @@ */ -# define UPDATE_HEADER 0 +# define UPDATE_HEADER 1 bool GetAstrometryFPA (pmConfig *config, pmFPAview *view) { @@ -174,12 +174,13 @@ pmChip *outChip = pmFPAviewThisChip(view, output->fpa); ///< Chip in the output -# if (UPDATE_HEADER) pmHDU *outHDU = pmFPAviewThisHDU (view, output->fpa); if (!outHDU) { - pmFPAAddSourceFromView(output->fpa, "name", view, output->format); + pmFPAAddSourceFromView(output->fpa, view, output->format); outHDU = pmFPAviewThisHDU (view, output->fpa); psAssert (outHDU, "failed to make HDU"); } -# endif + if (!outHDU->header) { + outHDU->header = psMetadataAlloc(); + } if (bilevelAstrometry) { @@ -188,10 +189,8 @@ continue; } -# if (UPDATE_HEADER) if (!pmAstromWriteBilevelChip(outHDU->header, outChip, WCS_NONLIN_TOL)) { psWarning("Unable to generate WCS header."); continue; } -# endif } else { // we use a default FPA pixel scale of 1.0 @@ -200,10 +199,8 @@ continue; } -# if (UPDATE_HEADER) - if (UPDATE_HEADER && !pmAstromWriteWCS(outHDU->header, output->fpa, outChip, WCS_NONLIN_TOL)) { + if (!pmAstromWriteWCS(outHDU->header, output->fpa, outChip, WCS_NONLIN_TOL)) { psWarning("Unable to generate WCS header."); continue; } -# endif } } @@ -225,16 +222,17 @@ psAssert (output, "missing file?"); -# if (UPDATE_HEADER) pmHDU *PHU = pmFPAviewThisPHU(view, output->fpa); if (!PHU) { - pmFPAAddSourceFromView(output->fpa, "name", view, output->format); + pmFPAAddSourceFromView(output->fpa, view, output->format); PHU = pmFPAviewThisPHU (view, output->fpa); psAssert (PHU, "failed to make PHU"); } + if (!PHU->header) { + PHU->header = psMetadataAlloc(); + } if (!pmAstromWriteBilevelMosaic(PHU->header, output->fpa, WCS_NONLIN_TOL)) { psWarning("Unable to generate WCS header."); } -# endif } Index: /trunk/psphot/src/psphotStackMatchPSFs.c =================================================================== --- /trunk/psphot/src/psphotStackMatchPSFs.c (revision 30623) +++ /trunk/psphot/src/psphotStackMatchPSFs.c (revision 30624) @@ -65,4 +65,13 @@ bool psphotStackMatchPSFsReadout (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index) { + psImageMaskType maskValue; + psImageMaskType markValue; + + // get the PSPHOT.MASK value from the config + if (!pmConfigMaskSetBits (&maskValue, &markValue, config)) { + psError (PS_ERR_UNKNOWN, true, "Unable to define the mask bit values"); + return false; + } + pmFPAfile *fileSrc = psphotStackGetConvolveSource(config, options, index); if (!fileSrc) { @@ -84,7 +93,6 @@ // set NAN pixels to 'SAT' - // XXX replace this is pmReadoutMaskInvalid? - psImageMaskType maskVal = pmConfigMaskGet("SAT", config); - if (!pmReadoutMaskNonfinite(readoutSrc, maskVal)) { + psImageMaskType maskSat = pmConfigMaskGet("SAT", config); + if (!pmReadoutMaskInvalid(readoutSrc, maskValue, maskSat)) { psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels in readout."); return false; @@ -95,27 +103,18 @@ matchKernel(config, readoutOut, readoutSrc, options, index); saveMatchData(readoutOut, options, index); - // renormKernel(readoutCnv, options, index); - } else { - // only match the flux (NO! not for multi-filter, at least!) - // XXX do not generate readoutCnv in this case? - // float norm = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation - // psBinaryOp(readoutRaw->image, readoutRaw->image, "*", psScalarAlloc(norm, PS_TYPE_F32)); - // psBinaryOp(readoutRaw->variance, readoutRaw->variance, "*", psScalarAlloc(PS_SQR(norm), PS_TYPE_F32)); } - rescaleData(readoutOut, config, options, index); - // dumpImage(readoutOut, readoutSrc, index, "convolved"); + // save the output fwhm values in the readout->analysis. we may have / will have multiple output PSF sizes, + // so we save this in a vector. if the vector is not yet defined, create it + bool mdok = false; + psVector *fwhmValues = psMetadataLookupVector(&mdok, readoutOut->analysis, "STACK.PSF.FWHM.VALUES"); + if (!fwhmValues) { + fwhmValues = psVectorAllocEmpty(10, PS_TYPE_F32); + psMetadataAddVector(readoutOut->analysis, PS_LIST_TAIL, "STACK.PSF.FWHM.VALUES", PS_META_REPLACE, "PSF sizes", fwhmValues); + psFree(fwhmValues); // drops the extra copy + } + psVectorAppend(fwhmValues, options->targetSeeing); return true; } - - -# if (0) -// Read previously produced kernel -if (psMetadataLookupBool(NULL, config->arguments, "PPSTACK.DEBUG.STACK")) { - loadKernel(config, readoutCnv, options, index); -} else { - matchKernel(config, readoutCnv, readoutRaw, options, index); -} -# endif Index: /trunk/psphot/src/psphotStackMatchPSFsNext.c =================================================================== --- /trunk/psphot/src/psphotStackMatchPSFsNext.c (revision 30624) +++ /trunk/psphot/src/psphotStackMatchPSFsNext.c (revision 30624) @@ -0,0 +1,151 @@ +# include "psphotInternal.h" + +// smooth the input image to match the next target PSF +// this function assumes the image has already been smoothed to match the first value (array element 0), +// and that the smoothing can use a 1D Gaussian kernel of width sqrt(TARGET^2 - CURRENT^2) +// each subsequent call +bool psphotStackMatchPSFsNext(bool *smoothAgain, pmConfig *config, const pmFPAview *view, const char *filerule, int lastSize) +{ + bool status = true; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + psVector *fwhmValues = psMetadataLookupVector(&status, recipe, "PSPHOT.STACK.TARGET.PSF.FWHM"); // Magnitude offsets + if (!status) { + // must not be a vector, only one value requested + *smoothAgain = false; + return true; + } + + if (lastSize + 1 >= fwhmValues->n) { + // all done with target FWHM values + *smoothAgain = false; + return true; + } + + int num = psphotFileruleCount(config, filerule); + + // smooth the image and variance map + bool oldThreads = psImageConvolveSetThreads(true); // Old value of threading in psImageConvolve + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotStackMatchPSFsNextReadout (config, view, filerule, i, recipe, fwhmValues, lastSize)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to smooth image %s (%d) to target PSF %f", filerule, i, fwhmValues->data.F32[lastSize+1]); + psImageConvolveSetThreads(oldThreads); + return false; + } + } + + psImageConvolveSetThreads(oldThreads); + return true; +} + +bool psphotStackMatchPSFsNextReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, psVector *fwhmValues, int lastSize) { + + bool status = false; + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + psLogMsg("psphot", PS_LOG_INFO, "smooth %d to next psf", index); + psphotVisualShowImage(readout); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + psAssert (maskVal, "missing mask value?"); + + float minGauss = psMetadataLookupF32(NULL, recipe, "PEAKS_MIN_GAUSS"); // Minimum valid fraction of kernel + if (!isfinite(minGauss)) { + psWarning("PEAKS_MIN_GAUSS is not set in recipe; using default value"); + minGauss = 0.5; + } + + float targetFWHM = fwhmValues->data.F32[lastSize + 1]; + float currentFWHM = fwhmValues->data.F32[lastSize]; + + if (targetFWHM <= currentFWHM) { + psError (PSPHOT_ERR_CONFIG, true, "target FWHM cannot be smaller than current FWHM"); + return false; + } + + float smoothFWHM = sqrt(PS_SQR(targetFWHM) - PS_SQR(currentFWHM)); + float SIGMA_SMTH = smoothFWHM / (2.0*sqrt(2.0*log(2.0))); + float NSIGMA_SMTH = 3.0; + + // record the actual smoothing sigma + // psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SIGMA_SMOOTH", PS_META_REPLACE, "Smoothing sigma for detections", SIGMA_SMTH); + + // smooth the image, applying the mask as we go + psImageSmoothMask_Threaded(readout->image, readout->image, readout->mask, maskVal, SIGMA_SMTH, NSIGMA_SMTH, minGauss); + psLogMsg("psphot", PS_LOG_MINUTIA, "smooth image: %f sec\n", psTimerMark("psphot.smooth")); + + // Smooth the variance, applying the mask as we go. The variance is smoothed by the PSF^2, + // renomalized to maintain the input level of the variance. We achieve this by smoothing + // with a Gaussian with sigma = SIGMA_SMTH/sqrt(2) with unity normalization. Note that + // this process yields a smoothed image with correlated errors. The pixel-to-pixel + // variations in smooth_im will be decreased by a factor of 4*pi*SIGMA_SMTH^2, but for + // measurements based on apertures comparable to or larger than the smoothing kernel, the + // effective per-pixel variance is maintained. + psImageSmoothMask_Threaded(readout->variance, readout->variance, readout->mask, maskVal, SIGMA_SMTH * M_SQRT1_2, NSIGMA_SMTH, minGauss); + psLogMsg("psphot", PS_LOG_MINUTIA, "smooth variance: %f sec\n", psTimerMark("psphot.smooth")); + + psLogMsg("psphot", PS_LOG_INFO, "smoothed"); + psphotVisualShowImage(readout); + + // optionally save example images under trace + if (psTraceGetLevel("psphot") > 5) { + static int pass = 0; + char name[64]; + sprintf (name, "stksm.v%d.fits", pass); + psphotSaveImage(NULL, readout->image, name); + sprintf (name, "stkwt.v%d.fits", pass); + psphotSaveImage(NULL, readout->variance, name); + pass ++; + } + + // XXX need to apply this to the radial apertures somehow. + // Calculate correction factor for the covariance produced by the (potentially multiple) smoothing + // psKernel *kernel = psImageSmoothKernel(SIGMA_SMTH, NSIGMA_SMTH); // Kernel used for smoothing + // double sum2 = 0.0; // Sum of kernel squared + // for (int y = kernel->yMin; y <= kernel->yMax; y++) { + // for (int x = kernel->xMin; x <= kernel->xMax; x++) { + // sum2 += PS_SQR(kernel->kernel[y][x]); + // } + // } + // float factor = 1.0 / (sum2 * psImageCovarianceCalculateFactor(kernel, readout->covariance)); + // psFree(kernel); + + // record the effective area and significance scaling factor + float effArea = 8.0 * M_PI * PS_SQR(SIGMA_SMTH); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "EFFECTIVE_AREA", PS_META_REPLACE, "Effective Area", effArea); + // psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SIGNIFICANCE_SCALE_FACTOR", PS_META_REPLACE, "Signicance scale factor", factor); + + // save the output fwhm values in the readout->analysis. we may have / will have multiple output PSF sizes, + // so we save this in a vector. if the vector is not yet defined, create it + + psVector *fwhmValuesOut = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES"); + psAssert(fwhmValuesOut, "should already exist.."); + psVectorAppend(fwhmValuesOut, targetFWHM); + + // do not generate a PSF if we already were supplied one + pmPSF *psfOld = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); + if (psfOld) { + // save PSF on readout->analysis + char psfEntry[64]; + snprintf (psfEntry, 64, "PSPHOT.PSF.V%d", lastSize); + if (!psMetadataAddPtr (readout->analysis, PS_LIST_TAIL, psfEntry, PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot psf model", psfOld)) { + psError (PSPHOT_ERR_UNKNOWN, false, "problem saving sources on readout"); + return false; + } + psMetadataRemoveKey(readout->analysis, "PSPHOT.PSF"); + } + + return true; +} Index: /trunk/psphot/src/psphotStackMatchPSFsUtils.c =================================================================== --- /trunk/psphot/src/psphotStackMatchPSFsUtils.c (revision 30623) +++ /trunk/psphot/src/psphotStackMatchPSFsUtils.c (revision 30624) @@ -365,8 +365,10 @@ pmSubtractionSetFWHMs(options->inputSeeing->data.F32[index], options->targetSeeing); - if (scale && !pmSubtractionParamsScale(&size, &footprint, widthsCopy, scaleRef, scaleMin, scaleMax)) { - psError(psErrorCodeLast(), false, "Unable to scale kernel parameters"); - goto escape; - } + pmSubtractionParamScaleOptions(scale, scaleRef, scaleMin, scaleMax); + + // if (scale && !pmSubtractionParamsScale(&size, &footprint, widthsCopy, scaleRef, scaleMin, scaleMax)) { + // psError(psErrorCodeLast(), false, "Unable to scale kernel parameters"); + // goto escape; + // } if (!pmSubtractionMatch(NULL, readoutOut, fake, readoutSrc, footprint, stride, regionSize, spacing, threshold, stampSources, stampsName, type, size, order, widthsCopy, orders, inner, ringsOrder, binning, penalty, optimum, optWidths, optOrder, optThresh, iter, rej, normFrac, sysError, skyErr, kernelError, covarFrac, maskVal, maskBad, maskPoor, poorFrac, badFrac, PM_SUBTRACTION_MODE_2)) { @@ -383,4 +385,10 @@ goto escape; } + + // save the PSF on the new readout->analysis: + // if (!psMetadataAddPtr (readoutOut->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot psf model", options->psf)) { + // psError (PSPHOT_ERR_UNKNOWN, false, "problem saving sources on readout"); + // return false; + // } // dumpImage(readoutOut, readoutSrc, index, "conv"); Index: /trunk/psphot/src/psphotStackPSF.c =================================================================== --- /trunk/psphot/src/psphotStackPSF.c (revision 30623) +++ /trunk/psphot/src/psphotStackPSF.c (revision 30624) @@ -12,12 +12,13 @@ bool autoPSF = psMetadataLookupBool (&mdok, psphotRecipe, "PSPHOT.STACK.TARGET.PSF.AUTO"); + // Get the recipe values + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSTACK"); // ppStack recipe + psAssert(recipe, "We've thrown an error on this before."); + + char *psfModel = psMetadataLookupStr(NULL, recipe, "PSF.MODEL"); // Model for PSF + if (autoPSF) { - // Get the recipe values - psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSTACK"); // ppStack recipe - psAssert(recipe, "We've thrown an error on this before."); - int psfInstances = psMetadataLookupS32(NULL, recipe, "PSF.INSTANCES"); // Number of instances for PSF float psfRadius = psMetadataLookupF32(NULL, recipe, "PSF.RADIUS"); // Radius for PSF - const char *psfModel = psMetadataLookupStr(NULL, recipe, "PSF.MODEL"); // Model for PSF int psfOrder = psMetadataLookupS32(NULL, recipe, "PSF.ORDER"); // Spatial order for PSF @@ -49,10 +50,14 @@ float targetFWHM = psMetadataLookupF32 (&mdok, psphotRecipe, "PSPHOT.STACK.TARGET.PSF.FWHM"); - psAssert (isfinite(targetFWHM), "missing psphot recipe value PSPHOT.STACK.TARGET.PSF.FWHM"); + if (!mdok) { + psVector *fwhmValues = psMetadataLookupVector(&mdok, psphotRecipe, "PSPHOT.STACK.TARGET.PSF.FWHM"); // Magnitude offsets + psAssert (mdok, "missing psphot recipe value PSPHOT.STACK.TARGET.PSF.FWHM"); + targetFWHM = fwhmValues->data.F32[0]; + } float Sxx = sqrt(2.0)*targetFWHM / 2.35; // XXX probably should make the model type (and par 7) optional from recipe - psf = pmPSFBuildSimple("PS_MODEL_PS1_V1", Sxx, Sxx, 0.0, 1.0); + psf = pmPSFBuildSimple(psfModel, Sxx, Sxx, 0.0, 1.0); if (!psf) { psError(PSPHOT_ERR_PSF, false, "Unable to build dummy PSF."); Index: /trunk/psphot/src/psphotStackParseCamera.c =================================================================== --- /trunk/psphot/src/psphotStackParseCamera.c (revision 30623) +++ /trunk/psphot/src/psphotStackParseCamera.c (revision 30624) @@ -15,4 +15,6 @@ } + int nRaw = 0; + int nCnv = 0; int nInputs = inputs->list->n; for (int i = 0; i < nInputs; i++) { @@ -57,4 +59,5 @@ } } + nRaw ++; } @@ -88,4 +91,5 @@ } } + nCnv ++; } @@ -95,6 +99,13 @@ } + // XXX what if they do not match in length + if (nCnv && nRaw) { + if (nCnv != nRaw) { + psError (PSPHOT_ERR_CONFIG, true, "if both RAW and CNV images are supplied, the number must match"); + return false; + } + } + psString sources = psMetadataLookupStr(&status, input, "SOURCES"); // Name of mask - // pmFPAfile *srcInputFile = rawInputFile ? rawInputFile : cnvInputFile; if (sources && strlen(sources) > 0) { if (!defineFile(config, NULL, "PSPHOT.STACK.SOURCES", sources, PM_FPA_FILE_CMF)) { @@ -107,4 +118,5 @@ // XXX output of these files should be optional { + // pmFPAfile *srcInputFile = rawInputFile ? rawInputFile : cnvInputFile; pmFPAfile *outputImage = pmFPAfileDefineOutput(config, NULL, "PSPHOT.STACK.OUTPUT.IMAGE"); if (!outputImage) { @@ -150,4 +162,7 @@ } psMetadataRemoveKey(config->arguments, "FILENAMES"); + psMetadataAddS32 (config->arguments, PS_LIST_TAIL, "PSPHOT.STACK.INPUT.RAW.NUM", PS_META_REPLACE, "number of inputs", nRaw); + psMetadataAddS32 (config->arguments, PS_LIST_TAIL, "PSPHOT.STACK.INPUT.CNV.NUM", PS_META_REPLACE, "number of inputs", nCnv); + psMetadataAddS32 (config->arguments, PS_LIST_TAIL, "PSPHOT.STACK.OUTPUT.IMAGE.NUM", PS_META_REPLACE, "number of inputs", nInputs); psMetadataAddS32 (config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.NUM", PS_META_REPLACE, "number of inputs", nInputs); Index: /trunk/psphot/src/psphotStackReadout.c =================================================================== --- /trunk/psphot/src/psphotStackReadout.c (revision 30623) +++ /trunk/psphot/src/psphotStackReadout.c (revision 30624) @@ -1,6 +1,44 @@ # include "psphotInternal.h" +// we have 3 possible real filesets: # define STACK_RAW "PSPHOT.STACK.INPUT.RAW" -# define STACK_OUT "PSPHOT.STACK.OUTPUT.IMAGE" +# define STACK_CNV "PSPHOT.STACK.INPUT.CNV" +# define STACK_OUT "PSPHOT.STACK.OUTPUT.IMAGE" /* the psf-matched image */ + +// we have 3 files on which we operate: +// DET (detection image) : nominally RAW (optionally CNV?) +// SRC (source analysis image) : nominally CNV (optionally RAW) +// OUT (psf-matched images) : always OUT + +bool psphotStackVisualFilerule(pmConfig *config, const pmFPAview *view, const char *filerule) { + + bool status = false; + + int num = psphotFileruleCount(config, filerule); + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); + + psphotVisualShowResidualImage (readout, true); + psphotVisualShowObjectRegions (readout, recipe, sources); + } + return true; +} bool psphotStackReadout (pmConfig *config, const pmFPAview *view) { @@ -20,8 +58,13 @@ PS_ASSERT_PTR_NON_NULL (breakPt, false); + // XXX or do I set OUT to be a pmFPAfile pointing at the input of interest? + bool useRaw = psMetadataLookupBool (NULL, recipe, "PSPHOT.STACK.USE.RAW"); + char *STACK_SRC = useRaw ? STACK_RAW : STACK_CNV; + char *STACK_DET = STACK_RAW; // XXX optionally allow this to be CNV? + // we have 3 relevant files: RAW, CNV, OUT // set the photcode for each image - if (!psphotAddPhotcode (config, view, STACK_OUT)) { + if (!psphotAddPhotcode (config, view, STACK_SRC)) { psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode"); return false; @@ -30,9 +73,9 @@ // Generate the mask and weight images // XXX this should be done before we perform the convolutions - if (!psphotSetMaskAndVariance (config, view, STACK_RAW)) { - return psphotReadoutCleanup (config, view, STACK_OUT); + if (!psphotSetMaskAndVariance (config, view, STACK_DET)) { + return psphotReadoutCleanup (config, view, STACK_SRC); } if (!strcasecmp (breakPt, "NOTHING")) { - return psphotReadoutCleanup (config, view, STACK_OUT); + return psphotReadoutCleanup (config, view, STACK_SRC); } @@ -40,57 +83,63 @@ // XXX I think this is not defined correctly for an array of images. // XXX probably need to subtract the model (same model?) for both RAW and OUT - if (!psphotModelBackground (config, view, STACK_RAW)) { - return psphotReadoutCleanup (config, view, STACK_OUT); - } - if (!psphotSubtractBackground (config, view, STACK_RAW)) { - return psphotReadoutCleanup (config, view, STACK_OUT); + if (!psphotModelBackground (config, view, STACK_DET)) { + return psphotReadoutCleanup (config, view, STACK_SRC); + } + if (!psphotSubtractBackground (config, view, STACK_DET)) { + return psphotReadoutCleanup (config, view, STACK_SRC); } if (!strcasecmp (breakPt, "BACKMDL")) { - return psphotReadoutCleanup (config, view, STACK_OUT); - } - - // load the psf model, if suppled. FWHM_X,FWHM_Y,etc are determined and saved on - // readout->analysis XXX this function currently only works with a single PSPHOT.INPUT - if (!psphotLoadPSF (config, view, STACK_RAW)) { - psError (PSPHOT_ERR_UNKNOWN, false, "error loading psf model"); - return psphotReadoutCleanup (config, view, STACK_OUT); - } - - if (!psphotStackChisqImage(config, view, STACK_RAW, STACK_OUT)) { + return psphotReadoutCleanup (config, view, STACK_SRC); + } + + if (!psphotStackChisqImage(config, view, STACK_DET, STACK_SRC)) { psError (PSPHOT_ERR_UNKNOWN, false, "failure to generate chisq image"); - return psphotReadoutCleanup (config, view, STACK_OUT); + return psphotReadoutCleanup (config, view, STACK_SRC); } if (!strcasecmp (breakPt, "CHISQ")) { - return psphotReadoutCleanup (config, view, STACK_OUT); + return psphotReadoutCleanup (config, view, STACK_SRC); } // find the detections (by peak and/or footprint) in the image. // This finds the detections on Chisq image as well as the individuals - if (!psphotFindDetections (config, view, STACK_RAW, true)) { // pass 1 + if (!psphotFindDetections (config, view, STACK_DET, true)) { // pass 1 // this only happens if we had an error in psphotFindDetections psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis"); - return psphotReadoutCleanup (config, view, STACK_OUT); - } - - // copy the detections from RAW to OUT - if (!psphotCopySources (config, view, STACK_OUT, STACK_RAW)) { - psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis"); - return psphotReadoutCleanup (config, view, STACK_OUT); + return psphotReadoutCleanup (config, view, STACK_SRC); + } + + // copy the detections from DET to SRC + if (strcmp(STACK_SRC, STACK_DET)) { + if (!psphotCopySources (config, view, STACK_SRC, STACK_DET)) { + psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis"); + return psphotReadoutCleanup (config, view, STACK_SRC); + } } // construct sources and measure basic stats (saved on detections->newSources) - // only run this on detections from the input images, not chisq image - if (!psphotSourceStats (config, view, STACK_OUT, true)) { // pass 1 + if (!psphotSourceStats (config, view, STACK_SRC, true)) { // pass 1 psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources"); - return psphotReadoutCleanup (config, view, STACK_OUT); + return psphotReadoutCleanup (config, view, STACK_SRC); + } + + if (!strcasecmp (breakPt, "TEST1")) { + return psphotReadoutCleanup (config, view, STACK_SRC); } // generate the objects (object unify the sources from the different images) - psArray *objects = psphotMatchSources (config, view, STACK_OUT); + // XXX this could just match the detections for the chisq image, and not bother measuring the + // source stats in that case... + psArray *objects = psphotMatchSources (config, view, STACK_SRC); + + if (!strcasecmp (breakPt, "TEST2")) { + psFree(objects); + return psphotReadoutCleanup (config, view, STACK_SRC); + } // construct sources for the newly-generated sources (from other images) - if (!psphotSourceStats (config, view, STACK_OUT, false)) { // pass 1 + if (!psphotSourceStats (config, view, STACK_SRC, false)) { // pass 1 + psFree(objects); psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources"); - return psphotReadoutCleanup (config, view, STACK_OUT); + return psphotReadoutCleanup (config, view, STACK_SRC); } @@ -98,5 +147,5 @@ // if (!psphotDeblendSatstars (config, view)) { // psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis"); - // return psphotReadoutCleanup (config, view, STACK_OUT); + // return psphotReadoutCleanup (config, view, STACK_SRC); // } @@ -104,66 +153,102 @@ // if (!psphotBasicDeblend (config, view)) { // psError (PSPHOT_ERR_UNKNOWN, false, "failed on deblend analysis"); - // return psphotReadoutCleanup (config, view, STACK_OUT); + // return psphotReadoutCleanup (config, view, STACK_SRC); // } // classify sources based on moments, brightness // only run this on detections from the input images, not chisq image - if (!psphotRoughClass (config, view, STACK_OUT)) { + if (!psphotRoughClass (config, view, STACK_SRC)) { + psFree(objects); psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications"); - return psphotReadoutCleanup (config, view, STACK_OUT); + return psphotReadoutCleanup (config, view, STACK_SRC); } // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources) // only run this on detections from the input images, not chisq image - if (!psphotImageQuality (config, view, STACK_OUT)) { // pass 1 + if (!psphotImageQuality (config, view, STACK_SRC)) { // pass 1 + psFree(objects); psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality"); - return psphotReadoutCleanup (config, view, STACK_OUT); + return psphotReadoutCleanup (config, view, STACK_SRC); } if (!strcasecmp (breakPt, "MOMENTS")) { - return psphotReadoutCleanup (config, view, STACK_OUT); + psFree(objects); + return psphotReadoutCleanup (config, view, STACK_SRC); } // use bright stellar objects to measure PSF if we were supplied a PSF for any input file, // this step is skipped - if (!psphotChoosePSF (config, view, STACK_OUT)) { // pass 1 + if (!psphotChoosePSF (config, view, STACK_SRC, true)) { // pass 1 + psFree(objects); psLogMsg ("psphot", 3, "failure to construct a psf model"); - return psphotReadoutCleanup (config, view, STACK_OUT); + return psphotReadoutCleanup (config, view, STACK_SRC); } if (!strcasecmp (breakPt, "PSFMODEL")) { - return psphotReadoutCleanup (config, view, STACK_OUT); + psFree(objects); + return psphotReadoutCleanup (config, view, STACK_SRC); } // construct an initial model for each object, set the radius to fitRadius, set circular fit mask - psphotGuessModels (config, view, STACK_OUT); + psphotGuessModels (config, view, STACK_SRC); // merge the newly selected sources into the existing list // NOTE: merge OLD and NEW - psphotMergeSources (config, view, STACK_OUT); + psphotMergeSources (config, view, STACK_SRC); // linear PSF fit to source peaks, subtract the models from the image (in PSF mask) psphotFitSourcesLinearStack (config, objects, FALSE); + psphotStackVisualFilerule(config, view, STACK_SRC); // identify CRs and extended sources - psphotSourceSize (config, view, STACK_OUT, TRUE); + psphotSourceSize (config, view, STACK_SRC, TRUE); + + // XXX do we want to do a preliminary (unconvolved) model fit here, and then + // do a second detection pass? (like standard psphot) // measure aperture photometry corrections - if (!psphotApResid (config, view, STACK_OUT)) { + if (!psphotApResid (config, view, STACK_SRC)) { psFree (objects); psLogMsg ("psphot", 3, "failed on psphotApResid"); - return psphotReadoutCleanup (config, view, STACK_OUT); + return psphotReadoutCleanup (config, view, STACK_SRC); } psphotStackObjectsUnifyPosition (objects); - // measure circular, radial apertures (objects sorted by S/N) - psphotRadialAperturesByObject (config, objects, view, STACK_OUT); - // measure elliptical apertures, petrosians (objects sorted by S/N) - psphotExtendedSourceAnalysisByObject (config, objects, view, STACK_OUT); // pass 1 (detections->allSources) + psphotExtendedSourceAnalysisByObject (config, objects, view, STACK_SRC); // pass 1 (detections->allSources) // measure non-linear extended source models (exponential, deVaucouleur, Sersic) (sources sorted by S/N) - psphotExtendedSourceFits (config, view, STACK_OUT); // pass 1 (detections->allSources) + psphotExtendedSourceFits (config, view, STACK_SRC); // pass 1 (detections->allSources) // calculate source magnitudes - psphotMagnitudes(config, view, STACK_OUT); + psphotMagnitudes(config, view, STACK_SRC); + + // create source children for the OUT filerule (for radial aperture photometry) + psArray *objectsRadial = psphotSourceChildrenByObject (config, view, STACK_OUT, objects); + if (!objectsRadial) { + psFree(objects); + psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis"); + return psphotReadoutCleanup (config, view, STACK_SRC); + } + + bool smoothAgain = true; + for (int nMatchedPSF = 0; smoothAgain; nMatchedPSF++) { + + // re-measure the PSF for the smoothed image (using entries in 'allSources') + psphotChoosePSF (config, view, STACK_OUT, false); + + // this is necessary to update the models based on the new PSF + psphotResetModels (config, view, STACK_OUT); + + // this is necessary to get the right normalization for the new models + psphotFitSourcesLinear (config, view, STACK_OUT, false); + + // measure circular, radial apertures (objects sorted by S/N) + psphotRadialAperturesByObject (config, objectsRadial, view, STACK_OUT, nMatchedPSF); + + // replace the flux in the image so it is returned to its original state + psphotReplaceAllSources (config, view, STACK_OUT); + + // smooth to the next FWHM, or set 'smoothAgain' to false if no more + psphotStackMatchPSFsNext(&smoothAgain, config, view, STACK_OUT, nMatchedPSF); + } if (0 && !psphotEfficiency(config, view, STACK_OUT)) { @@ -176,16 +261,116 @@ // replace background in residual image - psphotSkyReplace (config, view, STACK_RAW); + psphotSkyReplace (config, view, STACK_DET); // drop the references to the image pixels held by each source - psphotSourceFreePixels (config, view, STACK_OUT); - - // remove chisq image from config->file:PSPHOT.INPUT (why?) - psphotStackRemoveChisqFromInputs(config, STACK_RAW); + // psphotSourceFreePixels (config, view, STACK_OUT); + psphotSourceFreePixels (config, view, STACK_SRC); + + // remove chisq image from config->file:PSPHOT.INPUT + psphotStackRemoveChisqFromInputs(config, STACK_DET); + if (strcmp(STACK_SRC, STACK_DET)) { + psphotStackRemoveChisqFromInputs(config, STACK_SRC); + } psFree (objects); + psFree (objectsRadial); // create the exported-metadata and free local data - return psphotReadoutCleanup (config, view, STACK_OUT); + return psphotReadoutCleanup (config, view, STACK_SRC); } +/* here is the process: + + * we have three(*) images: + * RAW : unconvolved image stack + * CNV : input convolved image + + * OUT : psf-matched output image (there may be more than one of + * these. we will generate the first matched image by selecting the + * target PSF and doing a full psf-maching process (as used by ppStack + * and ppSub). But, additional target output files should use a + * simple gaussian convolution kernel determind from therms of the + * current and the target). + + * the output should be / could be one of the matched images, but not + * all. should we ensure the first gets written out, and ot save the + * others (or only optionally). + + * by default, we probably only sve the cmf ffile outputs. + + * load the RAW image (unconvolved stacks) + * add photcode to the output headers / readout->analysis + * generate mask and variance image (this is probably never needed in + practice: we always load an input mask & var. + * generate & subtract a model background for ?? (RAW? CNV? OUT? all?) + * load a PSF (probably not yet working) + + * generate the CHISQ image from the RAW input images (why save on OUT?) + + * find detections on RAW + + * copy detections to OUT + + * generate source stats (moments) for OUT + + * match sources across inputs (on OUT?) + + * generate source stats for the new constructions + + * rough class (star, galaxy, cosmic, etc) + + * Image quality + + * generate PSF + + * guess models + + * merge sources (new -> old) + + * linear fit to the psf + + * find ApResid + + * assign common positions + + * radial apertures (** this should be on the PSF-matched images + + * extended analysis (elliptical profile & petrosian) + + * extended fits (sersic, etc) + + * psphot magnitudes + + + ****** + + the above is all wrong: first, we should be doing the full + morphology analysis (ExtendedAnalysis & ExtendedFits) on the CNV or + RAW image (as desired optionally), etc. + + In the discussion below, 'BST' (best) means optionally RAW or CNV + + * detection : RAW & CHISQ (of RAW) + * moments : used by psf analysis & classification (BST) + * rough class : uses moments, not pixels + * image quality : uses moments as well + * generate PSF : (BST) + * guess models (BST) + * linear fit (BST) + * find ApResid (BST) -- uses sources not pixels + * extended analysis (BST) + * extended fits (BST) + * detection efficiency (BST) + + * somehow need to copy the sources so they point at the pixels on the + * OUT image + + * foreach target PSF + * radial aperture + * convolve to next target PSF + + * somehow need to organize the output file to have the values from + * the different PSFs in separate tables (with header info to + * specify the size of that PSF) + + */ Index: /trunk/psphot/src/psphotVisual.c =================================================================== --- /trunk/psphot/src/psphotVisual.c (revision 30623) +++ /trunk/psphot/src/psphotVisual.c (revision 30624) @@ -93,11 +93,4 @@ strcpy (coords.ctype, "RA---TAN"); - psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); - psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); - if (!psImageBackground(stats, NULL, inImage, NULL, 0, rng)) { - fprintf (stderr, "failed to get background values\n"); - return false; - } - image.Nx = inImage->numCols; image.Ny = inImage->numRows; @@ -125,7 +118,101 @@ free (image.data2d); + return true; +} + +bool psphotVisualShowObjectRegions (pmReadout *readout, psMetadata *recipe, psArray *sources) { + + KiiImage image; + KapaImageData data; + Coords coords; + + bool status = false; + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (maskVal); + + psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); // Mask value for bad pixels + assert (maskVal); + + maskVal |= markVal; + + if (!pmVisualTestLevel("psphot.image.objects", 2)) return true; + + int kapa = psphotKapaChannel (1); + if (kapa == -1) return false; + + strcpy (coords.ctype, "RA---TAN"); + + psImage *inImage = readout->image; + psImage *inMask = readout->mask; + image.Nx = inImage->numCols; + image.Ny = inImage->numRows; + + psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); + if (!psImageBackground(stats, NULL, inImage, inMask, 0xffff, rng)) { + fprintf (stderr, "failed to get background values\n"); + return false; + } + + ALLOCATE (image.data2d, float *, image.Ny); + for (int iy = 0; iy < image.Ny; iy++) { + ALLOCATE (image.data2d[iy], float, image.Nx); + for (int ix = 0; ix < image.Nx; ix++) { + image.data2d[iy][ix] = 0; + } + } + + // loop over sources and set unmasked pixels to 0 + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + if (source == NULL) continue; + + psImage *mask = source->maskObj; + if (mask == NULL) continue; + + for (int iy = 0; iy < mask->numRows; iy++) { + int jy = iy + mask->row0; + if (jy < 0) continue; + if (jy >= inImage->numRows) continue; + for (int ix = 0; ix < mask->numCols; ix++) { + int jx = ix + mask->col0; + if (jx < 0) continue; + if (jx >= inImage->numCols) continue; + + if (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) continue; + image.data2d[jy][jx] = 1; + } + } + } + + for (int iy = 0; iy < image.Ny; iy++) { + for (int ix = 0; ix < image.Nx; ix++) { + image.data2d[iy][ix] = (image.data2d[iy][ix] == 0.0) ? NAN : inImage->data.F32[iy][ix]; + } + } + + strcpy (data.name, "maskObj"); + strcpy (data.file, "maskObj"); + // data.zero = 0.0; + // data.range = 1.0; + data.zero = stats->robustMedian - stats->robustStdev; + data.range = 5*stats->robustStdev; + data.logflux = 0; + + KiiSetChannel (kapa, 2); + KiiNewPicture2D (kapa, &image, &data, &coords); + + for (int iy = 0; iy < image.Ny; iy++) { + free (image.data2d[iy]); + } + free (image.data2d); + psFree (stats); psFree (rng); - + + pmVisualAskUser(NULL); return true; } @@ -2401,5 +2488,6 @@ } -bool psphotVisualShowResidualImage (pmReadout *readout) { +// option to redo variance since in some cases we may have displayed a different image in the meanwhile +bool psphotVisualShowResidualImage (pmReadout *readout, bool reshow) { if (!pmVisualTestLevel("psphot.image.resid", 2)) return true; @@ -2408,5 +2496,9 @@ if (myKapa == -1) return false; - psphotVisualScaleImage (myKapa, readout->image, readout->mask, "resid", 1); + if (reshow) { + psphotVisualShowMask (myKapa, readout->mask, "mask", 2); + psphotVisualScaleImage (myKapa, readout->variance, readout->mask, "variance", 1); + } + psphotVisualScaleImage (myKapa, readout->image, readout->mask, "resid", 0); pmVisualAskUser(NULL); Index: /trunk/psphot/test/tap_psphot_deteff.pro =================================================================== --- /trunk/psphot/test/tap_psphot_deteff.pro (revision 30624) +++ /trunk/psphot/test/tap_psphot_deteff.pro (revision 30624) @@ -0,0 +1,396 @@ +#!/usr/bin/env mana +# -*-sh-*- + +# detection efficiency: +# * does it depend on the stellar density? +# * does it scale correctly with the depth? + +# config for ppImage to generate chip, mask, weight +$ppImageConfig = -recipe PPIMAGE PPIMAGE_N +$ppImageConfig = $ppImageConfig -Db BACKGROUND T +$ppImageConfig = $ppImageConfig -Db CHIP.FITS T +$ppImageConfig = $ppImageConfig -Db CHIP.MASK.FITS T +$ppImageConfig = $ppImageConfig -Db CHIP.VARIANCE.FITS T +$ppImageConfig = $ppImageConfig -Db BASE.FITS F +$ppImageConfig = $ppImageConfig -Db VARIANCE.BUILD T +$ppImageConfig = $ppImageConfig -Db PHOTOM T + +# basic options for the these images (filter, location, obstype) +$BaseOptions = -type OBJECT -filter r -skymags 20.86 -ra 270.70 -dec -23.70 -pa 0.0 +$BaseOptions = $BaseOptions -Df PSASTRO:DVO.GETSTAR.MAX.RHO 50000.0 +$BaseOptions = $BaseOptions -nx 2500 -ny 2500 + +# options for the reference image +$RefOptions = $BaseOptions -exptime 100.0 -seeing 1.0 +$RefOptions = $RefOptions -D PSF.MODEL PS_MODEL_GAUSS +$RefOptions = $RefOptions -Df STARS.SIGMA.LIM 0.5 +$RefOptions = $RefOptions -Db PSF.CONVOLVE T +# add the density in mkref +# $RefOptions = $RefOptions -Df STARS.DENSITY 10.0 + +# options for the repeated images +$FakeOptions = $BaseOptions -exptime 30.0 + +# basic config for ppSim with randomly distributed stars and gridded galaxies +$RealConfig = -camera SIMTEST -recipe PPSIM STACKTEST.RUN +$RealConfig = $RealConfig -Db STARS.FAKE F -Db MATCH.DENSITY F -Db PSF.CONVOLVE T +# we add the catdir in mkexp so we can change the refbase +# $RealConfig = $RealConfig -D PSASTRO:PSASTRO.CATDIR catdir.$refbase + +# add galaxies for this test or not? +if (1) + $RealConfig = $RealConfig -Db GALAXY.FAKE F +else + $RealConfig = $RealConfig -Db GALAXY.FAKE T -Df GALAXY.MAG 17.0 + $RealConfig = $RealConfig -Db GALAXY.GRID T -D GALAXY.MODEL PS_MODEL_SERSIC + $RealConfig = $RealConfig -Df GALAXY.ARATIO.MIN 0.5 + $RealConfig = $RealConfig -Df GALAXY.ARATIO.MAX 0.5 + $RealConfig = $RealConfig -Df GALAXY.THETA.MAX 180 + $RealConfig = $RealConfig -Df GALAXY.INDEX.MIN 1 + $RealConfig = $RealConfig -Df GALAXY.INDEX.MAX 1 + $RealConfig = $RealConfig -Di GALAXY.GRID.DX 120 + $RealConfig = $RealConfig -Di GALAXY.GRID.DY 120 +end + +$ppSimOptions = $FakeOptions -D PSF.MODEL PS_MODEL_GAUSS + +list fwhm + 1.0 + 1.1 + 1.2 + 1.5 +end + +macro fulltest + if ($0 != 2) + echo "USAGE: fulltest (Npass)" + break + end + + mkref + + local i Npass imname stname + + $Npass = $1 + + $imname = image.fg + $stname = stack.fg + + for i 0 $Npass + sprintf N "%02d" $i + onestack $N + end +end + +macro summary.stats + if ($0 != 2) + echo "USAGE: summary.stats (Npass)" + break + end + + local i Npass base stack imname stname + + $Npass = $1 + + exec /bin/rm -f psf.image.summary.stats + exec /bin/rm -f psf.stack.inputs.summary.stats + exec /bin/rm -f psf.stack.conv.summary.stats + exec /bin/rm -f psf.stack.unconv.summary.stats + + $imname = image.fg + $stname = stack.fg + + for i 0 $Npass + sprintf base "$imname.%02d/$imname.stats" $i + sprintf stack "$stname.%02d/$stname" $i + echo $i $base $stack + exec grep PSF_INST_MAG $base >> psf.image.summary.stats + exec grep PSF_INST_MAG $stack.conv.stats >> psf.stack.inputs.summary.stats + exec grep PSF_INST_MAG $stack.pht.stats >> psf.stack.conv.summary.stats + exec grep PSF_INST_MAG $stack.unconv.pht.stats >> psf.stack.unconv.summary.stats + end +end + +macro onestack + if ($0 != 2) + echo "USAGE: onestack (N)" + break + end + + mkstack image.$1/image stack.$1/stack + ckstack image.$1/image stack.$1/stack + rmstack image.$1/image stack.$1/stack +end + +# create a reference database of fake stars to be used by ppSim below +macro mkref + if ($0 != 3) + echo "mkref (refbase) (density)" + break + end + + local refbase + $refbase = $1 + + exec rm -rf $refbase.catdir + exec rm -f $refbase.fits + + $RefOptions = $RefOptions -Df STARS.DENSITY $2 + + # create an image with fake sources and insert the resulting cmf file into a dvodb + $RefConfig = -camera SIMTEST -recipe PPSIM STACKTEST.MAKE -D PSASTRO:PSASTRO.CATDIR $refbase.catdir + + exec ppSim $RefOptions $RefConfig $refbase + + file synth.photcodes found + if (not($found)) + echo "making photcodes file" + mkphotcodes synth.photcodes + end + + exec addstar -D CAMERA simtest -D CATDIR $refbase.catdir -accept-astrom -photcode SYNTH.r -D PHOTCODE_FILE synth.photcodes $refbase.cmf + exec relphot -averages -D CATDIR $refbase.catdir -update -region 260 280 -33 -13 +end + +# create a realistic distribution of fake stars, GAUSS PSF +macro mkexp + if ($0 != 4) + echo "USAGE: mkexp basename fwhm refbase" + break + end + + local fwhm basename refbase + $basename = $1 + $fwhm = $2 + $refbase = $3 + + $RealConfig = $RealConfig -D PSASTRO:PSASTRO.CATDIR $refbase.catdir + + # create the raw image + echo ppSim $ppSimOptions $RealConfig $basename -seeing $fwhm + exec ppSim $ppSimOptions $RealConfig $basename -seeing $fwhm + exec /bin/mv -f $basename.cmf $basename.in.cmf + + # create the chip output + echo ppImage $ppImageConfig -file $basename.fits $basename + exec ppImage $ppImageConfig -file $basename.fits $basename +end + +macro completeness + if ($0 != 5) + echo "USAGE: completeness (raw) (out) (output) (radius)" + break + end + + load.cmf $1 Chip.psf raw + load.cmf $2 Chip.psf out + set X_raw = int(X_PSF_raw) + 0.5 + set Y_raw = int(Y_PSF_raw) + 0.5 + match2d X_PSF_raw Y_PSF_raw X_PSF_out Y_PSF_out $4 -index1 index1 -index2 index2 -closest + + histogram PSF_INST_MAG_raw nMag -16.0 -3.0 0.25 -range dMag + set fMag = zero(dMag) + for i 0 {dMag[]-1} + set inrange = (PSF_INST_MAG_raw > dMag[$i]) && (PSF_INST_MAG_raw <= dMag[$i+1]) + subset all = index1 if (inrange) + subset got = index1 if (inrange) && (index1 >= 0) + if (all[] == 0) + fMag[$i] = 0 + else + fMag[$i] = got[] / all[] + end + end + + device -n complete + resize 1000 600 + + clear -s + + section default 0 0 1 1 + lim dMag fMag; clear; box -ypad 5 +ypad 5 -ticks 1110; plot -x 1 dMag fMag + label -x mag_inst -y det_frac + + set found = (index1 >= 0) + plot PSF_INST_MAG_raw found -c red + + set ldmag = log(PSF_INST_MAG_SIG_raw) + section overlay 0 0 1 1; lim dMag -5 1.2; box -ypad 5 +ypad 5 -ticks 1011 -labels 1001; plot PSF_INST_MAG_raw ldmag + label +y log(S/N) + + line -c red -20 {log(1/25.0)} to 20 {log(1/25.0)} + line -c blue -20 {log(1/5.0)} to 20 {log(1/5.0)} + + png -name $3 +end + +macro completestack + if ($0 != 6) + echo "USAGE: completeness (raw) (out) (output) (radius) (zeropt)" + break + end + + load.cmf $1 Chip.psf raw + load.cmf $2 SkyChip.psf out + set X_raw = int(X_PSF_raw) + 0.5 + set Y_raw = int(Y_PSF_raw) + 0.5 + set Mraw = PSF_INST_MAG_raw + $5 + match2d X_PSF_raw Y_PSF_raw X_PSF_out Y_PSF_out $4 -index1 index1 -index2 index2 -closest + + histogram Mraw nMag -16.0 -3.0 0.25 -range dMag + set fMag = zero(dMag) + for i 0 {dMag[]-1} + set inrange = (Mraw > dMag[$i]) && (Mraw <= dMag[$i+1]) + subset all = index1 if (inrange) + subset got = index1 if (inrange) && (index1 >= 0) + if (all[] == 0) + fMag[$i] = 0 + else + fMag[$i] = got[] / all[] + end + end + + device -n complete + resize 1000 600 + + clear -s + + section default 0 0 1 1 + lim dMag fMag; clear; box -ypad 5 +ypad 5 -ticks 1110; plot -x 1 dMag fMag + label -x mag_inst -y det_frac + + set found = (index1 >= 0) + plot Mraw found -c red + + set ldmag = log(PSF_INST_MAG_SIG_raw) + section overlay 0 0 1 1; lim dMag -5 1.2; box -ypad 5 +ypad 5 -ticks 1011 -labels 1001; plot Mraw ldmag + label +y log(S/N) + + line -c red -20 {log(1/25.0)} to 20 {log(1/25.0)} + line -c blue -20 {log(1/5.0)} to 20 {log(1/5.0)} + + png -name $3 +end + +macro falsepositive + if ($0 != 5) + echo "USAGE: falsepositive (raw) (out) (output) (radius)" + break + end + + load.cmf $1 Chip.psf raw + load.cmf $2 Chip.psf out + set X_raw = int(X_PSF_raw) + 0.5 + set Y_raw = int(Y_PSF_raw) + 0.5 + match2d X_PSF_raw Y_PSF_raw X_PSF_out Y_PSF_out $4 -index1 index1 -index2 index2 -closest + + histogram PSF_INST_MAG_out nMag -16.0 -3.0 0.25 -range dMag + set fMag = zero(dMag) + set fAll = zero(dMag) + for i 0 {dMag[]-1} + set inrange = (PSF_INST_MAG_out > dMag[$i]) && (PSF_INST_MAG_out <= dMag[$i+1]) + subset all = index2 if (inrange) + subset got = index2 if (inrange) && (index2 == -1) + if (all[] == 0) + fMag[$i] = 0 + else + # fraction of detections that are false + fMag[$i] = got[] / all[] + end + # fraction of detections from this bin that are fake + fAll[$i] = got[] / PSF_INST_MAG_out[] + end + + device -n complete + resize 1000 600 + + clear -s + + section default 0 0 1 1 + # lim dMag fMag; clear; box -ypad 5 +ypad 5 -ticks 1110; plot -x 1 dMag fMag + lim dMag fAll; box -ypad 5 +ypad 5 -ticks 1110; plot -x 1 dMag fAll + label -x mag_inst -y false_frac + + section overlay1 0 0 1 1 + lim dMag fMag; plot -x 1 dMag fMag -c blue + + set found = (index2 >= 0) + plot PSF_INST_MAG_out found -c red + + set ldmag = log(PSF_INST_MAG_SIG_out) + section overlay 0 0 1 1; lim dMag -5 1.2; box -ypad 5 +ypad 5 -ticks 1011 -labels 1001; plot PSF_INST_MAG_out ldmag + label +y log(S/N) + + line -c red -20 {log(1/25.0)} to 20 {log(1/25.0)} + line -c blue -20 {log(1/5.0)} to 20 {log(1/5.0)} + + png -name $3 +end + +macro load.cmf + if ($0 != 4) + echo "load.cmf (filename) (ext) (label)" + break + end + + data $1 + + # create the list of fields to load + delete -q myFields + for i 0 $fields:n + if ($?myFields) + $myFields = $myFields $fields:$i + else + $myFields = $fields:$i + end + end + + read -fits $2 $myFields + + # rename the loaded vectors appending the supplied lable + for i 0 $fields:n + set $fields:$i\_$3 = $fields:$i + delete $fields:$i + end +end + +# this list defines the fields to be loaded from file +# for completeness, we only need a few of the fields +list fields + X_PSF + Y_PSF + X_PSF_SIG + Y_PSF_SIG + PSF_INST_MAG + PSF_INST_MAG_SIG +# PSF_INST_FLUX +# PSF_INST_FLUX_SIG +end + +# if we run this test as a stand-alone program somewhere, we may need to create a local copy of the photcode file: +macro mkphotcodes + if ($0 != 2) + echo "USAGE: mkphotcodes (filename)" + break + end + + exec /bin/rm -f $1 + output $1 + echo "# airmass color astrometry mag photom astrom mask photom mask" + echo "# code name type zero slope offset c1 c2 slope zero equiv sys scale scale sys poor bad poor bad" + echo " 1 g_SYNTH sec 0.000 0.000 0.000 1 3 0.0000 0 21 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 2 r_SYNTH sec 0.000 0.000 0.000 2 3 0.0000 0 22 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3 i_SYNTH sec 0.000 0.000 0.000 2 3 0.0000 0 23 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 4 z_SYNTH sec 0.000 0.000 0.000 3 4 0.0000 0 24 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 5 y_SYNTH sec 0.000 0.000 0.000 4 5 0.0000 0 25 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3001 SYNTH.g ref 0.000 0.000 0.000 - - 0.0000 0 1 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3002 SYNTH.r ref 0.000 0.000 0.000 - - 0.0000 0 2 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3003 SYNTH.i ref 0.000 0.000 0.000 - - 0.0000 0 3 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3004 SYNTH.z ref 0.000 0.000 0.000 - - 0.0000 0 4 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3005 SYNTH.y ref 0.000 0.000 0.000 - - 0.0000 0 5 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + output stdout +end + +if ($SCRIPT) + fulltest 4 + exit 0 +end Index: /trunk/psphot/test/tap_psphot_diff.pro =================================================================== --- /trunk/psphot/test/tap_psphot_diff.pro (revision 30624) +++ /trunk/psphot/test/tap_psphot_diff.pro (revision 30624) @@ -0,0 +1,576 @@ +#!/usr/bin/env mana +# -*-sh-*- + +# config for ppImage to generate chip, mask, weight +$ppImageConfig = -recipe PPIMAGE PPIMAGE_N +$ppImageConfig = $ppImageConfig -Db BACKGROUND T +$ppImageConfig = $ppImageConfig -Db CHIP.FITS T +$ppImageConfig = $ppImageConfig -Db CHIP.MASK.FITS T +$ppImageConfig = $ppImageConfig -Db CHIP.VARIANCE.FITS T +$ppImageConfig = $ppImageConfig -Db BASE.FITS F +$ppImageConfig = $ppImageConfig -Db VARIANCE.BUILD T +$ppImageConfig = $ppImageConfig -Db PHOTOM T + +# basic options for the these images (filter, location, obstype) +$BaseOptions = -type OBJECT -filter r -skymags 20.86 -ra 270.70 -dec -23.70 -pa 0.0 +$BaseOptions = $BaseOptions -Df PSASTRO:DVO.GETSTAR.MAX.RHO 50000.0 + +# options for the reference image +$RefOptions = $BaseOptions -exptime 100.0 -seeing 1.0 -D PSF.MODEL PS_MODEL_GAUSS -Df STARS.DENSITY 10.0 -Df STARS.SIGMA.LIM 0.5 + +# options for the repeated images +$FakeOptions = $BaseOptions -exptime 30.0 + +# basic config for ppSim with randomly distributed stars and gridded galaxies +$RealConfig = -camera SIMTEST -recipe PPSIM STACKTEST.RUN -D PSASTRO:PSASTRO.CATDIR catdir.ref +$RealConfig = $RealConfig -Db STARS.FAKE F -Db MATCH.DENSITY F -Db PSF.CONVOLVE F + +if (1) + $RealConfig = $RealConfig -Db GALAXY.FAKE F +else + $RealConfig = $RealConfig -Db GALAXY.FAKE T -Df GALAXY.MAG 17.0 + $RealConfig = $RealConfig -Db GALAXY.GRID T -D GALAXY.MODEL PS_MODEL_SERSIC + $RealConfig = $RealConfig -Df GALAXY.ARATIO.MIN 0.5 + $RealConfig = $RealConfig -Df GALAXY.ARATIO.MAX 0.5 + $RealConfig = $RealConfig -Df GALAXY.THETA.MAX 180 + $RealConfig = $RealConfig -Df GALAXY.INDEX.MIN 1 + $RealConfig = $RealConfig -Df GALAXY.INDEX.MAX 1 + $RealConfig = $RealConfig -Di GALAXY.GRID.DX 120 + $RealConfig = $RealConfig -Di GALAXY.GRID.DY 120 +end + +# sample alternate options: +# $ppSimOptions = $FakeOptions -D PSF.MODEL PS_MODEL_PS1_V1 +# $ppSimOptions = $FakeOptions -Df PSF.ARATIO 1.2 +# $ppSimOptions = $FakeOptions -Df PSF.THETA +30.0 +$ppSimOptions = $FakeOptions -D PSF.MODEL PS_MODEL_GAUSS + +list fwhm + 1.00 + 1.25 +end + +list fwhmMax + 1.25 + 1.00 +end + +list aratio + 1.0 + 1.0 +end + +list thetas + 0.0 + 0.0 +end + +# create a reference database of fake stars to be used by ppSim below +macro mkref + exec rm -rf catdir.ref + exec rm -f refimage.fits + + # create an image with fake sources and insert the resulting cmf file into a dvodb + $RefConfig = -camera SIMTEST -recipe PPSIM STACKTEST.MAKE -D PSASTRO:PSASTRO.CATDIR catdir.ref + + exec ppSim $RefOptions $RefConfig refimage + + file synth.photcodes found + if (not($found)) + echo "making photcodes file" + mkphotcodes synth.photcodes + end + + exec addstar -D CAMERA simtest -D CATDIR catdir.ref -accept-astrom -photcode SYNTH.r -D PHOTCODE_FILE synth.photcodes refimage.cmf + exec relphot -averages -D CATDIR catdir.ref -update -region 260 280 -33 -13 +end + +# create a realistic distribution of fake stars, GAUSS PSF +macro mkdiff + if ($0 != 5) + echo "USAGE: mkdiff (rawbase) (diffbase) (init) (order)" + break + end + + # XXX check for local reference db and run mkref if not found + file catdir.ref found + if (not($found)) + echo "make reference catalog" + mkref + end + + local i base diff + + $base = $1 + $diff = $2 + + dirname $base -var dir + mkdir $dir + + dirname $diff -var dir + mkdir $dir + + $command = ppSub + + for i 0 $fwhm:n + if ($3 == 1) + mkexp $base.$i $fwhm:$i $aratio:$i $thetas:$i $fwhmMax:$i + end + + if ($i == 0) + $command = $command -inimage $base.$i.wrp.fits + $command = $command -inmask $base.$i.wrp.mask.fits + $command = $command -invariance $base.$i.wrp.wt.fits + $command = $command -insources $base.$i.wrp.cmf + else + $command = $command -refimage $base.$i.wrp.fits + $command = $command -refmask $base.$i.wrp.mask.fits + $command = $command -refvariance $base.$i.wrp.wt.fits + $command = $command -refsources $base.$i.wrp.cmf + end + end + +# $command = $command -recipe PSPHOT DIFF +# $command = $command -F PSPHOT.PSF.SAVE PSPHOT.PSF.SKY.SAVE +# $command = $command -F PSPHOT.OUTPUT PSPHOT.OUT.CMF.MEF +# $command = $command -F PSPHOT.BACKMDL PSPHOT.BACKMDL.MEF -D PSPHOT:OUTPUT.FORMAT PS1_DV2 + $command = $command -Db DUAL F -Db ADD.NOT.SUBTRACT T + $command = $command -Db INVERSE T -convolve 1 -Di SPATIAL.ORDER $4 + $command = $command -save-inconv -save-refconv + $command = $command -threads 4 + $command = $command -image_id 1 -source_id 1 + + echo $command $diff + # exec $command $diff +end + +# create a realistic distribution of fake stars, GAUSS PSF +macro mkexp + if ($0 != 6) + echo "USAGE: mkexp basename fwhm aratio theta fwhmMax" + break + end + + local fwhm basename + $basename = $1 + $fwhm = $2 + $aratio = $3 + $theta = $4 + $fwhmMax = $5 + + echo ppSim $ppSimOptions $RealConfig $basename -seeing $fwhm -seeing-max $fwhmMax -seeing-ramp -Df PSF.ARATIO $aratio -Df PSF.THETA $theta + exec ppSim $ppSimOptions $RealConfig $basename -seeing $fwhm -seeing-max $fwhmMax -seeing-ramp -Df PSF.ARATIO $aratio -Df PSF.THETA $theta + exec /bin/mv -f $basename.cmf $basename.in.cmf + + # create the chip output + echo ppImage $ppImageConfig -file $basename.fits $basename + exec ppImage $ppImageConfig -file $basename.fits $basename + + # XXX pswarp is using all 7k+ sources to measure the PSF : can we reduce this? + echo pswarp -threads 4 -astrom $basename.cmf -file $basename.ch.fits -mask $basename.ch.mk.fits -variance $basename.ch.wt.fits $basename.wrp refimage.fits + exec pswarp -threads 4 -astrom $basename.cmf -file $basename.ch.fits -mask $basename.ch.mk.fits -variance $basename.ch.wt.fits $basename.wrp refimage.fits -D PSPHOT:OUTPUT.FORMAT PS1_V2 +end + +# compare chip to warp +macro ckwarp + if ($0 != 5) + echo "USAGE: ckwarp (raw) (out) (output) (Next)" + break + end + + load.cmf $1 Chip.psf raw + load.cmf $2 SkyChip.psf out + set X_raw = int(X_PSF_raw) + 0.5 + set Y_raw = int(Y_PSF_raw) + 0.5 + + # M_raw is the predicted mag based on the raw input mags (need to add the zero point) + + set M_raw = PSF_INST_MAG_raw - 2.5*log($4) + set dM_raw = PSF_INST_MAG_SIG_raw / sqrt($4) + set K_out = -2.5*log(KRON_FLUX_out) + match2d X_PSF_out Y_PSF_out X_PSF_raw Y_PSF_raw 1.0 -index1 index1 -index2 index2 + + local i nx ny NX NY N + + device -n compare + resize 1000 1000 + + # plot trends as a function of mag + $NX = 2 + $NY = 5 + $nx = 0 + $ny = 0 + $N = 0 + clear -s + + for i 0 $pairs:n + section a$nx\$ny {$nx/$NX} {$ny/$NY} {1/$NX} {1/$NY} + show.pair $i + $ny ++ + if ($ny == $NY) + $ny = 0 + $nx ++ + end + if ($nx == $NX) + png -name $3.$N.png + clear -s + $nx = 0 + $ny = 0 + $N ++ + end + end +end + +# compare two cmf files with extname Chip.psf +# things to compare: +# * completeness (which sources in (1) are not detected in (2) +# * positions (X_PSF, Y_PSF) +# * instrumental psf mags +# * position errors (no input errors; use a model?) +# * measured FWHM? +# * kron mags (fluxes) +# * etc, etc +macro ckchip + if ($0 != 5) + echo "USAGE: ckchip (raw) (out) (output) (zpt_off)" + break + end + + load.cmf $1 Chip.psf raw + load.cmf $2 Chip.psf out + + # images generated with convolution will not have the right output positions + set X_raw = int(X_PSF_raw) + 0.5 + set Y_raw = int(Y_PSF_raw) + 0.5 + set M_raw = PSF_INST_MAG_raw + $4 + set K_out = -2.5*log(KRON_FLUX_out) + match2d X_PSF_out Y_PSF_out X_PSF_raw Y_PSF_raw 1.5 -index1 index1 -index2 index2 + + local i NX NY nx ny N + + device -n compare + resize 1000 1000 + + # plot trends as a function of mag + $NX = 2 + $NY = 5 + $nx = 0 + $ny = 0 + $N = 0 + clear -s + for i 0 $pairs:n + section a$nx\$ny {$nx/$NX} {$ny/$NY} {1/$NX} {1/$NY} + show.pair $i + $ny ++ + if ($ny == $NY) + $ny = 0 + $nx ++ + end + if ($nx == $NX) + png -name $3.$N.png + clear -s + $nx = 0 + $ny = 0 + $N ++ + end + end + + # plot (input - output) vs mag +end + +macro stats.pair + if ($0 != 3) + echo "USAGE: stats.pair (N) (output)" + break + end + + list word -split $spairs:$1 + if ($word:n != 8) + echo "invalid pair $1" + break + end + + $Nr = $word:3 + + reindex v1 = $word:0 using index1 + reindex v2 = $word:1 using index2 + reindex rv = $word:2 using index$Nr + + set delta = v1 - v2 + subset d1 = delta if ($word:4 < rv) && (rv < $word:5) && (abs(delta) < $word:7) + subset d2 = delta if ($word:5 < rv) && (rv < $word:6) && (abs(delta) < $word:7) + + vstats -q d1 -sigma-clip 3.0 + $M1 = $MEAN + $S1 = $SIGMA + vstats -q d2 -sigma-clip 3.0 + $M2 = $MEAN + $S2 = $SIGMA + + output $2 + fprintf "%-18s %7.4f %7.4f %7.4f %7.4f" $word:0 $M1 $S1 $M2 $S2 + output stdout +end + +macro show.pair + if ($0 != 2) + echo "USAGE: show.pair (N)" + break + end + + list word -split $pairs:$1 + if ($word:n != 7) + echo "invalid pair $1" + break + end + + $Nr = $word:3 + + reindex v1 = $word:0 using index1 + if ("$word:6" == "V") + reindex v2 = $word:1 using index2 + end + if ("$word:6" == "S") + set v2 = $word:1 + zero(index1) + end + reindex rv = $word:2 using index$Nr + + set delta = v1 - v2 + if (("$word:4" == "def") || ("$word:5" == "def")) + lim rv delta; box; plot rv delta + else + lim rv $word:4 $word:5; box; plot rv delta + end + label -y '$word:0' -x '$word:2' +end + +# This list is used to compare a pair of vectors (sans error) or a +# vector and an expected (constant) value. The last field defines a +# vector or constant for the comparison. It is assumed that the +# vector sets have been loaded and matched with match2d to generate +# index vectors 'index1' and index2'. The macro 'show.pair' generates +# a plot of the range vector vs (v1 - v2). The indices for v1, v2 are +# index1 and 2 respectively. The index for the range vector is defined +# by the integer following that vector. the y-limits of the plot are +# given by the last two numbers +list pairs + X_PSF_out X_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + Y_PSF_out Y_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + X_PSF_out X_PSF_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + Y_PSF_out Y_PSF_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + X_PSF_SIG_out X_PSF_SIG_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + Y_PSF_SIG_out Y_PSF_SIG_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + #PSF_INST_MAG_out PSF_INST_MAG_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + PSF_INST_MAG_out M_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + PSF_INST_MAG_SIG_out dM_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + #PSF_INST_FLUX_out PSF_INST_FLUX_raw PSF_INST_MAG_raw 2 def def V + #PSF_INST_FLUX_SIG_out PSF_INST_FLUX_SIG_raw PSF_INST_MAG_raw 2 def def V + AP_MAG_out M_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + AP_MAG_RAW_out M_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + AP_MAG_RADIUS_out 0.0 PSF_INST_MAG_raw 2 -0.01 20.1 S + SKY_out 0.0 PSF_INST_MAG_raw 2 def def S + SKY_SIGMA_out 0.0 PSF_INST_MAG_raw 2 def def S + PSF_CHISQ_out 1.0 PSF_INST_MAG_raw 2 def def S + CR_NSIGMA_out 0.0 PSF_INST_MAG_raw 2 def def S + EXT_NSIGMA_out 0.0 PSF_INST_MAG_raw 2 -5.01 5.01 S + PSF_MAJOR_out 0.0 PSF_INST_MAG_raw 2 -0.01 5.01 S + PSF_MINOR_out 0.0 PSF_INST_MAG_raw 2 -0.01 5.01 S + PSF_THETA_out 0.0 PSF_INST_MAG_raw 2 -1.61 1.61 S + PSF_QF_out 0.0 PSF_INST_MAG_raw 2 -0.10 1.10 S + PSF_QF_PERFECT_out 0.0 PSF_INST_MAG_raw 2 -0.10 1.10 S + PSF_NDOF_out 0.0 PSF_INST_MAG_raw 2 def def S + PSF_NPIX_out 0.0 PSF_INST_MAG_raw 2 def def S + MOMENTS_XX_out 0.0 PSF_INST_MAG_raw 2 -0.01 3.01 S + MOMENTS_XY_out 0.0 PSF_INST_MAG_raw 2 -3.01 3.01 S + MOMENTS_YY_out 0.0 PSF_INST_MAG_raw 2 -0.01 3.01 S + MOMENTS_M3C_out 0.0 PSF_INST_MAG_raw 2 -3.01 3.01 S + MOMENTS_M3S_out 0.0 PSF_INST_MAG_raw 2 -3.01 3.01 S + MOMENTS_M4C_out 0.0 PSF_INST_MAG_raw 2 -2.01 2.01 S + MOMENTS_M4S_out 0.0 PSF_INST_MAG_raw 2 -2.01 2.01 S + MOMENTS_R1_out 0.0 PSF_INST_MAG_raw 2 -5.01 5.01 S + MOMENTS_RH_out 0.0 PSF_INST_MAG_raw 2 -5.01 5.01 S + K_out M_raw PSF_INST_MAG_raw 2 def def V + KRON_FLUX_ERR_out 0.0 PSF_INST_MAG_raw 2 def def S + KRON_FLUX_INNER_out 0.0 PSF_INST_MAG_raw 2 def def S + KRON_FLUX_OUTER_out 0.0 PSF_INST_MAG_raw 2 def def S +# CAL_PSF_MAG CAL_PSF_MAG none Mraw 2 -1.01 1.01 V +# CAL_PSF_MAG_SIG CAL_PSF_MAG_SIG none Mraw 2 -1.01 1.01 V +# RA_PSF RA_PSF none Mraw 2 -1.01 1.01 V +# DEC_PSF DEC_PSF none Mraw 2 -1.01 1.01 V +# PEAK_FLUX_AS_MAG PEAK_FLUX_AS_MAG none Mraw 2 -1.01 1.01 V +# FLAGS FLAGS 0.0 Mraw 2 -1.01 1.01 S +# FLAGS2 FLAGS2 0.0 Mraw 2 -1.01 1.01 S +end + +macro load.cmf + if ($0 != 4) + echo "load.cmf (filename) (ext) (label)" + break + end + + data $1 + + # create the list of fields to load + delete -q myFields + for i 0 $fields:n + if ($?myFields) + $myFields = $myFields $fields:$i + else + $myFields = $fields:$i + end + end + + read -fits $2 $myFields + + # rename the loaded vectors appending the supplied lable + for i 0 $fields:n + set $fields:$i\_$3 = $fields:$i + delete $fields:$i + end +end + +# this list defines the fields to be loaded from file +list fields + X_PSF + Y_PSF + X_PSF_SIG + Y_PSF_SIG + PSF_INST_MAG + PSF_INST_MAG_SIG + PSF_INST_FLUX + PSF_INST_FLUX_SIG + AP_MAG + AP_MAG_RAW + AP_MAG_RADIUS + SKY + SKY_SIGMA + PSF_CHISQ + CR_NSIGMA + EXT_NSIGMA + PSF_MAJOR + PSF_MINOR + PSF_THETA + PSF_QF + PSF_QF_PERFECT + PSF_NDOF + PSF_NPIX + MOMENTS_XX + MOMENTS_XY + MOMENTS_YY + MOMENTS_M3C + MOMENTS_M3S + MOMENTS_M4C + MOMENTS_M4S + MOMENTS_R1 + MOMENTS_RH + KRON_FLUX + KRON_FLUX_ERR + KRON_FLUX_INNER + KRON_FLUX_OUTER +# CAL_PSF_MAG +# CAL_PSF_MAG_SIG +# RA_PSF +# DEC_PSF +# PEAK_FLUX_AS_MAG +# FLAGS +# FLAGS2 +end + +# use these cmf entries to measure average stats of the given pairs +list spairs + X_PSF_out X_raw PSF_INST_MAG_raw 2 -15 -10 -8.5 0.1 + Y_PSF_out Y_raw PSF_INST_MAG_raw 2 -15 -10 -8.5 0.1 + PSF_INST_MAG_out M_raw PSF_INST_MAG_raw 2 -15 -10 -8.5 0.1 + AP_MAG_out M_raw PSF_INST_MAG_raw 2 -15 -10 -8.5 0.1 +end + +# # XXX this is a hack: the cmf file created by ppSim is not +# # compatible with the image for warp (because one is a chip-mosaic +# # and the other is not) +# exec rm -f fix.hdr +# output fix.hdr +# echo "PSMOSAIC= 'CHIP ' / Mosaicked level" +# output stdout +# fix the header to be compatible with the chip file (to avoid running psphot) +# exec fits_insert $basename.cmf fix.hdr + +# if we run this test as a stand-alone program somewhere, we may need to create a local copy of the photcode file: +macro mkphotcodes + if ($0 != 2) + echo "USAGE: mkphotcodes (filename)" + break + end + + exec /bin/rm -f $1 + output $1 + echo "# airmass color astrometry mag photom astrom mask photom mask" + echo "# code name type zero slope offset c1 c2 slope zero equiv sys scale scale sys poor bad poor bad" + echo " 1 g_SYNTH sec 0.000 0.000 0.000 1 3 0.0000 0 21 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 2 r_SYNTH sec 0.000 0.000 0.000 2 3 0.0000 0 22 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3 i_SYNTH sec 0.000 0.000 0.000 2 3 0.0000 0 23 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 4 z_SYNTH sec 0.000 0.000 0.000 3 4 0.0000 0 24 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 5 y_SYNTH sec 0.000 0.000 0.000 4 5 0.0000 0 25 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3001 SYNTH.g ref 0.000 0.000 0.000 - - 0.0000 0 1 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3002 SYNTH.r ref 0.000 0.000 0.000 - - 0.0000 0 2 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3003 SYNTH.i ref 0.000 0.000 0.000 - - 0.0000 0 3 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3004 SYNTH.z ref 0.000 0.000 0.000 - - 0.0000 0 4 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3005 SYNTH.y ref 0.000 0.000 0.000 - - 0.0000 0 5 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + output stdout +end + +if ($SCRIPT) + fulltest 4 + exit 0 +end + +macro ckvar + if ($0 != 2) + echo "ckvar (buffer)" + break + end + + cursor + stats $1 {$X1-10} {$Y1-10} 21 21 + $K2sum0 = $TOTAL + stats $1 {$X2-10} {$Y2-10} 21 21 + $K2sum1 = $TOTAL + echo {$K2sum0 * 193 + $K2sum1 * 193} +end + +macro ckresid + + rd res resid.fits + rd con1 conv1.fits + rd con2 conv2.fits + dev -n im; tvch 1; tv res -50 100; tvch 2; tv con1 -20 100; tvch 3; tv con2 -20 100 +end + +macro ckchisq + data chisq2.dat + read flux2 3 Mxx 9 Myy 11 chisq 13 npix 15 + set chisqv = chisq / npix + set flux = sqrt(flux2) + + set lflux = log(flux) + set lchi = log(chisqv) + set lMxx = log(Mxx) + + set dev = sqrt(chisq) + set ndev = sqrt(chisq)/npix + set lndev = log(ndev) + + lim -n 0 lchi lMxx; clear; box; plot lchi lMxx + lim -n 1 lflux lMxx; clear; box; plot lflux lMxx + lim -n 2 lflux lchi; clear; box; plot lflux lchi + lim -n 3 lflux lndev; clear; box; plot lflux lndev + + # set flux2 = flux + # set flux = sqrt(flux2) + # set dev = sqrt(chisq) + # set ndev = sqrt(chisq)/npix + # set lndev = log(ndev) + # set lflux = log(flux) + # lim -n 0 lflux lndev; clear; box; plot lflux lndev + # lim -n 1 lndev lMxx; clear; box; plot lndev lMxx + +end Index: /trunk/psphot/test/tap_psphot_diffsuite.pro =================================================================== --- /trunk/psphot/test/tap_psphot_diffsuite.pro (revision 30624) +++ /trunk/psphot/test/tap_psphot_diffsuite.pro (revision 30624) @@ -0,0 +1,629 @@ +#!/usr/bin/env mana +# -*-sh-*- + +if (not($?TESTING)) set TESTING = 0 + +# config for ppImage to generate chip, mask, weight +$ppImageConfig = -recipe PPIMAGE PPIMAGE_N +$ppImageConfig = $ppImageConfig -Db BACKGROUND T +$ppImageConfig = $ppImageConfig -Db CHIP.FITS T +$ppImageConfig = $ppImageConfig -Db CHIP.MASK.FITS T +$ppImageConfig = $ppImageConfig -Db CHIP.VARIANCE.FITS T +$ppImageConfig = $ppImageConfig -Db BASE.FITS F +$ppImageConfig = $ppImageConfig -Db VARIANCE.BUILD T +$ppImageConfig = $ppImageConfig -Db PHOTOM T + +# basic options for the these images (filter, location, obstype) +$BaseOptions = -type OBJECT -filter r -skymags 20.86 -ra 270.70 -dec -23.70 -pa 0.0 +$BaseOptions = $BaseOptions -Df PSASTRO:DVO.GETSTAR.MAX.RHO 50000.0 + +# options for the reference image +$RefOptions = $BaseOptions -exptime 100.0 -seeing 1.0 -D PSF.MODEL PS_MODEL_GAUSS -Df STARS.DENSITY 10.0 -Df STARS.SIGMA.LIM 0.5 + +# options for the repeated images +$FakeOptions = $BaseOptions -exptime 30.0 + +# basic config for ppSim with randomly distributed stars and gridded galaxies +$RealConfig = -camera SIMTEST -recipe PPSIM STACKTEST.RUN -D PSASTRO:PSASTRO.CATDIR catdir.ref +$RealConfig = $RealConfig -Db STARS.FAKE F -Db MATCH.DENSITY F -Db PSF.CONVOLVE F + +if (1) + $RealConfig = $RealConfig -Db GALAXY.FAKE F +else + $RealConfig = $RealConfig -Db GALAXY.FAKE T -Df GALAXY.MAG 17.0 + $RealConfig = $RealConfig -Db GALAXY.GRID T -D GALAXY.MODEL PS_MODEL_SERSIC + $RealConfig = $RealConfig -Df GALAXY.ARATIO.MIN 0.5 + $RealConfig = $RealConfig -Df GALAXY.ARATIO.MAX 0.5 + $RealConfig = $RealConfig -Df GALAXY.THETA.MAX 180 + $RealConfig = $RealConfig -Df GALAXY.INDEX.MIN 1 + $RealConfig = $RealConfig -Df GALAXY.INDEX.MAX 1 + $RealConfig = $RealConfig -Di GALAXY.GRID.DX 120 + $RealConfig = $RealConfig -Di GALAXY.GRID.DY 120 +end + +# sample alternate options: +# $ppSimOptions = $FakeOptions -D PSF.MODEL PS_MODEL_PS1_V1 +# $ppSimOptions = $FakeOptions -Df PSF.ARATIO 1.2 +# $ppSimOptions = $FakeOptions -Df PSF.THETA +30.0 +$ppSimOptions = $FakeOptions -D PSF.MODEL PS_MODEL_GAUSS + +list params + 1.00 1.00 1.0 0.0 + 1.10 1.10 1.0 0.0 + 1.25 1.25 1.0 0.0 + 1.40 1.40 1.0 0.0 + 1.00 1.00 1.3 30.0 + 1.10 1.10 1.3 30.0 + 1.25 1.25 1.3 30.0 + 1.40 1.40 1.3 30.0 + 1.00 1.00 1.3 -30.0 + 1.10 1.10 1.3 -30.0 + 1.25 1.25 1.3 -30.0 + 1.40 1.40 1.3 -30.0 + 1.00 1.25 1.0 0.0 + 1.10 1.25 1.0 0.0 + 1.25 1.25 1.0 0.0 + 1.40 1.25 1.0 0.0 + 1.00 1.25 1.3 30.0 + 1.10 1.25 1.3 30.0 + 1.25 1.25 1.3 30.0 + 1.40 1.25 1.3 30.0 + 1.00 1.25 1.3 -30.0 + 1.10 1.25 1.3 -30.0 + 1.25 1.25 1.3 -30.0 + 1.40 1.25 1.3 -30.0 +end + +# create a reference database of fake stars to be used by ppSim below +macro mkref + exec rm -rf catdir.ref + exec rm -f refimage.fits + + # create an image with fake sources and insert the resulting cmf file into a dvodb + $RefConfig = -camera SIMTEST -recipe PPSIM STACKTEST.MAKE -D PSASTRO:PSASTRO.CATDIR catdir.ref + + exec ppSim $RefOptions $RefConfig refimage + + file synth.photcodes found + if (not($found)) + echo "making photcodes file" + mkphotcodes synth.photcodes + end + + exec addstar -D CAMERA simtest -D CATDIR catdir.ref -accept-astrom -photcode SYNTH.r -D PHOTCODE_FILE synth.photcodes refimage.cmf + exec relphot -averages -D CATDIR catdir.ref -update -region 260 280 -33 -13 +end + +# create a set of input exposures using the range of parameters given in the lists +macro mkdiffsuite + if ($0 != 3) + echo "USAGE: mkdiffsuite (rawbase) (diffbase)" + break + end + + local i1 i2 rawbase diffbase + $rawbase = $1 + $diffbase = $2 + + for i1 0 $params:n + for i2 {$i1 + 1} $params:n + mkdiff $rawbase $diffbase.$i1.$i2 $i1 $i2 + end + end +end + +# create diff for a single pair of input +macro mkdiff + if ($0 != 5) + echo "USAGE: mkdiff (rawbase) (diffbase) (in1) (in2)" + break + end + + # XXX check for local reference db and run mkref if not found + + local in1 in2 diff + + $base = $1 + + $diff = $2 + dirname $diff -var dir + mkdir $dir + + $in1 = $3 + $in2 = $4 + + $command = ppSub + + $command = $command -inimage $base.$in1.wrp.fits + $command = $command -inmask $base.$in1.wrp.mask.fits + $command = $command -invariance $base.$in1.wrp.wt.fits + $command = $command -insources $base.$in1.wrp.cmf + + $command = $command -refimage $base.$in2.wrp.fits + $command = $command -refmask $base.$in2.wrp.mask.fits + $command = $command -refvariance $base.$in2.wrp.wt.fits + $command = $command -refsources $base.$in2.wrp.cmf + +# no photometry for now +# $command = $command -recipe PSPHOT DIFF +# $command = $command -F PSPHOT.PSF.SAVE PSPHOT.PSF.SKY.SAVE +# $command = $command -F PSPHOT.OUTPUT PSPHOT.OUT.CMF.MEF +# $command = $command -F PSPHOT.BACKMDL PSPHOT.BACKMDL.MEF -D PSPHOT:OUTPUT.FORMAT PS1_DV2 + + $command = $command -Db DUAL T -Db INVERSE T -convolve 1 -Df PENALTY 1000000 -Di SPATIAL.ORDER 1 + # $command = $command -Db DUAL T -Db INVERSE T -convolve 1 -Df PENALTY 3000000 -Di SPATIAL.ORDER 1 + # $command = $command -Db DUAL T -Db INVERSE T -convolve 1 -Df PENALTY 3000000000 -Di SPATIAL.ORDER 1 + $command = $command -Di KERNEL.SIZE 5 -Di STAMP.FOOTPRINT 5 + + $command = $command -save-inconv -save-refconv + $command = $command -threads 4 + $command = $command -image_id 1 -source_id 1 + + echo $command $diff + if ($TESTING) + exec $command $diff -visual ppsub 10 + else + exec $command $diff >& $diff.log + end +end + +# create a set of input exposures using the range of parameters given in the lists +macro mkexpsuite + if ($0 != 2) + echo "USAGE: mkexpsuite (rawbase)" + break + end + + local i base + + $base = $1 + + dirname $base -var dir + mkdir $dir + + for i 0 $params:n + list words -split $params:$i + mkexp $base.$i $words:0 $words:1 $words:2 $words:3 + end +end + +# create a realistic distribution of fake stars, GAUSS PSF +macro mkexp + if ($0 != 6) + echo "USAGE: mkexp basename fwhmMin fwhmMax aratio theta" + break + end + + local fwhm basename + $basename = $1 + $fwhmMin = $2 + $fwhmMax = $3 + $aratio = $4 + $theta = $5 + + echo ppSim $ppSimOptions $RealConfig $basename -seeing $fwhmMin -seeing-max $fwhmMax -seeing-ramp -Df PSF.ARATIO $aratio -Df PSF.THETA $theta + exec ppSim $ppSimOptions $RealConfig $basename -seeing $fwhmMin -seeing-max $fwhmMax -seeing-ramp -Df PSF.ARATIO $aratio -Df PSF.THETA $theta + exec /bin/mv -f $basename.cmf $basename.in.cmf + + # create the chip output + echo ppImage $ppImageConfig -file $basename.fits $basename + exec ppImage $ppImageConfig -file $basename.fits $basename + + # XXX pswarp is using all 7k+ sources to measure the PSF : can we reduce this? + echo pswarp -threads 4 -astrom $basename.cmf -file $basename.ch.fits -mask $basename.ch.mk.fits -variance $basename.ch.wt.fits $basename.wrp refimage.fits + exec pswarp -threads 4 -astrom $basename.cmf -file $basename.ch.fits -mask $basename.ch.mk.fits -variance $basename.ch.wt.fits $basename.wrp refimage.fits -D PSPHOT:OUTPUT.FORMAT PS1_V2 +end + +macro stats.pair + if ($0 != 3) + echo "USAGE: stats.pair (N) (output)" + break + end + + list word -split $spairs:$1 + if ($word:n != 8) + echo "invalid pair $1" + break + end + + $Nr = $word:3 + + reindex v1 = $word:0 using index1 + reindex v2 = $word:1 using index2 + reindex rv = $word:2 using index$Nr + + set delta = v1 - v2 + subset d1 = delta if ($word:4 < rv) && (rv < $word:5) && (abs(delta) < $word:7) + subset d2 = delta if ($word:5 < rv) && (rv < $word:6) && (abs(delta) < $word:7) + + vstats -q d1 -sigma-clip 3.0 + $M1 = $MEAN + $S1 = $SIGMA + vstats -q d2 -sigma-clip 3.0 + $M2 = $MEAN + $S2 = $SIGMA + + output $2 + fprintf "%-18s %7.4f %7.4f %7.4f %7.4f" $word:0 $M1 $S1 $M2 $S2 + output stdout +end + +macro show.pair + if ($0 != 2) + echo "USAGE: show.pair (N)" + break + end + + list word -split $pairs:$1 + if ($word:n != 7) + echo "invalid pair $1" + break + end + + $Nr = $word:3 + + reindex v1 = $word:0 using index1 + if ("$word:6" == "V") + reindex v2 = $word:1 using index2 + end + if ("$word:6" == "S") + set v2 = $word:1 + zero(index1) + end + reindex rv = $word:2 using index$Nr + + set delta = v1 - v2 + if (("$word:4" == "def") || ("$word:5" == "def")) + lim rv delta; box; plot rv delta + else + lim rv $word:4 $word:5; box; plot rv delta + end + label -y '$word:0' -x '$word:2' +end + +# This list is used to compare a pair of vectors (sans error) or a +# vector and an expected (constant) value. The last field defines a +# vector or constant for the comparison. It is assumed that the +# vector sets have been loaded and matched with match2d to generate +# index vectors 'index1' and index2'. The macro 'show.pair' generates +# a plot of the range vector vs (v1 - v2). The indices for v1, v2 are +# index1 and 2 respectively. The index for the range vector is defined +# by the integer following that vector. the y-limits of the plot are +# given by the last two numbers +list pairs + X_PSF_out X_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + Y_PSF_out Y_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + X_PSF_out X_PSF_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + Y_PSF_out Y_PSF_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + X_PSF_SIG_out X_PSF_SIG_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + Y_PSF_SIG_out Y_PSF_SIG_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + #PSF_INST_MAG_out PSF_INST_MAG_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + PSF_INST_MAG_out M_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + PSF_INST_MAG_SIG_out PSF_INST_MAG_SIG_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + #PSF_INST_FLUX_out PSF_INST_FLUX_raw PSF_INST_MAG_raw 2 def def V + #PSF_INST_FLUX_SIG_out PSF_INST_FLUX_SIG_raw PSF_INST_MAG_raw 2 def def V + AP_MAG_out M_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + AP_MAG_RAW_out M_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + AP_MAG_RADIUS_out 0.0 PSF_INST_MAG_raw 2 -0.01 20.1 S + SKY_out 0.0 PSF_INST_MAG_raw 2 def def S + SKY_SIGMA_out 0.0 PSF_INST_MAG_raw 2 def def S + PSF_CHISQ_out 1.0 PSF_INST_MAG_raw 2 def def S + CR_NSIGMA_out 0.0 PSF_INST_MAG_raw 2 def def S + EXT_NSIGMA_out 0.0 PSF_INST_MAG_raw 2 -5.01 5.01 S + PSF_MAJOR_out 0.0 PSF_INST_MAG_raw 2 -0.01 5.01 S + PSF_MINOR_out 0.0 PSF_INST_MAG_raw 2 -0.01 5.01 S + PSF_THETA_out 0.0 PSF_INST_MAG_raw 2 -1.61 1.61 S + PSF_QF_out 0.0 PSF_INST_MAG_raw 2 -0.10 1.10 S + PSF_QF_PERFECT_out 0.0 PSF_INST_MAG_raw 2 -0.10 1.10 S + PSF_NDOF_out 0.0 PSF_INST_MAG_raw 2 def def S + PSF_NPIX_out 0.0 PSF_INST_MAG_raw 2 def def S + MOMENTS_XX_out 0.0 PSF_INST_MAG_raw 2 -0.01 3.01 S + MOMENTS_XY_out 0.0 PSF_INST_MAG_raw 2 -3.01 3.01 S + MOMENTS_YY_out 0.0 PSF_INST_MAG_raw 2 -0.01 3.01 S + MOMENTS_M3C_out 0.0 PSF_INST_MAG_raw 2 -3.01 3.01 S + MOMENTS_M3S_out 0.0 PSF_INST_MAG_raw 2 -3.01 3.01 S + MOMENTS_M4C_out 0.0 PSF_INST_MAG_raw 2 -2.01 2.01 S + MOMENTS_M4S_out 0.0 PSF_INST_MAG_raw 2 -2.01 2.01 S + MOMENTS_R1_out 0.0 PSF_INST_MAG_raw 2 -5.01 5.01 S + MOMENTS_RH_out 0.0 PSF_INST_MAG_raw 2 -5.01 5.01 S + K_out M_raw PSF_INST_MAG_raw 2 def def V + KRON_FLUX_ERR_out 0.0 PSF_INST_MAG_raw 2 def def S + KRON_FLUX_INNER_out 0.0 PSF_INST_MAG_raw 2 def def S + KRON_FLUX_OUTER_out 0.0 PSF_INST_MAG_raw 2 def def S +# CAL_PSF_MAG CAL_PSF_MAG none Mraw 2 -1.01 1.01 V +# CAL_PSF_MAG_SIG CAL_PSF_MAG_SIG none Mraw 2 -1.01 1.01 V +# RA_PSF RA_PSF none Mraw 2 -1.01 1.01 V +# DEC_PSF DEC_PSF none Mraw 2 -1.01 1.01 V +# PEAK_FLUX_AS_MAG PEAK_FLUX_AS_MAG none Mraw 2 -1.01 1.01 V +# FLAGS FLAGS 0.0 Mraw 2 -1.01 1.01 S +# FLAGS2 FLAGS2 0.0 Mraw 2 -1.01 1.01 S +end + +macro load.cmf + if ($0 != 4) + echo "load.cmf (filename) (ext) (label)" + break + end + + data $1 + + # create the list of fields to load + delete -q myFields + for i 0 $fields:n + if ($?myFields) + $myFields = $myFields $fields:$i + else + $myFields = $fields:$i + end + end + + read -fits $2 $myFields + + # rename the loaded vectors appending the supplied lable + for i 0 $fields:n + set $fields:$i\_$3 = $fields:$i + delete $fields:$i + end +end + +# this list defines the fields to be loaded from file +list fields + X_PSF + Y_PSF + X_PSF_SIG + Y_PSF_SIG + PSF_INST_MAG + PSF_INST_MAG_SIG + PSF_INST_FLUX + PSF_INST_FLUX_SIG + AP_MAG + AP_MAG_RAW + AP_MAG_RADIUS + SKY + SKY_SIGMA + PSF_CHISQ + CR_NSIGMA + EXT_NSIGMA + PSF_MAJOR + PSF_MINOR + PSF_THETA + PSF_QF + PSF_QF_PERFECT + PSF_NDOF + PSF_NPIX + MOMENTS_XX + MOMENTS_XY + MOMENTS_YY + MOMENTS_M3C + MOMENTS_M3S + MOMENTS_M4C + MOMENTS_M4S + MOMENTS_R1 + MOMENTS_RH + KRON_FLUX + KRON_FLUX_ERR + KRON_FLUX_INNER + KRON_FLUX_OUTER +# CAL_PSF_MAG +# CAL_PSF_MAG_SIG +# RA_PSF +# DEC_PSF +# PEAK_FLUX_AS_MAG +# FLAGS +# FLAGS2 +end + +# use these cmf entries to measure average stats of the given pairs +list spairs + X_PSF_out X_raw PSF_INST_MAG_raw 2 -15 -10 -8.5 0.1 + Y_PSF_out Y_raw PSF_INST_MAG_raw 2 -15 -10 -8.5 0.1 + PSF_INST_MAG_out M_raw PSF_INST_MAG_raw 2 -15 -10 -8.5 0.1 + AP_MAG_out M_raw PSF_INST_MAG_raw 2 -15 -10 -8.5 0.1 +end + +# # XXX this is a hack: the cmf file created by ppSim is not +# # compatible with the image for warp (because one is a chip-mosaic +# # and the other is not) +# exec rm -f fix.hdr +# output fix.hdr +# echo "PSMOSAIC= 'CHIP ' / Mosaicked level" +# output stdout +# fix the header to be compatible with the chip file (to avoid running psphot) +# exec fits_insert $basename.cmf fix.hdr + +# if we run this test as a stand-alone program somewhere, we may need to create a local copy of the photcode file: +macro mkphotcodes + if ($0 != 2) + echo "USAGE: mkphotcodes (filename)" + break + end + + exec /bin/rm -f $1 + output $1 + echo "# airmass color astrometry mag photom astrom mask photom mask" + echo "# code name type zero slope offset c1 c2 slope zero equiv sys scale scale sys poor bad poor bad" + echo " 1 g_SYNTH sec 0.000 0.000 0.000 1 3 0.0000 0 21 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 2 r_SYNTH sec 0.000 0.000 0.000 2 3 0.0000 0 22 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3 i_SYNTH sec 0.000 0.000 0.000 2 3 0.0000 0 23 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 4 z_SYNTH sec 0.000 0.000 0.000 3 4 0.0000 0 24 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 5 y_SYNTH sec 0.000 0.000 0.000 4 5 0.0000 0 25 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3001 SYNTH.g ref 0.000 0.000 0.000 - - 0.0000 0 1 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3002 SYNTH.r ref 0.000 0.000 0.000 - - 0.0000 0 2 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3003 SYNTH.i ref 0.000 0.000 0.000 - - 0.0000 0 3 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3004 SYNTH.z ref 0.000 0.000 0.000 - - 0.0000 0 4 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3005 SYNTH.y ref 0.000 0.000 0.000 - - 0.0000 0 5 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + output stdout +end + +macro ckvar + if ($0 != 2) + echo "ckvar (buffer)" + break + end + + cursor + stats $1 {$X1-10} {$Y1-10} 21 21 + $K2sum0 = $TOTAL + stats $1 {$X2-10} {$Y2-10} 21 21 + $K2sum1 = $TOTAL + echo {$K2sum0 * 193 + $K2sum1 * 193} +end + +macro ckresid + + rd res resid.fits + rd con1 conv1.fits + rd con2 conv2.fits + dev -n im; tvch 1; tv res -50 100; tvch 2; tv con1 -20 100; tvch 3; tv con2 -20 100 +end + +macro ckchisq + data chisq2.dat + read flux2 3 Mxx 9 Myy 11 chisq 13 npix 15 + set chisqv = chisq / npix + set flux = sqrt(flux2) + + set lflux = log(flux) + set lchi = log(chisqv) + set lMxx = log(Mxx) + + set dev = sqrt(chisq) + set ndev = sqrt(chisq)/npix + set lndev = log(ndev) + + lim -n 0 lchi lMxx; clear; box; plot lchi lMxx + lim -n 1 lflux lMxx; clear; box; plot lflux lMxx + lim -n 2 lflux lchi; clear; box; plot lflux lchi + lim -n 3 lflux lndev; clear; box; plot lflux lndev + + # set flux2 = flux + # set flux = sqrt(flux2) + # set dev = sqrt(chisq) + # set ndev = sqrt(chisq)/npix + # set lndev = log(ndev) + # set lflux = log(flux) + # lim -n 0 lflux lndev; clear; box; plot lflux lndev + # lim -n 1 lndev lMxx; clear; box; plot lndev lMxx + +end + +macro ckscore + if ($0 != 3) + echo "USAGE: ckscore (diffbase) (output)" + break + end + + local i j diffbase + $diffbase = $1 + + output $2 + for i 0 $params:n + for j {$i+1} $params:n + $line = `grep applying $diffbase.$i.$j.log` + list word -split $line + $order = $word:2 + $mode = $word:4 + $val = `grep chisq: $diffbase.$i.$j.log | head -n 1` + sprintf chisq1 "%s %6.3f %s %6.3f %s %6.3f %s %6.3f %s %6.3f" $val + $val = `grep chisq: $diffbase.$i.$j.log | head -n 2 | tail -n 1` + sprintf chisq2 "%s %6.3f %s %6.3f %s %6.3f %s %6.3f %s %6.3f" $val + $val = `grep chisq: $diffbase.$i.$j.log | head -n 3 | tail -n 1` + sprintf chisq3 "%s %6.3f %s %6.3f %s %6.3f %s %6.3f %s %6.3f" $val + echo "$order $mode $chisq1 $chisq2 $chisq3 $params:$i $params:$j" + end + end + output stdout +end + +macro ckslope + if ($0 != 3) + echo "USAGE: ckslope (diffbase) (output)" + break + end + + local i j diffbase + $diffbase = $1 + + output $2 + for i 0 $params:n + for j {$i+1} $params:n + $line = `grep applying $diffbase.$i.$j.log` + list word -split $line + $order = $word:2 + $mode = $word:4 + if ($mode == 1) + $val = `grep chisq: $diffbase.$i.$j.log | head -n 1` + sprintf chisq1 "%s %6.3f %s %6.3f %s %6.3f %s %6.3f %s %6.3f" $val + $val = `grep chisq: $diffbase.$i.$j.log | head -n 4 | tail -n 1` + sprintf chisq2 "%s %6.3f %s %6.3f %s %6.3f %s %6.3f %s %6.3f" $val + end + if ($mode == 2) + $val = `grep chisq: $diffbase.$i.$j.log | head -n 2 | tail -n 1` + sprintf chisq1 "%s %6.3f %s %6.3f %s %6.3f %s %6.3f %s %6.3f" $val + $val = `grep chisq: $diffbase.$i.$j.log | head -n 5 | tail -n 1` + sprintf chisq2 "%s %6.3f %s %6.3f %s %6.3f %s %6.3f %s %6.3f" $val + end + if ($mode == 4) + $val = `grep chisq: $diffbase.$i.$j.log | head -n 3 | tail -n 1` + sprintf chisq1 "%s %6.3f %s %6.3f %s %6.3f %s %6.3f %s %6.3f" $val + $val = `grep chisq: $diffbase.$i.$j.log | head -n 6 | tail -n 1` + sprintf chisq2 "%s %6.3f %s %6.3f %s %6.3f %s %6.3f %s %6.3f" $val + end + echo "$order $mode $chisq1 $chisq2 $params:$i $params:$j" + end + end + output stdout +end + +macro plotscore + if ($0 != 2) + echo "USAGE: plotscore (input)" + break + end + + data $1 + read x1 4 f1 10 x2 14 f2 20 xd 24 fd1 30 fd2 32 + set X1 = 2 * x1 / (1 + f1) + set X2 = 2 * x2 / (1 + f2) + set Xd = 2 * xd / (fd1 + fd2) + + subset X1_1 = X1 if (X1 < X2) && (X1 < Xd) + subset X2_1 = X2 if (X1 < X2) && (X1 < Xd) + subset Xd_1 = Xd if (X1 < X2) && (X1 < Xd) + + subset X1_2 = X1 if (X1 > X2) && (X2 < Xd) + subset X2_2 = X2 if (X1 > X2) && (X2 < Xd) + subset Xd_2 = Xd if (X1 > X2) && (X2 < Xd) + + subset X1_d = X1 if (Xd < X2) && (X1 > Xd) + subset X2_d = X2 if (Xd < X2) && (X1 > Xd) + subset Xd_d = Xd if (Xd < X2) && (X1 > Xd) + + lim -n 0 0.3 1.1 0.3 1.1; clear; box + line 0 0 to 2 2 + + plot X1_1 X2_1 -pt 0 -c blue -sz 0.5 + plot X1_2 X2_2 -pt 2 -c red + plot X1_d X2_d -pt 7 -c black + + lim -n 1 0.3 1.1 0.3 1.1; clear; box + line 0 0 to 2 2 + + plot X1_1 Xd_1 -pt 0 -c blue -sz 0.5 + plot X1_2 Xd_2 -pt 2 -c red + plot X1_d Xd_d -pt 7 -c black + + lim -n 2 0.3 1.1 0.3 1.1; clear; box + line 0 0 to 2 2 + + plot X2_1 Xd_1 -pt 0 -c blue -sz 0.5 + plot X2_2 Xd_2 -pt 2 -c red + plot X2_d Xd_d -pt 7 -c black + +end + +if ($SCRIPT) + # mkexpsuite image/image + mkdiffsuite image/image diff.wt4/diff + exit 0 +end + Index: /trunk/psphot/test/tap_psphot_stack.pro =================================================================== --- /trunk/psphot/test/tap_psphot_stack.pro (revision 30624) +++ /trunk/psphot/test/tap_psphot_stack.pro (revision 30624) @@ -0,0 +1,768 @@ +#!/usr/bin/env mana +# -*-sh-*- + +# config for ppImage to generate chip, mask, weight +$ppImageConfig = -recipe PPIMAGE PPIMAGE_N +$ppImageConfig = $ppImageConfig -Db BACKGROUND T +$ppImageConfig = $ppImageConfig -Db CHIP.FITS T +$ppImageConfig = $ppImageConfig -Db CHIP.MASK.FITS T +$ppImageConfig = $ppImageConfig -Db CHIP.VARIANCE.FITS T +$ppImageConfig = $ppImageConfig -Db BASE.FITS F +$ppImageConfig = $ppImageConfig -Db VARIANCE.BUILD T +$ppImageConfig = $ppImageConfig -Db PHOTOM T + +# basic options for the these images (filter, location, obstype) +$BaseOptions = -type OBJECT -filter r -skymags 20.86 -ra 270.70 -dec -23.70 -pa 0.0 +$BaseOptions = $BaseOptions -Df PSASTRO:DVO.GETSTAR.MAX.RHO 50000.0 + +# options for the reference image +$RefOptions = $BaseOptions -exptime 100.0 -seeing 1.0 -D PSF.MODEL PS_MODEL_GAUSS -Df STARS.DENSITY 10.0 -Df STARS.SIGMA.LIM 0.5 + +# options for the repeated images +$FakeOptions = $BaseOptions -exptime 30.0 + +# basic config for ppSim with randomly distributed stars and gridded galaxies +$RealConfig = -camera SIMTEST -recipe PPSIM STACKTEST.RUN -D PSASTRO:PSASTRO.CATDIR catdir.ref +$RealConfig = $RealConfig -Db STARS.FAKE F -Db MATCH.DENSITY F -Db PSF.CONVOLVE T +$RealConfig = $RealConfig -Db GALAXY.FAKE T -Df GALAXY.MAG 17.0 +$RealConfig = $RealConfig -Db GALAXY.GRID T -D GALAXY.MODEL PS_MODEL_SERSIC +$RealConfig = $RealConfig -Df GALAXY.ARATIO.MIN 0.5 +$RealConfig = $RealConfig -Df GALAXY.ARATIO.MAX 0.5 +$RealConfig = $RealConfig -Df GALAXY.THETA.MAX 180 +$RealConfig = $RealConfig -Df GALAXY.INDEX.MIN 1 +$RealConfig = $RealConfig -Df GALAXY.INDEX.MAX 1 +$RealConfig = $RealConfig -Di GALAXY.GRID.DX 120 +$RealConfig = $RealConfig -Di GALAXY.GRID.DY 120 + +# sample alternate options: +# $ppSimOptions = $FakeOptions -D PSF.MODEL PS_MODEL_PS1_V1 +# $ppSimOptions = $FakeOptions -Df PSF.ARATIO 1.2 +# $ppSimOptions = $FakeOptions -Df PSF.THETA +30.0 +$ppSimOptions = $FakeOptions -D PSF.MODEL PS_MODEL_GAUSS + +list fwhm + 1.0 + 1.1 + 1.2 + 1.5 +end + +macro fulltest + if ($0 != 2) + echo "USAGE: fulltest (Npass)" + break + end + + mkref + + local i Npass imname stname + + $Npass = $1 + + $imname = image.fg + $stname = stack.fg + + for i 0 $Npass + sprintf N "%02d" $i + onestack $N + end +end + +macro summary.stats + if ($0 != 2) + echo "USAGE: summary.stats (Npass)" + break + end + + local i Npass base stack imname stname + + $Npass = $1 + + exec /bin/rm -f psf.image.summary.stats + exec /bin/rm -f psf.stack.inputs.summary.stats + exec /bin/rm -f psf.stack.conv.summary.stats + exec /bin/rm -f psf.stack.unconv.summary.stats + + $imname = image.fg + $stname = stack.fg + + for i 0 $Npass + sprintf base "$imname.%02d/$imname.stats" $i + sprintf stack "$stname.%02d/$stname" $i + echo $i $base $stack + exec grep PSF_INST_MAG $base >> psf.image.summary.stats + exec grep PSF_INST_MAG $stack.conv.stats >> psf.stack.inputs.summary.stats + exec grep PSF_INST_MAG $stack.pht.stats >> psf.stack.conv.summary.stats + exec grep PSF_INST_MAG $stack.unconv.pht.stats >> psf.stack.unconv.summary.stats + end +end + +# create a reference database of fake stars to be used by ppSim below +macro mkref + exec rm -rf catdir.ref + exec rm -f refimage.fits + + # create an image with fake sources and insert the resulting cmf file into a dvodb + $RefConfig = -camera SIMTEST -recipe PPSIM STACKTEST.MAKE -D PSASTRO:PSASTRO.CATDIR catdir.ref + + exec ppSim $RefOptions $RefConfig refimage + + file synth.photcodes found + if (not($found)) + echo "making photcodes file" + mkphotcodes synth.photcodes + end + + exec addstar -D CAMERA simtest -D CATDIR catdir.ref -accept-astrom -photcode SYNTH.r -D PHOTCODE_FILE synth.photcodes refimage.cmf + exec relphot -averages -D CATDIR catdir.ref -update -region 260 280 -33 -13 +end + +macro onestack + if ($0 != 2) + echo "USAGE: onestack (N)" + break + end + + mkstack image.$1/image stack.$1/stack + ckstack image.$1/image stack.$1/stack + rmstack image.$1/image stack.$1/stack +end + +# create a realistic distribution of fake stars, GAUSS PSF +macro mkstack + if ($0 != 3) + echo "USAGE: mkstack (rawbase) (stackbase)" + break + end + + local i base stack + + $base = $1 + $stack = $2 + + dirname $base -var dir + mkdir $dir + + dirname $stack -var dir + mkdir $dir + + exec echo "INPUT MULTI" > $stack.mdc + exec echo "" >> $stack.mdc + + for i 0 $fwhm:n + mkexp $base.$i $fwhm:$i + + exec echo "INPUT METADATA" >> $stack.mdc + exec echo "IMAGE STR $base.$i.wrp.fits" >> $stack.mdc + exec echo "MASK STR $base.$i.wrp.mask.fits" >> $stack.mdc + exec echo "VARIANCE STR $base.$i.wrp.wt.fits" >> $stack.mdc + exec echo "PSF STR $base.$i.wrp.psf" >> $stack.mdc + exec echo "SOURCES STR $base.$i.wrp.cmf" >> $stack.mdc + exec echo "END" >> $stack.mdc + exec echo "" >> $stack.mdc + end + + exec ppStack -threads 4 -input $stack.mdc $stack -Db TEMP.DELETE F + + basename $stack -var stackbase + for i 0 $fwhm:n + exec /bin/cp -f /tmp/$stackbase.$i.conv.im.fits $stack.$i.conv.im.fits + exec /bin/cp -f /tmp/$stackbase.$i.conv.var.fits $stack.$i.conv.var.fits + exec /bin/cp -f /tmp/$stackbase.$i.conv.mk.fits $stack.$i.conv.mk.fits + end + + # XXX note that the output convolved variance is missing the covariance + for i 0 $fwhm:n + exec psphot -file $stack.$i.conv.im.fits -mask $stack.$i.conv.mk.fits $stack.$i.conv + end + + # basic photometry for the stack and unconvolved stack + exec psphot -file $stack.fits -mask $stack.mask.fits -variance $stack.weight.fits $stack.pht + exec psphot -file $stack.unconv.fits -mask $stack.unconv.mask.fits -variance $stack.unconv.wt.fits $stack.unconv.pht +end + +# run this on the files created with mkstack +macro ckstack + if ($0 != 3) + echo "USAGE: mkstack (rawbase) (stackbase)" + break + end + + local i ZPT base stack + + $base = $1 + $stack = $2 + + # generate plots and stats for the per-image chip analysis + # and the psf-matched convolved images + $ZPT = -2.5*log($fwhm:n) + for i 0 $fwhm:n + ckchip $base.$i.in.cmf $base.$i.cmf $base.$i 0.0 + ckwarp $base.$i.in.cmf $stack.$i.conv.cmf $stack.$i.conv $ZPT + + stchip $base.$i.in.cmf $base.$i.cmf $base.stats 0.0 + stwarp $base.$i.in.cmf $stack.$i.conv.cmf $stack.conv.stats $ZPT + + completeness $base.$i.in.cmf $base.$i.cmf $base.$i.complete.png + end + + ckwarp $base.$i.in.cmf $stack.pht.cmf $stack.pht $ZPT + ckwarp $base.$i.in.cmf $stack.unconv.pht.cmf $stack.unconv.pht $ZPT + + stwarp $base.$i.in.cmf $stack.pht.cmf $stack.pht.stats $ZPT + stwarp $base.$i.in.cmf $stack.unconv.pht.cmf $stack.unconv.pht.stats $ZPT +end + +macro rmstack + if ($0 != 3) + echo "USAGE: mkstack (rawbase) (stackbase)" + break + end + + local base stack + $base = $1 + $stack = $2 + + for i 0 $fwhm:n + exec rm -f $base.$i.fits + exec rm -f $base.$i.ch.fits + exec rm -f $base.$i.ch.mk.fits + exec rm -f $base.$i.ch.wt.fits + + exec rm -f $base.$i.wrp.fits + exec rm -f $base.$i.wrp.mask.fits + exec rm -f $base.$i.wrp.wt.fits + + exec rm -f $stack.$i.conv.im.fits + exec rm -f $stack.$i.conv.mk.fits + exec rm -f $stack.$i.conv.wt.fits + end + + exec rm -f $stack.fits + exec rm -f $stack.mask.fits + exec rm -f $stack.weight.fits + + exec rm -f $stack.exp.fits + exec rm -f $stack.num.fits + exec rm -f $stack.expwt.fits + + exec rm -f $stack.unconv.fits + exec rm -f $stack.unconv.mask.fits + exec rm -f $stack.unconv.wt.fits + + exec rm -f $stack.unconv.exp.fits + exec rm -f $stack.unconv.num.fits + exec rm -f $stack.unconv.expwt.fits +end + +# create a realistic distribution of fake stars, GAUSS PSF +macro mkexp + if ($0 != 3) + echo "USAGE: mkexp basename fwhm" + break + end + + local fwhm basename + $basename = $1 + $fwhm = $2 + + # create the raw image + echo ppSim $ppSimOptions $RealConfig $basename -seeing $fwhm + exec ppSim $ppSimOptions $RealConfig $basename -seeing $fwhm + exec /bin/mv -f $basename.cmf $basename.in.cmf + + # create the chip output + echo ppImage $ppImageConfig -file $basename.fits $basename + exec ppImage $ppImageConfig -file $basename.fits $basename + + # XXX pswarp is using all 7k+ sources to measure the PSF : can we reduce this? + echo pswarp -threads 4 -astrom $basename.cmf -file $basename.ch.fits -mask $basename.ch.mk.fits -variance $basename.ch.wt.fits $basename.wrp refimage.fits + exec pswarp -threads 4 -astrom $basename.cmf -file $basename.ch.fits -mask $basename.ch.mk.fits -variance $basename.ch.wt.fits $basename.wrp refimage.fits + + echo psphot -file $basename.wrp.fits -mask $basename.wrp.mask.fits -variance $basename.wrp.wt.fits $basename.wrp.pht + exec psphot -file $basename.wrp.fits -mask $basename.wrp.mask.fits -variance $basename.wrp.wt.fits $basename.wrp.pht +end + +# compare two cmf files with extname Chip.psf +# things to compare: +# * completeness (which sources in (1) are not detected in (2) +# * positions (X_PSF, Y_PSF) +# * instrumental psf mags +# * position errors (no input errors; use a model?) +# * measured FWHM? +# * kron mags (fluxes) +# * etc, etc +macro ckchip + if ($0 != 5) + echo "USAGE: ckchip (raw) (out) (output) (zpt_off)" + break + end + + load.cmf $1 Chip.psf raw + load.cmf $2 Chip.psf out + + # images generated with convolution will not have the right output positions + set X_raw = int(X_PSF_raw) + 0.5 + set Y_raw = int(Y_PSF_raw) + 0.5 + set M_raw = PSF_INST_MAG_raw + $4 + set K_out = -2.5*log(KRON_FLUX_out) + match2d X_PSF_out Y_PSF_out X_PSF_raw Y_PSF_raw 1.5 -index1 index1 -index2 index2 + + local i NX NY nx ny N + + device -n compare + resize 1000 1000 + + # plot trends as a function of mag + $NX = 2 + $NY = 5 + $nx = 0 + $ny = 0 + $N = 0 + clear -s + for i 0 $pairs:n + section a$nx\$ny {$nx/$NX} {$ny/$NY} {1/$NX} {1/$NY} + show.pair $i + $ny ++ + if ($ny == $NY) + $ny = 0 + $nx ++ + end + if ($nx == $NX) + png -name $3.$N.png + clear -s + $nx = 0 + $ny = 0 + $N ++ + end + end + + # plot (input - output) vs mag +end + +macro stchip + if ($0 != 5) + echo "USAGE: stchip (raw) (out) (output) (zpt_off)" + break + end + + load.cmf $1 Chip.psf raw + load.cmf $2 Chip.psf out + + # images generated with convolution will not have the right output positions + set X_raw = int(X_PSF_raw) + 0.5 + set Y_raw = int(Y_PSF_raw) + 0.5 + set M_raw = PSF_INST_MAG_raw + $4 + set K_out = -2.5*log(KRON_FLUX_out) + match2d X_PSF_out Y_PSF_out X_PSF_raw Y_PSF_raw 1.5 -index1 index1 -index2 index2 + + local i + + for i 0 $spairs:n + stats.pair $i $3 + end +end + +# compare chip to warp +macro ckwarp + if ($0 != 5) + echo "USAGE: ckwarp (raw) (out) (output) (zpt_off)" + break + end + + load.cmf $1 Chip.psf raw + load.cmf $2 SkyChip.psf out + set X_raw = int(X_PSF_raw) + 0.5 + set Y_raw = int(Y_PSF_raw) + 0.5 + set M_raw = PSF_INST_MAG_raw + $4 + set K_out = -2.5*log(KRON_FLUX_out) + match2d X_PSF_out Y_PSF_out X_PSF_raw Y_PSF_raw 1.0 -index1 index1 -index2 index2 + + local i nx ny NX NY N + + device -n compare + resize 1000 1000 + + # plot trends as a function of mag + $NX = 2 + $NY = 5 + $nx = 0 + $ny = 0 + $N = 0 + clear -s + + for i 0 $pairs:n + section a$nx\$ny {$nx/$NX} {$ny/$NY} {1/$NX} {1/$NY} + show.pair $i + $ny ++ + if ($ny == $NY) + $ny = 0 + $nx ++ + end + if ($nx == $NX) + png -name $3.$N.png + clear -s + $nx = 0 + $ny = 0 + $N ++ + end + end +end + +macro stwarp + if ($0 != 5) + echo "USAGE: stwarp (raw) (out) (output) (zpt_off)" + break + end + + load.cmf $1 Chip.psf raw + load.cmf $2 SkyChip.psf out + + # images generated with convolution will not have the right output positions + set X_raw = int(X_PSF_raw) + 0.5 + set Y_raw = int(Y_PSF_raw) + 0.5 + set M_raw = PSF_INST_MAG_raw + $4 + set K_out = -2.5*log(KRON_FLUX_out) + match2d X_PSF_out Y_PSF_out X_PSF_raw Y_PSF_raw 1.5 -index1 index1 -index2 index2 + + local i + + for i 0 $spairs:n + stats.pair $i $3 + end +end + +macro stats.pair + if ($0 != 3) + echo "USAGE: stats.pair (N) (output)" + break + end + + list word -split $spairs:$1 + if ($word:n != 8) + echo "invalid pair $1" + break + end + + $Nr = $word:3 + + reindex v1 = $word:0 using index1 + reindex v2 = $word:1 using index2 + reindex rv = $word:2 using index$Nr + + set delta = v1 - v2 + subset d1 = delta if ($word:4 < rv) && (rv < $word:5) && (abs(delta) < $word:7) + subset d2 = delta if ($word:5 < rv) && (rv < $word:6) && (abs(delta) < $word:7) + + vstats -q d1 -sigma-clip 3.0 + $M1 = $MEAN + $S1 = $SIGMA + vstats -q d2 -sigma-clip 3.0 + $M2 = $MEAN + $S2 = $SIGMA + + output $2 + fprintf "%-18s %7.4f %7.4f %7.4f %7.4f" $word:0 $M1 $S1 $M2 $S2 + output stdout +end + +macro show.pair + if ($0 != 2) + echo "USAGE: show.pair (N)" + break + end + + list word -split $pairs:$1 + if ($word:n != 7) + echo "invalid pair $1" + break + end + + $Nr = $word:3 + + reindex v1 = $word:0 using index1 + if ("$word:6" == "V") + reindex v2 = $word:1 using index2 + end + if ("$word:6" == "S") + set v2 = $word:1 + zero(index1) + end + reindex rv = $word:2 using index$Nr + + set delta = v1 - v2 + if (("$word:4" == "def") || ("$word:5" == "def")) + lim rv delta; box; plot rv delta + else + lim rv $word:4 $word:5; box; plot rv delta + end + label -y '$word:0' -x '$word:2' +end + +# This list is used to compare a pair of vectors (sans error) or a +# vector and an expected (constant) value. The last field defines a +# vector or constant for the comparison. It is assumed that the +# vector sets have been loaded and matched with match2d to generate +# index vectors 'index1' and index2'. The macro 'show.pair' generates +# a plot of the range vector vs (v1 - v2). The indices for v1, v2 are +# index1 and 2 respectively. The index for the range vector is defined +# by the integer following that vector. the y-limits of the plot are +# given by the last two numbers +list pairs + X_PSF_out X_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + Y_PSF_out Y_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + X_PSF_out X_PSF_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + Y_PSF_out Y_PSF_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + X_PSF_SIG_out X_PSF_SIG_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + Y_PSF_SIG_out Y_PSF_SIG_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + #PSF_INST_MAG_out PSF_INST_MAG_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + PSF_INST_MAG_out M_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + PSF_INST_MAG_SIG_out PSF_INST_MAG_SIG_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + #PSF_INST_FLUX_out PSF_INST_FLUX_raw PSF_INST_MAG_raw 2 def def V + #PSF_INST_FLUX_SIG_out PSF_INST_FLUX_SIG_raw PSF_INST_MAG_raw 2 def def V + AP_MAG_out M_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + AP_MAG_RAW_out M_raw PSF_INST_MAG_raw 2 -1.01 1.01 V + AP_MAG_RADIUS_out 0.0 PSF_INST_MAG_raw 2 -0.01 20.1 S + SKY_out 0.0 PSF_INST_MAG_raw 2 def def S + SKY_SIGMA_out 0.0 PSF_INST_MAG_raw 2 def def S + PSF_CHISQ_out 1.0 PSF_INST_MAG_raw 2 def def S + CR_NSIGMA_out 0.0 PSF_INST_MAG_raw 2 def def S + EXT_NSIGMA_out 0.0 PSF_INST_MAG_raw 2 -5.01 5.01 S + PSF_MAJOR_out 0.0 PSF_INST_MAG_raw 2 -0.01 5.01 S + PSF_MINOR_out 0.0 PSF_INST_MAG_raw 2 -0.01 5.01 S + PSF_THETA_out 0.0 PSF_INST_MAG_raw 2 -1.61 1.61 S + PSF_QF_out 0.0 PSF_INST_MAG_raw 2 -0.10 1.10 S + PSF_QF_PERFECT_out 0.0 PSF_INST_MAG_raw 2 -0.10 1.10 S + PSF_NDOF_out 0.0 PSF_INST_MAG_raw 2 def def S + PSF_NPIX_out 0.0 PSF_INST_MAG_raw 2 def def S + MOMENTS_XX_out 0.0 PSF_INST_MAG_raw 2 -0.01 3.01 S + MOMENTS_XY_out 0.0 PSF_INST_MAG_raw 2 -3.01 3.01 S + MOMENTS_YY_out 0.0 PSF_INST_MAG_raw 2 -0.01 3.01 S + MOMENTS_M3C_out 0.0 PSF_INST_MAG_raw 2 -3.01 3.01 S + MOMENTS_M3S_out 0.0 PSF_INST_MAG_raw 2 -3.01 3.01 S + MOMENTS_M4C_out 0.0 PSF_INST_MAG_raw 2 -2.01 2.01 S + MOMENTS_M4S_out 0.0 PSF_INST_MAG_raw 2 -2.01 2.01 S + MOMENTS_R1_out 0.0 PSF_INST_MAG_raw 2 -5.01 5.01 S + MOMENTS_RH_out 0.0 PSF_INST_MAG_raw 2 -5.01 5.01 S + K_out M_raw PSF_INST_MAG_raw 2 def def V + KRON_FLUX_ERR_out 0.0 PSF_INST_MAG_raw 2 def def S + KRON_FLUX_INNER_out 0.0 PSF_INST_MAG_raw 2 def def S + KRON_FLUX_OUTER_out 0.0 PSF_INST_MAG_raw 2 def def S +# CAL_PSF_MAG CAL_PSF_MAG none Mraw 2 -1.01 1.01 V +# CAL_PSF_MAG_SIG CAL_PSF_MAG_SIG none Mraw 2 -1.01 1.01 V +# RA_PSF RA_PSF none Mraw 2 -1.01 1.01 V +# DEC_PSF DEC_PSF none Mraw 2 -1.01 1.01 V +# PEAK_FLUX_AS_MAG PEAK_FLUX_AS_MAG none Mraw 2 -1.01 1.01 V +# FLAGS FLAGS 0.0 Mraw 2 -1.01 1.01 S +# FLAGS2 FLAGS2 0.0 Mraw 2 -1.01 1.01 S +end + +macro load.cmf + if ($0 != 4) + echo "load.cmf (filename) (ext) (label)" + break + end + + data $1 + + # create the list of fields to load + delete -q myFields + for i 0 $fields:n + if ($?myFields) + $myFields = $myFields $fields:$i + else + $myFields = $fields:$i + end + end + + read -fits $2 $myFields + + # rename the loaded vectors appending the supplied lable + for i 0 $fields:n + set $fields:$i\_$3 = $fields:$i + delete $fields:$i + end +end + +# this list defines the fields to be loaded from file +list fields + X_PSF + Y_PSF + X_PSF_SIG + Y_PSF_SIG + PSF_INST_MAG + PSF_INST_MAG_SIG + PSF_INST_FLUX + PSF_INST_FLUX_SIG + AP_MAG + AP_MAG_RAW + AP_MAG_RADIUS + SKY + SKY_SIGMA + PSF_CHISQ + CR_NSIGMA + EXT_NSIGMA + PSF_MAJOR + PSF_MINOR + PSF_THETA + PSF_QF + PSF_QF_PERFECT + PSF_NDOF + PSF_NPIX + MOMENTS_XX + MOMENTS_XY + MOMENTS_YY + MOMENTS_M3C + MOMENTS_M3S + MOMENTS_M4C + MOMENTS_M4S + MOMENTS_R1 + MOMENTS_RH + KRON_FLUX + KRON_FLUX_ERR + KRON_FLUX_INNER + KRON_FLUX_OUTER +# CAL_PSF_MAG +# CAL_PSF_MAG_SIG +# RA_PSF +# DEC_PSF +# PEAK_FLUX_AS_MAG +# FLAGS +# FLAGS2 +end + +# use these cmf entries to measure average stats of the given pairs +list spairs + X_PSF_out X_raw PSF_INST_MAG_raw 2 -15 -10 -8.5 0.1 + Y_PSF_out Y_raw PSF_INST_MAG_raw 2 -15 -10 -8.5 0.1 + PSF_INST_MAG_out M_raw PSF_INST_MAG_raw 2 -15 -10 -8.5 0.1 + AP_MAG_out M_raw PSF_INST_MAG_raw 2 -15 -10 -8.5 0.1 +end + +# # XXX this is a hack: the cmf file created by ppSim is not +# # compatible with the image for warp (because one is a chip-mosaic +# # and the other is not) +# exec rm -f fix.hdr +# output fix.hdr +# echo "PSMOSAIC= 'CHIP ' / Mosaicked level" +# output stdout +# fix the header to be compatible with the chip file (to avoid running psphot) +# exec fits_insert $basename.cmf fix.hdr + +macro completeness + if ($0 != 4) + echo "USAGE: completeness (raw) (out) (output)" + break + end + + load.cmf $1 Chip.psf raw + load.cmf $2 Chip.psf out + set X_raw = int(X_PSF_raw) + 0.5 + set Y_raw = int(Y_PSF_raw) + 0.5 + match2d X_PSF_raw Y_PSF_raw X_PSF_out Y_PSF_out 1.5 -index1 index1 -index2 index2 -closest + + histogram PSF_INST_MAG_raw nMag -16.0 -3.0 0.25 -range dMag + set fMag = zero(dMag) + for i 0 {dMag[]-1} + set inrange = (PSF_INST_MAG_raw > dMag[$i]) && (PSF_INST_MAG_raw <= dMag[$i+1]) + subset all = index1 if (inrange) + subset got = index1 if (inrange) && (index1 >= 0) + if (all[] == 0) + fMag[$i] = 0 + else + fMag[$i] = got[] / all[] + end + end + + device -n complete + resize 1000 600 + + clear -s + + section default 0 0 1 1 + lim dMag fMag; clear; box -ypad 5 +ypad 5 -ticks 1110; plot -x 1 dMag fMag + label -x mag_inst -y det_frac + + set found = (index1 >= 0) + plot PSF_INST_MAG_raw found -c red + + set ldmag = log(PSF_INST_MAG_SIG_raw) + section overlay 0 0 1 1; lim dMag -5 1.2; box -ypad 5 +ypad 5 -ticks 1011 -labels 1001; plot PSF_INST_MAG_raw ldmag + label +y log(S/N) + + png -name $3 +end + +macro show.dpair + if ($0 != 2) + echo "USAGE: show.dpair (N)" + break + end + + list word -split $dpairs:$1 + if ($word:n != 7) + echo "invalid dpair $1" + break + end + + $Nr = $word:4 + + reindex v1 = $word:0 using index1 + reindex dv = $word:1 using index1 + reindex v2 = $word:2 using index2 + reindex rv = $word:3 using index$Nr + + set delta = (v1 - v2) / dv + if (("$word:5" == "def") || ("$word:6" == "def")) + lim rv delta; box; plot rv delta + else + lim rv $word:5 $word:6; box; plot rv delta + end + label -y '$word:0' -x '$word:3' +end + +# this list is used to compare a pair of vectors with an error it is +# assumed that the vector sets have been loaded and matched with +# match2d to generate index vectors 'index1' and index2'. the macro +# show.dpair generates a plot of the range vector vs (v1 - v2) / dv1. +# The indices for v1, dv1, and v2 are index1,1, and 2 respectively. The +# index for the range vector is defined by the integer following that +# vector. The y-limits of the plot are given by the last two numbers +# (use 'def') for the full default range of the delta vector +list dpairs + # v1 dv v2 range + X_PSF X_PSF_SIG Xraw Mraw 2 -10.0 10.0 + Y_PSF Y_PSF_SIG Yraw Mraw 2 -10.0 10.0 + PSF_INST_MAG PSF_INST_MAG_SIG Mraw Mraw 2 -10.0 10.0 + PSF_INST_FLUX PSF_INST_FLUX_SIG Fraw Mraw 2 -10.0 10.0 + AP_MAG PSF_INST_MAG_SIG Mraw Mraw 2 -10.0 10.0 + AP_MAG_RAW PSF_INST_MAG_SIG Mraw Mraw 2 -10.0 10.0 +end + +# if we run this test as a stand-alone program somewhere, we may need to create a local copy of the photcode file: +macro mkphotcodes + if ($0 != 2) + echo "USAGE: mkphotcodes (filename)" + break + end + + exec /bin/rm -f $1 + output $1 + echo "# airmass color astrometry mag photom astrom mask photom mask" + echo "# code name type zero slope offset c1 c2 slope zero equiv sys scale scale sys poor bad poor bad" + echo " 1 g_SYNTH sec 0.000 0.000 0.000 1 3 0.0000 0 21 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 2 r_SYNTH sec 0.000 0.000 0.000 2 3 0.0000 0 22 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3 i_SYNTH sec 0.000 0.000 0.000 2 3 0.0000 0 23 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 4 z_SYNTH sec 0.000 0.000 0.000 3 4 0.0000 0 24 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 5 y_SYNTH sec 0.000 0.000 0.000 4 5 0.0000 0 25 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3001 SYNTH.g ref 0.000 0.000 0.000 - - 0.0000 0 1 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3002 SYNTH.r ref 0.000 0.000 0.000 - - 0.0000 0 2 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3003 SYNTH.i ref 0.000 0.000 0.000 - - 0.0000 0 3 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3004 SYNTH.z ref 0.000 0.000 0.000 - - 0.0000 0 4 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + echo " 3005 SYNTH.y ref 0.000 0.000 0.000 - - 0.0000 0 5 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000" + output stdout +end + +if ($SCRIPT) + fulltest 4 + exit 0 +end