Index: /trunk/psphot/doc/stack.txt =================================================================== --- /trunk/psphot/doc/stack.txt (revision 31153) +++ /trunk/psphot/doc/stack.txt (revision 31154) @@ -1,2 +1,71 @@ + +20110303 + + psphotStack sample timing: + + * convolve images & stamps: 78 sec + * standard photometry : 250? + * petrosians : 6 sec + * ext fits : 350 sec + * radial aps : 140 sec + +20110302 + + * total number of images generated by this process: + + load: + raw : Nfilter x (image + variance + mask) (20M) + cnv : Nfilter x (image + variance + mask) (20M) + (8 * 4B) + (4 * 2B) + + match: (added) + match : Nfilter x (image + variance + mask) (20M) +(4 * 4B) + (2 * 2B) + + phot: (added) + chisq : 1 x (image + variance + mask) (10M) (freed at end of psphotStackReadout?) + backgnd : 1 x (image) (10M) or Nfilter (20M) [unsure] + + load: + 1089174 4000000 psFitsImage.c:467 + 1092246 4000000 psFitsImage.c:467 + 1095902 4000000 psFitsImage.c:467 + 1098974 4000000 psFitsImage.c:467 + 1102619 4000000 psFitsImage.c:467 + 1105669 4000000 psFitsImage.c:467 + 1109303 4000000 psFitsImage.c:467 + 1112353 4000000 psFitsImage.c:467 + 1090766 2000000 psFitsImage.c:467 + 1097494 2000000 psFitsImage.c:467 + 1104200 2000000 psFitsImage.c:467 + 1110884 2000000 psFitsImage.c:467 + + match: + 1128976 4000000 pmSubtractionMatch.c:178 + 1128979 4000000 pmSubtractionMatch.c:183 + 1652553 4000000 pmSubtractionMatch.c:178 + 1652556 4000000 pmSubtractionMatch.c:183 + 1128982 2000000 pmSubtractionMatch.c:189 + 1652559 2000000 pmSubtractionMatch.c:189 + + chisq: + 2133004 4000000 pmFPACopy.c:71 + 2133010 4000000 pmFPACopy.c:71 + 2133007 2000000 pmFPACopy.c:71 + + backmdl: + 2132413 4000000 pmFPAfileDefine.c:1275 + + ??? + 8614548 2000000 psBinaryOp.c:502 + 8617313 2000000 psBinaryOp.c:502 + + pmSource (modelFlx): + 2775 * 6k = 16.9M + + pmSource (maskObj): + 2775 * 3k = 8.4M + + (span + footprints account for ~5 - 10 M, rest in little things) 20101221 Index: /trunk/psphot/src/Makefile.am =================================================================== --- /trunk/psphot/src/Makefile.am (revision 31153) +++ /trunk/psphot/src/Makefile.am (revision 31154) @@ -69,7 +69,6 @@ psphotForced.c \ psphotForcedArguments.c \ - psphotForcedImageLoop.c \ - psphotForcedReadout.c \ psphotParseCamera.c \ + psphotImageLoop.c \ psphotMosaicChip.c \ psphotCleanup.c @@ -79,7 +78,6 @@ psphotMakePSF.c \ psphotMakePSFArguments.c \ - psphotMakePSFImageLoop.c \ - psphotMakePSFReadout.c \ psphotParseCamera.c \ + psphotImageLoop.c \ psphotMosaicChip.c \ psphotCleanup.c @@ -128,5 +126,8 @@ psphotReadoutForcedKnownSources.c \ psphotReadoutMinimal.c \ + psphotForcedReadout.c \ + psphotMakePSFReadout.c \ psphotModelBackground.c \ + psphotMaskBackground.c \ psphotSubtractBackground.c \ psphotFindDetections.c \ Index: /trunk/psphot/src/models/pmModel_STRAIL.c =================================================================== --- /trunk/psphot/src/models/pmModel_STRAIL.c (revision 31153) +++ /trunk/psphot/src/models/pmModel_STRAIL.c (revision 31154) @@ -496,5 +496,5 @@ params[PM_PAR_SKY] = Smoments->Sky; - params[PM_PAR_I0] = peak->flux; + params[PM_PAR_I0] = peak->rawFlux; params[PM_PAR_XPOS] = peak->xf; params[PM_PAR_YPOS] = peak->yf; Index: /trunk/psphot/src/models/pmModel_TEST1.c =================================================================== --- /trunk/psphot/src/models/pmModel_TEST1.c (revision 31153) +++ /trunk/psphot/src/models/pmModel_TEST1.c (revision 31154) @@ -139,5 +139,5 @@ PAR[PM_PAR_SKY] = moments->Sky; - PAR[PM_PAR_I0] = peak->flux; + PAR[PM_PAR_I0] = peak->rawFlux; PAR[PM_PAR_XPOS] = peak->xf; PAR[PM_PAR_YPOS] = peak->yf; Index: /trunk/psphot/src/psphot.c =================================================================== --- /trunk/psphot/src/psphot.c (revision 31153) +++ /trunk/psphot/src/psphot.c (revision 31154) @@ -1,3 +1,4 @@ # include "psphotStandAlone.h" +# define FORCED_PHOTOMETRY 0 int main (int argc, char **argv) { @@ -20,5 +21,5 @@ // call psphot for each readout - if (!psphotImageLoop (config)) { + if (!psphotImageLoop (config, PSPHOT_SINGLE)) { psErrorStackPrint(stderr, "Error in the psphot image loop\n"); exit (psphotGetExitStatus()); Index: /trunk/psphot/src/psphot.h =================================================================== --- /trunk/psphot/src/psphot.h (revision 31153) +++ /trunk/psphot/src/psphot.h (revision 31154) @@ -12,4 +12,12 @@ #define PSPHOT_RECIPE_PSF_FAKE_ALLOW "PSF.FAKE.ALLOW" // Name for recipe component permitting fake PSFs + +# define READOUT_OR_INTERNAL(VIEW,FILE)((FILE)->mode == PM_FPA_MODE_INTERNAL) ? (FILE)->readout : pmFPAviewThisReadout((VIEW), (FILE)->fpa) + +typedef enum { + PSPHOT_SINGLE, + PSPHOT_FORCED, + PSPHOT_MAKE_PSF, +} psphotImageLoopMode; // top-level psphot functions @@ -22,4 +30,6 @@ void psphotVersionPrint(void); +bool psphotImageLoop (pmConfig *config, psphotImageLoopMode mode); + bool psphotModelTest (pmConfig *config, const pmFPAview *view, psMetadata *recipe); bool psphotInit (void); @@ -53,4 +63,7 @@ bool psphotModelBackgroundReadoutFileIndex (pmConfig *config, const pmFPAview *view, const char *filerule, int index); +bool psphotMaskBackground (pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotMaskBackgroundReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); + bool psphotSubtractBackground (pmConfig *config, const pmFPAview *view, const char *filerule); bool psphotSubtractBackgroundReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); @@ -168,8 +181,8 @@ // used by psphotFindDetections -psImage *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, psImageMaskType maskVal); -psArray *psphotFindPeaks (psImage *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax); -bool psphotFindFootprints (pmDetections *detections, psImage *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int pass, psImageMaskType maskVal); -psErrorCode psphotCullPeaks(const pmReadout *readout, const psMetadata *recipe, psArray *footprints); +pmReadout *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, psImageMaskType maskVal); +psArray *psphotFindPeaks (pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax); +bool psphotFindFootprints (pmDetections *detections, pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int pass, psImageMaskType maskVal); +psErrorCode psphotCullPeaks(const pmReadout *readout, const pmReadout *signifRO, const psMetadata *recipe, psArray *footprints); // in psphotApResid.c: @@ -248,4 +261,5 @@ bool psphotVisualShowLogSignificance (psImage *image, float min, float max); bool psphotVisualShowPeaks (pmDetections *detections); +bool psphotVisualShowSources (psArray *sources); bool psphotVisualShowFootprints (pmDetections *detections); bool psphotVisualShowMoments (psArray *sources); @@ -305,9 +319,7 @@ pmConfig *psphotForcedArguments(int argc, char **argv); -bool psphotForcedImageLoop (pmConfig *config); bool psphotForcedReadout(pmConfig *config, const pmFPAview *view, const char *filerule); pmConfig *psphotMakePSFArguments(int argc, char **argv); -bool psphotMakePSFImageLoop (pmConfig *config); bool psphotMakePSFReadout(pmConfig *config, const pmFPAview *view, const char *filerule); @@ -430,5 +442,5 @@ 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, int entry); +bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, psImageMaskType maskVal, const psVector *radMax, int entry); bool psphotExtendedSourceAnalysisByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule); Index: /trunk/psphot/src/psphotApResid.c =================================================================== --- /trunk/psphot/src/psphotApResid.c (revision 31153) +++ /trunk/psphot/src/psphotApResid.c (revision 31154) @@ -1,3 +1,4 @@ # include "psphotInternal.h" +// # define DEBUG # define SKIPSTAR(MSG) { psTrace ("psphot", 3, "invalid : %s", MSG); continue; } @@ -116,5 +117,5 @@ // set limits on the aperture magnitudes - pmSourceMagnitudesInit (recipe); + pmSourceMagnitudesInit (config, recipe); // threaded measurement of the source magnitudes @@ -465,5 +466,5 @@ psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", markVal); - bool status = pmSourceMagnitudes (source, psf, photMode, maskVal, markVal); + bool status = pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, source->apRadius); // clear the mask bit @@ -481,5 +482,5 @@ if (!isfinite(source->apMag) || !isfinite(source->psfMag)) { Nfail ++; - psTrace ("psphot", 3, "fail : nan mags : %f %f", source->apMag, source->psfMag); + psTrace ("psphot", 3, "fail : %f, %f : nan mags : %f %f", source->peak->xf, source->peak->yf, source->apMag, source->psfMag); continue; } Index: /trunk/psphot/src/psphotArguments.c =================================================================== --- /trunk/psphot/src/psphotArguments.c (revision 31153) +++ /trunk/psphot/src/psphotArguments.c (revision 31154) @@ -109,12 +109,13 @@ } - PSARGUMENTS_INSTANTIATE_GENERICS( psphot, config, argc, argv ); + // generic arguments (version, dumpconfig) + PS_ARGUMENTS_GENERIC( psphot, config, argc, argv ); + + // thread arguments + PS_ARGUMENTS_THREADS( psphot, config, argc, argv ) // save the following additional recipe values based on command-line options // these options override the PSPHOT recipe values loaded from recipe files psMetadata *options = pmConfigRecipeOptions (config, PSPHOT_RECIPE); - - // Number of threads is handled - PSARGUMENTS_INSTANTIATE_THREADSARG( psphot, config, argc, argv ) // run the test model (requires X,Y coordinate) Index: /trunk/psphot/src/psphotBasicDeblend.c =================================================================== --- /trunk/psphot/src/psphotBasicDeblend.c (revision 31153) +++ /trunk/psphot/src/psphotBasicDeblend.c (revision 31154) @@ -63,5 +63,5 @@ for (int i = 0; i < SN->n; i++) { source = sources->data[i]; - SN->data.F32[i] = source->peak->SN; + SN->data.F32[i] = source->peak->rawFlux; } psVector *index = psVectorSortIndex (NULL, SN); @@ -115,5 +115,5 @@ // threshold is fraction of the source peak flux // image is background subtracted; source->moments->Sky should always be 0.0 - threshold = FRACTION * source->peak->SN; + threshold = FRACTION * sqrt(source->peak->detValue); // threshold is no less than NSIGMA threshold = PS_MAX (threshold, NSIGMA); @@ -138,5 +138,5 @@ testSource = overlap->data[k]; if (testSource->mode & PM_SOURCE_MODE_BLEND) continue; - if (testSource->peak->value > source->peak->value) continue; + if (testSource->peak->rawFlux > source->peak->rawFlux) continue; for (int j = 0; j < xv->n; j+=2) { if (fabs(yv->data.F32[j] - testSource->peak->y) > 0.5) continue; Index: /trunk/psphot/src/psphotBlendFit.c =================================================================== --- /trunk/psphot/src/psphotBlendFit.c (revision 31153) +++ /trunk/psphot/src/psphotBlendFit.c (revision 31154) @@ -94,4 +94,5 @@ fitOptions->weight = PS_SQR(skySig); fitOptions->mode = PM_SOURCE_FIT_PSF; + fitOptions->covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix psphotInitLimitsPSF (recipe, readout); @@ -103,5 +104,5 @@ // source analysis is done in S/N order (brightest first) - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); if (!sources->n) { psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend"); @@ -176,5 +177,5 @@ } psFree(job); - } + } } psFree (cellGroups); @@ -242,5 +243,5 @@ // limit selection to some SN limit - if (source->peak->SN < FIT_SN_LIM) continue; + if (sqrt(source->peak->detValue) < FIT_SN_LIM) continue; // exclude sources outside optional analysis region Index: /trunk/psphot/src/psphotChoosePSF.c =================================================================== --- /trunk/psphot/src/psphotChoosePSF.c (revision 31153) +++ /trunk/psphot/src/psphotChoosePSF.c (revision 31154) @@ -70,5 +70,5 @@ // examine PSF sources in S/N order (brightest first) - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); // structure to store user options defining the psf @@ -160,4 +160,6 @@ options->fitOptions->weight = PS_SQR(SKY_SIG); options->fitOptions->mode = PM_SOURCE_FIT_PSF; + options->fitOptions->covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix + psArray *stars = psArrayAllocEmpty (sources->n); @@ -178,5 +180,5 @@ psphotCheckStarDistribution (stars, sources, options); - psLogMsg ("psphot.pspsf", PS_LOG_DETAIL, "selected candidate %ld PSF objects\n", stars->n); + psLogMsg ("psphot.pspsf", PS_LOG_DETAIL, "selected %ld candidate PSF objects\n", stars->n); if (stars->n < 50) { @@ -336,4 +338,7 @@ } + // XXX does this work here? + psphotVisualShowPSFStars (recipe, try->psf, try->sources); + // build the flux-to-magnitude conversion information if (!psphotMakeFluxScale (readout->image, recipe, try->psf)) { @@ -443,4 +448,7 @@ } psFree (modelPSF); + + // float fwhmtest = pmPSFtoFWHM(psf, xc, yc); + // fprintf (stderr, "fwhm: %f, %f : %f\n", FWHM_MAJOR, FWHM_MINOR, fwhmtest); } } Index: /trunk/psphot/src/psphotCullPeaks.c =================================================================== --- /trunk/psphot/src/psphotCullPeaks.c (revision 31153) +++ /trunk/psphot/src/psphotCullPeaks.c (revision 31154) @@ -6,4 +6,5 @@ psErrorCode psphotCullPeaks(const pmReadout *readout, // the image wherein lives the footprint + const pmReadout *signifR, // smoothed image and significance const psMetadata *recipe, psArray *footprints) { // array of pmFootprints @@ -26,16 +27,74 @@ psAssert (status, "cannot find SKY_STDEV in readout->analysis"); - const float min_threshold = nsigma_min*skyStdev; + // the sky has been smoothed, so we need to scale down the raw sky stdev by this amount: + float scaleFactor = psMetadataLookupF32(&status, readout->analysis, "SIGNIFICANCE_SCALE_FACTOR"); + if (!status) scaleFactor = 1.0; + + const float MIN_THRESHOLD = nsigma_min*skyStdev / sqrt(scaleFactor); + // for saturated stars, we should be somewhat more agressive about culling: instead of + // letting the height of the intervening col (saddle point) be set by the S/N of the image + // pixel, it should be set to a fraction of the saturation value. example for a + // near-saturation pixel: + // flux = 40k, sigma = 200 + // nsigma_delta = 4.0, nsigma_min = 1.0, fPadding = 0.01, + // fStdev = 0.005, stdev_pad = 0.011*40k = 400 + // threshold = 40k - 4*400 = 38400 + // this gives too tight of a tolerance on the bright stars + // in practice, we should make the threshold much lower. + // below I am using 0.05 * saturation (eg, 2000 DN above sky in a GPC1 image) + + float SATURATION = NAN; + + // do not completely trust the values in the header... + float CELL_SATURATION = psMetadataLookupF32 (&status, readout->parent->concepts, "CELL.SATURATION"); + float MIN_SATURATION = psMetadataLookupF32 (&status, recipe, "DEBLEND_MIN_SATURATION"); + if (!status || !isfinite(MIN_SATURATION)) { + MIN_SATURATION = 40000.0; + } + if (!isfinite(CELL_SATURATION)) { + SATURATION = MIN_SATURATION; + } else { + SATURATION = PS_MAX(MIN_SATURATION, CELL_SATURATION); + } + float SAT_TEST_LEVEL = 0.50*SATURATION; + float SAT_THRESHOLD = 0.05*SATURATION; + +# if (PM_PEAKS_CULL_WITH_SMOOTHED_IMAGE) + psLogMsg ("psphot", PS_LOG_INFO, "Culling peaks from footprints using the smoothed image"); +# else + psLogMsg ("psphot", PS_LOG_INFO, "Culling peaks from footprints using the raw (unsmoothed) image"); +# endif + for (int i = 0; i < footprints->n; i++) { - // if (i % 50 == 0) fprintf (stderr, "cull %d\n", i); pmFootprint *fp = footprints->data[i]; + if (fp->peaks == NULL) continue; + if (fp->peaks->n == 0) continue; + if (fp->npix > 30000) { fprintf (stderr, "big footprint: %f %f to %f %f (%d pix)\n", fp->bbox.x0, fp->bbox.y0, fp->bbox.x1, fp->bbox.y1, fp->npix); } psTimerStart ("psphot.cull.footprints"); - if (pmFootprintCullPeaks(readout->image, readout->variance, fp, nsigma_delta, fPadding, min_threshold) != PS_ERR_NONE) { + + pmPeak *brightPeak = fp->peaks->data[0]; + float max_threshold = SAT_TEST_LEVEL; + if (brightPeak->rawFlux > SAT_TEST_LEVEL) { + max_threshold = SAT_THRESHOLD; + brightPeak->type = PM_PEAK_SUSPECT_SATURATION; + } + +# if (PM_PEAKS_CULL_WITH_SMOOTHED_IMAGE) + // New (post r30869) style of culling using the smoothed image and variance (S/N) + // if we cull using the significance image, then the definition of variance is different (thus the bool in arg 8) + if (pmFootprintCullPeaks(signifR->image, signifR->variance, fp, nsigma_delta, fPadding, MIN_THRESHOLD, max_threshold, true) != PS_ERR_NONE) { return psError(PS_ERR_UNKNOWN, false, "Culling pmFootprint %d", fp->id); } +# else + // Old (pre r30869) style of culling using the raw image and variance + if (pmFootprintCullPeaks(readout->image, readout->variance, fp, nsigma_delta, fPadding, MIN_THRESHOLD, max_threshold, false) != PS_ERR_NONE) { + return psError(PS_ERR_UNKNOWN, false, "Culling pmFootprint %d", fp->id); + } +# endif + float dtime = psTimerMark ("psphot.cull.footprints"); if (dtime > 1.0) { Index: /trunk/psphot/src/psphotDeblendSatstars.c =================================================================== --- /trunk/psphot/src/psphotDeblendSatstars.c (revision 31153) +++ /trunk/psphot/src/psphotDeblendSatstars.c (revision 31154) @@ -74,5 +74,5 @@ for (int i = 0; i < SN->n; i++) { source = sources->data[i]; - SN->data.F32[i] = source->peak->SN; + SN->data.F32[i] = source->peak->rawFlux; } psVector *index = psVectorSortIndex (NULL, SN); @@ -86,5 +86,5 @@ // XXX filter? if (source->mode & PM_SOURCE_MODE_SATSTAR) continue; if (source->mode & PM_SOURCE_MODE_BLEND) continue; - if (source->peak->flux < SATURATION) continue; + if (source->peak->rawFlux < SAT_TEST_LEVEL) continue; // save these for reference below Index: /trunk/psphot/src/psphotEllipticalProfile.c =================================================================== --- /trunk/psphot/src/psphotEllipticalProfile.c (revision 31153) +++ /trunk/psphot/src/psphotEllipticalProfile.c (revision 31154) @@ -82,5 +82,5 @@ // } - psphotPetrosianVisualProfileRadii (radius, flux, radiusRaw, fluxRaw, source->peak->flux, 0.0); + psphotPetrosianVisualProfileRadii (radius, flux, radiusRaw, fluxRaw, source->peak->rawFlux, 0.0); // psphotPetrosianVisualProfileByAngle (radius, flux); Index: /trunk/psphot/src/psphotExtendedSourceAnalysis.c =================================================================== --- /trunk/psphot/src/psphotExtendedSourceAnalysis.c (revision 31153) +++ /trunk/psphot/src/psphotExtendedSourceAnalysis.c (revision 31154) @@ -1,9 +1,6 @@ # include "psphotInternal.h" -// ?? these cannot happen here --> we would need to do this in psphotExtendedSourceAnalysis -// XXX option to choose a consistent position -// XXX option to choose a consistent elliptical contour -// XXX SDSS uses the r-band petrosian radius to measure petrosian fluxes in all bands -// XXX consistent choice of extendedness... +// measure the elliptical radial profile and use this to measure the petrosian parameters for the sources +// XXX this function needs to be threaded // for now, let's store the detections on the readout->analysis for each readout @@ -40,7 +37,5 @@ int Next = 0; int Npetro = 0; - int Nisophot = 0; int Nannuli = 0; - int Nkron = 0; psTimerStart ("psphot.extended"); @@ -68,8 +63,10 @@ assert (maskVal); - // XXX require petrosian analysis for non-linear fits? - - // XXX temporary user-supplied systematic sky noise measurement (derive from background model) - float skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE"); + // get the sky noise from the background analysis; if this is missing, get the user-supplied value + float skynoise = psMetadataLookupF32 (&status, readout->analysis, "SKY_STDEV"); + if (!status) { + skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE"); + psWarning ("failed to get sky noise level from background analysis; defaulting to user supplied value of %f\n", skynoise); + } // S/N limit to perform full non-linear fits @@ -78,17 +75,9 @@ // which extended source analyses should we perform? bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); - bool doIsophotal = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL"); bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); - bool doKron = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON"); - -# if (0) - // if backModel or backStdev are missing, the values of sky and/or skyErr will be set to NAN - // XXX use this to set skynoise - pmReadout *backModel = psphotSelectBackground (config, view); - pmReadout *backStdev = psphotSelectBackgroundStdev (config, view); -# endif + bool doPetroStars = psMetadataLookupBool (&status, recipe, "PETROSIAN_FOR_STARS"); // source analysis is done in S/N order (brightest first) - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); // option to limit analysis to a specific region @@ -103,14 +92,18 @@ // skip PSF-like and non-astronomical objects - if (source->type == PM_SOURCE_TYPE_STAR) continue; if (source->type == PM_SOURCE_TYPE_DEFECT) continue; if (source->type == PM_SOURCE_TYPE_SATURATED) continue; if (source->mode & PM_SOURCE_MODE_DEFECT) continue; if (source->mode & PM_SOURCE_MODE_SATSTAR) continue; - if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue; + + // optionally allow non-extended objects to get petrosians as well + if (!doPetroStars) { + if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue; + if (source->type == PM_SOURCE_TYPE_STAR) continue; + } // limit selection to some SN limit assert (source->peak); // how can a source not have a peak? - if (source->peak->SN < SN_LIM) continue; + if (sqrt(source->peak->detValue) < SN_LIM) continue; // limit selection by analysis region @@ -119,6 +112,4 @@ if (source->peak->x > AnalysisRegion.x1) continue; if (source->peak->y > AnalysisRegion.y1) continue; - - // fprintf (stderr, "xsrc: %f, %f : %f\n", source->peak->xf, source->peak->yf, source->peak->SN); // replace object in image @@ -136,5 +127,5 @@ // if we request any of these measurements, we require the radial profile - if (doPetrosian || doIsophotal || doAnnuli || doKron) { + if (doPetrosian || doAnnuli) { if (!psphotRadialProfile (source, recipe, skynoise, maskVal)) { // all measurements below require the radial profile; skip them all @@ -144,4 +135,5 @@ continue; } + Nannuli ++; source->mode |= PM_SOURCE_MODE_RADIAL_FLUX; } @@ -169,7 +161,5 @@ psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next); psLogMsg ("psphot", PS_LOG_INFO, " %d petrosian\n", Npetro); - psLogMsg ("psphot", PS_LOG_INFO, " %d isophotal\n", Nisophot); psLogMsg ("psphot", PS_LOG_INFO, " %d annuli\n", Nannuli); - psLogMsg ("psphot", PS_LOG_INFO, " %d kron\n", Nkron); psphotVisualShowResidualImage (readout, false); Index: /trunk/psphot/src/psphotExtendedSourceAnalysisByObject.c =================================================================== --- /trunk/psphot/src/psphotExtendedSourceAnalysisByObject.c (revision 31153) +++ /trunk/psphot/src/psphotExtendedSourceAnalysisByObject.c (revision 31154) @@ -61,5 +61,5 @@ // source analysis is done in S/N order (brightest first) - objects = psArraySort (objects, pmPhotObjSortBySN); + objects = psArraySort (objects, pmPhotObjSortByFlux); // process the objects in order. @@ -88,5 +88,5 @@ // limit selection to some SN limit assert (source->peak); // how can a source not have a peak? - if (source->peak->SN < SN_LIM) continue; + if (sqrt(source->peak->detValue) < SN_LIM) continue; measureSource = true; } @@ -148,4 +148,5 @@ psFree(source->extpars->radFlux); psFree(source->extpars->ellipticalFlux); + psFree(source->extpars->petProfile); } } Index: /trunk/psphot/src/psphotExtendedSourceFits.c =================================================================== --- /trunk/psphot/src/psphotExtendedSourceFits.c (revision 31153) +++ /trunk/psphot/src/psphotExtendedSourceFits.c (revision 31154) @@ -43,4 +43,5 @@ int NplainPass = 0; int Nfaint = 0; + int Nfail = 0; psTimerStart ("psphot.extended"); @@ -145,5 +146,5 @@ // source analysis is done in S/N order (brightest first) - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts) @@ -176,32 +177,38 @@ 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)) { + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfaint + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail + +// set this to 0 to run without threading +# if (1) + if (!psThreadJobAddPending(job)) { psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); psFree(AnalysisRegion); return false; - } else { - if (!psphotExtendedSourceFits_Threaded(job)) { - psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - psFree(AnalysisRegion); - return false; - } - psScalar *scalar = NULL; - scalar = job->args->data[7]; - Next += scalar->data.S32; - scalar = job->args->data[8]; - Nconvolve += scalar->data.S32; - scalar = job->args->data[9]; - NconvolvePass += scalar->data.S32; - scalar = job->args->data[10]; - Nplain += scalar->data.S32; - scalar = job->args->data[11]; - NplainPass += scalar->data.S32; - scalar = job->args->data[12]; - Nfaint += scalar->data.S32; - psFree(job); + } +# else + if (!psphotExtendedSourceFits_Threaded(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + psFree(AnalysisRegion); + return false; } - } + psScalar *scalar = NULL; + scalar = job->args->data[7]; + Next += scalar->data.S32; + scalar = job->args->data[8]; + Nconvolve += scalar->data.S32; + scalar = job->args->data[9]; + NconvolvePass += scalar->data.S32; + scalar = job->args->data[10]; + Nplain += scalar->data.S32; + scalar = job->args->data[11]; + NplainPass += scalar->data.S32; + scalar = job->args->data[12]; + Nfaint += scalar->data.S32; + scalar = job->args->data[13]; + Nfail += scalar->data.S32; + psFree(job); +# endif + } // wait for the threads to finish and manage results @@ -231,4 +238,6 @@ scalar = job->args->data[12]; Nfaint += scalar->data.S32; + scalar = job->args->data[13]; + Nfail += scalar->data.S32; } psFree(job); @@ -238,8 +247,8 @@ psFree(AnalysisRegion); - psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next); + psLogMsg ("psphot", PS_LOG_INFO, "extended source model fits: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next); 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); + psLogMsg ("psphot", PS_LOG_INFO, " %d too faint to fit, %d failed\n", Nfaint, Nfail); return true; } @@ -250,8 +259,9 @@ bool status; int Next = 0; + int Nfaint = 0; + int Nfail = 0; int Nconvolve = 0; int NconvolvePass = 0; int Nplain = 0; - int Nfaint = 0; int NplainPass = 0; bool savePics = false; @@ -268,7 +278,12 @@ psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA); + // pthread_t tid = pthread_self(); // Thread identifier + // Define source fitting parameters for extended source fits pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc(); - fitOptions->mode = PM_SOURCE_FIT_EXT; + fitOptions->mode = PM_SOURCE_FIT_EXT; + fitOptions->saveCovariance = true; // XXX make this a user option? + fitOptions->covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix + // XXX for now, use the defaults for the rest: // fitOptions->nIter = fitIter; @@ -296,4 +311,6 @@ // XXX use the parameters guessed from moments // if (source->modelEXT == NULL) continue; + + // fprintf (stderr, "fit %d,%d in thread %d\n", source->peak->x, source->peak->y, (int) tid); // replace object in image @@ -366,5 +383,5 @@ // limit selection to some SN limit assert (source->peak); // how can a source not have a peak? - if (source->peak->SN < SNlim) { + if (sqrt(source->peak->detValue) < SNlim) { Nfaint ++; continue; @@ -388,4 +405,5 @@ if (!modelFit) { psTrace ("psphot", 5, "failed to fit psf-conv model for object at %f, %f", source->moments->Mx, source->moments->My); + Nfail ++; continue; } @@ -402,4 +420,5 @@ if (!modelFit) { psTrace ("psphot", 5, "failed to fit plain model for object at %f, %f", source->moments->Mx, source->moments->My); + Nfail ++; continue; } @@ -519,4 +538,7 @@ scalar->data.S32 = Nfaint; + scalar = job->args->data[13]; + scalar->data.S32 = Nfail; + return true; } Index: /trunk/psphot/src/psphotFindDetections.c =================================================================== --- /trunk/psphot/src/psphotFindDetections.c (revision 31153) +++ /trunk/psphot/src/psphotFindDetections.c (revision 31154) @@ -84,11 +84,11 @@ // generate the smoothed significance image - psImage *significance = psphotSignificanceImage (readout, recipe, maskVal); + pmReadout *significance = psphotSignificanceImage (readout, recipe, maskVal); // display the log significance image - psphotVisualShowLogSignificance (significance, 0.0, 4.5); + psphotVisualShowLogSignificance (significance->variance, 0.0, 4.5); // display the significance image - psphotVisualShowSignificance (significance, 0.98*threshold, 1.02*threshold); + psphotVisualShowSignificance (significance->variance, 0.98*threshold, 1.02*threshold); // detect the peaks in the significance image Index: /trunk/psphot/src/psphotFindFootprints.c =================================================================== --- /trunk/psphot/src/psphotFindFootprints.c (revision 31153) +++ /trunk/psphot/src/psphotFindFootprints.c (revision 31154) @@ -1,5 +1,5 @@ # include "psphotInternal.h" -bool psphotFindFootprints (pmDetections *detections, psImage *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int pass, psImageMaskType maskVal) { +bool psphotFindFootprints (pmDetections *detections, pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int pass, psImageMaskType maskVal) { bool status; @@ -18,6 +18,6 @@ PS_ASSERT (status, false); - // find the raw footprints & assign the peaks to those footprints - psArray *footprints = pmFootprintsFind (significance, threshold, npixMin); + // find the raw footprints in the smoothed significance image & assign the peaks to those footprints + psArray *footprints = pmFootprintsFind (significance->variance, threshold, npixMin); if (pmFootprintsAssignPeaks(footprints, detections->peaks) != PS_ERR_NONE) { @@ -56,5 +56,5 @@ } - psphotCullPeaks(readout, recipe, detections->footprints); + psphotCullPeaks(readout, significance, recipe, detections->footprints); detections->peaks = pmFootprintArrayToPeaks(detections->footprints); psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks, %ld total footprints: %f sec\n", detections->peaks->n, detections->footprints->n, psTimerMark ("psphot.footprints")); Index: /trunk/psphot/src/psphotFindPeaks.c =================================================================== --- /trunk/psphot/src/psphotFindPeaks.c (revision 31153) +++ /trunk/psphot/src/psphotFindPeaks.c (revision 31154) @@ -4,5 +4,5 @@ // image must be constructed to represent (S/N)^2. If nMax is non-zero, only return a maximum // of nMax peaks -psArray *psphotFindPeaks (psImage *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax) { +psArray *psphotFindPeaks (pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax) { bool status = false; @@ -11,5 +11,6 @@ // find the peaks in the smoothed image - psArray *peaks = pmPeaksInImage (significance, threshold); + // NOTE : significance->variance actually carries the detection S/N image + psArray *peaks = pmPeaksInImage (significance->variance, threshold); if (peaks == NULL) { // we only get a NULL peaks array due to a programming or config error. @@ -34,12 +35,13 @@ for (int i = 0; i < peaks->n; i++) { pmPeak *peak = peaks->data[i]; - peak->SN = sqrt(peak->value); - peak->flux = readout->image->data.F32[peak->y-row0][peak->x-col0]; - // if (peak->flux / peak->value > 5.0/12.0) { - // psWarning ("odd peak levels (1)"); - // } - // if (peak->value / peak->flux > 5*12.0) { - // psWarning ("odd peak levels (2)"); - // } + peak->rawFlux = readout->image->data.F32[peak->y-row0][peak->x-col0]; + peak->rawFluxStdev = sqrt(readout->variance->data.F32[peak->y-row0][peak->x-col0]); + peak->smoothFlux = significance->image->data.F32[peak->y-row0][peak->x-col0]; + peak->smoothFluxStdev = peak->smoothFlux / sqrt(significance->variance->data.F32[peak->y-row0][peak->x-col0]); + // NOTE smoothFluxStdev is actually (sqrt(variance) / covar_factor) + + // do we need this or not? + // peak->SN = sqrt(peak->detValue); + if (readout->variance && isfinite (peak->dx)) { peak->dx *= sqrt(readout->variance->data.F32[peak->y-row0][peak->x-col0]); @@ -51,5 +53,5 @@ // limit the total number of returned peaks as specified - psArraySort (peaks, pmPeakSortBySN); + psArraySort (peaks, pmPeaksSortByRawFluxDescend); if (nMax && (peaks->n > nMax)) { psArray *tmpPeaks = psArrayAllocEmpty (nMax); Index: /trunk/psphot/src/psphotFitSet.c =================================================================== --- /trunk/psphot/src/psphotFitSet.c (revision 31153) +++ /trunk/psphot/src/psphotFitSet.c (revision 31154) @@ -27,4 +27,5 @@ pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc(); fitOptions->mode = PM_SOURCE_FIT_EXT; + fitOptions->covarFactor = 1.0; // XXX for now, use the defaults for the rest: // fitOptions->nIter = fitIter; Index: /trunk/psphot/src/psphotFitSourcesLinear.c =================================================================== --- /trunk/psphot/src/psphotFitSourcesLinear.c (revision 31153) +++ /trunk/psphot/src/psphotFitSourcesLinear.c (revision 31154) @@ -121,4 +121,6 @@ // covarFactor = 1.0; + int Nsat = 0; + // select the sources which will be used for the fitting analysis for (int i = 0; i < sources->n; i++) { @@ -134,4 +136,9 @@ // do not include CRs in the full ensemble fit if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue; + + // XXX count saturated stars + if (source->mode & PM_SOURCE_MODE_SATSTAR) { + Nsat ++; + } if (final) { @@ -167,5 +174,5 @@ 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); + source->peak->xf, source->peak->yf, modelSum, source->peak->rawFlux); } @@ -179,4 +186,6 @@ } psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "built fitSources: %f sec (%ld objects)\n", psTimerMark ("psphot.linear"), sources->n); + + // fprintf (stderr, "****** Nsat : %d ********\n", Nsat); if (fitSources->n == 0) { @@ -303,7 +312,9 @@ for (int i = 0; i < fitSources->n; i++) { pmSource *source = fitSources->data[i]; - if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) continue; pmModel *model = pmSourceGetModel (NULL, source); - pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal, covarFactor); + if (!(source->mode & PM_SOURCE_MODE_NONLINEAR_FIT)) { + model->nPar = 1; // LINEAR-only sources have 1 parameter; NONLINEAR sources have their original value + } + pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal); } psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "get chisqs: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem); @@ -324,5 +335,5 @@ // We have to place this visualization here because the models are not realized until // psphotGuessModels or fitted until psphotFitSourcesLinear. - psphotVisualShowPSFStars (recipe, psf, sources); + // psphotVisualShowPSFStars (recipe, psf, sources); return true; Index: /trunk/psphot/src/psphotFitSourcesLinearStack.c =================================================================== --- /trunk/psphot/src/psphotFitSourcesLinearStack.c (revision 31153) +++ /trunk/psphot/src/psphotFitSourcesLinearStack.c (revision 31154) @@ -163,7 +163,9 @@ for (int i = 0; i < fitSources->n; i++) { pmSource *source = fitSources->data[i]; - if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) continue; pmModel *model = pmSourceGetModel (NULL, source); - pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal, COVAR_FACTOR); + if (!(source->mode & PM_SOURCE_MODE_NONLINEAR_FIT)) { + model->nPar = 1; // LINEAR-only sources have 1 parameter; NONLINEAR sources have their original value + } + pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal); } psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "get chisqs: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem); Index: /trunk/psphot/src/psphotForced.c =================================================================== --- /trunk/psphot/src/psphotForced.c (revision 31153) +++ /trunk/psphot/src/psphotForced.c (revision 31154) @@ -1,3 +1,4 @@ # include "psphotStandAlone.h" +# define FORCED_PHOTOMETRY 1 int main (int argc, char **argv) { @@ -20,5 +21,5 @@ // call psphot for each readout - if (!psphotForcedImageLoop (config)) { + if (!psphotImageLoop (config, PSPHOT_FORCED)) { psErrorStackPrint(stderr, "Error in the psphot image loop\n"); exit (psphotGetExitStatus()); Index: /trunk/psphot/src/psphotForcedArguments.c =================================================================== --- /trunk/psphot/src/psphotForcedArguments.c (revision 31153) +++ /trunk/psphot/src/psphotForcedArguments.c (revision 31154) @@ -103,5 +103,5 @@ } - PSARGUMENTS_INSTANTIATE_GENERICS( psphot, config, argc, argv ); + PS_ARGUMENTS_GENERIC( psphot, config, argc, argv ); // save the following additional recipe values based on command-line options @@ -110,6 +110,10 @@ // Number of threads is handled - PSARGUMENTS_INSTANTIATE_THREADSARG( psphot, config, argc, argv ) + PS_ARGUMENTS_THREADS( psphot, config, argc, argv ) + if (psArgumentGet(argc, argv, "-help") || + psArgumentGet(argc, argv, "-h")) + writeHelpInfo(argv[0], config, stdout); + // visual : interactive display mode if ((N = psArgumentGet (argc, argv, "-visual"))) { @@ -176,8 +180,4 @@ } - if (psArgumentGet(argc, argv, "-help") || - psArgumentGet(argc, argv, "-h")) - writeHelpInfo(argv[0], config, stdout); - // the input file is a required argument; if not found, we will exit status = pmConfigFileSetsMD (config->arguments, &argc, argv, "INPUT", "-file", "-list"); Index: /trunk/psphot/src/psphotForcedImageLoop.c =================================================================== --- /trunk/psphot/src/psphotForcedImageLoop.c (revision 31153) +++ /trunk/psphot/src/psphotForcedImageLoop.c (revision 31154) @@ -1,9 +1,9 @@ # include "psphotStandAlone.h" -# define ESCAPE(MESSAGE) { \ - psError(PSPHOT_ERR_DATA, false, MESSAGE); \ - psFree (view); \ - return false; \ -} +# define ESCAPE(MESSAGE) { \ + psError(PSPHOT_ERR_DATA, false, MESSAGE); \ + psFree (view); \ + return false; \ + } bool psphotForcedImageLoop (pmConfig *config) { @@ -31,4 +31,10 @@ if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa in psphot."); + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + psImageMaskType maskTest = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); + // for psphot, we force data to be read at the chip level while ((chip = pmFPAviewNextChip (view, load->fpa, 1)) != NULL) { @@ -45,4 +51,21 @@ // mosaic the cells of a chip into a single contiguous (trimmed) chip if (!psphotMosaicChip(config, view, "PSPHOT.INPUT", "PSPHOT.LOAD")) ESCAPE ("Unable to mosaic chip."); + + // Read WCS if easy. + // XXX Since we're mosaicking cells, we ignore the case where the WCS is defined for a cell. + { + pmChip *inChip = pmFPAviewThisChip(view, input->fpa); // Mosaicked chip + pmHDU *hduLow = pmHDUGetLowest(input->fpa, inChip, NULL); + if (hduLow && !pmAstromReadWCS(input->fpa, inChip, hduLow->header, 1.0)) { + psWarning("Unable to read WCS astrometry from header."); + psErrorClear(); + pmHDU *hduHigh = pmHDUGetHighest(input->fpa, inChip, NULL); + if (hduHigh && hduHigh != hduLow && + !pmAstromReadWCS(input->fpa, chip, hduHigh->header, 1.0)) { + psWarning("Unable to read WCS astrometry from primary header."); + psErrorClear(); + } + } + } // try to load other supporting data (PSF, SRC, etc). @@ -76,5 +99,5 @@ if (readout->mask) { psImageMaskType maskSat = pmConfigMaskGet("SAT", config); // Mask value for saturated pixels - if (!pmReadoutMaskNonfinite(readout, maskSat)) { + if (!pmReadoutMaskInvalid(readout, maskTest, maskSat)) { psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels."); psFree(view); @@ -91,4 +114,5 @@ } + // drop all versions of the internal files status = true; status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL"); Index: /trunk/psphot/src/psphotGuessModels.c =================================================================== --- /trunk/psphot/src/psphotGuessModels.c (revision 31153) +++ /trunk/psphot/src/psphotGuessModels.c (revision 31154) @@ -160,4 +160,8 @@ pmSource *source = sources->data[i]; + if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) { + fprintf (stderr, "moment failure\n"); + } + // this is used to mark sources for which the model is measured. We check later that // all are used. @@ -172,5 +176,8 @@ // the guess central intensity comes from the peak: - float Io = source->peak->flux; + float Io = source->peak->rawFlux; + if (!isfinite(Io) && source->moments) { + Io = source->moments->Peak; + } // We have two options to get a guess for the object position: the position from the @@ -178,14 +185,5 @@ // moments -# if (1) - bool useMoments = true; -# else - bool useMoments = false; - useMoments = (source->mode & PM_SOURCE_MODE_SATSTAR); // we only want to try if SATSTAR is set, but.. -# endif - - useMoments = (useMoments && source->moments); // can't if there are no moments - useMoments = (useMoments && source->moments->nPixels); // can't if the moments were not measured - useMoments = (useMoments && !(source->mode && PM_SOURCE_MODE_MOMENTS_FAILURE)); // can't if the moments failed... + bool useMoments = pmSourcePositionUseMoments(source); float Xo, Yo; @@ -198,4 +196,11 @@ } + if (source->mode & PM_SOURCE_MODE_SATSTAR) { + fprintf (stderr, "satstar: %f,%f vs %f,%f : %c\n", + source->moments->Mx, source->moments->My, + source->peak->xf, source->peak->yf, + (useMoments ? 'T' : 'F')); + } + // set PSF parameters for this model (apply 2D shape model to coordinates Xo, Yo) pmModel *modelPSF = pmModelFromPSFforXY(psf, Xo, Yo, Io); Index: /trunk/psphot/src/psphotImageLoop.c =================================================================== --- /trunk/psphot/src/psphotImageLoop.c (revision 31153) +++ /trunk/psphot/src/psphotImageLoop.c (revision 31154) @@ -1,11 +1,11 @@ # include "psphotStandAlone.h" -# define ESCAPE(MESSAGE) { \ - psError(PSPHOT_ERR_DATA, false, MESSAGE); \ - psFree (view); \ - return false; \ -} +# define ESCAPE(MESSAGE) { \ + psError(PSPHOT_ERR_DATA, false, MESSAGE); \ + psFree (view); \ + return false; \ + } -bool psphotImageLoop (pmConfig *config) { +bool psphotImageLoop (pmConfig *config, psphotImageLoopMode mode) { bool status; @@ -90,10 +90,8 @@ // Update the header - { - pmHDU *hdu = pmHDUGetHighest(input->fpa, chip, cell); - if (hdu && hdu != lastHDU) { - psphotVersionHeaderFull(hdu->header); - lastHDU = hdu; - } + pmHDU *hdu = pmHDUGetHighest(input->fpa, chip, cell); + if (hdu && hdu != lastHDU) { + psphotVersionHeaderFull(hdu->header); + lastHDU = hdu; } @@ -109,9 +107,27 @@ // run the actual photometry analysis on this chip/cell/readout - if (!psphotReadout (config, view, "PSPHOT.INPUT")) { - psError(psErrorCodeLast(), false, "failure in psphotReadout for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout); - psFree (view); - return false; - } + switch (mode) { + case PSPHOT_SINGLE: + if (!psphotReadout (config, view, "PSPHOT.INPUT")) { + psError(psErrorCodeLast(), false, "failure in psphotReadout for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout); + psFree (view); + return false; + } + break; + case PSPHOT_FORCED: + if (!psphotForcedReadout (config, view, "PSPHOT.INPUT")) { + psError(psErrorCodeLast(), false, "failure in psphotReadout for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout); + psFree (view); + return false; + } + break; + case PSPHOT_MAKE_PSF: + if (!psphotMakePSFReadout (config, view, "PSPHOT.INPUT")) { + psError(psErrorCodeLast(), false, "failure in psphotReadout for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout); + psFree (view); + return false; + } + break; + } } Index: /trunk/psphot/src/psphotLoadSRCTEXT.c =================================================================== --- /trunk/psphot/src/psphotLoadSRCTEXT.c (revision 31153) +++ /trunk/psphot/src/psphotLoadSRCTEXT.c (revision 31154) @@ -78,5 +78,4 @@ source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE); - source->peak->flux = peakFlux; source->peak->dx = dPAR[PM_PAR_XPOS]; source->peak->dy = dPAR[PM_PAR_YPOS]; Index: /trunk/psphot/src/psphotMagnitudes.c =================================================================== --- /trunk/psphot/src/psphotMagnitudes.c (revision 31153) +++ /trunk/psphot/src/psphotMagnitudes.c (revision 31154) @@ -76,5 +76,5 @@ maskVal |= markVal; - pmSourceMagnitudesInit (recipe); + pmSourceMagnitudesInit (config, recipe); // the binning details are saved on the analysis metadata @@ -176,5 +176,5 @@ psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", markVal); - status = pmSourceMagnitudes (source, psf, photMode, maskVal, markVal); + status = pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, source->apRadius); if (status && isfinite(source->apMag)) Nap ++; @@ -295,5 +295,4 @@ psArray *sources = job->args->data[0]; psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[1],PS_TYPE_IMAGE_MASK_DATA); - psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[2],PS_TYPE_IMAGE_MASK_DATA); for (int i = 0; i < sources->n; i++) { @@ -309,5 +308,5 @@ } - status = pmSourcePixelWeight (&source->pixWeightNotBad, &source->pixWeightNotPoor, model, source->maskObj, maskVal, markVal); + status = pmSourcePixelWeight (source, model, source->maskObj, maskVal, source->apRadius); if (!status) { psTrace ("psphot", 3, "fail to measure pixel weight"); Index: /trunk/psphot/src/psphotMakeGrowthCurve.c =================================================================== --- /trunk/psphot/src/psphotMakeGrowthCurve.c (revision 31153) +++ /trunk/psphot/src/psphotMakeGrowthCurve.c (revision 31154) @@ -8,5 +8,5 @@ // set limits on the aperture magnitudes - pmSourceMagnitudesInit (recipe); + pmSourceMagnitudesInit (NULL, recipe); // bit-masks to test for good/bad pixels Index: /trunk/psphot/src/psphotMakePSF.c =================================================================== --- /trunk/psphot/src/psphotMakePSF.c (revision 31153) +++ /trunk/psphot/src/psphotMakePSF.c (revision 31154) @@ -20,5 +20,5 @@ // call psphot for each readout - if (!psphotMakePSFImageLoop (config)) { + if (!psphotImageLoop (config, PSPHOT_MAKE_PSF)) { psErrorStackPrint(stderr, "Error in the psphot image loop\n"); exit (psphotGetExitStatus()); Index: /trunk/psphot/src/psphotMakePSFArguments.c =================================================================== --- /trunk/psphot/src/psphotMakePSFArguments.c (revision 31153) +++ /trunk/psphot/src/psphotMakePSFArguments.c (revision 31154) @@ -103,5 +103,5 @@ } - PSARGUMENTS_INSTANTIATE_GENERICS( psphot, config, argc, argv ); + PS_ARGUMENTS_GENERIC( psphot, config, argc, argv ); // save the following additional recipe values based on command-line options @@ -110,5 +110,5 @@ // Number of threads is handled - PSARGUMENTS_INSTANTIATE_THREADSARG( psphot, config, argc, argv ) + PS_ARGUMENTS_THREADS( psphot, config, argc, argv ) // visual : interactive display mode Index: /trunk/psphot/src/psphotMaskBackground.c =================================================================== --- /trunk/psphot/src/psphotMaskBackground.c (revision 31154) +++ /trunk/psphot/src/psphotMaskBackground.c (revision 31154) @@ -0,0 +1,102 @@ +# include "psphotInternal.h" +# define NSIGMA 2.0 + +// mask pixel in the input image above model + N*stdev (using mark) +bool psphotMaskBackgroundReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) +{ + bool status = true; + + psTimerStart ("psphot.background"); + + // 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?"); + + // find the currently selected readout (XXX note that the model is saved on PSPHOT.BACKMDL regardless of 'filename' + pmFPAfile *modelFile = pmFPAfileSelectSingle(config->files, "PSPHOT.BACKMDL", index); // File of interest + assert (modelFile); + + pmFPAfile *stdevFile = pmFPAfileSelectSingle(config->files, "PSPHOT.BACKMDL.STDEV", index); + assert (stdevFile); + + float skyMean = psMetadataLookupF32(&status, readout->analysis, "SKY_MEAN"); + float skyStdv = psMetadataLookupF32(&status, readout->analysis, "SKY_STDEV"); + + pmReadout *model = READOUT_OR_INTERNAL(view, modelFile); + psAssert (model, "this must exist"); + + pmReadout *stdev = READOUT_OR_INTERNAL(view, stdevFile); + psAssert (stdev, "this must exist"); + + // 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); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); // Mask value for bad pixels + assert (markVal); + + psImageBinning *binning = psMetadataLookupPtr(&status, model->analysis, "PSPHOT.BACKGROUND.BINNING"); + assert (binning); + + // create a binned version of the threshold image (= model + NSIGMA*stdev) + psImage *threshBinned = psImageCopy (NULL, model->image, PS_TYPE_F32); + for (int iy = 0; iy < model->image->numRows; iy++) { + for (int ix = 0; ix < model->image->numCols; ix++) { + // threshBinned->data.F32[iy][ix] = model->image->data.F32[iy][ix] + NSIGMA*stdev->image->data.F32[iy][ix]; + threshBinned->data.F32[iy][ix] = skyMean + NSIGMA*skyStdv; + } + } + psphotSaveImage (NULL, threshBinned, "threshbin.fits"); + + // create a threshold image + psImage *threshold = psImageCopy (NULL, readout->image, PS_TYPE_F32); + + // linear interpolation to full-scale + if (!psImageUnbin (threshold, threshBinned, binning)) { + psError (PSPHOT_ERR_PROG, true, "inconsistent sizes for unbinning"); + return false; + } + psphotSaveImage (NULL, threshold, "threshold.fits"); + + psLogMsg ("psphot", PS_LOG_MINUTIA, "build threshold image: %f sec\n", psTimerMark ("psphot.background")); + + // raise the 'markVal' mask bit for pixels above threshold + for (int iy = 0; iy < readout->image->numRows; iy++) { + for (int ix = 0; ix < readout->image->numCols; ix++) { + if (readout->mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) continue; + if (readout->image->data.F32[iy][ix] < threshold->data.F32[iy][ix]) continue; + readout->mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= markVal; + } + } + psphotSaveImage (NULL, readout->mask, "newmask.fits"); + psLogMsg ("psphot", PS_LOG_INFO, "masked image based on 1st pass background model: %f sec\n", psTimerMark ("psphot.background")); + + psFree (threshold); + psFree (threshBinned); + + return true; +} + +bool psphotMaskBackground (pmConfig *config, const pmFPAview *view, const char *filerule) +{ + bool status = false; + + // 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 (!psphotMaskBackgroundReadout (config, view, filerule, i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to subtract background for %s entry %d", filerule, i); + return false; + } + } + return true; +} Index: /trunk/psphot/src/psphotMergeSources.c =================================================================== --- /trunk/psphot/src/psphotMergeSources.c (revision 31153) +++ /trunk/psphot/src/psphotMergeSources.c (revision 31154) @@ -1,7 +1,8 @@ # include "psphotInternal.h" +// Mask to apply for PSF sources : only exclude bad sources -- we will re-test for extendedness #define PSF_SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_SATSTAR | PM_SOURCE_MODE_BLEND | \ PM_SOURCE_MODE_BADPSF | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | \ - PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply for PSF sources + PM_SOURCE_MODE_CR_LIMIT) // for now, let's store the detections on the readout->analysis for each readout @@ -110,6 +111,7 @@ // the supplied peak flux needs to be re-normalized - source->peak->flux = 1.0; - source->peak->value = 1.0; + source->peak->rawFlux = 1.0; + source->peak->smoothFlux = 1.0; + source->peak->detValue = 1.0; // drop the loaded source modelPSF @@ -281,5 +283,6 @@ for (int i = 0; i < sources->n; i++) { pmSource *source = sources->data[i]; - source->peak->flux = source->moments->Peak; + source->peak->rawFlux = source->moments->Peak; + source->peak->smoothFlux = source->moments->Peak; } @@ -331,14 +334,7 @@ float ypos = model->params->data.F32[PM_PAR_YPOS]; - pmPeak *peak = pmPeakAlloc(xpos, ypos, 1.0, PM_PEAK_LONE); + pmPeak *peak = pmPeakAlloc(xpos, ypos, flux, PM_PEAK_LONE); peak->xf = xpos; peak->yf = ypos; - peak->flux = flux; // this are being set wrong, but does it matter? - - if (isfinite (source->errMag) && (source->errMag > 0.0)) { - peak->SN = 1.0 / source->errMag; - } else { - peak->SN = 0.0; - } psArrayAdd (detections->peaks, 100, peak); @@ -347,4 +343,6 @@ psLogMsg ("psphot", 3, "%ld PSF sources loaded", detections->peaks->n); + psphotVisualShowSources (sources); + psphotVisualShowPeaks (detections); // save detections on the readout->analysis @@ -354,8 +352,10 @@ } psFree (detections); + return true; } // generate the detection structure for the supplied array of sources +// XXX this function is currently unused bool psphotSetSourceParams (pmConfig *config, psArray *sources, pmPSF *psf) { @@ -378,14 +378,9 @@ float ypos = model->params->data.F32[PM_PAR_YPOS]; - pmPeak *peak = pmPeakAlloc(xpos, ypos, 1.0, PM_PEAK_LONE); + pmPeak *peak = pmPeakAlloc(xpos, ypos, flux, PM_PEAK_LONE); peak->xf = xpos; peak->yf = ypos; - peak->flux = flux; // this are being set wrong, but does it matter? - - if (isfinite (source->errMag) && (source->errMag > 0.0)) { - peak->SN = 1.0 / source->errMag; - } else { - peak->SN = 0.0; - } + peak->rawFlux = flux; // this are being set wrong, but does it matter? + peak->smoothFlux = flux; // this are being set wrong, but does it matter? source->peak = peak; @@ -613,5 +608,7 @@ // allocate image, weight, mask for the new image for each peak - pmSourceRedefinePixels (sourceOut, readoutOut, sourceOut->peak->x, sourceOut->peak->y, sourceOut->modelPSF->fitRadius); + if (sourceOut->modelPSF) { + 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 @@ -679,7 +676,7 @@ objectsOut->data[k] = objectOut; - objectOut->SN = objectSrc->SN; - objectOut->x = objectSrc->x; - objectOut->y = objectSrc->y; + objectOut->flux = objectSrc->flux; + objectOut->x = objectSrc->x; + objectOut->y = objectSrc->y; objectOut->sources = psArrayAlloc(objectSrc->sources->n); @@ -697,5 +694,6 @@ sourceOut->parent = sourceSrc; - // keep the original source flags + // keep the original source flags and sequence ID (if set) + sourceOut->seq = sourceSrc->seq; sourceOut->type = sourceSrc->type; sourceOut->mode = sourceSrc->mode; @@ -723,5 +721,7 @@ // allocate image, weight, mask for the new image for each peak - pmSourceRedefinePixels (sourceOut, readout, sourceOut->peak->x, sourceOut->peak->y, sourceOut->modelPSF->fitRadius); + if (sourceOut->modelPSF) { + 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 Index: /trunk/psphot/src/psphotModelBackground.c =================================================================== --- /trunk/psphot/src/psphotModelBackground.c (revision 31153) +++ /trunk/psphot/src/psphotModelBackground.c (revision 31154) @@ -59,5 +59,16 @@ assert (maskVal); + // mark pixel may be used for optional iteration -- is not required by all processes + psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); // Mask value for bad pixels + + // apply both MASK and MARK to make background model + maskVal |= markVal; + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); + + float dXsample = psMetadataLookupF32(&status, recipe, "BACKGROUND.XSAMPLE"); + if (!status) dXsample = 1.0; + float dYsample = psMetadataLookupF32(&status, recipe, "BACKGROUND.YSAMPLE"); + if (!status) dYsample = 1.0; // subtract this amount extra from the sky @@ -78,8 +89,5 @@ PS_STAT_ROBUST_QUARTILE | PS_STAT_CLIPPED_MEAN | - PS_STAT_FITTED_MEAN | - PS_STAT_FITTED_MEAN_V2 | - PS_STAT_FITTED_MEAN_V3 | - PS_STAT_FITTED_MEAN_V4))) { + PS_STAT_FITTED_MEAN))) { statsOptionLocation = PS_STAT_FITTED_MEAN; } @@ -89,19 +97,9 @@ statsOptionWidth = PS_STAT_SAMPLE_STDEV; } else if (statsOptionLocation & (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_QUARTILE)) { -#if 1 - statsOptionWidth = PS_STAT_ROBUST_STDEV; // not set; => NaN -#else - statsOptionWidth = PS_STAT_FITTED_STDEV; -#endif + statsOptionWidth = PS_STAT_ROBUST_STDEV; } else if (statsOptionLocation & PS_STAT_FITTED_MEAN) { statsOptionWidth = PS_STAT_FITTED_STDEV; } else if (statsOptionLocation & PS_STAT_CLIPPED_MEAN) { statsOptionWidth = PS_STAT_CLIPPED_STDEV; - } else if (statsOptionLocation & PS_STAT_FITTED_MEAN_V2) { - statsOptionWidth = PS_STAT_FITTED_STDEV_V2; - } else if (statsOptionLocation & PS_STAT_FITTED_MEAN_V3) { - statsOptionWidth = PS_STAT_FITTED_STDEV_V3; - } else if (statsOptionLocation & PS_STAT_FITTED_MEAN_V4) { - statsOptionWidth = PS_STAT_FITTED_STDEV_V4; } else { psAbort("Unable to estimate variance of selected statsOptionLocations 0x%x", statsOptionLocation); @@ -132,5 +130,5 @@ stats->clipSigma = psMetadataLookupF32 (&status, recipe, "SKY_CLIP_SIGMA"); if (!status) { - if ((stats->options & PS_STAT_FITTED_MEAN) || (stats->options & PS_STAT_FITTED_MEAN_V2)) { + if (stats->options & PS_STAT_FITTED_MEAN) { stats->clipSigma = 1.0; } else { @@ -154,4 +152,8 @@ int nFailures = 0; psImageBackgroundInit(); + + // we have Nx * Ny model points, but we can use a window which is larger (or smaller) than + // 1 superpixel. If we have a window of size dXsample * dYsample, then the regions run from: + // (ix + 0.5) - dX to (ix + 0.5) + dX // measure clipped median for subimages @@ -162,6 +164,6 @@ // convert the ruff grid cell to the equivalent fine grid cell - // XXX we need to watch out for row0,col0 - ruffRegion = psRegionSet (ix, ix + 1, iy, iy + 1); + ruffRegion = psRegionSet (ix + 0.5 - 0.5*dXsample, ix + 0.5 + 0.5*dXsample, + iy + 0.5 - 0.5*dYsample, iy + 0.5 + 0.5*dYsample); fineRegion = psImageBinningSetFineRegion (binning, ruffRegion); fineRegion = psRegionForImage (image, fineRegion); @@ -199,8 +201,8 @@ if (gotX && gotY) { - (void) psTraceSetLevel ("psLib.math.vectorFittedStats_v4", 6); + (void) psTraceSetLevel ("psLib.math.vectorFittedStats", 6); (void) psTraceSetLevel ("psLib.math.vectorRobustStats", 6); } else { - (void) psTraceSetLevel ("psLib.math.vectorFittedStats_v4", 0); + (void) psTraceSetLevel ("psLib.math.vectorFittedStats", 0); (void) psTraceSetLevel ("psLib.math.vectorRobustStats", 0); } @@ -264,5 +266,5 @@ for (int iy = 0; iy < model->numRows; iy++) { for (int ix = 0; ix < model->numCols; ix++) { - if (!isnan(modelData[iy][ix])) { + if (isfinite(modelData[iy][ix])) { Value += modelData[iy][ix]; ValueStdev += modelStdevData[iy][ix]; @@ -280,9 +282,13 @@ if (jx < 0) continue; if (jx >= model->numCols) continue; + if (!isfinite(modelData[jy][jx])) continue; value += modelData[jy][jx]; count += 1.0; } } - if (count > 0) modelData[iy][ix] = value / count; + if (count > 0) { + psLogMsg ("psphot", PS_LOG_DETAIL, "patching background %d, %d: %f (%d pts)\n", ix, iy, (value / count), (int) count); + modelData[iy][ix] = value / count; + } } } @@ -306,4 +312,12 @@ modelStdevData[iy][ix] = ValueStdev; } + } + + if (psTraceGetLevel("psphot") > 5) { + char name[256]; + sprintf (name, "backfill.%02d.fits", npass); + psphotSaveImage (NULL, model, name); + sprintf (name, "backfist.%02d.fits", npass); + psphotSaveImage (NULL, modelStdev, name); } Index: /trunk/psphot/src/psphotMosaicSubimage.c =================================================================== --- /trunk/psphot/src/psphotMosaicSubimage.c (revision 31153) +++ /trunk/psphot/src/psphotMosaicSubimage.c (revision 31154) @@ -29,5 +29,5 @@ inRegion = psRegionForImage (inImage, inRegion); - float peak = source->peak->flux; + float peak = source->peak->rawFlux; psImage *subImage = psImageSubset (inImage, inRegion); Index: /trunk/psphot/src/psphotOutput.c =================================================================== --- /trunk/psphot/src/psphotOutput.c (revision 31153) +++ /trunk/psphot/src/psphotOutput.c (revision 31154) @@ -29,33 +29,24 @@ } -pmReadout *psphotSelectBackground (pmConfig *config, - const pmFPAview *view) { +pmReadout *psphotSelectBackground (pmConfig *config, const pmFPAview *view) { bool status; - pmReadout *background; + pmFPAfile *file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKMDL"); if (!file) return NULL; - if (file->mode == PM_FPA_MODE_INTERNAL) { - background = file->readout; - } else { - background = pmFPAviewThisReadout (view, file->fpa); - } + + pmReadout *background = READOUT_OR_INTERNAL(view, file); return background; } -pmReadout *psphotSelectBackgroundStdev (pmConfig *config, - const pmFPAview *view) { +pmReadout *psphotSelectBackgroundStdev (pmConfig *config, const pmFPAview *view) { bool status; - pmReadout *background; pmFPAfile *file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKMDL.STDEV"); if (!file) return NULL; - if (file->mode == PM_FPA_MODE_INTERNAL) { - background = file->readout; - } else { - background = pmFPAviewThisReadout (view, file->fpa); - } + + pmReadout *background = READOUT_OR_INTERNAL(view, file); return background; } @@ -79,5 +70,5 @@ // float mpeak = model ? model->params->data.F32[PM_PAR_I0] : NAN; // bool subtracted = source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED; - // fprintf (stderr, "%d %d : %d : %f %f : %f %f\n", source->peak->x, source->peak->y, subtracted, source->peak->flux, source->pixels->data.F32[yc][xc], mcore, mpeak); + // fprintf (stderr, "%d %d : %d : %f %f : %f %f\n", source->peak->x, source->peak->y, subtracted, source->peak->rawFlux, source->pixels->data.F32[yc][xc], mcore, mpeak); fprintf (f, "%6.1f %6.1f : %6.1f %6.1f : %8.3f %8.3f %8.3f : %f : %f %f %f : %f\n", @@ -445,5 +436,5 @@ if (!source->moments) continue; - float Io = source->peak->flux; + float Io = source->peak->rawFlux; fprintf (f, "%f %f : %f %f :: %f\n", Index: /trunk/psphot/src/psphotPetrosianAnalysis.c =================================================================== --- /trunk/psphot/src/psphotPetrosianAnalysis.c (revision 31153) +++ /trunk/psphot/src/psphotPetrosianAnalysis.c (revision 31154) @@ -29,5 +29,5 @@ // source analysis is done in S/N order (brightest first) - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); // choose the sources of interest Index: /trunk/psphot/src/psphotPetrosianRadialBins.c =================================================================== --- /trunk/psphot/src/psphotPetrosianRadialBins.c (revision 31153) +++ /trunk/psphot/src/psphotPetrosianRadialBins.c (revision 31154) @@ -51,4 +51,5 @@ psVector *binRad = psVectorAllocEmpty(nMax, PS_TYPE_F32); // mean radius of radial bin psVector *binArea = psVectorAllocEmpty(nMax, PS_TYPE_F32); // area of radial bin (contiguous, non-overlapping) + psVector *binFill = psVectorAllocEmpty(nMax, PS_TYPE_F32); // fraction of radial bin with valid pixels psVectorInit (binSB, 0.0); @@ -144,4 +145,5 @@ binSB->data.F32[nOut] = value; binSBstdev->data.F32[nOut] = sqrt(PS_SQR(dvalue) / values->n + skyModelErrorSQ); + binFill->data.F32[nOut] = values->n / binArea->data.F32[nOut]; // error in the SB is the stdev per bin / sqrt (number of pixels) @@ -176,4 +178,5 @@ psFree(binRad); psFree(binArea); + psFree(binFill); psFree(radMin); psFree(radMax); @@ -202,8 +205,10 @@ psFree(profile->radialBins); psFree(profile->area); + psFree(profile->binFill); // save the vectors profile->radialBins = binRad; profile->area = binArea; + profile->binFill = binFill; profile->binSB = binSB; profile->binSBstdev = binSBstdev; Index: /trunk/psphot/src/psphotPetrosianStats.c =================================================================== --- /trunk/psphot/src/psphotPetrosianStats.c (revision 31153) +++ /trunk/psphot/src/psphotPetrosianStats.c (revision 31154) @@ -6,5 +6,7 @@ // generate the Petrosian radius and flux from the mean surface brightness (r_i) -float InterpolateValues (float X0, float Y0, float X1, float Y1, float X); +float InterpolateValues (float X0, float Y0, float X1, float Y1, float X); +float InterpolateValuesErrX (float X0, float Y0, float X1, float Y1, float X, float dX0, float dX1); +float InterpolateValuesErrY (float X0, float Y0, float X1, float Y1, float X, float dY0, float dY1); bool psphotPetrosianStats (pmSource *source) { @@ -25,4 +27,5 @@ psVector *binRad = profile->radialBins; psVector *area = profile->area; + psVector *binFill = profile->binFill; psVector *fluxSum = psVectorAllocEmpty(binSB->n, PS_TYPE_F32); @@ -33,7 +36,12 @@ psVector *meanSB = psVectorAllocEmpty(binSB->n, PS_TYPE_F32); psVector *areaSum = psVectorAllocEmpty(binSB->n, PS_TYPE_F32); + psVector *apixSum = psVectorAllocEmpty(binSB->n, PS_TYPE_F32); float petRadius = NAN; + float petRadiusErr = NAN; float petFlux = NAN; + float petFluxErr = NAN; + float petArea = NAN; + float petApix = NAN; bool anyPetro = false; @@ -41,4 +49,5 @@ bool above = true; float Asum = 0.0; + float Psum = 0.0; float Fsum = 0.0; float dFsum2 = 0.0; @@ -56,7 +65,8 @@ float Area = area->data.F32[i]; - Asum += Area; + Asum += Area; // Asum is the cumulative area interior to this bin Fsum += binSB->data.F32[i] * Area; dFsum2 += PS_SQR(binSBstdev->data.F32[i] * Area); + Psum += Area*binFill->data.F32[i]; // Psum is the cumulative number of pixels interior to this bin float areaInner = 0.5 * Area; @@ -87,14 +97,15 @@ psVectorAppend(areaSum, Asum); + psVectorAppend(apixSum, Psum); psVectorAppend(fluxSum, Fsum); psVectorAppend(fluxSumErr2, dFsum2); psVectorAppend(refRadius, binRad->data.F32[i]); - psTrace ("psphot", 4, "%3d : %5.2f : %5.3f %5.3f : %5.3f %5.3f : %5.3f %5.3f : %5.3f %5.3f : %5.1f %5.1f\n", + psTrace ("psphot", 4, "%3d : %5.2f : %5.3f %5.3f : %5.3f %5.3f : %5.3f %5.3f : %5.3f %5.3f : %5.1f %5.1f %5.1f\n", i, refRadius->data.F32[nOut], binSB->data.F32[i], binSBstdev->data.F32[i], meanSB->data.F32[nOut], meanSBerr, petRatio->data.F32[nOut], petRatioErr->data.F32[nOut], - fluxSum->data.F32[nOut], sqrt(fluxSumErr2->data.F32[nOut]), areaSum->data.F32[nOut], areaInner); + fluxSum->data.F32[nOut], sqrt(fluxSumErr2->data.F32[nOut]), areaSum->data.F32[nOut], apixSum->data.F32[nOut], areaInner); // anytime we transition below the PETROSIAN_RATIO, calculate the radius and flux @@ -104,7 +115,9 @@ if (i == 0) { // assume Fmax @ R = 0.0 - petRadius = InterpolateValues (1.0, 0.0, petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO); - } else { - petRadius = InterpolateValues (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO); + petRadius = InterpolateValues (1.0, 0.0, petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO); + petRadiusErr = InterpolateValuesErrX (1.0, 0.0, petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO, 0.0, petRatioErr->data.F32[nOut]); + } else { + petRadius = InterpolateValues (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO); + petRadiusErr = InterpolateValuesErrX (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO, petRatioErr->data.F32[nOut-1], petRatioErr->data.F32[nOut]); } above = false; @@ -128,11 +141,17 @@ } + // if we failed to reach the PETROSIAN_RATIO, use the lowest significant ratio instead (flag this!) if (!anyPetro) { // interpolate Rvec between i-1 and i to PETROSIAN_RATIO to get flux (Fvec) and radius (rvec) if (lowestSignificantRadius == 0) { // assume Fmax @ R = 0.0 - petRadius = InterpolateValues (1.0, 0.0, petRatio->data.F32[lowestSignificantRadius], refRadius->data.F32[lowestSignificantRadius], PETROSIAN_RATIO); + petRadius = InterpolateValues (1.0, 0.0, petRatio->data.F32[lowestSignificantRadius], refRadius->data.F32[lowestSignificantRadius], PETROSIAN_RATIO); + petRadiusErr = InterpolateValuesErrX (1.0, 0.0, petRatio->data.F32[lowestSignificantRadius], refRadius->data.F32[lowestSignificantRadius], PETROSIAN_RATIO, 0.0, petRatioErr->data.F32[lowestSignificantRadius]); + } else { - petRadius = InterpolateValues (petRatio->data.F32[lowestSignificantRadius-1], refRadius->data.F32[lowestSignificantRadius-1], petRatio->data.F32[lowestSignificantRadius], refRadius->data.F32[lowestSignificantRadius], PETROSIAN_RATIO); + int n0 = lowestSignificantRadius-1; + int n1 = lowestSignificantRadius; + petRadius = InterpolateValues (petRatio->data.F32[n0], refRadius->data.F32[n0], petRatio->data.F32[n1], refRadius->data.F32[n1], PETROSIAN_RATIO); + petRadiusErr = InterpolateValuesErrX (petRatio->data.F32[n0], refRadius->data.F32[n0], petRatio->data.F32[n1], refRadius->data.F32[n1], PETROSIAN_RATIO, petRatioErr->data.F32[n0], petRatioErr->data.F32[n1]); } } @@ -147,5 +166,8 @@ continue; } else { - petFlux = InterpolateValues (refRadius->data.F32[i-1], fluxSum->data.F32[i-1], refRadius->data.F32[i], fluxSum->data.F32[i], apRadius); + petFlux = InterpolateValues (refRadius->data.F32[i-1], fluxSum->data.F32[i-1], refRadius->data.F32[i], fluxSum->data.F32[i], apRadius); + petFluxErr = InterpolateValuesErrY (refRadius->data.F32[i-1], fluxSum->data.F32[i-1], refRadius->data.F32[i], fluxSum->data.F32[i], apRadius, sqrt(fluxSumErr2->data.F32[i-1]), sqrt(fluxSumErr2->data.F32[i])); + petArea = InterpolateValues (refRadius->data.F32[i-1], areaSum->data.F32[i-1], refRadius->data.F32[i], areaSum->data.F32[i], apRadius); + petApix = InterpolateValues (refRadius->data.F32[i-1], apixSum->data.F32[i-1], refRadius->data.F32[i], apixSum->data.F32[i], apRadius); break; } @@ -158,4 +180,6 @@ float R50 = NAN; float R90 = NAN; + float R50err = NAN; + float R90err = NAN; bool found50 = false; bool found90 = false; @@ -167,5 +191,6 @@ continue; } else { - R50 = InterpolateValues (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux50); + R50 = InterpolateValues (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux50); + R50err = InterpolateValuesErrX (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux50, sqrt(fluxSumErr2->data.F32[i-1]), sqrt(fluxSumErr2->data.F32[i])); found50 = true; } @@ -176,5 +201,6 @@ continue; } else { - R90 = InterpolateValues (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux90); + R90 = InterpolateValues (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux90); + R90err = InterpolateValuesErrX (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux90, sqrt(fluxSumErr2->data.F32[i-1]), sqrt(fluxSumErr2->data.F32[i])); found90 = true; } @@ -188,10 +214,11 @@ source->extpars->petrosianR50 = R50; source->extpars->petrosianR90 = R90; + source->extpars->petrosianFill = petApix / petArea; // XXX add the errors - source->extpars->petrosianRadiusErr = NAN; - source->extpars->petrosianFluxErr = NAN; - source->extpars->petrosianR50Err = NAN; - source->extpars->petrosianR90Err = NAN; + source->extpars->petrosianRadiusErr = petRadiusErr; + source->extpars->petrosianFluxErr = petFluxErr; + source->extpars->petrosianR50Err = R50err; + source->extpars->petrosianR90Err = R90err; // fprintf (stderr, "source @ %f,%f\n", source->peak->xf, source->peak->yf); @@ -205,4 +232,5 @@ psFree(meanSB); psFree(areaSum); + psFree(apixSum); return true; @@ -210,6 +238,24 @@ float InterpolateValues (float X0, float Y0, float X1, float Y1, float X) { - float Y = Y0 + (Y1 - Y0) * (X - X0) / (X1 - X0); + float dydx = (Y1 - Y0) / (X1 - X0); + float Y = Y0 + dydx * (X - X0); return Y; } +float InterpolateValuesErrX (float X0, float Y0, float X1, float Y1, float X, float dX0, float dX1) { + + float dydx = (Y1 - Y0) / (X1 - X0); + float dxdx = (X - X0) / (X1 - X0); + + float dY = sqrt(PS_SQR(dX1*dydx*dxdx) + PS_SQR(dX0*dydx*(dxdx - 1.0))); + return dY; +} + +float InterpolateValuesErrY (float X0, float Y0, float X1, float Y1, float X, float dY0, float dY1) { + + float dxdx = (X - X0) / (X1 - X0); + + float dY = sqrt(PS_SQR(dY1*dxdx) + PS_SQR(dY0*(1.0 - dxdx))); + return dY; +} + Index: /trunk/psphot/src/psphotRadialApertures.c =================================================================== --- /trunk/psphot/src/psphotRadialApertures.c (revision 31153) +++ /trunk/psphot/src/psphotRadialApertures.c (revision 31154) @@ -37,4 +37,10 @@ int Nradial = 0; + // perform full non-linear fits / extended source analysis? + if (!psMetadataLookupBool (&status, recipe, "RADIAL_APERTURES")) { + psLogMsg ("psphot", PS_LOG_INFO, "skipping radial apertures\n"); + return true; + } + psTimerStart ("psphot.radial"); @@ -62,4 +68,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) @@ -67,7 +74,4 @@ assert (maskVal); - // XXX temporary user-supplied systematic sky noise measurement (derive from background model) - float skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE"); - // S/N limit to perform full non-linear fits float SN_LIM = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM"); @@ -75,5 +79,5 @@ // source analysis is done in S/N order (brightest first) // XXX are we getting the objects out of order? does it matter? - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); // option to limit analysis to a specific region @@ -95,5 +99,5 @@ // limit selection to some SN limit assert (source->peak); // how can a source not have a peak? - if (source->peak->SN < SN_LIM) continue; + if (sqrt(source->peak->detValue) < SN_LIM) continue; // limit selection by analysis region @@ -112,14 +116,18 @@ 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 ++; // 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); - - if (!psphotRadialApertureSource (source, recipe, skynoise, maskVal, radMax, 0)) { + pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, outerRadius + 2); + + if (!psphotRadialApertureSource (source, recipe, 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 { @@ -127,4 +135,9 @@ } + 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); @@ -135,5 +148,5 @@ } -bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *aperRadii, int entry) { +bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, psImageMaskType maskVal, const psVector *aperRadii, int entry) { // if we are a child source, save the results to the parent source radial aperture array @@ -163,7 +176,12 @@ } - // 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 + float xCM = NAN, yCM = NAN; + if (pmSourcePositionUseMoments(source)) { + xCM = source->moments->Mx - 0.5 - source->pixels->col0; // coord of peak in subimage + yCM = source->moments->My - 0.5 - source->pixels->row0; // coord of peak in subimage + } else { + xCM = source->peak->xf - 0.5 - source->pixels->col0; // coord of peak in subimage + yCM = source->peak->yf - 0.5 - source->pixels->row0; // coord of peak in subimage + } // one pass through the pixels to select the valid pixels and calculate R^2 @@ -203,7 +221,9 @@ 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 *fluxStd = 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 (fluxStd, 0.0); psVectorInit (fluxErr, 0.0); psVectorInit (fill, 0.0); @@ -212,22 +232,40 @@ for (int i = 0; i < pixRadius2->n; i++, rPix2++) { + int j = 0; float *aRad2 = aperRadii2->data.F32; - for (int j = 0; (*rPix2 < *aRad2) && (j < aperRadii2->n); j++, aRad2++) { + for (; (*aRad2 < *rPix2) && (j < aperRadii2->n); j++, aRad2++); + for (; j < aperRadii2->n; j++, aRad2++) { flux->data.F32[j] += pixFlux->data.F32[i]; + fluxStd->data.F32[j] += PS_SQR(pixFlux->data.F32[i]); fluxErr->data.F32[j] += pixVar->data.F32[i]; fill->data.F32[j] += 1.0; } } + + /* for each radial bin, R(i), we measure: + 1) the flux within that aperture: F(i) = \sum_{r_jn; 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); + + int nPix = fill->data.F32[i]; + float SBmean = flux->data.F32[i] / nPix; + float SBstdv = sqrt((fluxStd->data.F32[i] / nPix) - PS_SQR(SBmean)); + + // flux->data.F32[i] = SBmean * Area; + fluxStd->data.F32[i] = SBstdv * Area; + fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]) * Area / nPix; + + // fill->data.F32[i] /= Area; + psTrace ("psphot", 5, "radial bins: %3d %5.1f : %8.1f +/- %7.2f : %8.1f +/- %8.1f : %4.2f %6.1f\n", + i, aperRadii->data.F32[i], flux->data.F32[i], fluxErr->data.F32[i], SBmean, SBstdv, fill->data.F32[i], Area); } radialAper->flux = flux; + radialAper->fluxStdev = fluxStd; radialAper->fluxErr = fluxErr; radialAper->fill = fill; @@ -245,5 +283,5 @@ static int nPix = 0; -bool psphotRadialApertureSource_With_Sort (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *radMax, int entry) { +bool psphotRadialApertureSource_With_Sort (pmSource *source, psMetadata *recipe, psImageMaskType maskVal, const psVector *radMax, int entry) { psAssert(source->radialAper->data[entry] == NULL, "why is this already defined?"); Index: /trunk/psphot/src/psphotRadialAperturesByObject.c =================================================================== --- /trunk/psphot/src/psphotRadialAperturesByObject.c (revision 31153) +++ /trunk/psphot/src/psphotRadialAperturesByObject.c (revision 31154) @@ -43,7 +43,4 @@ assert (maskVal); - // XXX temporary user-supplied systematic sky noise measurement (derive from background model) - float skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE"); - // S/N limit to perform full non-linear fits float SN_LIM = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM"); @@ -63,5 +60,5 @@ // source analysis is done in S/N order (brightest first) - objects = psArraySort (objects, pmPhotObjSortBySN); + objects = psArraySort (objects, pmPhotObjSortByFlux); // generate look-up arrays for readouts @@ -109,5 +106,5 @@ // limit selection to some SN limit assert (source->peak); // how can a source not have a peak? - if (source->peak->SN < SN_LIM) continue; + if (sqrt(source->peak->detValue) < SN_LIM) continue; int index = source->imageID; @@ -147,14 +144,13 @@ // 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, outerRadius + 2); - if (!psphotRadialApertureSource (source, recipe, skynoise, maskVal, radMax, nMatchedPSF)) { + if (!psphotRadialApertureSource (source, recipe, 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 { source->mode |= PM_SOURCE_MODE_RADIAL_FLUX; + if (source->parent) { + source->parent->mode |= PM_SOURCE_MODE_RADIAL_FLUX; + } } Index: /trunk/psphot/src/psphotRadialProfile.c =================================================================== --- /trunk/psphot/src/psphotRadialProfile.c (revision 31153) +++ /trunk/psphot/src/psphotRadialProfile.c (revision 31154) @@ -14,5 +14,5 @@ float Rmax = 200; float fluxMin = 0.0; - float fluxMax = source->peak->flux; + float fluxMax = source->peak->rawFlux; bool RAW_RADIUS = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_RAW_RADIUS"); @@ -25,4 +25,5 @@ return false; } + // allocate: extpars->radFlux->radii,fluxes,theta // use the radial profiles to determine the radius of a given isophote. this isophote @@ -33,4 +34,6 @@ return false; } + // allocate : extpars->radFlux->isophotalRadii (use profile->radii,fluxes) + // convert the isophotal radius vs angle measurements to an elliptical contour @@ -39,5 +42,6 @@ return false; } - + // use extpars->radFlux->isophotalRadii,theta (result in extpars->axes) + // generate a single, normalized radial profile following the elliptical contours. // the radius is normalized by the axis ratio so that on the major axis, 1 pixel = 1 pixel @@ -46,4 +50,5 @@ return false; } + // allocate extpars->ellipticalFlux->radiusElliptical,fluxElliptical (use axes to scale raw pixels) // generated profile in averaged bins @@ -52,4 +57,6 @@ return false; } + // allocate extpars->radProfile->binSB, binSBstdv, binSum, binFill, radialBins, area (small lengths) + // use radiusElliptical, fluxElliptical, return true; Index: /trunk/psphot/src/psphotReadout.c =================================================================== --- /trunk/psphot/src/psphotReadout.c (revision 31153) +++ /trunk/psphot/src/psphotReadout.c (revision 31154) @@ -9,4 +9,55 @@ } +// for now, let's store the detections on the readout->analysis for each readout +bool psphotDumpChisqs (pmConfig *config, const pmFPAview *view, const char *filerule) +{ + static int npass = 0; + char filename[64]; + + return true; + + bool status = true; + + int num = psphotFileruleCount(config, filerule); + + snprintf (filename, 64, "chisq.%02d.dat", npass); + FILE *f = fopen (filename, "w"); + + // 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?"); + + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + if (!source) continue; + + pmModel *model = pmSourceGetModel (NULL, source); + if (!model) continue; + + if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) { + fprintf (f, "%f %f %f %d %d %f 1 NONLINEAR\n", model->mag, model->params->data.F32[1], model->chisq, model->nDOF, model->nPix, model->chisqNorm); + } else { + fprintf (f, "%f %f %f %d %d %f 0 LINEAR\n", model->mag, model->params->data.F32[1], model->chisq, model->nDOF, model->nPix, model->chisqNorm); + } + } + } + fclose (f); + npass ++; + + return true; +} + bool psphotReadout(pmConfig *config, const pmFPAview *view, const char *filerule) { @@ -52,4 +103,17 @@ } +# if (0) + // XXX test to mask outliers, not very helpful + // mask the high values in the image (with MARK) + if (!psphotMaskBackground (config, view, filerule)) { + return psphotReadoutCleanup (config, view, filerule); + } + + // generate a background model (median, smoothed image) + if (!psphotModelBackground (config, view, filerule)) { + return psphotReadoutCleanup (config, view, filerule); + } +# endif + if (!psphotSubtractBackground (config, view, filerule)) { return psphotReadoutCleanup (config, view, filerule); @@ -83,8 +147,8 @@ // find blended neighbors of very saturated stars (detections->newSources) - if (!psphotDeblendSatstars (config, view, filerule)) { - psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis"); - return psphotReadoutCleanup (config, view, filerule); - } + // if (!psphotDeblendSatstars (config, view, filerule)) { + // psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis"); + // return psphotReadoutCleanup (config, view, filerule); + // } // mark blended peaks PS_SOURCE_BLEND (detections->newSources) @@ -134,4 +198,5 @@ // linear PSF fit to source peaks, subtract the models from the image (in PSF mask) psphotFitSourcesLinear (config, view, filerule, false); // pass 1 (detections->allSources) + psphotDumpChisqs (config, view, filerule); // identify CRs and extended sources (only unmeasured sources are measured) @@ -144,4 +209,5 @@ // replace model flux, adjust mask as needed, fit, subtract the models (full stamp) psphotBlendFit (config, view, filerule); // pass 1 (detections->allSources) + psphotDumpChisqs (config, view, filerule); // replace all sources @@ -151,4 +217,5 @@ // NOTE : apply to ALL sources (extended + psf) psphotFitSourcesLinear (config, view, filerule, true); // pass 2 (detections->allSources) + psphotDumpChisqs (config, view, filerule); // if we only do one pass, skip to extended source analysis @@ -157,41 +224,86 @@ // NOTE: possibly re-measure background model here with objects subtracted / or masked - // add noise for subtracted objects - psphotAddNoise (config, view, filerule); // pass 1 (detections->allSources) - - // find fainter sources - // NOTE: finds new peaks and new footprints, OLD and FULL set are saved on detections - psphotFindDetections (config, view, filerule, false); // pass 2 (detections->peaks, detections->footprints) - - // remove noise for subtracted objects (ie, return to normal noise level) - // NOTE: this needs to operate only on the OLD sources - psphotSubNoise (config, view, filerule); // pass 1 (detections->allSources) - - // define new sources based on only the new peaks - // NOTE: new sources are saved on detections->newSources - psphotSourceStats (config, view, filerule, false); // pass 2 (detections->newSources) - - // set source type - // NOTE: apply only to detections->newSources - if (!psphotRoughClass (config, view, filerule)) { // pass 2 (detections->newSources) - psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image"); - return psphotReadoutCleanup (config, view, filerule); - } - - // create full input models, set the radius to fitRadius, set circular fit mask - // NOTE: apply only to detections->newSources - psphotGuessModels (config, view, filerule); // pass 2 (detections->newSources) - - // replace all sources so fit below applies to all at once - // NOTE: apply only to OLD sources (which have been subtracted) - psphotReplaceAllSources (config, view, filerule); // pass 2 - - // merge the newly selected sources into the existing list - // NOTE: merge OLD and NEW - // XXX check on free of sources... - psphotMergeSources (config, view, filerule); // (detections->newSources + detections->allSources -> detections->allSources) - - // NOTE: apply to ALL sources - psphotFitSourcesLinear (config, view, filerule, true); // pass 3 (detections->allSources) + // NOTE: this block performs the 2nd pass low-significance PSF detection stage + { + // add noise for subtracted objects + psphotAddNoise (config, view, filerule); // pass 1 (detections->allSources) + + // find fainter sources + // NOTE: finds new peaks and new footprints, OLD and FULL set are saved on detections + psphotFindDetections (config, view, filerule, false); // pass 2 (detections->peaks, detections->footprints) + + // remove noise for subtracted objects (ie, return to normal noise level) + // NOTE: this needs to operate only on the OLD sources + psphotSubNoise (config, view, filerule); // pass 1 (detections->allSources) + + // define new sources based on only the new peaks + // NOTE: new sources are saved on detections->newSources + psphotSourceStats (config, view, filerule, false); // pass 2 (detections->newSources) + + // set source type + // NOTE: apply only to detections->newSources + if (!psphotRoughClass (config, view, filerule)) { // pass 2 (detections->newSources) + psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image"); + return psphotReadoutCleanup (config, view, filerule); + } + + // create full input models, set the radius to fitRadius, set circular fit mask + // NOTE: apply only to detections->newSources + psphotGuessModels (config, view, filerule); // pass 2 (detections->newSources) + + // replace all sources so fit below applies to all at once + // NOTE: apply only to OLD sources (which have been subtracted) + psphotReplaceAllSources (config, view, filerule); // pass 2 + + // merge the newly selected sources into the existing list + // NOTE: merge OLD and NEW + // XXX check on free of sources... + psphotMergeSources (config, view, filerule); // (detections->newSources + detections->allSources -> detections->allSources) + + // NOTE: apply to ALL sources + psphotFitSourcesLinear (config, view, filerule, true); // pass 3 (detections->allSources) + psphotDumpChisqs (config, view, filerule); + } + + // NOTE: this block performs the 2nd pass low-significance EXT detection stage (smooth or rebin by NxN times PSF size) + if (0) { + // add noise for subtracted objects + psphotAddNoise (config, view, filerule); // pass 1 (detections->allSources) + + // find fainter sources + // NOTE: finds new peaks and new footprints, OLD and FULL set are saved on detections + psphotFindDetections (config, view, filerule, false); // pass 2 (detections->peaks, detections->footprints) + + // remove noise for subtracted objects (ie, return to normal noise level) + // NOTE: this needs to operate only on the OLD sources + psphotSubNoise (config, view, filerule); // pass 1 (detections->allSources) + + // define new sources based on only the new peaks + // NOTE: new sources are saved on detections->newSources + psphotSourceStats (config, view, filerule, false); // pass 2 (detections->newSources) + + // set source type + // NOTE: apply only to detections->newSources + if (!psphotRoughClass (config, view, filerule)) { // pass 2 (detections->newSources) + psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image"); + return psphotReadoutCleanup (config, view, filerule); + } + + // create full input models, set the radius to fitRadius, set circular fit mask + // NOTE: apply only to detections->newSources + psphotGuessModels (config, view, filerule); // pass 2 (detections->newSources) + + // replace all sources so fit below applies to all at once + // NOTE: apply only to OLD sources (which have been subtracted) + psphotReplaceAllSources (config, view, filerule); // pass 2 + + // merge the newly selected sources into the existing list + // NOTE: merge OLD and NEW + // XXX check on free of sources... + psphotMergeSources (config, view, filerule); // (detections->newSources + detections->allSources -> detections->allSources) + + // NOTE: apply to ALL sources + psphotFitSourcesLinear (config, view, filerule, true); // pass 3 (detections->allSources) + } pass1finish: @@ -203,4 +315,5 @@ psphotExtendedSourceAnalysis (config, view, filerule); // pass 1 (detections->allSources) psphotExtendedSourceFits (config, view, filerule); // pass 1 (detections->allSources) + psphotRadialApertures(config, view, filerule); finish: Index: /trunk/psphot/src/psphotReplaceUnfit.c =================================================================== --- /trunk/psphot/src/psphotReplaceUnfit.c (revision 31153) +++ /trunk/psphot/src/psphotReplaceUnfit.c (revision 31154) @@ -155,4 +155,5 @@ // sources have not yet been subtracted in this image (but this flag may be raised) source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED; + if (!source->modelPSF) continue; float Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS]; @@ -231,4 +232,6 @@ bool isPSF = false; pmModel *model = pmSourceGetModel(&isPSF, source); + if (!model) continue; + float radius = model->fitRadius; // save for future use below @@ -236,5 +239,5 @@ if (isPSF || model->isPCM) { // the guess central intensity comes from the peak: - float Io = source->peak->flux; + float Io = source->peak->rawFlux; float Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS]; float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS]; Index: /trunk/psphot/src/psphotRoughClass.c =================================================================== --- /trunk/psphot/src/psphotRoughClass.c (revision 31153) +++ /trunk/psphot/src/psphotRoughClass.c (revision 31154) @@ -63,13 +63,39 @@ } + int NstarsInClump = psMetadataLookupS32 (&status, readout->analysis, "PSF_CLUMP_NSTARS"); + // if NstarsInClump is not defined, use the user-selected option: + if (!status) { + NstarsInClump = 1000; + } + + int ScaleForClump = 1; + if (NstarsInClump >= 12) ScaleForClump = 2; // 4 cells + if (NstarsInClump >= 27) ScaleForClump = 3; // 9 cells + if (NstarsInClump >= 48) ScaleForClump = 4; // 16 cells + if (NstarsInClump > 75) ScaleForClump = 5; // 25 cells + // we make this measurement on a NxM grid of regions across the readout int NX = psMetadataLookupS32 (&status, recipe, "PSF_CLUMP_NX"); CHECK_STATUS (status, "PSF_CLUMP_NX"); int NY = psMetadataLookupS32 (&status, recipe, "PSF_CLUMP_NY"); CHECK_STATUS (status, "PSF_CLUMP_NY"); - int dX = readout->image->numCols / NX; - int dY = readout->image->numRows / NY; + + int NXuse, NYuse; + + ScaleForClump = PS_MIN(ScaleForClump, PS_MAX(NX, NY)); + if (NX > NY) { + NXuse = ScaleForClump; + NYuse = (int) (ScaleForClump * (NY / NX) + 0.5); + } else { + NYuse = ScaleForClump; + NXuse = (int) (ScaleForClump * (NY / NX) + 0.5); + } + + psLogMsg ("psphot", 4, "With %d stars, using %d x %d grid for PSF clump\n", NstarsInClump, NXuse, NYuse); + + int dX = readout->image->numCols / NXuse; + int dY = readout->image->numRows / NYuse; int nRegion = 0; - for (int ix = 0; ix < NX; ix ++) { - for (int iy = 0; iy < NY; iy ++) { + for (int ix = 0; ix < NXuse; ix ++) { + for (int iy = 0; iy < NYuse; iy ++) { psRegion *region = psRegionAlloc (ix*dX, (ix + 1)*dX, iy*dY, (iy + 1)*dY); @@ -181,6 +207,5 @@ return false; } - psLogMsg ("psphot", 3, "psf clump X, Y: %f, %f\n", psfClump.X, psfClump.Y); - psLogMsg ("psphot", 3, "psf clump DX, DY: %f, %f\n", psfClump.dX, psfClump.dY); + psLogMsg ("psphot", 3, "psf clump X, Y: %f, %f : DX, DY: %f, %f : nStars %d of %d\n", psfClump.X, psfClump.Y, psfClump.dX, psfClump.dY, psfClump.nStars, psfClump.nTotal); // get basic parameters, or set defaults Index: /trunk/psphot/src/psphotSavePSFStars.c =================================================================== --- /trunk/psphot/src/psphotSavePSFStars.c (revision 31153) +++ /trunk/psphot/src/psphotSavePSFStars.c (revision 31154) @@ -17,5 +17,5 @@ // examine PSF sources in S/N order (brightest first) - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); // counters to track the size of the image and area used in a row Index: /trunk/psphot/src/psphotSetThreads.c =================================================================== --- /trunk/psphot/src/psphotSetThreads.c (revision 31153) +++ /trunk/psphot/src/psphotSetThreads.c (revision 31154) @@ -35,5 +35,5 @@ psFree(task); - task = psThreadTaskAlloc("PSPHOT_EXTENDED_FIT", 13); + task = psThreadTaskAlloc("PSPHOT_EXTENDED_FIT", 14); task->function = &psphotExtendedSourceFits_Threaded; psThreadTaskAdd(task); Index: /trunk/psphot/src/psphotSignificanceImage.c =================================================================== --- /trunk/psphot/src/psphotSignificanceImage.c (revision 31153) +++ /trunk/psphot/src/psphotSignificanceImage.c (revision 31154) @@ -4,5 +4,5 @@ // (S/N)^2. If FWMH_X,Y have been recorded, use them, otherwise use PEAKS_SMOOTH_SIGMA for the // smoothing kernel. -psImage *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, psImageMaskType maskVal) { +pmReadout *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, psImageMaskType maskVal) { float SIGMA_SMTH, NSIGMA_SMTH; @@ -21,4 +21,6 @@ minGauss = 0.5; } + + // NOTE: for a faint extended-source detection pass, we over-smooth by SOMETHING // if we have already determined the PSF model, then we have a better idea how to smooth this image @@ -92,13 +94,18 @@ psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SIGNIFICANCE_SCALE_FACTOR", PS_META_REPLACE, "Signicance scale factor", factor); - // XXX multithread this if needed + // we are going to return both the image and the weight here: the image contains the signal + // while the 'weight' will contain the significance (NOTE the deviation from the usual + // definition) + + // save the smoothed significance image in the weight array for (int j = 0; j < smooth_im->numRows; j++) { for (int i = 0; i < smooth_im->numCols; i++) { float value = smooth_im->data.F32[j][i]; if (value < 0 || smooth_wt->data.F32[j][i] <= 0 || (mask->data.PS_TYPE_IMAGE_MASK_DATA[j][i] & maskVal)) { - smooth_im->data.F32[j][i] = 0.0; + smooth_wt->data.F32[j][i] = 0.0; } else { - float v2 = value + PS_SQR(value/1000.0); - smooth_im->data.F32[j][i] = factor * PS_SQR(v2) / smooth_wt->data.F32[j][i]; + // XXX the value of 100 here (or 1000 before) must depend on the FWHM of the smoothing kernel, right?? + float v2 = value + PS_SQR(value/100.0); + smooth_wt->data.F32[j][i] = factor * PS_SQR(v2) / smooth_wt->data.F32[j][i]; } } @@ -112,13 +119,14 @@ static int pass = 0; sprintf (name, "snsmooth.v%d.fits", pass); - psphotSaveImage (NULL, smooth_im, name); + psphotSaveImage (NULL, smooth_wt, name); pass ++; } - - psFree(smooth_wt); - psImageConvolveSetThreads(oldThreads); - return smooth_im; + pmReadout *significanceRO = pmReadoutAlloc(NULL); + significanceRO->variance = smooth_wt; + significanceRO->image = smooth_im; + + return significanceRO; } Index: /trunk/psphot/src/psphotSourceFits.c =================================================================== --- /trunk/psphot/src/psphotSourceFits.c (revision 31153) +++ /trunk/psphot/src/psphotSourceFits.c (revision 31154) @@ -51,5 +51,5 @@ psArrayAdd (sourceSet, 16, source); - psTrace ("psphot", 4, "fitting blended source at %f %f : %f\n", source->peak->xf, source->peak->yf, source->peak->flux); + psTrace ("psphot", 4, "fitting blended source at %f %f : %f\n", source->peak->xf, source->peak->yf, source->peak->rawFlux); // we need to include all blends in the fit (unless primary is saturated?) @@ -69,5 +69,5 @@ // XXX assume local sky is 0.0? - model->params->data.F32[PM_PAR_I0] = blend->peak->flux; + model->params->data.F32[PM_PAR_I0] = blend->peak->rawFlux; model->params->data.F32[PM_PAR_XPOS] = blend->peak->xf; model->params->data.F32[PM_PAR_YPOS] = blend->peak->yf; @@ -83,5 +83,5 @@ psArrayAdd (sourceSet, 16, blend); - psTrace ("psphot", 5, "adding source at %f %f : %f\n", blend->peak->xf, blend->peak->yf, blend->peak->flux); + psTrace ("psphot", 5, "adding source at %f %f : %f\n", blend->peak->xf, blend->peak->yf, blend->peak->rawFlux); // free to avoid double counting model @@ -100,8 +100,4 @@ if (!isfinite(PSF->params->data.F32[PM_PAR_I0])) psAbort("nan in fit"); - - // correct model chisq for flux trend - double chiTrend = psPolynomial1DEval (psf->ChiTrend, PSF->params->data.F32[PM_PAR_I0]); - PSF->chisqNorm = PSF->chisq / chiTrend; // evaluate the blend objects, subtract if good, free otherwise @@ -111,8 +107,4 @@ if (!isfinite(model->params->data.F32[PM_PAR_I0])) psAbort("nan in fit"); - - // correct model chisq for flux trend - chiTrend = psPolynomial1DEval (psf->ChiTrend, model->params->data.F32[PM_PAR_I0]); - model->chisqNorm = model->chisq / chiTrend; // if this one failed, skip it @@ -159,5 +151,4 @@ bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) { - double chiTrend; pmSourceFitOptions options = *fitOptions; @@ -182,8 +173,4 @@ // clear the circular mask psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); - - // correct model chisq for flux trend - chiTrend = psPolynomial1DEval (psf->ChiTrend, PSF->params->data.F32[PM_PAR_I0]); - PSF->chisqNorm = PSF->chisq / chiTrend; // does the PSF model succeed? @@ -224,6 +211,4 @@ float radius; bool okEXT, okDBL; - float chiEXT, chiDBL; - double chiTrend; pmModel *ONE = NULL; pmSource *tmpSrc = NULL; @@ -258,4 +243,7 @@ maskVal |= markVal; + float chiEXT = NAN; + float chiDBL = NAN; + // this temporary source is used as a place-holder by the psphotEval functions below tmpSrc = pmSourceAlloc (); @@ -271,11 +259,8 @@ // correct first model chisqs for flux trend - chiDBL = NAN; ONE = DBL->data[0]; if (ONE) { if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit"); - chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]); - ONE->chisqNorm = ONE->chisq / chiTrend; - chiDBL = ONE->chisq / ONE->nDOF; // save chisq for double-star/galaxy comparison + chiDBL = ONE->chisqNorm; // save chisq for double-star/galaxy comparison ONE->fitRadius = radius; } @@ -285,6 +270,4 @@ if (ONE) { if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit"); - chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]); - ONE->chisqNorm = ONE->chisq / chiTrend; ONE->fitRadius = radius; } @@ -298,5 +281,5 @@ okEXT = psphotEvalEXT (tmpSrc, EXT); - chiEXT = EXT ? EXT->chisq / EXT->nDOF : NAN; + chiEXT = EXT ? EXT->chisqNorm : NAN; } Index: /trunk/psphot/src/psphotSourceMatch.c =================================================================== --- /trunk/psphot/src/psphotSourceMatch.c (revision 31153) +++ /trunk/psphot/src/psphotSourceMatch.c (revision 31154) @@ -86,10 +86,10 @@ if (!src) NEXT1; if (!src->peak) NEXT1; - if (!finite(src->peak->xf)) NEXT1; - if (!finite(src->peak->yf)) NEXT1; + if (!isfinite(src->peak->xf)) NEXT1; + if (!isfinite(src->peak->yf)) NEXT1; if (!obj) NEXT2; - if (!finite(obj->x)) NEXT2; - if (!finite(obj->y)) NEXT2; + if (!isfinite(obj->x)) NEXT2; + if (!isfinite(obj->y)) NEXT2; dx = src->peak->xf - obj->x; @@ -229,10 +229,7 @@ int col0 = readout->image->col0; - // XXX the peak type is not really used in psphot - // PM_PEAK_LONE is certainly not true, but irrelevant + // The peak type is not used in psphot. PM_PEAK_LONE may be wrong, but irrelevant float peakFlux = readout->image->data.F32[(int)(obj->y-row0-0.5)][(int)(obj->x-col0-0.5)]; pmPeak *peak = pmPeakAlloc(obj->x, obj->y, peakFlux, PM_PEAK_LONE); - peak->flux = peakFlux; - peak->SN = 1.0; peak->xf = obj->x; peak->yf = obj->y; Index: /trunk/psphot/src/psphotSourcePlots.c =================================================================== --- /trunk/psphot/src/psphotSourcePlots.c (revision 31153) +++ /trunk/psphot/src/psphotSourcePlots.c (revision 31154) @@ -18,5 +18,5 @@ // examine PSF sources in S/N order (brightest first) - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); // counters to track the size of the image and area used in a row Index: /trunk/psphot/src/psphotSourceSize.c =================================================================== --- /trunk/psphot/src/psphotSourceSize.c (revision 31153) +++ /trunk/psphot/src/psphotSourceSize.c (revision 31154) @@ -207,5 +207,5 @@ num++; - pmSourceMagnitudes (source, psf, photMode, maskVal, markVal); + pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, source->apRadius); float kMag = -2.5*log10(source->moments->KronFlux); @@ -327,5 +327,5 @@ // XXX can we test if psfMag is set and calculate only if needed? - pmSourceMagnitudes (source, psf, photMode, maskVal, markVal); + pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, source->apRadius); // convert to Mmaj, Mmin: @@ -501,5 +501,5 @@ // 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); + pmSourceMagnitudes (source, psf, photMode, maskVal, markVal, source->apRadius); // clear the mask bit @@ -1268,5 +1268,5 @@ // Soften variances (add systematic error) - float softening = options->soft * PS_SQR(source->peak->flux); // Softening for variances + float softening = options->soft * PS_SQR(source->peak->rawFlux); // Softening for variances // Across the middle: y = 0 Index: /trunk/psphot/src/psphotSourceStats.c =================================================================== --- /trunk/psphot/src/psphotSourceStats.c (revision 31153) +++ /trunk/psphot/src/psphotSourceStats.c (revision 31154) @@ -122,4 +122,8 @@ // allocate space for moments source->moments = pmMomentsAlloc(); + + if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) { + fprintf (stderr, "moment failure\n"); + } // allocate image, weight, mask arrays for each peak (square of radius OUTER) @@ -361,5 +365,4 @@ bool status = false; - float BIG_RADIUS; psScalar *scalar = NULL; @@ -373,4 +376,9 @@ psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA); psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA); + + // if no valid pixels, massive swing or very large Mrf, likely saturated source, + // try a much larger box + float BIG_RADIUS = 3.0*RADIUS; + float BIG_SIGMA = 3.0*SIGMA; // maskVal is used to test for rejected pixels, and must include markVal @@ -390,5 +398,5 @@ // skip faint sources for moments measurement - if (source->peak->SN < MIN_SN) { + if (sqrt(source->peak->detValue) < MIN_SN) { source->mode |= PM_SOURCE_MODE_BELOW_MOMENTS_SN; Nfaint++; @@ -416,39 +424,40 @@ } - // measure basic source moments (no S/N clipping on input pixels) - status = pmSourceMoments (source, RADIUS, SIGMA, 0.0, maskVal); - // XXX moments / aperture test: - if (0) { - // clear the mask bit and set the circular mask pixels - psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, 4.0*source->moments->Mrf, "OR", markVal); - - float apMag = NAN; - pmSourcePhotometryAper (&apMag, NULL, NULL, NULL, source->pixels, source->variance, source->maskObj, maskVal); - fprintf (stderr, "apMag: %f, kronMag: %f\n", apMag, -2.5*log10(source->moments->KronFlux)); - - // clear the mask bit - psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); - } - if (status) { - Nmoments ++; - continue; - } - - // if no valid pixels, or massive swing, likely saturated source, - // try a much larger box - BIG_RADIUS = PS_MIN (INNER, 3*RADIUS); - psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y); - status = pmSourceMoments (source, BIG_RADIUS, 3.0*SIGMA, 0.0, maskVal); - if (status) { - source->mode |= PM_SOURCE_MODE_BIG_RADIUS; - Nmoments ++; - continue; - } - - source->mode |= PM_SOURCE_MODE_MOMENTS_FAILURE; - Nfail ++; - psErrorClear(); - continue; + // XXX how can this be set: this is raised in pmSourceRoughClass, called after this function? + if (source->mode & PM_SOURCE_MODE_SATSTAR) { + fprintf (stderr, "satstar: %f,%f\n", source->peak->xf, source->peak->yf); + } + + if (!(source->peak->type == PM_PEAK_SUSPECT_SATURATION)) { + // measure basic source moments (no S/N clipping on input pixels) + status = pmSourceMoments (source, RADIUS, SIGMA, 0.0, maskVal); + } else { + // For saturated stars, choose a much larger box NOTE this is slightly sleazy, but + // only slightly: pmSourceRedefinePixels uses the readout to pass the pointers to + // the parent image data. I guess the API could be simplified: we could recover + // this from the source in the function + + pmReadout tmpReadout; + tmpReadout.image = (psImage *)source->pixels->parent; + tmpReadout.mask = (psImage *)source->maskView->parent; + tmpReadout.variance = (psImage *)source->variance->parent; + + // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS + pmSourceRedefinePixels (source, &tmpReadout, source->peak->x, source->peak->y, BIG_RADIUS + 2); + + psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y); + status = pmSourceMoments (source, BIG_RADIUS, BIG_SIGMA, 0.0, maskVal); + source->mode |= PM_SOURCE_MODE_BIG_RADIUS; + } + if (!status) { + source->mode |= PM_SOURCE_MODE_MOMENTS_FAILURE; + Nfail ++; + psErrorClear(); + continue; + } + Nmoments ++; + + // re-try big sources or not?? + // if (source->moments->Mrf < 2.0*SIGMA) } @@ -473,6 +482,9 @@ float MIN_SN = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN"); psAssert (status, "missing MOMENTS_SN_MIN"); - float PSF_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); psAssert (status, "missing PSF_SN_LIM"); psF32 MOMENTS_AR_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_AR_MAX"); psAssert (status, "missing MOMENTS_AR_MAX"); + + // when we set the window, we are not attempting to measure spatial variations; we can use a somewhat higher S/N limit + // since we are using all sources (true?) + float PSF_SN_LIM = 2.0*psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); psAssert (status, "missing PSF_SN_LIM"); // XXX move this to a config file? @@ -480,7 +492,8 @@ float sigma[NSIGMA] = {1.0, 2.0, 3.0, 4.5, 6.0, 9.0, 12.0, 18.0}; float Sout[NSIGMA]; - - // this sorts by peak->SN - sources = psArraySort (sources, pmSourceSortBySN); + int Nout[NSIGMA]; // number of stars found in clump : use this to control the number of regions measured by psphotRoughClass + + // this sorts by peak->rawFlux + sources = psArraySort (sources, pmSourceSortByFlux); // loop over radii: @@ -495,5 +508,5 @@ // skip faint sources for moments measurement - if (source->peak->SN < MIN_SN) { + if (sqrt(source->peak->detValue) < MIN_SN) { continue; } @@ -510,5 +523,5 @@ // determine the PSF parameters from the source moment values pmPSFClump psfClump = pmSourcePSFClump (NULL, NULL, sources, PSF_SN_LIM, PSF_CLUMP_GRID_SCALE, MOMENTS_SX_MAX, MOMENTS_SY_MAX, MOMENTS_AR_MAX); - psLogMsg ("psphot", 3, "radius %.1f, nStars: %d, nSigma: %5.2f, X, Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]); + psLogMsg ("psphot", 3, "radius %.1f, nStars: %d of %d in clump, nSigma: %5.2f, X, Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nTotal, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]); #if 0 @@ -531,7 +544,9 @@ Sout[i] = sqrt(0.5*(psfClump.X + psfClump.Y)) / sigma[i]; + Nout[i] = psfClump.nStars; } // we are looking for sigma for which Sout = 0.65 (or some other value) + int Nstars = 0; float Sigma = NAN; float minS = Sout[0]; @@ -549,4 +564,5 @@ if ((Sout[i] < 0.65) && (Sout[i+1] < 0.65)) continue; Sigma = sigma[i] + (0.65 - Sout[i])*(sigma[i+1] - sigma[i])/(Sout[i+1] - Sout[i]); + Nstars = 0.5*(Nout[i] + Nout[i+1]); } psAssert (isfinite(Sigma), "did we miss a case?"); @@ -558,4 +574,5 @@ psMetadataAddF32(analysis, PS_LIST_TAIL, "MOMENTS_GAUSS_SIGMA", PS_META_REPLACE, "moments limit", Sigma); psMetadataAddF32(analysis, PS_LIST_TAIL, "PSF_MOMENTS_RADIUS", PS_META_REPLACE, "moments limit", 4.0*Sigma); + psMetadataAddF32(analysis, PS_LIST_TAIL, "PSF_CLUMP_NSTARS", PS_META_REPLACE, "number of stars in clump", Nstars); psLogMsg ("psphot", 3, "using window function with sigma = %f\n", Sigma); Index: /trunk/psphot/src/psphotStack.c =================================================================== --- /trunk/psphot/src/psphotStack.c (revision 31153) +++ /trunk/psphot/src/psphotStack.c (revision 31154) @@ -2,4 +2,7 @@ int main (int argc, char **argv) { + + // uncomment to turn on memory dumps (move this to an option) + // psMemDumpSetState(true); psTimerStart ("complete"); @@ -12,4 +15,6 @@ psphotVersionPrint(); + + psMemDump("start"); // load input data (config and images (signal, noise, mask) Index: /trunk/psphot/src/psphotStackArguments.c =================================================================== --- /trunk/psphot/src/psphotStackArguments.c (revision 31153) +++ /trunk/psphot/src/psphotStackArguments.c (revision 31154) @@ -22,6 +22,9 @@ } - // -version and -dumpconfig arguments - PSARGUMENTS_INSTANTIATE_GENERICS( psphot, config, argc, argv ); + // generic arguments (version, dumpconfig) + PS_ARGUMENTS_GENERIC( psphot, config, argc, argv ); + + // thread arguments + PS_ARGUMENTS_THREADS( psphot, config, argc, argv ) // save the following additional recipe values based on command-line options @@ -29,11 +32,14 @@ psMetadata *options = pmConfigRecipeOptions (config, PSPHOT_RECIPE); - // Number of threads is handled - PSARGUMENTS_INSTANTIATE_THREADSARG( psphot, config, argc, argv ) - // visual : interactive display mode if ((N = psArgumentGet (argc, argv, "-visual"))) { psArgumentRemove (N, &argc, argv); pmVisualSetVisual(true); + } + + // memdump : enable memory spot checks + if ((N = psArgumentGet (argc, argv, "-memdump"))) { + psArgumentRemove (N, &argc, argv); + psMemDumpSetState(true); } Index: /trunk/psphot/src/psphotStackChisqImage.c =================================================================== --- /trunk/psphot/src/psphotStackChisqImage.c (revision 31153) +++ /trunk/psphot/src/psphotStackChisqImage.c (revision 31154) @@ -101,7 +101,9 @@ for (int ix = 0; ix < inImage->numCols; ix++) { if (inMask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) continue; - chiImage->data.F32[iy][ix] += PS_SQR(inImage->data.F32[iy][ix]) / inVariance->data.F32[iy][ix]; + // XXX TEST chiImage->data.F32[iy][ix] += PS_SQR(inImage->data.F32[iy][ix]) / inVariance->data.F32[iy][ix]; + chiImage->data.F32[iy][ix] = 0.0; chiVariance->data.F32[iy][ix] = 1.0; // ?? what is the right value? just init to this? - chiMask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] = 0x00; // we have valid data so unmask this pixel + // XXX TEST chiMask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] = 0x00; // we have valid data so unmask this pixel + chiMask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] = 0x01; // we have valid data so unmask this pixel } } Index: /trunk/psphot/src/psphotStackImageLoop.c =================================================================== --- /trunk/psphot/src/psphotStackImageLoop.c (revision 31153) +++ /trunk/psphot/src/psphotStackImageLoop.c (revision 31154) @@ -1,4 +1,3 @@ # include "psphotStandAlone.h" -#define WCS_NONLIN_TOL 0.001 // Non-linear tolerance for header WCS # define ESCAPE(MESSAGE) { \ @@ -8,7 +7,7 @@ } -bool GetAstrometryFPA (pmConfig *config, pmFPAview *view); -bool SetAstrometryFPA (pmConfig *config, pmFPAview *view, bool bilevelAstrometry); -bool GetAstrometryChip (pmConfig *config, pmFPAview *view, bool bilevelAstrometry); +// XXX this implementation is not smart about multi-level astrometry headers +bool UpdateHeadersForFPA (pmConfig *config, pmFPAview *view); +bool UpdateHeadersForChip (pmConfig *config, pmFPAview *view); bool psphotStackImageLoop (pmConfig *config) { @@ -18,4 +17,6 @@ pmCell *cell; pmReadout *readout; + + psMemDump("startloop"); pmFPAfile *inputRaw = psMetadataLookupPtr (&status, config->files, "PSPHOT.STACK.INPUT.RAW"); @@ -32,6 +33,4 @@ // XXX for now, just load the full set of images up front if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa in psphot."); - - bool bilevelAstrometry = GetAstrometryFPA(config, view); // for psphot, we force data to be read at the chip level @@ -51,4 +50,6 @@ if (! readout->data_exists) { continue; } + psMemDump("load"); + // PSF matching if (!psphotStackMatchPSFs (config, view)) { @@ -57,4 +58,5 @@ return false; } + psMemDump("stackmatch"); // XXX for now, we assume there is only a single chip in the PHU: @@ -64,5 +66,5 @@ return false; } - + psMemDump("psphot"); } // drop all versions of the internal files @@ -78,12 +80,13 @@ } - GetAstrometryChip(config, view, bilevelAstrometry); + UpdateHeadersForChip(config, view); // save output which is saved at the chip level if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed output for Chip in psphot."); } + psMemDump("doneloop"); - SetAstrometryFPA(config, view, bilevelAstrometry); - + UpdateHeadersForFPA(config, view); + // save output which is saved at the fpa level if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed ouput for FPA in psphot."); @@ -93,18 +96,10 @@ psFree (view); + + psMemDump("doneoutput"); return true; } -/* - the easiest way to implement this is to assume we can pre-load the full set of images up front. - with 5 filters and 6000^2 (image, mask, var = 10 byte per pixel), we need 1.8GB, which is not too bad. -*/ - -# define UPDATE_HEADER 1 - -bool GetAstrometryFPA (pmConfig *config, pmFPAview *view) { - - bool foundAstrometry = false; - bool bilevelAstrometry = false; +bool UpdateHeadersForChip (pmConfig *config, pmFPAview *view) { int num = psphotFileruleCount(config, "PSPHOT.INPUT"); @@ -122,36 +117,25 @@ psAssert (input, "missing input file"); - // find the FPA phu - pmHDU *phu = pmFPAviewThisPHU(view, input->fpa); - if (!phu) { - psWarning("Unable to read bilevel mosaic astrometry for input FPA entry %d", i); - continue; + // just copy the input headers to the output headers, then update version info + pmChip *inChip = pmFPAviewThisChip(view, input->fpa); ///< Chip in the input + pmHDU *inHDU = pmFPAviewThisHDU (view, input->fpa); + + pmChip *outChip = pmFPAviewThisChip(view, output->fpa); ///< Chip in the output + pmHDU *outHDU = pmFPAviewThisHDU (view, output->fpa); + if (!outHDU) { + pmFPAAddSourceFromView(output->fpa, view, output->format); + outHDU = pmFPAviewThisHDU (view, output->fpa); + psAssert (outHDU, "failed to make HDU"); } - - char *ctype = psMetadataLookupStr(NULL, phu->header, "CTYPE1"); - if (!ctype) { - psWarning("Error in WCS keywords for input FPA entry %d", i); - continue; + if (!outHDU->header) { + outHDU->header = psMetadataCopy(NULL, inHDU->header); } - - if (!foundAstrometry) { - bilevelAstrometry = !strcmp (&ctype[4], "-DIS"); - } else { - if (bilevelAstrometry != !strcmp (&ctype[4], "-DIS")) { - psAbort("astrometry type mis-match"); - } - } - - if (bilevelAstrometry) { - // update the output structures - if (!pmAstromReadBilevelMosaic(output->fpa, phu->header)) { - psWarning("Unable to read bilevel mosaic astrometry for input FPA."); - } - } + outChip->toFPA = psMemIncrRefCounter(inChip->toFPA); + outChip->fromFPA = psMemIncrRefCounter(inChip->fromFPA); } - return bilevelAstrometry; + return true; } -bool GetAstrometryChip (pmConfig *config, pmFPAview *view, bool bilevelAstrometry) { +bool UpdateHeadersForFPA (pmConfig *config, pmFPAview *view) { int num = psphotFileruleCount(config, "PSPHOT.INPUT"); @@ -169,72 +153,23 @@ psAssert (input, "missing input file"); - // Need to update the output for astrometry. Read WCS data from the corresponding - // header and save in the output fpa - pmHDU *inHDU = pmFPAviewThisHDU (view, input->fpa); - pmChip *outChip = pmFPAviewThisChip(view, output->fpa); ///< Chip in the output + output->fpa->toTPA = psMemIncrRefCounter(input->fpa->toTPA); + output->fpa->fromTPA = psMemIncrRefCounter(input->fpa->fromTPA); + output->fpa->toSky = psMemIncrRefCounter(input->fpa->toSky); - pmHDU *outHDU = pmFPAviewThisHDU (view, output->fpa); - if (!outHDU) { - pmFPAAddSourceFromView(output->fpa, view, output->format); - outHDU = pmFPAviewThisHDU (view, output->fpa); - psAssert (outHDU, "failed to make HDU"); - } - if (!outHDU->header) { - outHDU->header = psMetadataAlloc(); - } + pmConceptsCopyFPA(output->fpa, input->fpa, true, true); - if (bilevelAstrometry) { - if (!pmAstromReadBilevelChip (outChip, inHDU->header)) { - psWarning("Unable to read bilevel chip astrometry for input FPA."); - continue; - } - if (!pmAstromWriteBilevelChip(outHDU->header, outChip, WCS_NONLIN_TOL)) { - psWarning("Unable to generate WCS header."); - continue; - } - } else { - // we use a default FPA pixel scale of 1.0 - if (!pmAstromReadWCS (output->fpa, outChip, inHDU->header, 1.0)) { - psWarning("Unable to read WCS astrometry for input FPA."); - continue; - } - if (!pmAstromWriteWCS(outHDU->header, output->fpa, outChip, WCS_NONLIN_TOL)) { - psWarning("Unable to generate WCS header."); - continue; - } - } + // XXX TEST + // pmFPASetFileStatus(output->fpa, true); + // pmFPASetDataStatus(output->fpa, true); + // pmChip *chip = output->fpa->chips->data[0]; + // pmCell *cell = chip->cells->data[0]; + // pmReadout *readout = cell->readouts->data[0]; } - return true; } -bool SetAstrometryFPA (pmConfig *config, pmFPAview *view, bool bilevelAstrometry) { +/* + the easiest way to implement this is to assume we can pre-load the full set of images up front. + with 5 filters and 6000^2 (image, mask, var = 10 byte per pixel), we need 1.8GB, which is not too bad. +*/ - if (!bilevelAstrometry) return true; - - int num = psphotFileruleCount(config, "PSPHOT.INPUT"); - - // loop over the available readouts - for (int i = 0; i < num; i++) { - - // find the currently selected readout - pmFPAfile *output = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.OUTPUT.IMAGE", i); // File of interest - psAssert (output, "missing file?"); - - pmHDU *PHU = pmFPAviewThisPHU(view, output->fpa); - if (!PHU) { - 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."); - } - } - - return true; -} Index: /trunk/psphot/src/psphotStackMatchPSFs.c =================================================================== --- /trunk/psphot/src/psphotStackMatchPSFs.c (revision 31153) +++ /trunk/psphot/src/psphotStackMatchPSFs.c (revision 31154) @@ -92,4 +92,9 @@ } + // copy the header data from Src to Out + // pmHDU *hduSrc = pmHDUFromReadout(readoutSrc); + // psAssert (hduSrc, "input missing hdu?"); + + // set NAN pixels to 'SAT' psImageMaskType maskSat = pmConfigMaskGet("SAT", config); Index: /trunk/psphot/src/psphotStackMatchPSFsUtils.c =================================================================== --- /trunk/psphot/src/psphotStackMatchPSFsUtils.c (revision 31153) +++ /trunk/psphot/src/psphotStackMatchPSFsUtils.c (revision 31154) @@ -64,5 +64,5 @@ } if (!source->peak) continue; - if (source->peak->SN < SN_MIN) continue; + if (sqrt(source->peak->detValue) < SN_MIN) continue; coordsFromSource(&x->data.F32[numGood], &y->data.F32[numGood], source); numGood++; @@ -81,5 +81,5 @@ } if (!source->peak) continue; - if (source->peak->SN < SN_MIN) continue; + if (sqrt(source->peak->detValue) < SN_MIN) continue; float xSource, ySource; // Coordinates of source coordsFromSource(&xSource, &ySource, source); @@ -287,5 +287,5 @@ float penalty = psMetadataLookupF32(NULL, subRecipe, "PENALTY"); // Penalty for wideness - int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads + int threads = psMetadataLookupS32(NULL, config->arguments, "NTHREADS"); // Number of threads int order = psMetadataLookupS32(NULL, subRecipe, "SPATIAL.ORDER"); // Spatial polynomial order @@ -363,5 +363,5 @@ // we need to register the FWHM values for use downstream - pmSubtractionSetFWHMs(options->inputSeeing->data.F32[index], options->targetSeeing); + pmSubtractionSetFWHMs(options->targetSeeing, options->inputSeeing->data.F32[index]); pmSubtractionParamScaleOptions(scale, scaleRef, scaleMin, scaleMax); Index: /trunk/psphot/src/psphotStackPSF.c =================================================================== --- /trunk/psphot/src/psphotStackPSF.c (revision 31153) +++ /trunk/psphot/src/psphotStackPSF.c (revision 31154) @@ -56,8 +56,34 @@ } - float Sxx = sqrt(2.0)*targetFWHM / 2.35; + // measured scale factors (fwhm = Sxx * 2.35 * scaleFactor / sqrt(2.0)) + // GAUSS : 1.000 + // PGAUSS : 1.006 + // QGAUSS : 1.151 + // RGAUSS : 0.883 + // PS1_V1 : 1.134 + + float scaleFactor = NAN; + if (!strcmp(psfModel, "PS_MODEL_GAUSS")) { + scaleFactor = 1.000; + } + if (!strcmp(psfModel, "PS_MODEL_PGAUSS")) { + scaleFactor = 1.0006; + } + if (!strcmp(psfModel, "PS_MODEL_QGAUSS")) { + scaleFactor = 1.151; + } + if (!strcmp(psfModel, "PS_MODEL_RGAUSS")) { + scaleFactor = 0.883; + } + if (!strcmp(psfModel, "PS_MODEL_PS1_V1")) { + scaleFactor = 1.134; + } + psAssert (isfinite(scaleFactor), "invalid model for PSF"); + + float Sxx = sqrt(2.0)*targetFWHM / 2.35 / scaleFactor; // XXX probably should make the model type (and par 7) optional from recipe - psf = pmPSFBuildSimple(psfModel, Sxx, Sxx, 0.0, 1.0); + // psf = pmPSFBuildSimple(psfModel, Sxx, Sxx, 0.0, 1.0); + psf = pmPSFBuildSimple(psfModel, Sxx, Sxx, 0.0, 0.2); if (!psf) { psError(PSPHOT_ERR_PSF, false, "Unable to build dummy PSF."); Index: /trunk/psphot/src/psphotStackParseCamera.c =================================================================== --- /trunk/psphot/src/psphotStackParseCamera.c (revision 31153) +++ /trunk/psphot/src/psphotStackParseCamera.c (revision 31154) @@ -172,4 +172,23 @@ } + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + bool saveCnv = psMetadataLookupBool(&status, recipe, "SAVE.CNV"); + bool saveChisq = psMetadataLookupBool(&status, recipe, "SAVE.CHISQ"); + + // loop over the available readouts + for (int i = 0; i < nInputs; i++) { + pmFPAfile *file = NULL; + + file = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.OUTPUT.IMAGE", i); + file->save = saveCnv; + + file = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.OUTPUT.MASK", i); + file->save = saveCnv; + + file = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.OUTPUT.VARIANCE", i); + file->save = saveCnv; + } + // generate an pmFPAimage for the chisqImage // XXX output of these files should be optional @@ -180,5 +199,5 @@ return false; } - chisqImage->save = true; + chisqImage->save = saveChisq; pmFPAfile *chisqMask = pmFPAfileDefineOutput(config, chisqImage->fpa, "PSPHOT.CHISQ.MASK"); @@ -191,5 +210,5 @@ return NULL; } - chisqMask->save = true; + chisqMask->save = saveChisq; pmFPAfile *chisqVariance = pmFPAfileDefineOutput(config, chisqImage->fpa, "PSPHOT.CHISQ.VARIANCE"); @@ -202,5 +221,5 @@ return NULL; } - chisqVariance->save = true; + chisqVariance->save = saveChisq; } Index: /trunk/psphot/src/psphotStackReadout.c =================================================================== --- /trunk/psphot/src/psphotStackReadout.c (revision 31153) +++ /trunk/psphot/src/psphotStackReadout.c (revision 31154) @@ -127,4 +127,6 @@ } + psMemDump("sourcestats"); + // generate the objects (object unify the sources from the different images) // XXX this could just match the detections for the chisq image, and not bother measuring the @@ -132,4 +134,6 @@ psArray *objects = psphotMatchSources (config, view, STACK_SRC); + psMemDump("matchsources"); + if (!strcasecmp (breakPt, "TEST2")) { psFree(objects); @@ -143,4 +147,6 @@ return psphotReadoutCleanup (config, view, STACK_SRC); } + + psMemDump("sourcestats"); // find blended neighbors of very saturated stars (detections->newSources) @@ -211,4 +217,6 @@ } + psMemDump("psfstats"); + psphotStackObjectsUnifyPosition (objects); @@ -230,4 +238,6 @@ } + psMemDump("extmeas"); + bool smoothAgain = true; for (int nMatchedPSF = 0; smoothAgain; nMatchedPSF++) { @@ -250,4 +260,5 @@ // smooth to the next FWHM, or set 'smoothAgain' to false if no more psphotStackMatchPSFsNext(&smoothAgain, config, view, STACK_OUT, nMatchedPSF); + psMemDump("matched"); } @@ -264,5 +275,5 @@ // drop the references to the image pixels held by each source - // psphotSourceFreePixels (config, view, STACK_OUT); + psphotSourceFreePixels (config, view, STACK_OUT); psphotSourceFreePixels (config, view, STACK_SRC); Index: /trunk/psphot/src/psphotStandAlone.h =================================================================== --- /trunk/psphot/src/psphotStandAlone.h (revision 31153) +++ /trunk/psphot/src/psphotStandAlone.h (revision 31154) @@ -14,5 +14,4 @@ pmConfig *psphotArguments (int argc, char **argv); bool psphotParseCamera (pmConfig *config); -bool psphotImageLoop (pmConfig *config); bool psphotMosaicChip(pmConfig *config, const pmFPAview *view, char *outFile, char *inFile); void psphotCleanup (pmConfig *config); Index: /trunk/psphot/src/psphotSubtractBackground.c =================================================================== --- /trunk/psphot/src/psphotSubtractBackground.c (revision 31153) +++ /trunk/psphot/src/psphotSubtractBackground.c (revision 31154) @@ -24,10 +24,5 @@ assert (modelFile); - pmReadout *model = NULL; - if (modelFile->mode == PM_FPA_MODE_INTERNAL) { - model = modelFile->readout; - } else { - model = pmFPAviewThisReadout (view, modelFile->fpa); - } + pmReadout *model = READOUT_OR_INTERNAL(view, modelFile); assert (model); @@ -44,9 +39,5 @@ if (file) { // we are using PSPHOT.BACKGND as an I/O file: select readout or create - if (file->mode == PM_FPA_MODE_INTERNAL) { - background = file->readout; - } else { - background = pmFPAviewThisReadout (view, file->fpa); - } + background = READOUT_OR_INTERNAL(view, file); if (background == NULL) { // readout does not yet exist: create from input Index: /trunk/psphot/src/psphotTestPSF.c =================================================================== --- /trunk/psphot/src/psphotTestPSF.c (revision 31153) +++ /trunk/psphot/src/psphotTestPSF.c (revision 31154) @@ -17,5 +17,5 @@ // examine PSF sources in S/N order (brightest first) - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); // array to store candidate PSF stars Index: /trunk/psphot/src/psphotVisual.c =================================================================== --- /trunk/psphot/src/psphotVisual.c (revision 31153) +++ /trunk/psphot/src/psphotVisual.c (revision 31154) @@ -39,29 +39,11 @@ switch (channel) { case 1: - if (kapa1 == -1) { - kapa1 = KapaOpenNamedSocket ("kapa", "psphot:images"); - if (kapa1 == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - pmVisualSetVisual(false); - } - } - return kapa1; + pmVisualInitWindow (&kapa1, "psphot:images"); + return kapa1; case 2: - if (kapa2 == -1) { - kapa2 = KapaOpenNamedSocket ("kapa", "psphot:plots"); - if (kapa2 == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - pmVisualSetVisual(false); - } - } + pmVisualInitWindow (&kapa2, "psphot:plots"); return kapa2; case 3: - if (kapa3 == -1) { - kapa3 = KapaOpenNamedSocket ("kapa", "psphot:stamps"); - if (kapa3 == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - pmVisualSetVisual(false); - } - } + pmVisualInitWindow (&kapa3, "psphot:stamps"); return kapa3; default: @@ -74,4 +56,6 @@ int myKapa = psphotKapaChannel (channel); + if (myKapa == -1) return false; + if (!(strcasecmp (overlay, "all"))) { KiiEraseOverlay (myKapa, "red"); @@ -218,5 +202,5 @@ } -bool psphotVisualScaleImage (int kapaFD, psImage *inImage, psImage *inMask, const char *name, int channel) { +bool psphotVisualScaleImage (int kapaFD, psImage *inImage, psImage *inMask, const char *name, float factor, int channel) { KiiImage image; @@ -239,6 +223,8 @@ strcpy (data.name, name); strcpy (data.file, name); - data.zero = stats->robustMedian - stats->robustStdev; - data.range = 5*stats->robustStdev; + data.zero = stats->robustMedian - factor*stats->robustStdev; + + // XXX I we have a smoothed image, this make a much-too-tight display range + data.range = 5*factor*stats->robustStdev; data.logflux = 0; @@ -283,7 +269,12 @@ if (kapa == -1) return false; + float factor = 1.0; + if (readout->covariance) { + factor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); + } + psphotVisualShowMask (kapa, readout->mask, "mask", 2); - psphotVisualScaleImage (kapa, readout->variance, readout->mask, "variance", 1); - psphotVisualScaleImage (kapa, readout->image, readout->mask, "image", 0); + psphotVisualScaleImage (kapa, readout->variance, readout->mask, "variance", 1.0, 1); + psphotVisualScaleImage (kapa, readout->image, readout->mask, "image", sqrt(factor), 0); pmVisualAskUser(NULL); @@ -292,6 +283,4 @@ bool psphotVisualShowBackground (pmConfig *config, const pmFPAview *view, pmReadout *readout) { - - pmReadout *backgnd; if (!pmVisualTestLevel("psphot.image.backgnd", 2)) return true; @@ -303,12 +292,13 @@ pmFPAfile *file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKGND"); - if (file->mode == PM_FPA_MODE_INTERNAL) { - backgnd = file->readout; - } else { - backgnd = pmFPAviewThisReadout (view, file->fpa); - } - - psphotVisualScaleImage (kapa, backgnd->image, readout->mask, "backgnd", 2); - psphotVisualScaleImage (kapa, readout->image, readout->mask, "backsub", 0); + pmReadout *backgnd = READOUT_OR_INTERNAL(view, file); + + float factor = 1.0; + if (readout->covariance) { + factor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); + } + + psphotVisualScaleImage (kapa, backgnd->image, readout->mask, "backgnd", 1.0, 2); + psphotVisualScaleImage (kapa, readout->image, readout->mask, "backsub", sqrt(factor), 0); pmVisualAskUser(NULL); @@ -339,4 +329,45 @@ psphotVisualRangeImage (kapa, lsig, "log-signif", 2, min, max); psFree (lsig); + + pmVisualAskUser(NULL); + return true; +} + +// requires psphotVisualShowImage +bool psphotVisualShowSources (psArray *sources) { + + int Noverlay; + KiiOverlay *overlay; + + if (!pmVisualTestLevel("psphot.objects.sources", 1)) return true; + + int kapa = psphotKapaChannel (1); + if (kapa == -1) return false; + + // note: this uses the Ohana allocation tools: + // ALLOCATE (overlay, KiiOverlay, 3*peaks->n + 1); + ALLOCATE (overlay, KiiOverlay, sources->n + 2); + + Noverlay = 0; + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + if (!source) continue; + + pmPeak *peak = source->peak; + if (!peak) continue; + + overlay[Noverlay].type = KII_OVERLAY_BOX; + overlay[Noverlay].x = peak->xf; + overlay[Noverlay].y = peak->yf; + overlay[Noverlay].dx = 4.0; + overlay[Noverlay].dy = 4.0; + overlay[Noverlay].angle = 0.0; + overlay[Noverlay].text = NULL; + Noverlay ++; + } + + KiiLoadOverlay (kapa, overlay, Noverlay, "blue"); + FREE (overlay); pmVisualAskUser(NULL); @@ -839,5 +870,13 @@ if (source->type != type) continue; - if (mode && !(source->mode & mode)) continue; + + if (mode == PM_SOURCE_MODE_PSFSTAR) { + bool keep = false; + keep |= (source->tmpFlags & PM_SOURCE_TMPF_CANDIDATE_PSFSTAR); + keep |= (source->mode & PM_SOURCE_MODE_PSFSTAR); + if (!keep) continue; + } else { + if (mode && !(source->mode & mode)) continue; + } pmMoments *moments = source->moments; @@ -974,5 +1013,5 @@ // examine PSF sources in S/N order (brightest first) - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); // counters to track the size of the image and area used in a row @@ -1126,5 +1165,5 @@ // examine PSF sources in S/N order (brightest first) - sources = psArraySort (sources, pmSourceSortBySN); + sources = psArraySort (sources, pmSourceSortByFlux); // counters to track the size of the image and area used in a row @@ -1227,5 +1266,5 @@ } - psphotVisualScaleImage (myKapa, outsat, NULL, "satstar", 2); + psphotVisualScaleImage (myKapa, outsat, NULL, "satstar", 1.0, 2); pmVisualAskUser(NULL); @@ -1355,5 +1394,5 @@ graphdata.xmax = +30.05; graphdata.ymin = -0.05; - graphdata.ymax = +5.05; + graphdata.ymax = +8.05; KapaSetLimits (myKapa, &graphdata); @@ -1415,5 +1454,5 @@ graphdata.xmax = +1.51; graphdata.ymin = -0.05; - graphdata.ymax = +5.05; + graphdata.ymax = +8.05; graphdata.color = KapaColorByName ("black"); KapaSetLimits (myKapa, &graphdata); @@ -2496,9 +2535,14 @@ if (myKapa == -1) return false; + float factor = 1.0; + if (readout->covariance) { + factor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); + } + if (reshow) { psphotVisualShowMask (myKapa, readout->mask, "mask", 2); - psphotVisualScaleImage (myKapa, readout->variance, readout->mask, "variance", 1); - } - psphotVisualScaleImage (myKapa, readout->image, readout->mask, "resid", 0); + psphotVisualScaleImage (myKapa, readout->variance, readout->mask, "variance", 1.0, 1); + } + psphotVisualScaleImage (myKapa, readout->image, readout->mask, "resid", sqrt(factor), 0); pmVisualAskUser(NULL); @@ -2653,6 +2697,4 @@ graphdata.ymax = -32.0; - FILE *f = fopen ("chisq.dat", "w"); - // construct the plot vectors int n = 0; @@ -2672,11 +2714,7 @@ graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[n]); graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[n]); - - fprintf (f, "%d %d %f %f\n", i, n, x->data.F32[n], y->data.F32[n]); - n++; } x->n = y->n = n; - fclose (f); float range; Index: /trunk/psphot/test/tap_psphot_forced.pro =================================================================== --- /trunk/psphot/test/tap_psphot_forced.pro (revision 31154) +++ /trunk/psphot/test/tap_psphot_forced.pro (revision 31154) @@ -0,0 +1,264 @@ +#!/usr/bin/env mana +# -*-sh-*- + +# forced photometry + +# 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 + + local i Npass + $Npass = $1 + + # use a limiting stellar density of 1 / deg^2 + mkref ref 1.0 + + # use a FWHM of 1.0 arcsec + mkexp image/image 1.0 ref + + for i 0 $Npass + mkforce image/image v$i + forcedCompare image/image.cmf image/image.v$i.frc.cmf image/image.v$i.frc.png 1.0 + end +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 + + $dirname = `dirname $basename` + exec mkdir -p $dirname + + $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 -seed 12345 + exec ppImage $ppImageConfig -file $basename.fits $basename -seed 12345 +end + +macro mkforce + if ($0 != 3) + echo "USAGE: mkforce (basename) (version) " + break + end + + local basename + $basename = $1 + echo $basename + + data $basename.cmf + read -fits Chip.psf X_PSF Y_PSF + write $basename.$2.frc.dat X_PSF Y_PSF + echo "here" + + $forcedOpt = -file $basename.ch.fits + $forcedOpt = $forcedOpt -mask $basename.ch.mk.fits + $forcedOpt = $forcedOpt -variance $basename.ch.wt.fits + $forcedOpt = $forcedOpt -psf $basename.psf + $forcedOpt = $forcedOpt -srctext $basename.$2.frc.dat + exec psphotForced $forcedOpt $basename.$2.frc -seed 12345 +end + +macro forcedCompare + if ($0 != 5) + echo "USAGE: forcedCompare (raw) (out) (output) (radius)" + break + end + + load.cmf $1 Chip.psf raw + load.cmf $2 Chip.psf out + + # if we compare to the pre-convolution deltas, we need to do this: + # 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 + + set m1 = PSF_INST_MAG_raw + set m1err = PSF_INST_MAG_SIG_raw + set m1err10 = PSF_INST_MAG_SIG_raw / 10.0 + reindex m2 = PSF_INST_MAG_out using index1 + set dm = m1 - m2 + + resize 1000 600 + label -fn courier 14 + lim m1 dm; clear; box; plot m1 dm -dy m1err10 + + set chi = dm^2 / m1err^2 + vstat -q chi + echo "median offset: $MEDIAN, chisq : $TOTAL, reduced chisq : {$TOTAL / $NPTS}" + + sprintf line "median offset: %6.4f, chisq : %6.4f, reduced chisq : %6.4f" $MEDIAN $TOTAL {$TOTAL / $NPTS} + label -x "standard inst mag" -y "standard mag - forced mag" -ur "$line" + label -ul "plotted errors are poisson / 10.0 for clarity" + 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_stackphot.pro =================================================================== --- /trunk/psphot/test/tap_psphot_stackphot.pro (revision 31154) +++ /trunk/psphot/test/tap_psphot_stackphot.pro (revision 31154) @@ -0,0 +1,791 @@ +#!/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 -Db PSF.CONVOLVE T + + 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 -R PPSTACK.OUTPUT.VARIANCE FITS.TYPE NONE + + 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 + +macro ckradialflux + data $1 + read -fits Chip.xrad X_APER Y_APER PSF_FWHM APER_FLUX APER_FLUX_ERR APER_FLUX_STDEV APER_FILL + read -fits Chip.psf X_PSF Y_PSF PSF_INST_MAG + set mag = PSF_INST_MAG + # XXX include nradii in header + for i 0 11 + set mA$i = -2.5*log(APER_FLUX:$i) + set dM$i = mA$i - PSF_INST_MAG + vstat dM$i + $DM$i = $MEDIAN + end + for i 1 11 + $j = $i - 1 + echo {$DM$i - $DM$j} + end + + lim mag -2 2; clear; box + for i 0 11 + plot mag dM$i -c black + end +end + +if ($SCRIPT) + fulltest 4 + exit 0 +end