Index: /branches/eam_branches/20090715/psphot/src/Makefile.am =================================================================== --- /branches/eam_branches/20090715/psphot/src/Makefile.am (revision 25274) +++ /branches/eam_branches/20090715/psphot/src/Makefile.am (revision 25275) @@ -25,6 +25,7 @@ libpsphot_la_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) -# bin_PROGRAMS = psphot psphotTest psphotMomentsStudy psphotPetrosianStudy -bin_PROGRAMS = psphotPetrosianStudy +bin_PROGRAMS = psphot psphotTest psphotMomentsStudy psphotPetrosianStudy +# bin_PROGRAMS = psphotPetrosianStudy +# bin_PROGRAMS = psphot psphot_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) @@ -68,13 +69,5 @@ psphotPetrosianStudy_SOURCES = \ - psphotPetrosianStudy.c \ - psphotPetrosianProfile.c \ - psphotRadialProfileByAngle.c \ - psphotRadiiFromProfiles.c \ - psphotEllipticalContour.c \ - psphotEllipticalProfile.c \ - psphotPetrosianRadialBins.c \ - psphotPetrosianVisual.c \ - psphotPetrosianStats.c + psphotPetrosianStudy.c libpsphot_la_SOURCES = \ @@ -145,5 +138,15 @@ psphotThreadTools.c \ psphotAddNoise.c \ + psphotPetrosianProfile.c \ + psphotRadialProfileByAngle.c \ + psphotRadiiFromProfiles.c \ + psphotEllipticalContour.c \ + psphotEllipticalProfile.c \ + psphotPetrosianRadialBins.c \ + psphotPetrosianVisual.c \ + psphotPetrosianStats.c \ + psphotPetrosianAnalysis.c \ pmPetrosian.c + # dropped? psphotGrowthCurve.c Index: /branches/eam_branches/20090715/psphot/src/pmPetrosian.h =================================================================== --- /branches/eam_branches/20090715/psphot/src/pmPetrosian.h (revision 25274) +++ /branches/eam_branches/20090715/psphot/src/pmPetrosian.h (revision 25275) @@ -35,5 +35,5 @@ bool psphotPetrosianFreeVectors(pmPetrosian *petrosian); -bool psphotPetrosianProfile (pmSource *source); +bool psphotPetrosianProfile (pmSource *source, float skynoise); bool psphotRadialProfilesByAngles (pmPetrosian *petro, pmSource *source, int Nsec, float Rmax); float psphotRadiusFromProfile (psVector *radius, psVector *flux, float fluxMin, float fluxMax); @@ -41,5 +41,5 @@ bool psphotEllipticalProfile (pmSource *source, pmPetrosian *petrosian); bool psphotEllipticalContour (pmPetrosian *petrosian); -bool psphotPetrosianRadialBins (pmSource *source, pmPetrosian *petrosian, float radiusMax); +bool psphotPetrosianRadialBins (pmSource *source, pmPetrosian *petrosian, float radiusMax, float skynoise); bool psphotPetrosianStats (pmPetrosian *petrosian); @@ -48,5 +48,5 @@ bool psphotPetrosianVisualProfileByAngle (psVector *radius, psVector *flux); -bool psphotPetrosianVisualProfileRadii (psVector *radius, psVector *flux, psVector *radiusBin, psVector *fluxBin, float RadiusRef); +bool psphotPetrosianVisualProfileRadii (psVector *radius, psVector *flux, psVector *radiusBin, psVector *fluxBin, float peakFlux, float RadiusRef); bool psphotPetrosianVisualEllipticalContour (pmPetrosian *petrosian); @@ -57,6 +57,4 @@ float petFlux, float radiusForFlux); -bool pmVisualLimitsFromVectors (Graphdata *graphdata, psVector *xVec, psVector *yVec); - /// @} Index: /branches/eam_branches/20090715/psphot/src/psphot.h =================================================================== --- /branches/eam_branches/20090715/psphot/src/psphot.h (revision 25274) +++ /branches/eam_branches/20090715/psphot/src/psphot.h (revision 25275) @@ -8,4 +8,5 @@ #include #include "psphotErrorCodes.h" +#include "pmPetrosian.h" #define PSPHOT_RECIPE "PSPHOT" // Name of the recipe to use @@ -192,4 +193,7 @@ bool psphotVisualPlotSourceSize (psArray *sources); +bool psphotVisualShowPetrosian (pmSource *source, pmPetrosian *petrosian); +bool psphotPetrosianAnalysis (pmReadout *readout, psArray *sources, psMetadata *recipe); + bool psphotImageQuality (psMetadata *recipe, psArray *sources); Index: /branches/eam_branches/20090715/psphot/src/psphotInternal.h =================================================================== --- /branches/eam_branches/20090715/psphot/src/psphotInternal.h (revision 25274) +++ /branches/eam_branches/20090715/psphot/src/psphotInternal.h (revision 25275) @@ -14,5 +14,4 @@ #include #include "psphot.h" -#include "pmPetrosian.h" #endif Index: /branches/eam_branches/20090715/psphot/src/psphotPetrosianAnalysis.c =================================================================== --- /branches/eam_branches/20090715/psphot/src/psphotPetrosianAnalysis.c (revision 25275) +++ /branches/eam_branches/20090715/psphot/src/psphotPetrosianAnalysis.c (revision 25275) @@ -0,0 +1,64 @@ +# include "psphotInternal.h" + +// aperture-like measurements for extended sources +bool psphotPetrosianAnalysis (pmReadout *readout, psArray *sources, psMetadata *recipe) { + + bool status; + + // 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); + + // 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, "EXTENDED_SOURCE_SN_LIM"); + + // option to limit analysis to a specific region + char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION"); + psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region)); + if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined"); + + // source analysis is done in S/N order (brightest first) + sources = psArraySort (sources, pmSourceSortBySN); + + // choose the sources of interest + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + + // 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; + + // 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 (source->extNsigma < 10.0) continue; // XXX this should not be hardwired + + // limit selection by analysis region + if (source->peak->x < AnalysisRegion.x0) continue; + if (source->peak->y < AnalysisRegion.y0) continue; + if (source->peak->x > AnalysisRegion.x1) continue; + if (source->peak->y > AnalysisRegion.y1) continue; + + // replace object in image + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + } + + psphotPetrosianProfile (source, skynoise); + } + + // pause and wait for user input: + // continue, save (provide name), ?? + char key[10]; + fprintf (stdout, "[c]ontinue? "); + if (!fgets(key, 8, stdin)) { + psWarning("Unable to read option"); + } + + return true; +} Index: /branches/eam_branches/20090715/psphot/src/psphotPetrosianProfile.c =================================================================== --- /branches/eam_branches/20090715/psphot/src/psphotPetrosianProfile.c (revision 25274) +++ /branches/eam_branches/20090715/psphot/src/psphotPetrosianProfile.c (revision 25275) @@ -7,5 +7,5 @@ // structure to something the pmRadialProfile -bool psphotPetrosianProfile (pmSource *source) { +bool psphotPetrosianProfile (pmSource *source, float skynoise) { // container to hold results from the radial profile analysis @@ -48,5 +48,5 @@ // integrate the radial profile for radial bins defined for the petrosian measurement: // SB_i (r_i) where \alpha r_i < r < \beta r_i - if (!psphotPetrosianRadialBins (source, petrosian, Rmax)) { + if (!psphotPetrosianRadialBins (source, petrosian, Rmax, skynoise)) { psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile"); return false; @@ -59,4 +59,8 @@ } + // XXX this will only work in the psphot context, not the psphotPetrosianStudy... + // XXX add the petrosian to the pmSource structure... + psphotVisualShowPetrosian (source, petrosian); + psphotPetrosianFreeVectors(petrosian); Index: /branches/eam_branches/20090715/psphot/src/psphotPetrosianRadialBins.c =================================================================== --- /branches/eam_branches/20090715/psphot/src/psphotPetrosianRadialBins.c (revision 25274) +++ /branches/eam_branches/20090715/psphot/src/psphotPetrosianRadialBins.c (revision 25275) @@ -13,5 +13,15 @@ // track the non-overlapping radius values. -bool psphotPetrosianRadialBins (pmSource *source, pmPetrosian *petrosian, float radiusMax) { +bool psphotPetrosianRadialBins (pmSource *source, pmPetrosian *petrosian, float radiusMax, float skynoise) { + + // XXX for testing, let's just set this to a value + + float skyModelErrorSQ = PS_SQR(skynoise); + +# if (0) + // if backModel or backStdev are missing, the values of sky and/or skyErr will be set to NAN + pmReadout *backModel = psphotSelectBackground (config, view); + pmReadout *backStdev = psphotSelectBackgroundStdev (config, view); +# endif psVector *radius = petrosian->radiusElliptical; @@ -104,4 +114,5 @@ psVector *values = psVectorAllocEmpty (flux->n, PS_TYPE_F32); psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); + // psStats *stats = psStatsAlloc(PS_STAT_FITTED_MEAN_V4 | PS_STAT_FITTED_STDEV_V4); // integrate flux, radius for each of these bins. since flux is sorted by radius, @@ -123,5 +134,11 @@ // binSB->data.F32[nOut] = stats->sampleMedian; binSB->data.F32[nOut] = stats->robustMedian; - binSBstdev->data.F32[nOut] = stats->robustStdev / sqrt(values->n); + binSBstdev->data.F32[nOut] = sqrt(PS_SQR(stats->robustStdev) / values->n + skyModelErrorSQ); + // binSB->data.F32[nOut] = stats->fittedMean; + // binSBstdev->data.F32[nOut] = sqrt(PS_SQR(stats->fittedStdev) / values->n + skyModelErrorSQ); + + // error in the SB is the stdev per bin / sqrt (number of pixels) + // added in quadrature to a fraction of the local sky (not the + // residual flux, but the sky from the sky model) if (1) { @@ -153,5 +170,5 @@ petrosian->binSBstdev = binSBstdev; - psphotPetrosianVisualProfileRadii (radius, flux, binRad, binSB, 0.0); + psphotPetrosianVisualProfileRadii (radius, flux, binRad, binSB, source->peak->flux, 0.0); psFree(radMin); Index: /branches/eam_branches/20090715/psphot/src/psphotPetrosianStats.c =================================================================== --- /branches/eam_branches/20090715/psphot/src/psphotPetrosianStats.c (revision 25274) +++ /branches/eam_branches/20090715/psphot/src/psphotPetrosianStats.c (revision 25275) @@ -32,4 +32,8 @@ float dFsum2 = 0.0; + float nSigma = 3.0; + int lowestSignificantRadius = 0; + float lowestSignificantRatio = 1.0; + int nOut = 0; for (int i = 0; i < binSB->n; i++) { @@ -48,5 +52,5 @@ areaInner += areaSum->data.F32[nOut-1]; fluxInner += fluxSum->data.F32[nOut-1]; - fluxInnerErr2 = PS_SQR(fluxSumErr2->data.F32[nOut-1] * areaSum->data.F32[nOut-1]); + fluxInnerErr2 += fluxSumErr2->data.F32[nOut-1]; } @@ -64,5 +68,5 @@ float meanSBerr = sqrt(fluxInnerErr2) / areaInner; - float ratioErr = ratio * sqrt(PS_SQR(binSBstdev->data.F32[i]/binSB->data.F32[i]) + PS_SQR(meanSBerr/meanSB->data.F32[nOut])); + float ratioErr = fabs(ratio) * sqrt(PS_SQR(binSBstdev->data.F32[i]/binSB->data.F32[i]) + PS_SQR(meanSBerr/meanSB->data.F32[nOut])); psVectorAppend(petRatioErr, ratioErr); @@ -74,11 +78,15 @@ if (1) { - fprintf (stderr, "%3d : %5.2f : %5.3f %5.3f : %5.3f %5.3f : %5.1f %5.1f\n", - i, refRadius->data.F32[nOut], binSB->data.F32[i], meanSB->data.F32[nOut], petRatio->data.F32[nOut], petRatioErr->data.F32[nOut], fluxSum->data.F32[nOut], areaSum->data.F32[nOut]); + fprintf (stderr, "%3d : %5.2f : %5.3f %5.3f : %5.3f %5.3f : %5.3f %5.3f : %5.3f %5.3f : %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); } // anytime we transition below the PETROSIAN_RATIO, calculate the radius and flux // we will keep and report the last (largest radius) value - if (above && (petRatio->data.F32[nOut] < PETROSIAN_RATIO)) { + if (above && (petRatio->data.F32[nOut] < PETROSIAN_RATIO) && (petRatio->data.F32[nOut] > nSigma*petRatioErr->data.F32[nOut])) { // interpolate Rvec between i-1 and i to PETROSIAN_RATIO to get flux (Fvec) and radius (rvec) if (i == 0) { @@ -93,4 +101,12 @@ } + // anytime we transition below the PETROSIAN_RATIO, calculate the radius and flux + // we will keep and report the last (largest radius) value + // find the last signficant measurement of the petrosian ratio + if (above && (petRatio->data.F32[nOut] < lowestSignificantRatio) && (petRatio->data.F32[nOut] > nSigma*petRatioErr->data.F32[nOut])) { + lowestSignificantRadius = nOut; + lowestSignificantRatio = petRatio->data.F32[nOut]; + } + // reset on transitions up, but do not re-calculate rad_90, flux_90 if (!above && (petRatio->data.F32[nOut] >= PETROSIAN_RATIO)) { @@ -101,5 +117,11 @@ if (!anyPetro) { - // set default radius: + // 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); + } else { + petRadius = InterpolateValues (petRatio->data.F32[lowestSignificantRadius-1], refRadius->data.F32[lowestSignificantRadius-1], petRatio->data.F32[lowestSignificantRadius], refRadius->data.F32[lowestSignificantRadius], PETROSIAN_RATIO); + } } @@ -110,5 +132,6 @@ if (refRadius->data.F32[i] > apRadius) { if (i == 0) { - psAbort ("does this case make any sense?"); + psWarning ("does this case make any sense? (refRadius[0] > apRadius)"); + continue; } else { petFlux = InterpolateValues (refRadius->data.F32[i-1], fluxSum->data.F32[i-1], refRadius->data.F32[i], fluxSum->data.F32[i], apRadius); @@ -119,4 +142,5 @@ // save petRadius, petFlux + // XXX save flags (anyPetro, manyPetro) petrosian->petrosianRadius = petRadius; petrosian->petrosianFlux = petFlux; Index: /branches/eam_branches/20090715/psphot/src/psphotPetrosianStudy.c =================================================================== --- /branches/eam_branches/20090715/psphot/src/psphotPetrosianStudy.c (revision 25274) +++ /branches/eam_branches/20090715/psphot/src/psphotPetrosianStudy.c (revision 25275) @@ -136,5 +136,5 @@ } - psphotPetrosianProfile (source); + psphotPetrosianProfile (source, skynoise); psFree (source); Index: /branches/eam_branches/20090715/psphot/src/psphotPetrosianVisual.c =================================================================== --- /branches/eam_branches/20090715/psphot/src/psphotPetrosianVisual.c (revision 25274) +++ /branches/eam_branches/20090715/psphot/src/psphotPetrosianVisual.c (revision 25275) @@ -93,5 +93,5 @@ } -bool psphotPetrosianVisualProfileRadii (psVector *radius, psVector *flux, psVector *radiusBin, psVector *fluxBin, float RadiusRef) { +bool psphotPetrosianVisualProfileRadii (psVector *radius, psVector *flux, psVector *radiusBin, psVector *fluxBin, float peakFlux, float RadiusRef) { float FluxRef = 500.0; @@ -115,5 +115,7 @@ KapaSetFont (kapa, "courier", 14); - pmVisualLimitsFromVectors (&graphdata, radius, flux); + graphdata.ymax = +1.05*peakFlux; + graphdata.ymin = -0.05*peakFlux; + pmVisualLimitsFromVectors (&graphdata, radius, NULL); KapaSetLimits (kapa, &graphdata); Index: /branches/eam_branches/20090715/psphot/src/psphotRadiiFromProfiles.c =================================================================== --- /branches/eam_branches/20090715/psphot/src/psphotRadiiFromProfiles.c (revision 25274) +++ /branches/eam_branches/20090715/psphot/src/psphotRadiiFromProfiles.c (revision 25275) @@ -56,4 +56,5 @@ } } + Rbin = 3; psVector *fluxBinned = NULL; @@ -117,5 +118,6 @@ // XXX is there a macro in psLib that does this interpolation? if (i == 0) { - psAbort ("impossible condition f[0] < Fo"); + psWarning ("impossible condition f[0] < Fo"); + continue; } else { Ro = radiusBinned->data.F32[i-1] + (radiusBinned->data.F32[i] - radiusBinned->data.F32[i-1]) * (Fo - fluxBinned->data.F32[i-1]) / (fluxBinned->data.F32[i] - fluxBinned->data.F32[i-1]);