Index: branches/pap/psphot/doc/stack.txt =================================================================== --- branches/pap/psphot/doc/stack.txt (revision 27726) +++ branches/pap/psphot/doc/stack.txt (revision 28003) @@ -1,2 +1,12 @@ + +20100503: + + psphotStack : outline + + * load N images + * make chisq detection image + * detect sources on all (N+1) images + * convolve N input images (to match, to target) + * measure source properties on convolved images 20100126: @@ -4,8 +14,4 @@ * watch out for psphotSetMomentsWindow & MOMENTS_SX_MAX,etc * watch out for psphotSignificanceImage.c:,psphotEfficiency using the FWHM_MAJ from psphotChoosePSF - * ppSimDetections.c : SIGMA_SMOOTH -ppSmooth/src/ppSmoothReadout.c: psMetadataAddF32(recipe, PS_LIST_TAIL, "EFFECTIVE_AREA", PS_META_REPLACE, "Effective Area", effArea); -ppSmooth/src/ppSmoothReadout.c: psMetadataAddF32(recipe, PS_LIST_TAIL, "SIGNIFICANCE_SCALE_FACTOR", PS_META_REPLACE, "Signicance scale factor", factor); - 20100120 : more stack processing mods: Index: branches/pap/psphot/src/Makefile.am =================================================================== --- branches/pap/psphot/src/Makefile.am (revision 27726) +++ branches/pap/psphot/src/Makefile.am (revision 28003) @@ -185,7 +185,9 @@ psphotEllipticalContour.c \ psphotEllipticalProfile.c \ + psphotRadialBins.c \ psphotPetrosian.c \ psphotPetrosianRadialBins.c \ psphotPetrosianStats.c \ + psphotPetrosianVisual.c \ psphotEfficiency.c @@ -197,5 +199,4 @@ # psphotAnnuli.c \ # psphotKron.c \ -# psphotPetrosianVisual.c \ # Index: branches/pap/psphot/src/psphot.h =================================================================== --- branches/pap/psphot/src/psphot.h (revision 27726) +++ branches/pap/psphot/src/psphot.h (revision 28003) @@ -245,5 +245,5 @@ float psphotRadiusFromProfile (pmSource *source, psVector *radius, psVector *flux, float fluxMin, float fluxMax); bool psphotRadiiFromProfiles (pmSource *source, float fluxMin, float fluxMax); -bool psphotEllipticalProfile (pmSource *source); +bool psphotEllipticalProfile (pmSource *source, bool RAW_RADIUS); bool psphotEllipticalContour (pmSource *source); @@ -281,12 +281,14 @@ // XXX visualization functions related to radial profiles (disabled) -// bool psphotPetrosianVisualProfileByAngle (psVector *radius, psVector *flux); -// bool psphotPetrosianVisualProfileRadii (psVector *radius, psVector *flux, psVector *radiusBin, psVector *fluxBin, float peakFlux, float RadiusRef); -// bool psphotPetrosianVisualEllipticalContour (pmPetrosian *petrosian); -// bool psphotPetrosianVisualStats (psVector *radBin, psVector *fluxBin, -// psVector *refRadius, psVector *meanSB, -// psVector *petRatio, psVector *petRatioErr, psVector *fluxSum, -// float petRadius, float ratioForRadius, -// float petFlux, float radiusForFlux); +bool psphotPetrosianVisualProfileByAngle (psVector *radius, psVector *flux); +bool psphotPetrosianVisualProfileRadii (psVector *radius, psVector *flux, psVector *radiusBin, psVector *fluxBin, float peakFlux, float RadiusRef); +bool psphotPetrosianVisualEllipticalContour (pmSourceRadialFlux *radFlux, pmSourceExtendedPars *extpars); +bool psphotPetrosianVisualStats (psVector *radBin, psVector *fluxBin, + psVector *refRadius, psVector *meanSB, + psVector *petRatio, psVector *petRatioErr, psVector *fluxSum, + float petRadius, float ratioForRadius, + float petFlux, float radiusForFlux); + +bool psphotRadialBins (psMetadata *recipe, pmSource *source, float radiusMax, float skynoise); // structures & functions to support psf-convolved model fitting Index: branches/pap/psphot/src/psphotCleanup.c =================================================================== --- branches/pap/psphot/src/psphotCleanup.c (revision 27726) +++ branches/pap/psphot/src/psphotCleanup.c (revision 28003) @@ -14,9 +14,6 @@ psFree (config); - psThreadPoolFinalize (); - psTimerStop (); psMemCheckCorruption (stderr, true); pmModelClassCleanup (); - psTimeFinalize (); pmSourceFitSetDone (); pmConceptsDone (); Index: branches/pap/psphot/src/psphotEllipticalContour.c =================================================================== --- branches/pap/psphot/src/psphotEllipticalContour.c (revision 27726) +++ branches/pap/psphot/src/psphotEllipticalContour.c (revision 28003) @@ -7,5 +7,10 @@ bool psphotEllipticalContour (pmSource *source) { - pmSourceRadialProfile *profile = source->extpars->profile; + psAssert (source, "missing source"); + psAssert (source->extpars, "missing extpars"); + psAssert (source->extpars->radFlux, "missing radFlux"); + + pmSourceRadialFlux *profile = source->extpars->radFlux; + pmSourceExtendedPars *extpars = source->extpars; // use LMM to fit theta vs radius to an ellipse @@ -85,20 +90,20 @@ /// XXX rationalize? if epsilon > 1, flip major and minor axes (rotate by 90 degrees) if (params->data.F32[PAR_EPSILON] < 1.0) { - profile->axes.major = params->data.F32[PAR_RMIN] / params->data.F32[PAR_EPSILON]; - profile->axes.minor = params->data.F32[PAR_RMIN]; - profile->axes.theta = params->data.F32[PAR_PHI]; + extpars->axes.major = params->data.F32[PAR_RMIN] / params->data.F32[PAR_EPSILON]; + extpars->axes.minor = params->data.F32[PAR_RMIN]; + extpars->axes.theta = params->data.F32[PAR_PHI]; } else { - profile->axes.major = params->data.F32[PAR_RMIN]; - profile->axes.minor = params->data.F32[PAR_RMIN] / params->data.F32[PAR_EPSILON]; - profile->axes.theta = params->data.F32[PAR_PHI] + 0.5*M_PI; + extpars->axes.major = params->data.F32[PAR_RMIN]; + extpars->axes.minor = params->data.F32[PAR_RMIN] / params->data.F32[PAR_EPSILON]; + extpars->axes.theta = params->data.F32[PAR_PHI] + 0.5*M_PI; } psTrace ("psphot", 4, "# fitted values:\n"); - psTrace ("psphot", 4, "Phi: %f\n", profile->axes.theta*PS_DEG_RAD); - psTrace ("psphot", 4, "Rmaj: %f\n", profile->axes.major); - psTrace ("psphot", 4, "Rmin: %f\n", profile->axes.minor); + psTrace ("psphot", 4, "Phi: %f\n", extpars->axes.theta*PS_DEG_RAD); + psTrace ("psphot", 4, "Rmaj: %f\n", extpars->axes.major); + psTrace ("psphot", 4, "Rmin: %f\n", extpars->axes.minor); // show the results - // psphotPetrosianVisualEllipticalContour (petrosian); + // psphotPetrosianVisualEllipticalContour (profile, extpars); psFree (x); Index: branches/pap/psphot/src/psphotEllipticalProfile.c =================================================================== --- branches/pap/psphot/src/psphotEllipticalProfile.c (revision 27726) +++ branches/pap/psphot/src/psphotEllipticalProfile.c (revision 28003) @@ -1,7 +1,16 @@ # include "psphotInternal.h" -bool psphotEllipticalProfile (pmSource *source) { +bool psphotEllipticalProfile (pmSource *source, bool RAW_RADIUS) { - pmSourceRadialProfile *profile = source->extpars->profile; + psAssert (source, "missing source"); + psAssert (source->extpars, "missing extpars"); + psAssert (source->pixels, "missing pixels"); + + pmSourceExtendedPars *extpars = source->extpars; + + if (!source->extpars->ellipticalFlux) { + source->extpars->ellipticalFlux = pmSourceEllipticalFluxAlloc(); + } + pmSourceEllipticalFlux *profile = source->extpars->ellipticalFlux; profile->radiusElliptical = psVectorAllocEmpty(100, PS_TYPE_F32); @@ -21,11 +30,17 @@ psEllipseAxes axes; - axes.major = M_SQRT1_2; - axes.minor = M_SQRT1_2 * (profile->axes.minor / profile->axes.major); + if (RAW_RADIUS) { + // force circular profile + axes.major = M_SQRT1_2; + axes.minor = M_SQRT1_2; + } else { + axes.major = M_SQRT1_2; + axes.minor = M_SQRT1_2 * (extpars->axes.minor / extpars->axes.major); + } // axes.major = 1.0; - // axes.minor = profile->axes.minor / profile->axes.major; + // axes.minor = extpars->axes.minor / extpars->axes.major; - axes.theta = profile->axes.theta; + axes.theta = extpars->axes.theta; psEllipseShape shape = psEllipseAxesToShape (axes); @@ -46,9 +61,9 @@ float r2 = 0.5*PS_SQR(x/Sxx) + 0.5*PS_SQR(y/Syy) + x*y*Sxy; + float Rraw = hypot(x, y); psVectorAppend(radius, sqrt(r2)); psVectorAppend(flux, source->pixels->data.F32[iy][ix]); - float Rraw = hypot(x, y); psVectorAppend(radiusRaw, Rraw); psVectorAppend(fluxRaw, source->pixels->data.F32[iy][ix]); @@ -67,5 +82,5 @@ // } - // psphotPetrosianVisualProfileRadii (radius, flux, radiusRaw, fluxRaw, 0.0); + psphotPetrosianVisualProfileRadii (radius, flux, radiusRaw, fluxRaw, source->peak->flux, 0.0); // psphotPetrosianVisualProfileByAngle (radius, flux); Index: branches/pap/psphot/src/psphotExtendedSourceAnalysis.c =================================================================== --- branches/pap/psphot/src/psphotExtendedSourceAnalysis.c (revision 27726) +++ branches/pap/psphot/src/psphotExtendedSourceAnalysis.c (revision 28003) @@ -39,4 +39,6 @@ int Nkron = 0; + psTimerStart ("psphot.extended"); + // find the currently selected readout pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest @@ -66,11 +68,4 @@ float skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE"); -# 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 - // S/N limit to perform full non-linear fits float SN_LIM = psMetadataLookupF32 (&status, recipe, "EXTENDED_SOURCE_SN_LIM"); @@ -81,4 +76,11 @@ 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 // source analysis is done in S/N order (brightest first) @@ -119,4 +121,11 @@ Next ++; + // 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 we request any of these measurements, we require the radial profile if (doPetrosian || doIsophotal || doAnnuli || doKron) { @@ -134,5 +143,5 @@ if (doPetrosian) { if (!psphotPetrosian (source, recipe, skynoise, maskVal)) { - psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); + psTrace ("psphot", 5, "FAILED petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); } else { psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); @@ -169,9 +178,10 @@ if (source->extpars) { - pmSourceRadialProfileFreeVectors(source->extpars->profile); + psFree(source->extpars->radFlux); + psFree(source->extpars->ellipticalFlux); } } - psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot"), Next); + 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); Index: branches/pap/psphot/src/psphotMosaicChip.c =================================================================== --- branches/pap/psphot/src/psphotMosaicChip.c (revision 27726) +++ branches/pap/psphot/src/psphotMosaicChip.c (revision 28003) @@ -21,6 +21,5 @@ pmChip *inChip = pmFPAviewThisChip(view, in->fpa); if (!outChip->hdu && !outChip->parent->hdu) { - const char *name = psMetadataLookupStr(&status, in->fpa->concepts, "FPA.OBS"); // Name of FPA - pmFPAAddSourceFromView(out->fpa, name, view, out->format); + pmFPAAddSourceFromView(out->fpa, view, out->format); } Index: branches/pap/psphot/src/psphotPetrosian.c =================================================================== --- branches/pap/psphot/src/psphotPetrosian.c (revision 27726) +++ branches/pap/psphot/src/psphotPetrosian.c (revision 28003) @@ -7,5 +7,5 @@ psAssert (source->extpars, "need to run psphotRadialProfile first"); - psAssert (source->extpars->profile, "need to run psphotRadialProfile first"); + psAssert (source->extpars->ellipticalFlux, "need to run psphotRadialProfile first"); // integrate the radial profile for radial bins defined for the petrosian measurement: @@ -24,8 +24,8 @@ psTrace ("psphot", 3, "source at %f,%f: petrosian radius: %f, flux: %f, axis ratio: %f, angle: %f", source->peak->xf, source->peak->yf, - source->extpars->petrosian_80->radius, - source->extpars->petrosian_80->flux, - source->extpars->profile->axes.minor/source->extpars->profile->axes.major, - source->extpars->profile->axes.theta*PS_DEG_RAD); + source->extpars->petrosianRadius, + source->extpars->petrosianFlux, + source->extpars->axes.minor/source->extpars->axes.major, + source->extpars->axes.theta*PS_DEG_RAD); return true; Index: branches/pap/psphot/src/psphotPetrosianRadialBins.c =================================================================== --- branches/pap/psphot/src/psphotPetrosianRadialBins.c (revision 27726) +++ branches/pap/psphot/src/psphotPetrosianRadialBins.c (revision 28003) @@ -1,3 +1,4 @@ # include "psphotInternal.h" +float InterpolateValues (float X0, float Y0, float X1, float Y1, float X); // convert the flux vs elliptical radius to annular bins @@ -5,5 +6,5 @@ // we are guaranteed to be limited by either the seeing (1 - few pixels) or by the pixels // themselves. this function does not attempt to measure the radial profiles accurately -// for radii that are smaller than a minimum (currently 1.0 pixels). +// for radii that are smaller than ~2 pixels // for small radii, we are measuring the mean surface brightness in non-overlapping radial @@ -13,16 +14,27 @@ // track the non-overlapping radius values. +// Photo interpolates the image of interest to place the peak on the center of the central +// pixel, and then uses the exact fractions of the pixels in each of the first few annuli. +// Seems like a reasonable thing, but is there any significance to the difference? + // XXX move the resulting elements from profile to extpars->petrosian? bool psphotPetrosianRadialBins (pmSource *source, float radiusMax, float skynoise) { - pmSourceRadialProfile *profile = source->extpars->profile; - - float skyModelErrorSQ = PS_SQR(skynoise); - - psVector *radius = profile->radiusElliptical; - psVector *flux = profile->fluxElliptical; + psAssert (source, "missing source"); + psAssert (source->extpars, "missing extpars"); + psAssert (source->extpars->ellipticalFlux, "missing ellipticalFlux"); + + psVector *radius = source->extpars->ellipticalFlux->radiusElliptical; + psVector *flux = source->extpars->ellipticalFlux->fluxElliptical; // sort incoming vectors by radius pmSourceRadialProfileSortPair (radius, flux); + + if (!source->extpars->petProfile) { + source->extpars->petProfile = pmSourceRadialProfileAlloc(); + } + pmSourceRadialProfile *profile = source->extpars->petProfile; + + float skyModelErrorSQ = PS_SQR(skynoise); int nMax = radiusMax; @@ -107,8 +119,4 @@ 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, - // we can do this fairly quickly bool done = false; @@ -133,9 +141,7 @@ dvalue = NAN; } - // binSB->data.F32[nOut] = stats->sampleMedian; + binSB->data.F32[nOut] = value; binSBstdev->data.F32[nOut] = sqrt(PS_SQR(dvalue) / 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) @@ -143,6 +149,5 @@ // residual flux, but the sky from the sky model) - psTrace ("psphot", 5, "%3d %5.1f %5.1f : %5.1f %5.2f\n", - nOut, radAlp->data.F32[nOut], radBet->data.F32[nOut], binSB->data.F32[nOut], binSBstdev->data.F32[nOut]); + psTrace ("psphot", 5, "%3d %5.1f %5.1f : %5.1f %5.2f\n", nOut, radAlp->data.F32[nOut], radBet->data.F32[nOut], binSB->data.F32[nOut], binSBstdev->data.F32[nOut]); nOut ++; @@ -163,4 +168,39 @@ // XXX I think this misses the last radial bin -- do we care? + // interpolate any bins that were empty (extrapolate to center if needed) + if (!isfinite(binSB->data.F32[0]) && !isfinite(binSB->data.F32[1])) { + psWarning ("center 2 bins of source at %f, %f are NAN, skipping this source", source->peak->xf, source->peak->yf); + // XXX raise a flag + psFree(binSB); + psFree(binSBstdev); + psFree(binRad); + psFree(binArea); + psFree(radMin); + psFree(radMax); + psFree(radAlp); + psFree(radBet); + psFree(values); + psFree(stats); + return false; + } + + // if center bin is empty assume same SB as next radius (probably true due to PSF) + if (!isfinite(binSB->data.F32[0])) { + binSB->data.F32[0] = binSB->data.F32[1]; + binSBstdev->data.F32[0] = binSBstdev->data.F32[1]; + } + + // interpolate any bins that were empty (if center if needed) + for (int i = 1; i < binSB->n - 1; i++) { + if (isfinite(binSB->data.F32[i])) continue; + binSB->data.F32[i] = InterpolateValues (binRad->data.F32[i-1], binSB->data.F32[i-1], binRad->data.F32[i+1], binSB->data.F32[i+1], binRad->data.F32[i]); + binSBstdev->data.F32[i] = InterpolateValues (binRad->data.F32[i-1], binSBstdev->data.F32[i-1], binRad->data.F32[i+1], binSBstdev->data.F32[i+1], binRad->data.F32[i]); + } + + psFree(profile->binSB); + psFree(profile->binSBstdev); + psFree(profile->radialBins); + psFree(profile->area); + // save the vectors profile->radialBins = binRad; Index: branches/pap/psphot/src/psphotPetrosianStats.c =================================================================== --- branches/pap/psphot/src/psphotPetrosianStats.c (revision 27726) +++ branches/pap/psphot/src/psphotPetrosianStats.c (revision 28003) @@ -10,8 +10,9 @@ bool psphotPetrosianStats (pmSource *source) { - pmSourceRadialProfile *profile = source->extpars->profile; - - float petRadius = NAN; - float petFlux = NAN; + psAssert (source, "missing source"); + psAssert (source->extpars, "missing extpars"); + psAssert (source->extpars->petProfile, "missing petProfile"); + + pmSourceRadialProfile *profile = source->extpars->petProfile; psVector *binSB = profile->binSB; @@ -28,4 +29,7 @@ psVector *areaSum = psVectorAllocEmpty(binSB->n, PS_TYPE_F32); + float petRadius = NAN; + float petFlux = NAN; + bool anyPetro = false; bool manyPetro = false; @@ -38,4 +42,6 @@ int lowestSignificantRadius = 0; float lowestSignificantRatio = 1.0; + + // find the Petrosian Radius and Petrosian Flux int nOut = 0; @@ -142,14 +148,48 @@ } - if (!source->extpars->petrosian_80) { - source->extpars->petrosian_80 = pmSourceExtendedFluxAlloc (); - } - pmSourceExtendedFlux *petrosian = source->extpars->petrosian_80; + // now measure the radii R90 and R50 where flux = 0.9 (or 0.5) * petFlux; + float flux90 = 0.9 * petFlux; + float flux50 = 0.5 * petFlux; + float R50 = NAN; + float R90 = NAN; + bool found50 = false; + bool found90 = false; + // XXX use bisection to do this faster: + for (int i = 0; !(found50 && found90) && i < refRadius->n; i++) { + if (!found50 && (fluxSum->data.F32[i] > flux50)) { + if (i == 0) { + psWarning ("does this case make any sense? (fluxSum[0] > flux50)"); + continue; + } else { + R50 = InterpolateValues (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux50); + found50 = true; + } + } + if (!found90 && (fluxSum->data.F32[i] > flux90)) { + if (i == 0) { + psWarning ("does this case make any sense? (fluxSum[0] > flux90)"); + continue; + } else { + R90 = InterpolateValues (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux90); + found90 = true; + } + } + } + // XXX save flags (anyPetro, manyPetro) - petrosian->radius = petRadius; - petrosian->flux = petFlux; - - // psphotPetrosianVisualStats (binRad, binSB, refRadius, meanSB, petRatio, petRatioErr, fluxSum, petRadius, PETROSIAN_RATIO, petFlux, apRadius); + source->extpars->petrosianRadius = petRadius; + source->extpars->petrosianFlux = petFlux; + source->extpars->petrosianR50 = R50; + source->extpars->petrosianR90 = R90; + + // XXX add the errors + source->extpars->petrosianRadiusErr = NAN; + source->extpars->petrosianFluxErr = NAN; + source->extpars->petrosianR50Err = NAN; + source->extpars->petrosianR90Err = NAN; + + fprintf (stderr, "source @ %f,%f\n", source->peak->xf, source->peak->yf); + psphotPetrosianVisualStats (binRad, binSB, refRadius, meanSB, petRatio, petRatioErr, fluxSum, petRadius, PETROSIAN_RATIO, petFlux, apRadius); psFree(fluxSum); Index: branches/pap/psphot/src/psphotPetrosianVisual.c =================================================================== --- branches/pap/psphot/src/psphotPetrosianVisual.c (revision 27726) +++ branches/pap/psphot/src/psphotPetrosianVisual.c (revision 28003) @@ -1,3 +1,4 @@ # include "psphotInternal.h" +# define FORCE_VISUAL 0 // this function displays representative images as the psphot analysis progresses: @@ -54,5 +55,5 @@ // return true; - if (!pmVisualIsVisual()) return true; + if (!FORCE_VISUAL && !pmVisualIsVisual()) return true; if (kapa2 == -1) { @@ -101,5 +102,5 @@ // return true; - if (!pmVisualIsVisual()) return true; + if (!FORCE_VISUAL && !pmVisualIsVisual()) return true; if (kapa == -1) { @@ -172,5 +173,143 @@ KapaSection section; - if (!pmVisualIsVisual()) return true; + if (!FORCE_VISUAL && !pmVisualIsVisual()) return true; + + if (kapa2 == -1) { + kapa2 = KapaOpenNamedSocket ("kapa", "psphot:stats"); + if (kapa2 == -1) { + fprintf (stderr, "failure to open kapa; visual mode disabled\n"); + pmVisualSetVisual(false); + return false; + } + } + + KapaClearPlots (kapa2); + KapaInitGraph (&graphdata); + KapaSetFont (kapa2, "courier", 14); + + // radius vs flux + // radius vs mean SB + // radius vs petRatio + + // *** section 1: radius vs mean SB + section.dx = 1.00; + section.dy = 0.33; + section.x = 0.00; + section.y = 0.00; + section.name = psStringCopy ("meanSB"); + KapaSetSection (kapa2, §ion); + psFree (section.name); + + graphdata.color = KapaColorByName ("black"); + pmVisualLimitsFromVectors (&graphdata, radBin, fluxBin); + KapaSetLimits (kapa2, &graphdata); + + KapaBox (kapa2, &graphdata); + KapaSendLabel (kapa2, "radius", KAPA_LABEL_XM); + KapaSendLabel (kapa2, "mean SB", KAPA_LABEL_YM); + + graphdata.color = KapaColorByName ("black"); + graphdata.style = 2; + graphdata.ptype = 0; + graphdata.size = 1.0; + KapaPrepPlot (kapa2, radBin->n, &graphdata); + KapaPlotVector (kapa2, radBin->n, radBin->data.F32, "x"); + KapaPlotVector (kapa2, radBin->n, fluxBin->data.F32, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.style = 2; + graphdata.ptype = 1; + graphdata.size = 2.0; + KapaPrepPlot (kapa2, refRadius->n, &graphdata); + KapaPlotVector (kapa2, refRadius->n, refRadius->data.F32, "x"); + KapaPlotVector (kapa2, refRadius->n, meanSB->data.F32, "y"); + + // *** section 2: radius vs petrosian ratio + section.dx = 1.00; + section.dy = 0.33; + section.x = 0.00; + section.y = 0.33; + section.name = psStringCopy ("ratio"); + KapaSetSection (kapa2, §ion); + psFree (section.name); + + graphdata.color = KapaColorByName ("black"); + graphdata.ymax = +1.05; + graphdata.ymin = -0.05; + pmVisualLimitsFromVectors (&graphdata, radBin, NULL); + KapaSetLimits (kapa2, &graphdata); + + KapaBox (kapa2, &graphdata); + KapaSendLabel (kapa2, "radius", KAPA_LABEL_XM); + KapaSendLabel (kapa2, "ratio", KAPA_LABEL_YM); + + graphdata.color = KapaColorByName ("black"); + graphdata.style = 2; + graphdata.ptype = 0; + graphdata.size = 1.0; + graphdata.etype = 0x01; + KapaPrepPlot (kapa2, refRadius->n, &graphdata); + KapaPlotVector (kapa2, refRadius->n, refRadius->data.F32, "x"); + KapaPlotVector (kapa2, refRadius->n, petRatio->data.F32, "y"); + KapaPlotVector (kapa2, refRadius->n, petRatioErr->data.F32, "dym"); + KapaPlotVector (kapa2, refRadius->n, petRatioErr->data.F32, "dyp"); + graphdata.etype = 0; + + graphdata.color = KapaColorByName ("red"); + graphdata.style = 2; + graphdata.ptype = 2; + graphdata.size = 2.0; + KapaPrepPlot (kapa2, 1, &graphdata); + KapaPlotVector (kapa2, 1, &petRadius, "x"); + KapaPlotVector (kapa2, 1, &ratioForRadius, "y"); + + // *** section 3: radius vs integrated flux + section.dx = 1.00; + section.dy = 0.33; + section.x = 0.00; + section.y = 0.66; + section.name = psStringCopy ("flux"); + KapaSetSection (kapa2, §ion); + psFree (section.name); + + graphdata.color = KapaColorByName ("black"); + pmVisualLimitsFromVectors (&graphdata, radBin, fluxSum); + KapaSetLimits (kapa2, &graphdata); + + KapaBox (kapa2, &graphdata); + KapaSendLabel (kapa2, "radius", KAPA_LABEL_XM); + KapaSendLabel (kapa2, "integrated flux", KAPA_LABEL_YM); + + graphdata.color = KapaColorByName ("black"); + graphdata.style = 2; + graphdata.ptype = 0; + graphdata.size = 1.0; + KapaPrepPlot (kapa2, refRadius->n, &graphdata); + KapaPlotVector (kapa2, refRadius->n, refRadius->data.F32, "x"); + KapaPlotVector (kapa2, refRadius->n, fluxSum->data.F32, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 2; + graphdata.style = 2; + graphdata.size = 2.0; + KapaPrepPlot (kapa2, 1, &graphdata); + KapaPlotVector (kapa2, 1, &radiusForFlux, "x"); + KapaPlotVector (kapa2, 1, &petFlux, "y"); + + // 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; +} + +bool psphotPetrosianVisualEllipticalContour (pmSourceRadialFlux *radFlux, pmSourceExtendedPars *extpars) { + + Graphdata graphdata; + + if (!FORCE_VISUAL && !pmVisualIsVisual()) return true; if (kapa == -1) { @@ -187,112 +326,71 @@ KapaSetFont (kapa, "courier", 14); - // radius vs flux - // radius vs mean SB - // radius vs petRatio - - // *** section 1: radius vs mean SB - section.dx = 1.00; - section.dy = 0.33; - section.x = 0.00; - section.y = 0.00; - section.name = psStringCopy ("meanSB"); - KapaSetSection (kapa, §ion); - psFree (section.name); - - graphdata.color = KapaColorByName ("black"); - pmVisualLimitsFromVectors (&graphdata, radBin, fluxBin); + psVector *theta = radFlux->theta; + psVector *radius = radFlux->isophotalRadii; + + // find Rmin and Rmax for the initial guess + float Rmin = radius->data.F32[0]; + float Rmax = radius->data.F32[0]; + + psVector *Rx = psVectorAlloc(radius->n, PS_TYPE_F32); + psVector *Ry = psVectorAlloc(radius->n, PS_TYPE_F32); + + for (int i = 0; i < theta->n; i++) { + Rx->data.F32[i] = radius->data.F32[i]*cos(theta->data.F32[i]); + Ry->data.F32[i] = radius->data.F32[i]*sin(theta->data.F32[i]); + + // check the radius range + Rmin = MIN (Rmin, radius->data.F32[i]); + Rmax = MAX (Rmax, radius->data.F32[i]); + } + + psVector *rx = psVectorAlloc(361, PS_TYPE_F32); + psVector *ry = psVectorAlloc(361, PS_TYPE_F32); + + float epsilon = extpars->axes.minor / extpars->axes.major; + + for (int i = 0; i < 361; i++) { + + float alpha = PS_RAD_DEG * i; + + float cs_alpha = cos(alpha); + float sn_alpha = sin(alpha); + + float cs_phi = cos(alpha - extpars->axes.theta); + float sn_phi = sin(alpha - extpars->axes.theta); + + float r = 1.0 / sqrt(SQ(sn_phi) + SQ(epsilon*cs_phi)); + + // generate the model fit here + rx->data.F32[i] = extpars->axes.minor * cs_alpha * r; + ry->data.F32[i] = extpars->axes.minor * sn_alpha * r; + } + + graphdata.color = KapaColorByName ("black"); + graphdata.xmin = -1.1*Rmax; + graphdata.ymin = -1.1*Rmax; + graphdata.xmax = +1.1*Rmax; + graphdata.ymax = +1.1*Rmax; KapaSetLimits (kapa, &graphdata); KapaBox (kapa, &graphdata); - KapaSendLabel (kapa, "radius", KAPA_LABEL_XM); - KapaSendLabel (kapa, "mean SB", KAPA_LABEL_YM); - - graphdata.color = KapaColorByName ("black"); - graphdata.style = 2; - graphdata.ptype = 0; - graphdata.size = 1.0; - KapaPrepPlot (kapa, radBin->n, &graphdata); - KapaPlotVector (kapa, radBin->n, radBin->data.F32, "x"); - KapaPlotVector (kapa, radBin->n, fluxBin->data.F32, "y"); + KapaSendLabel (kapa, "R_x", KAPA_LABEL_XM); + KapaSendLabel (kapa, "R_y", KAPA_LABEL_YM); graphdata.color = KapaColorByName ("red"); graphdata.style = 2; - graphdata.ptype = 1; - graphdata.size = 2.0; - KapaPrepPlot (kapa, refRadius->n, &graphdata); - KapaPlotVector (kapa, refRadius->n, refRadius->data.F32, "x"); - KapaPlotVector (kapa, refRadius->n, meanSB->data.F32, "y"); - - // *** section 2: radius vs petrosian ratio - section.dx = 1.00; - section.dy = 0.33; - section.x = 0.00; - section.y = 0.33; - section.name = psStringCopy ("ratio"); - KapaSetSection (kapa, §ion); - psFree (section.name); - - graphdata.color = KapaColorByName ("black"); - graphdata.ymax = +1.05; - graphdata.ymin = -0.05; - pmVisualLimitsFromVectors (&graphdata, radBin, NULL); - KapaSetLimits (kapa, &graphdata); - - KapaBox (kapa, &graphdata); - KapaSendLabel (kapa, "radius", KAPA_LABEL_XM); - KapaSendLabel (kapa, "ratio", KAPA_LABEL_YM); - - graphdata.color = KapaColorByName ("black"); - graphdata.style = 2; - graphdata.ptype = 0; - graphdata.size = 1.0; - graphdata.etype = 0x01; - KapaPrepPlot (kapa, refRadius->n, &graphdata); - KapaPlotVector (kapa, refRadius->n, refRadius->data.F32, "x"); - KapaPlotVector (kapa, refRadius->n, petRatio->data.F32, "y"); - KapaPlotVector (kapa, refRadius->n, petRatioErr->data.F32, "dym"); - KapaPlotVector (kapa, refRadius->n, petRatioErr->data.F32, "dyp"); - graphdata.etype = 0; - - graphdata.color = KapaColorByName ("red"); - graphdata.style = 2; graphdata.ptype = 2; - graphdata.size = 2.0; - KapaPrepPlot (kapa, 1, &graphdata); - KapaPlotVector (kapa, 1, &petRadius, "x"); - KapaPlotVector (kapa, 1, &ratioForRadius, "y"); - - // *** section 3: radius vs integrated flux - section.dx = 1.00; - section.dy = 0.33; - section.x = 0.00; - section.y = 0.66; - section.name = psStringCopy ("flux"); - KapaSetSection (kapa, §ion); - psFree (section.name); - - graphdata.color = KapaColorByName ("black"); - pmVisualLimitsFromVectors (&graphdata, radBin, fluxSum); - KapaSetLimits (kapa, &graphdata); - - KapaBox (kapa, &graphdata); - KapaSendLabel (kapa, "radius", KAPA_LABEL_XM); - KapaSendLabel (kapa, "integrated flux", KAPA_LABEL_YM); - - graphdata.color = KapaColorByName ("black"); - graphdata.style = 2; - graphdata.ptype = 0; - graphdata.size = 1.0; - KapaPrepPlot (kapa, refRadius->n, &graphdata); - KapaPlotVector (kapa, refRadius->n, refRadius->data.F32, "x"); - KapaPlotVector (kapa, refRadius->n, fluxSum->data.F32, "y"); - - graphdata.color = KapaColorByName ("red"); - graphdata.ptype = 2; - graphdata.style = 2; - graphdata.size = 2.0; - KapaPrepPlot (kapa, 1, &graphdata); - KapaPlotVector (kapa, 1, &radiusForFlux, "x"); - KapaPlotVector (kapa, 1, &petFlux, "y"); + graphdata.size = 1.0; + KapaPrepPlot (kapa, Rx->n, &graphdata); + KapaPlotVector (kapa, Rx->n, Rx->data.F32, "x"); + KapaPlotVector (kapa, Rx->n, Ry->data.F32, "y"); + + graphdata.color = KapaColorByName ("black"); + graphdata.style = 0; + graphdata.ptype = 0; + graphdata.size = 1.0; + KapaPrepPlot (kapa, rx->n, &graphdata); + KapaPlotVector (kapa, rx->n, rx->data.F32, "x"); + KapaPlotVector (kapa, rx->n, ry->data.F32, "y"); // pause and wait for user input: @@ -306,100 +404,3 @@ } -bool psphotPetrosianVisualEllipticalContour (pmPetrosian *petrosian) { - - Graphdata graphdata; - - if (!pmVisualIsVisual()) return true; - - if (kapa == -1) { - kapa = KapaOpenNamedSocket ("kapa", "psphot:plots"); - if (kapa == -1) { - fprintf (stderr, "failure to open kapa; visual mode disabled\n"); - pmVisualSetVisual(false); - return false; - } - } - - KapaClearPlots (kapa); - KapaInitGraph (&graphdata); - KapaSetFont (kapa, "courier", 14); - - psVector *theta = petrosian->theta; - psVector *radius = petrosian->isophotalRadii; - - // find Rmin and Rmax for the initial guess - float Rmin = radius->data.F32[0]; - float Rmax = radius->data.F32[0]; - - psVector *Rx = psVectorAlloc(radius->n, PS_TYPE_F32); - psVector *Ry = psVectorAlloc(radius->n, PS_TYPE_F32); - - for (int i = 0; i < theta->n; i++) { - Rx->data.F32[i] = radius->data.F32[i]*cos(theta->data.F32[i]); - Ry->data.F32[i] = radius->data.F32[i]*sin(theta->data.F32[i]); - - // check the radius range - Rmin = MIN (Rmin, radius->data.F32[i]); - Rmax = MAX (Rmax, radius->data.F32[i]); - } - - psVector *rx = psVectorAlloc(361, PS_TYPE_F32); - psVector *ry = psVectorAlloc(361, PS_TYPE_F32); - - float epsilon = petrosian->axes.minor / petrosian->axes.major; - - for (int i = 0; i < 361; i++) { - - float alpha = PS_RAD_DEG * i; - - float cs_alpha = cos(alpha); - float sn_alpha = sin(alpha); - - float cs_phi = cos(alpha - petrosian->axes.theta); - float sn_phi = sin(alpha - petrosian->axes.theta); - - float r = 1.0 / sqrt(SQ(sn_phi) + SQ(epsilon*cs_phi)); - - // generate the model fit here - rx->data.F32[i] = petrosian->axes.minor * cs_alpha * r; - ry->data.F32[i] = petrosian->axes.minor * sn_alpha * r; - } - - graphdata.color = KapaColorByName ("black"); - graphdata.xmin = -1.1*Rmax; - graphdata.ymin = -1.1*Rmax; - graphdata.xmax = +1.1*Rmax; - graphdata.ymax = +1.1*Rmax; - KapaSetLimits (kapa, &graphdata); - - KapaBox (kapa, &graphdata); - KapaSendLabel (kapa, "R_x", KAPA_LABEL_XM); - KapaSendLabel (kapa, "R_y", KAPA_LABEL_YM); - - graphdata.color = KapaColorByName ("red"); - graphdata.style = 2; - graphdata.ptype = 2; - graphdata.size = 1.0; - KapaPrepPlot (kapa, Rx->n, &graphdata); - KapaPlotVector (kapa, Rx->n, Rx->data.F32, "x"); - KapaPlotVector (kapa, Rx->n, Ry->data.F32, "y"); - - graphdata.color = KapaColorByName ("black"); - graphdata.style = 0; - graphdata.ptype = 0; - graphdata.size = 1.0; - KapaPrepPlot (kapa, rx->n, &graphdata); - KapaPlotVector (kapa, rx->n, rx->data.F32, "x"); - KapaPlotVector (kapa, rx->n, ry->data.F32, "y"); - - // 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; -} - # endif Index: branches/pap/psphot/src/psphotRadialBins.c =================================================================== --- branches/pap/psphot/src/psphotRadialBins.c (revision 28003) +++ branches/pap/psphot/src/psphotRadialBins.c (revision 28003) @@ -0,0 +1,204 @@ +# include "psphotInternal.h" +float InterpolateValues (float X0, float Y0, float X1, float Y1, float X); + +// convert the flux vs elliptical radius to annular bins + +// we are guaranteed to be limited by either the seeing (1 - few pixels) or by the pixels +// themselves. this function does not attempt to measure the radial profiles accurately +// for radii that are smaller than ~2 pixels + +// for small radii, we are measuring the mean surface brightness in non-overlapping radial +// bins. for large radii (r > 2 pixels), we are measuring the surface brightness for a +// radius range \alpha r_i < i < \beta r_i, but performing this measurement for radii more +// finely spaced than r_{i+1} = r_i * \beta / \alpha. for the integration, we need to +// track the non-overlapping radius values. + +// Photo interpolates the image of interest to place the peak on the center of the central +// pixel, and then uses the exact fractions of the pixels in each of the first few annuli. +// Seems like a reasonable thing, but is there any significance to the difference? + +// XXX move the resulting elements from profile to extpars->petrosian? +bool psphotRadialBins (psMetadata *recipe, pmSource *source, float radiusMax, float skynoise) { + + psAssert (source, "missing source"); + psAssert (source->extpars, "missing extpars"); + psAssert (source->extpars->ellipticalFlux, "missing ellipticalFlux"); + + psVector *radius = source->extpars->ellipticalFlux->radiusElliptical; + psVector *flux = source->extpars->ellipticalFlux->fluxElliptical; + + // sort incoming vectors by radius + pmSourceRadialProfileSortPair (radius, flux); + + if (!source->extpars->radProfile) { + source->extpars->radProfile = pmSourceRadialProfileAlloc(); + } + pmSourceRadialProfile *profile = source->extpars->radProfile; + + float skyModelErrorSQ = PS_SQR(skynoise); + psEllipseAxes axes = source->extpars->axes; + float AxialRatio = axes.minor / axes.major; + + // radMin, radMax store the bounds of the annuli + bool status = false; + psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); + psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); + if (!radMin || !radMin->n) return false; + if (!radMax || !radMax->n) return false; + + psVector *binSB = psVectorAllocEmpty(radMin->n, PS_TYPE_F32); // surface brightness of radial bin + psVector *binSBstdev = psVectorAllocEmpty(radMin->n, PS_TYPE_F32); // surface brightness error of radial bin + psVector *binSum = psVectorAllocEmpty(radMin->n, PS_TYPE_F32); // surface brightness of radial bin + psVector *binFill = psVectorAllocEmpty(radMin->n, PS_TYPE_F32); // surface brightness of radial bin + psVector *binRad = psVectorAllocEmpty(radMin->n, PS_TYPE_F32); // mean radius of radial bin + psVector *binArea = psVectorAllocEmpty(radMin->n, PS_TYPE_F32); // area of radial bin (contiguous, non-overlapping) + + psVectorInit (binSB, 0.0); + psVectorInit (binSBstdev, 0.0); + psVectorInit (binSum, 0.0); + psVectorInit (binFill, 0.0); + + psVectorInit (binRad, 0.0); + + // generate radial area-weighted mean radius & non-overlapping areas + for (int i = 0; i < radMin->n; i++) { + float rMin2 = PS_SQR(radMin->data.F32[i]); + float rMax2 = PS_SQR(radMax->data.F32[i]); + + float rMin3 = rMin2*radMin->data.F32[i]; + float rMax3 = rMax2*radMax->data.F32[i]; + + float rBin = 2.0 * (rMax3 - rMin3) / (rMax2 - rMin2) / 3.0; + + // XXX calculate area-weighted radius rather than asserting? + binRad->data.F32[i] = rBin; + binArea->data.F32[i] = M_PI * (rMax2 - rMin2); + } + + // storage vector for stats + psVector *values = psVectorAllocEmpty (flux->n, PS_TYPE_F32); + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); + + float fluxSum = 0.0; + int nPixSum = 0; + + bool done = false; + int nOut = 0; + float Rmin = radMin->data.F32[nOut]; + float Rmax = radMax->data.F32[nOut]; + float Rnxt = radMin->data.F32[nOut+1]; // minimum radius for next range + int iNext = 0; + int iLast = -1; + for (int i = 0; !done && (i < radius->n); i++) { + if (radius->data.F32[i] < Rnxt) { + iNext = i; // on the next pass, we will start back here to ensure we hit all pixels in the next bin + } + if (radius->data.F32[i] > Rmax) { + // calculate the surface brightness for bin 'nOut' + float value, dvalue; + if (values->n > 0) { + psVectorStats (stats, values, NULL, NULL, 0); + value = stats->robustMedian; + dvalue = stats->robustStdev; + } else { + value = NAN; + dvalue = NAN; + } + + binSB->data.F32[nOut] = value; + binSBstdev->data.F32[nOut] = sqrt(PS_SQR(dvalue) / values->n + skyModelErrorSQ); + + // calculate the total flux for bin 'nOut' + float Area = M_PI*AxialRatio*PS_SQR(Rmax); + binSum->data.F32[nOut] = fluxSum; + binFill->data.F32[nOut] = nPixSum / Area; + + psTrace ("psphot", 5, "radial bins: %3d %5.1f %5.1f : %7.1f %6.2f : %8.1f %4.2f %6.1f\n", + nOut, radMin->data.F32[nOut], radMax->data.F32[nOut], + binSB->data.F32[nOut], binSBstdev->data.F32[nOut], + binSum->data.F32[nOut], binFill->data.F32[nOut], Area); + + nOut ++; + if (nOut >= radMin->n) break; + Rmin = radMin->data.F32[nOut]; + Rmax = radMax->data.F32[nOut]; + Rnxt = (nOut < radMin->n - 1) ? radMin->data.F32[nOut+1] : Rmax; // minimum radius for next range + values->n = 0; + psStatsInit(stats); + iLast = i; + i = iNext; + } + if (radius->data.F32[i] < Rmin) { + continue; + } + psVectorAppend (values, flux->data.F32[i]); + + if (i > iLast) { + fluxSum += flux->data.F32[i]; + nPixSum ++; + } + } + binSB->n = binSBstdev->n = binRad->n = binArea->n = nOut; + + // interpolate any bins that were empty (extrapolate to center if needed) + if (!isfinite(binSB->data.F32[0]) && !isfinite(binSB->data.F32[1])) { + psWarning ("center 2 bins of source at %f, %f are NAN, skipping this source", source->peak->xf, source->peak->yf); + // XXX raise a flag + psFree(binSB); + psFree(binSBstdev); + + psFree(binSum); + psFree(binFill); + + psFree(binRad); + psFree(binArea); + psFree(values); + psFree(stats); + return false; + } + + // if center bin is empty assume same SB as next radius (probably true due to PSF) + if (!isfinite(binSB->data.F32[0])) { + binSB->data.F32[0] = binSB->data.F32[1]; + binSBstdev->data.F32[0] = binSBstdev->data.F32[1]; + } + + // interpolate any bins that were empty (if center if needed) + for (int i = 1; i < binSB->n - 1; i++) { + if (isfinite(binSB->data.F32[i])) continue; + binSB->data.F32[i] = InterpolateValues (binRad->data.F32[i-1], binSB->data.F32[i-1], binRad->data.F32[i+1], binSB->data.F32[i+1], binRad->data.F32[i]); + binSBstdev->data.F32[i] = InterpolateValues (binRad->data.F32[i-1], binSBstdev->data.F32[i-1], binRad->data.F32[i+1], binSBstdev->data.F32[i+1], binRad->data.F32[i]); + } + + psFree(profile->binSB); + psFree(profile->binSBstdev); + + psFree(profile->binSum); + psFree(profile->binFill); + + psFree(profile->radialBins); + psFree(profile->area); + + // save the vectors + profile->binSB = binSB; + profile->binSBstdev = binSBstdev; + + profile->binSum = binSum; + profile->binFill = binFill; + + profile->radialBins = binRad; + profile->area = binArea; + + psFree(values); + psFree(stats); + + return true; +} + +// the area-weighted mean radius is given by: + +// integral r * 2 pi r dr / integral 2 pi r dr + +// = 2/3 pi (r_max^3 - r_min^3) / pi (r_max^2 - r_min^2) +// = 2/3 (r_max^3 - r_min^3) / (r_max^2 - r_min^2) + Index: branches/pap/psphot/src/psphotRadialProfile.c =================================================================== --- branches/pap/psphot/src/psphotRadialProfile.c (revision 27726) +++ branches/pap/psphot/src/psphotRadialProfile.c (revision 28003) @@ -3,11 +3,9 @@ bool psphotRadialProfile (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal) { + bool status; + // allocate pmSourceExtendedParameters, if not already defined if (!source->extpars) { source->extpars = pmSourceExtendedParsAlloc (); - } - - if (!source->extpars->profile) { - source->extpars->profile = pmSourceRadialProfileAlloc (); } @@ -17,4 +15,6 @@ float fluxMin = 0.0; float fluxMax = source->peak->flux; + + bool RAW_RADIUS = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_RAW_RADIUS"); // generate a series of radial profiles at Nsec evenly spaced angles. the profile flux @@ -28,5 +28,5 @@ // use the radial profiles to determine the radius of a given isophote. this isophote // is used to determine the elliptical shape of the object, so it has a relatively high - // value (nominally 50% of the peak) + // value (nominally 25% of the peak) if (!psphotRadiiFromProfiles (source, fluxMin, fluxMax)) { psError (PS_ERR_UNKNOWN, false, "failed to measure isophotal radii from profiles"); @@ -42,6 +42,12 @@ // 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 - if (!psphotEllipticalProfile (source)) { + if (!psphotEllipticalProfile (source, RAW_RADIUS)) { psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile"); + return false; + } + + // generated profile in averaged bins + if (!psphotRadialBins (recipe, source, Rmax, skynoise)) { + psError (PS_ERR_UNKNOWN, false, "failed to generate radial bins"); return false; } Index: branches/pap/psphot/src/psphotRadialProfileByAngles.c =================================================================== --- branches/pap/psphot/src/psphotRadialProfileByAngles.c (revision 27726) +++ branches/pap/psphot/src/psphotRadialProfileByAngles.c (revision 28003) @@ -12,9 +12,14 @@ bool psphotRadialProfilesByAngles (pmSource *source, int Nsec, float Rmax) { + psAssert (source->extpars, "define extpars"); + // we want to have an even number of sectors so we can do 180 deg symmetrizing Nsec = (Nsec % 2) ? Nsec + 1 : Nsec; float dtheta = 2.0*M_PI / Nsec; - pmSourceRadialProfile *profile = source->extpars->profile; + if (!source->extpars->radFlux) { + source->extpars->radFlux = pmSourceRadialFluxAlloc(); + } + pmSourceRadialFlux *profile = source->extpars->radFlux; psFree(profile->radii); psFree(profile->fluxes); @@ -24,5 +29,4 @@ profile->fluxes = psArrayAllocEmpty(Nsec); profile->theta = psVectorAllocEmpty(Nsec, PS_TYPE_F32); - for (int i = 0; i < Nsec; i++) { Index: branches/pap/psphot/src/psphotRadiiFromProfiles.c =================================================================== --- branches/pap/psphot/src/psphotRadiiFromProfiles.c (revision 27726) +++ branches/pap/psphot/src/psphotRadiiFromProfiles.c (revision 28003) @@ -5,5 +5,9 @@ bool psphotRadiiFromProfiles (pmSource *source, float fluxMin, float fluxMax) { - pmSourceRadialProfile *profile = source->extpars->profile; + psAssert (source, "missing source"); + psAssert (source->extpars, "missing extpars"); + psAssert (source->extpars->radFlux, "missing radFlux"); + + pmSourceRadialFlux *profile = source->extpars->radFlux; psFree(profile->isophotalRadii); @@ -146,5 +150,5 @@ // show the results - // psphotPetrosianVisualProfileRadii (radius, flux, radiusBinned, fluxBinned, Ro); + // psphotPetrosianVisualProfileRadii (radius, flux, radiusBinned, fluxBinned, fluxMax, Ro); psFree(fluxBinned); @@ -152,2 +156,5 @@ return Ro; } + + + Index: branches/pap/psphot/src/psphotStackImageLoop.c =================================================================== --- branches/pap/psphot/src/psphotStackImageLoop.c (revision 27726) +++ branches/pap/psphot/src/psphotStackImageLoop.c (revision 28003) @@ -41,7 +41,16 @@ if (! readout->data_exists) { continue; } +# if (0) + // uncomment to generate matched psfs + if (!psphotStackMatchPSFs (config, view)) { + psError(psErrorCodeLast(), false, "failure in psphotStackMatchPSFs for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout); + psFree (view); + return false; + } +# endif + // XXX for now, we assume there is only a single chip in the PHU: if (!psphotStackReadout (config, view)) { - psError(psErrorCodeLast(), false, "failure in psphotReadout for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout); + psError(psErrorCodeLast(), false, "failure in psphotStackReadout for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout); psFree (view); return false; Index: branches/pap/psphot/src/psphotStackMatchPSFs.c =================================================================== --- branches/pap/psphot/src/psphotStackMatchPSFs.c (revision 28003) +++ branches/pap/psphot/src/psphotStackMatchPSFs.c (revision 28003) @@ -0,0 +1,152 @@ +# include "psphotInternal.h" + +bool psphotStackMatchPSFs (pmConfig *config, const pmFPAview *view, bool firstPass) +{ + bool status = true; + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotStackMatchPSFsReadout (config, view, i)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to find initial detections for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +} + +// convolve the image to match desired PSF +bool psphotStackMatchPSFsReadout (pmConfig *config, const pmFPAview *view, int index) { + + bool status; + int pass; + float NSIGMA_PEAK = 25.0; + int NMAX = 0; + + // find the currently selected readout + pmFPAfile *fileRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", index); // File of interest + psAssert (fileRaw, "missing file?"); + + pmFPAfile *fileCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT.CONV", index); // File of interest + psAssert (fileCnv, "missing file?"); + + pmReadout *readoutRaw = pmFPAviewThisReadout(view, fileRaw->fpa); + psAssert (readoutRaw, "missing readout?"); + + pmReadout *readoutCnv = pmFPAviewThisReadout(view, fileCnv->fpa); + psAssert (readoutCnv, "missing readout?"); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + psAssert (maskVal, "missing mask value?"); + + /***** set up recipe options *****/ + + psMetadata *stackRecipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe + psAssert(stackRecipe, "We've thrown an error on this before."); + + // Look up appropriate values from the ppSub recipe + psMetadata *subRecipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe + psAssert(subRecipe, "recipe missing"); + + int size = psMetadataLookupS32(NULL, subRecipe, "KERNEL.SIZE"); // Kernel half-size + + float deconvLimit = psMetadataLookupF32(NULL, stackRecipe, "DECONV.LIMIT"); // Limit on deconvolution fraction + + psString maskValStr = psMetadataLookupStr(NULL, subRecipe, "MASK.VAL"); // Name of bits to mask going in + psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch + psString maskPoorStr = psMetadataLookupStr(NULL, stackRecipe, "MASK.POOR"); // Name of bits to mask for poor + psImageMaskType maskPoor = pmConfigMaskGet(maskPoorStr, config); // Bits to mask for poor pixels + psString maskBadStr = psMetadataLookupStr(NULL, stackRecipe, "MASK.BAD"); // Name of bits to mask for bad + psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels + + bool mdok; // Status of MD lookup + float penalty = psMetadataLookupF32(NULL, subRecipe, "PENALTY"); // Penalty for wideness + int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads + + if (!pmReadoutMaskNonfinite(readout, maskVal)) { + psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels in readout."); + return false; + } + + // Image Matching (PSFs or just flux) + if (options->convolve) { + // Full match of PSFs + pmReadout *conv = pmReadoutAlloc(NULL); // Conv readout, for holding results temporarily + + // Read previously produced kernel + if (psMetadataLookupBool(NULL, config->arguments, "PPSTACK.DEBUG.STACK")) { + loadKernel(); + } else { + matchKernel(); + } // !DEBUG.STACK + + saveMatchData(); + + saveChiSquare(); + + renormKernel(); + + // Reject image completely if the maximum deconvolution fraction exceeds the limit + float deconv = psMetadataLookupF32(NULL, conv->analysis, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Max deconvolution fraction + if (deconv > deconvLimit) { + psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting image %d\n", deconv, deconvLimit, index); + psFree(conv); + return NULL; + } + + readout->analysis = psMetadataCopy(readout->analysis, conv->analysis); + + psFree(conv); + } else { + // only match the flux + float norm = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation + psBinaryOp(readout->image, readout->image, "*", psScalarAlloc(norm, PS_TYPE_F32)); + psBinaryOp(readout->variance, readout->variance, "*", psScalarAlloc(PS_SQR(norm), PS_TYPE_F32)); + } + + // Ensure the background value is zero + psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for background + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator + if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) { + psWarning("Can't measure background for image."); + psErrorClear(); + } else { + if (!psMetadataLookupBool(NULL, config->arguments, "PPSTACK.SKIP.BG.SUB")) { + psLogMsg("ppStack", PS_LOG_INFO, "Correcting convolved image background by %lf (+/- %lf)", + psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), psStatsGetValue(bg, PS_STAT_ROBUST_STDEV)); + (void)psBinaryOp(readout->image, readout->image, "-", + psScalarAlloc(psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), PS_TYPE_F32)); + } + } + + if (!stackRenormaliseReadout(config, readout)) { + psFree(rng); + psFree(bg); + return false; + } + + // Measure the variance level for the weighting + if (psMetadataLookupBool(NULL, stackRecipe, "WEIGHTS")) { + if (!psImageBackground(bg, NULL, readout->variance, readout->mask, maskVal | maskBad, rng)) { + psError(PPSTACK_ERR_DATA, false, "Can't measure mean variance for image."); + psFree(rng); + psFree(bg); + return false; + } + options->weightings->data.F32[index] = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) * psImageCovarianceFactor(readout->covariance)); + } else { + options->weightings->data.F32[index] = 1.0; + } + psLogMsg("ppStack", PS_LOG_INFO, "Weighting for image %d is %f\n", + index, options->weightings->data.F32[index]); + + psFree(rng); + psFree(bg); + + dumpImage3(); + + return true; +} Index: branches/pap/psphot/src/psphotStackMatchPSFsUtils.c =================================================================== --- branches/pap/psphot/src/psphotStackMatchPSFsUtils.c (revision 28003) +++ branches/pap/psphot/src/psphotStackMatchPSFsUtils.c (revision 28003) @@ -0,0 +1,556 @@ +/***** defines *****/ + +#define ARRAY_BUFFER 16 // Number to add to array at a time +#define MAG_IGNORE 50 // Ignore magnitudes fainter than this --- they're not real! +#define FAKE_SIZE 1 // Size of fake convolution kernel +#define SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to input sources +#define NOISE_FRACTION 0.01 // Set minimum flux to this fraction of noise +#define COVAR_FRAC 0.01 // Truncation fraction for covariance matrix + +// XXX better name +bool readImage(psImage **target, // Target for image + const char *name, // Name of FITS file + const pmConfig *config // Configuration + ) +{ + psString resolved = pmConfigConvertFilename(name, config, false, false); // Resolved filename + psFits *fits = psFitsOpen(resolved, "r"); + psFree(resolved); + if (!fits) { + psError(PPSTACK_ERR_IO, false, "Unable to open previously produced image: %s", name); + return false; + } + psImage *image = psFitsReadImage(fits, psRegionSet(0,0,0,0), 0); // Image of interest + if (!image) { + psError(PPSTACK_ERR_IO, false, "Unable to read previously produced image: %s", name); + psFitsClose(fits); + return false; + } + psFitsClose(fits); + + psFree(*target); + *target = image; + + return true; +} + +// Get coordinates from a source +void coordsFromSource(float *x, float *y, const pmSource *source) +{ + assert(x && y); + assert(source); + + if (source->modelPSF) { + *x = source->modelPSF->params->data.F32[PM_PAR_XPOS]; + *y = source->modelPSF->params->data.F32[PM_PAR_YPOS]; + } else { + *x = source->peak->xf; + *y = source->peak->yf; + } + return; +} + +psArray *stackSourcesFilter(psArray *sources, // Source list to filter + int exclusion // Exclusion zone, pixels + ) +{ + psAssert(sources && sources->n > 0, "Require array of sources"); + if (exclusion <= 0) { + return psMemIncrRefCounter(sources); + } + + int num = sources->n; // Number of sources + psVector *x = psVectorAlloc(num, PS_TYPE_F32), *y = psVectorAlloc(num, PS_TYPE_F32); // Coordinates + int numGood = 0; // Number of good sources + for (int i = 0; i < num; i++) { + pmSource *source = sources->data[i]; // Source of interest + if (!source) { + continue; + } + coordsFromSource(&x->data.F32[numGood], &y->data.F32[numGood], source); + numGood++; + } + x->n = y->n = numGood; + + psTree *tree = psTreePlant(2, 2, PS_TREE_EUCLIDEAN, x, y); // kd tree + + psArray *filtered = psArrayAllocEmpty(numGood); // Filtered list of sources + psVector *coords = psVectorAlloc(2, PS_TYPE_F64); // Coordinates of source + int numFiltered = 0; // Number of filtered sources + for (int i = 0; i < num; i++) { + pmSource *source = sources->data[i]; // Source of interest + if (!source) { + continue; + } + float xSource, ySource; // Coordinates of source + coordsFromSource(&xSource, &ySource, source); + + coords->data.F64[0] = xSource; + coords->data.F64[1] = ySource; + + long numWithin = psTreeWithin(tree, coords, exclusion); // Number within exclusion zone + psTrace("ppStack", 9, "Source at %.0lf,%.0lf has %ld sources in exclusion zone", + coords->data.F64[0], coords->data.F64[1], numWithin); + if (numWithin == 1) { + // Only itself inside the exclusion zone + filtered = psArrayAdd(filtered, filtered->n, source); + } else { + numFiltered++; + } + } + psFree(coords); + psFree(tree); + psFree(x); + psFree(y); + + psLogMsg("ppStack", PS_LOG_INFO, "Filtered out %d of %d sources", numFiltered, numGood); + + return filtered; +} + +// Add background into the fake image +// Based on ppSubBackground() +static psImage *stackBackgroundModel(pmReadout *ro, // Readout for which to generate background model + const pmConfig *config // Configuration + ) +{ + psAssert(ro && ro->image, "Need readout image"); + psAssert(config, "Need configuration"); + + psImage *image = ro->image; // Image of interest + int numCols = image->numCols, numRows = image->numRows; // Size of image + + psMetadata *ppStackRecipe = psMetadataLookupPtr(NULL, config->recipes, PPSTACK_RECIPE); + psAssert(ppStackRecipe, "Need PPSTACK recipe"); + psMetadata *psphotRecipe = psMetadataLookupPtr(NULL, config->recipes, PSPHOT_RECIPE); + psAssert(psphotRecipe, "Need PSPHOT recipe"); + + psString maskBadStr = psMetadataLookupStr(NULL, ppStackRecipe, "MASK.BAD");// Name of bits to mask for bad + psMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psMetadataAddU8(psphotRecipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "user-defined mask", maskBad); + + psImage *binned = psphotModelBackgroundReadoutNoFile(ro, config); // Binned background model + psImageBinning *binning = psMetadataLookupPtr(NULL, ro->analysis, + "PSPHOT.BACKGROUND.BINNING"); // Binning for model + psAssert(binning, "Need binning parameters"); + psImage *unbinned = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Unbinned background model + if (!psImageUnbin(unbinned, binned, binning)) { + psError(PPSTACK_ERR_DATA, false, "Unable to unbin background model"); + psFree(binned); + psFree(unbinned); + return NULL; + } + psFree(binned); + + return unbinned; +} + +// Renormalise a readout's variance map +bool stackRenormaliseReadout(const pmConfig *config, // Configuration + pmReadout *readout // Readout to renormalise + ) +{ +#if 1 + bool mdok; // Status of metadata lookups + + psMetadata *recipe = psMetadataLookupPtr(NULL, config->recipes, PPSTACK_RECIPE); // Recipe for ppStack + psAssert(recipe, "Need PPSTACK recipe"); + + if (!psMetadataLookupBool(&mdok, recipe, "RENORM")) return true; + + int num = psMetadataLookupS32(&mdok, recipe, "RENORM.NUM"); + if (!mdok) { + psError(PPSTACK_ERR_CONFIG, true, "RENORM.NUM is not set in the recipe"); + return false; + } + float minValid = psMetadataLookupF32(&mdok, recipe, "RENORM.MIN"); + if (!mdok) { + psError(PPSTACK_ERR_CONFIG, true, "RENORM.MIN is not set in the recipe"); + return false; + } + float maxValid = psMetadataLookupF32(&mdok, recipe, "RENORM.MAX"); + if (!mdok) { + psError(PPSTACK_ERR_CONFIG, true, "RENORM.MAX is not set in the recipe"); + return false; + } + + psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask + + psImageCovarianceTransfer(readout->variance, readout->covariance); + return pmReadoutVarianceRenormalise(readout, maskBad, num, minValid, maxValid); +#else + return true; +#endif +} + +// This is a hack to use the temporary convolved images and kernel generated previously. +// This makes the 'matching' operation much faster, allowing debugging of the stack process easier. +// It implicitly assumes the output root name is the same between invocations. + +bool loadKernel () { + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.CONV.KERNEL", index); + psAssert(file, "Require file"); + + pmFPAview *view = pmFPAviewAlloc(0); // View to readout of interest + view->chip = view->cell = view->readout = 0; + psString filename = pmFPAfileNameFromRule(file->filerule, file, view); // Filename of interest + + // Read convolution kernel + psString resolved = pmConfigConvertFilename(filename, config, false, false); // Resolved filename + psFree(filename); + psFits *fits = psFitsOpen(resolved, "r"); // FITS file for subtraction kernel + psFree(resolved); + if (!fits || !pmReadoutReadSubtractionKernels(conv, fits)) { + psError(PPSTACK_ERR_IO, false, "Unable to read previously produced kernel"); + psFitsClose(fits); + return false; + } + psFitsClose(fits); + + if (!readImage(&readout->image, options->convImages->data[index], config) || + !readImage(&readout->mask, options->convMasks->data[index], config) || + !readImage(&readout->variance, options->convVariances->data[index], config)) { + psError(PPSTACK_ERR_IO, false, "Unable to read previously produced image."); + return false; + } + + psRegion *region = psMetadataLookupPtr(NULL, conv->analysis, + PM_SUBTRACTION_ANALYSIS_REGION); // Convolution region + pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, conv->analysis, + PM_SUBTRACTION_ANALYSIS_KERNEL); + + pmSubtractionAnalysis(conv->analysis, NULL, kernels, region, + readout->image->numCols, readout->image->numRows); + + psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel + bool oldThreads = psImageCovarianceSetThreads(true); // Old thread setting + psKernel *covar = psImageCovarianceCalculate(kernel, readout->covariance); // Covariance matrix + psImageCovarianceSetThreads(oldThreads); + psFree(readout->covariance); + readout->covariance = covar; + psFree(kernel); +} + +bool dumpImage() { + // XXX should be optional + { + pmHDU *hdu = pmHDUFromCell(readout->parent); + psString name = NULL; + psStringAppend(&name, "fake_%03d.fits", index); + pmStackVisualPlotTestImage(fake->image, name); + psFits *fits = psFitsOpen(name, "w"); + psFree(name); + psFitsWriteImage(fits, hdu->header, fake->image, 0, NULL); + psFitsClose(fits); + } + { + pmHDU *hdu = pmHDUFromCell(readout->parent); + psString name = NULL; + psStringAppend(&name, "real_%03d.fits", index); + pmStackVisualPlotTestImage(readout->image, name); + psFits *fits = psFitsOpen(name, "w"); + psFree(name); + psFitsWriteImage(fits, hdu->header, readout->image, 0, NULL); + psFitsClose(fits); + } +} + +bool dumpImage2() { + // XXX should be optional + + { + pmHDU *hdu = pmHDUFromCell(readout->parent); + psString name = NULL; + psStringAppend(&name, "conv_%03d.fits", index); + pmStackVisualPlotTestImage(conv->image, name); + psFits *fits = psFitsOpen(name, "w"); + psFree(name); + psFitsWriteImage(fits, hdu->header, conv->image, 0, NULL); + psFitsClose(fits); + } + { + pmHDU *hdu = pmHDUFromCell(readout->parent); + psString name = NULL; + psStringAppend(&name, "diff_%03d.fits", index); + pmStackVisualPlotTestImage(fake->image, name); + psFits *fits = psFitsOpen(name, "w"); + psFree(name); + psBinaryOp(fake->image, conv->image, "-", fake->image); + psFitsWriteImage(fits, hdu->header, fake->image, 0, NULL); + psFitsClose(fits); + } +} + +bool dumpImage3() + { + pmHDU *hdu = pmHDUFromCell(readout->parent); + psString name = NULL; + psStringAppend(&name, "convolved_%03d.fits", index); + pmStackVisualPlotTestImage(readout->image, name); + psFits *fits = psFitsOpen(name, "w"); + psFree(name); + psFitsWriteImage(fits, hdu->header, readout->image, 0, NULL); + psFitsClose(fits); + } + +bool matchKernel() { + // Normal operations here + psAssert(options->psf, "Require target PSF"); + psAssert(options->sourceLists && options->sourceLists->data[index], "Require source list"); + + int order = psMetadataLookupS32(NULL, subRecipe, "SPATIAL.ORDER"); // Spatial polynomial order + float regionSize = psMetadataLookupF32(NULL, subRecipe, "REGION.SIZE"); // Size of iso-kernel regs + float spacing = psMetadataLookupF32(NULL, subRecipe, "STAMP.SPACING"); // Typical stamp spacing + int footprint = psMetadataLookupS32(NULL, subRecipe, "STAMP.FOOTPRINT"); // Stamp half-size + float threshold = psMetadataLookupF32(NULL, subRecipe, "STAMP.THRESHOLD"); // Threshold for stmps + int stride = psMetadataLookupS32(NULL, subRecipe, "STRIDE"); // Size of convolution patches + int iter = psMetadataLookupS32(NULL, subRecipe, "ITER"); // Rejection iterations + float rej = psMetadataLookupF32(NULL, subRecipe, "REJ"); // Rejection threshold + float kernelError = psMetadataLookupF32(NULL, subRecipe, "KERNEL.ERR"); // Relative systematic error in kernel + float normFrac = psMetadataLookupF32(NULL, subRecipe, "NORM.FRAC"); // Fraction of window for normalisn windw + float sysError = psMetadataLookupF32(NULL, subRecipe, "SYS.ERR"); // Relative systematic error in images + float skyErr = psMetadataLookupF32(NULL, subRecipe, "SKY.ERR"); // Additional error in sky + float covarFrac = psMetadataLookupF32(NULL, subRecipe, "COVAR.FRAC"); // Fraction for covariance calculation + + const char *typeStr = psMetadataLookupStr(NULL, subRecipe, "KERNEL.TYPE"); // Kernel type + pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(typeStr); // Kernel type + psVector *widths = psMetadataLookupPtr(NULL, subRecipe, "ISIS.WIDTHS"); // ISIS Gaussian widths + psVector *orders = psMetadataLookupPtr(NULL, subRecipe, "ISIS.ORDERS"); // ISIS Polynomial orders + int inner = psMetadataLookupS32(NULL, subRecipe, "INNER"); // Inner radius + int ringsOrder = psMetadataLookupS32(NULL, subRecipe, "RINGS.ORDER"); // RINGS polynomial order + int binning = psMetadataLookupS32(NULL, subRecipe, "SPAM.BINNING"); // Binning for SPAM kernel + float badFrac = psMetadataLookupF32(NULL, subRecipe, "BADFRAC"); // Maximum bad fraction + bool optimum = psMetadataLookupBool(&mdok, subRecipe, "OPTIMUM"); // Derive optimum parameters? + float optMin = psMetadataLookupF32(&mdok, subRecipe, "OPTIMUM.MIN"); // Minimum width for search + float optMax = psMetadataLookupF32(&mdok, subRecipe, "OPTIMUM.MAX"); // Maximum width for search + float optStep = psMetadataLookupF32(&mdok, subRecipe, "OPTIMUM.STEP"); // Step for search + float optThresh = psMetadataLookupF32(&mdok, subRecipe, "OPTIMUM.TOL"); // Tolerance for search + int optOrder = psMetadataLookupS32(&mdok, subRecipe, "OPTIMUM.ORDER"); // Order for search + float poorFrac = psMetadataLookupF32(&mdok, subRecipe, "POOR.FRACTION"); // Fraction for "poor" + + bool scale = psMetadataLookupBool(NULL, subRecipe, "SCALE"); // Scale kernel parameters? + float scaleRef = psMetadataLookupF32(NULL, subRecipe, "SCALE.REF"); // Reference for scaling + float scaleMin = psMetadataLookupF32(NULL, subRecipe, "SCALE.MIN"); // Minimum for scaling + float scaleMax = psMetadataLookupF32(NULL, subRecipe, "SCALE.MAX"); // Maximum for scaling + if (!isfinite(scaleRef) || !isfinite(scaleMin) || !isfinite(scaleMax)) { + psError(PPSTACK_ERR_CONFIG, false, + "Scale parameters (SCALE.REF=%f, SCALE.MIN=%f, SCALE.MAX=%f) not set in PPSUB recipe.", + scaleRef, scaleMin, scaleMax); + return false; + } + + + // These values are specified specifically for stacking + const char *stampsName = psMetadataLookupStr(NULL, config->arguments, "STAMPS");// Stamps filename + + psVector *optWidths = NULL; // Vector with FWHMs for optimum search + if (optimum) { + optWidths = psVectorCreate(optWidths, optMin, optMax, optStep, PS_TYPE_F32); + } + + pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout with target PSF + + psStats *bg = psStatsAlloc(PS_STAT_ROBUST_STDEV); // Statistics for background + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator + if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) { + psError(PPSTACK_ERR_DATA, false, "Can't measure background for image."); + psFree(fake); + psFree(optWidths); + psFree(conv); + psFree(bg); + psFree(rng); + return false; + } + float minFlux = NOISE_FRACTION * bg->robustStdev; // Minimum flux level for fake image + psFree(rng); + psFree(bg); + + // For the sake of stamps, remove nearby sources + psArray *stampSources = stackSourcesFilter(options->sourceLists->data[index], + footprint); // Filtered list of sources + + bool oldThreads = pmReadoutFakeThreads(true); // Old threading state + if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows, + stampSources, SOURCE_MASK, NULL, NULL, options->psf, + minFlux, footprint + size, false, true)) { + psError(PPSTACK_ERR_DATA, false, "Unable to generate fake image with target PSF."); + psFree(fake); + psFree(optWidths); + psFree(conv); + return false; + } + pmReadoutFakeThreads(oldThreads); + + fake->mask = psImageCopy(NULL, readout->mask, PS_TYPE_IMAGE_MASK); + + // Add the background into the target image + psImage *bgImage = stackBackgroundModel(readout, config); // Image of background + psBinaryOp(fake->image, fake->image, "+", bgImage); + psFree(bgImage); + + dumpImage(); + + if (threads > 0) { + pmSubtractionThreadsInit(); + } + + // Do the image matching + pmSubtractionKernels *kernel = psMetadataLookupPtr(&mdok, readout->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL); // Conv kernel + if (kernel) { + if (!pmSubtractionMatchPrecalc(NULL, conv, fake, readout, readout->analysis, + stride, kernelError, covarFrac, maskVal, maskBad, maskPoor, + poorFrac, badFrac)) { + psError(psErrorCodeLast(), false, "Unable to convolve images."); + psFree(fake); + psFree(optWidths); + psFree(stampSources); + psFree(conv); + if (threads > 0) { + pmSubtractionThreadsFinalize(); + } + return false; + } + } else { + // Scale the input parameters + psVector *widthsCopy = psVectorCopy(NULL, widths, PS_TYPE_F32); // Copy of kernel widths + if (scale && !pmSubtractionParamsScale(&size, &footprint, widthsCopy, + options->inputSeeing->data.F32[index], + options->targetSeeing, scaleRef, scaleMin, scaleMax)) { + psError(psErrorCodeLast(), false, "Unable to scale kernel parameters"); + psFree(fake); + psFree(optWidths); + psFree(stampSources); + psFree(conv); + psFree(widthsCopy); + if (threads > 0) { + pmSubtractionThreadsFinalize(); + } + return false; + } + + if (!pmSubtractionMatch(NULL, conv, fake, readout, footprint, stride, regionSize, spacing, + threshold, stampSources, stampsName, type, size, order, widthsCopy, + orders, inner, ringsOrder, binning, penalty, + optimum, optWidths, optOrder, optThresh, iter, rej, normFrac, + sysError, skyErr, kernelError, covarFrac, maskVal, maskBad, maskPoor, + poorFrac, badFrac, PM_SUBTRACTION_MODE_2)) { + psError(psErrorCodeLast(), false, "Unable to match images."); + psFree(fake); + psFree(optWidths); + psFree(stampSources); + psFree(conv); + psFree(widthsCopy); + if (threads > 0) { + pmSubtractionThreadsFinalize(); + } + return false; + } + psFree(widthsCopy); + } + + dumpImage2(); + + psFree(fake); + psFree(optWidths); + psFree(stampSources); + + if (threads > 0) { + pmSubtractionThreadsFinalize(); + } + + // Replace original images with convolved + psFree(readout->image); + psFree(readout->mask); + psFree(readout->variance); + psFree(readout->covariance); + readout->image = psMemIncrRefCounter(conv->image); + readout->mask = psMemIncrRefCounter(conv->mask); + readout->variance = psMemIncrRefCounter(conv->variance); + readout->covariance = psImageCovarianceTruncate(conv->covariance, COVAR_FRAC); + +} + +bool saveMatchData () { + // Extract the regions and solutions used in the image matching + // This stops them from being freed when we iterate back up the FPA + psArray *regions = options->regions->data[index] = psArrayAllocEmpty(ARRAY_BUFFER); // Match regions + { + psString regex = NULL; // Regular expression + psStringAppend(®ex, "^%s$", PM_SUBTRACTION_ANALYSIS_REGION); + psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex); + psFree(regex); + psMetadataItem *item = NULL;// Item from iteration + while ((item = psMetadataGetAndIncrement(iter))) { + assert(item->type == PS_DATA_REGION); + regions = psArrayAdd(regions, ARRAY_BUFFER, item->data.V); + } + psFree(iter); + } + psArray *kernels = options->kernels->data[index] = psArrayAllocEmpty(ARRAY_BUFFER); // Match kernels + { + psString regex = NULL; // Regular expression + psStringAppend(®ex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL); + psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex); + psFree(regex); + psMetadataItem *item = NULL;// Item from iteration + while ((item = psMetadataGetAndIncrement(iter))) { + assert(item->type == PS_DATA_UNKNOWN); + pmSubtractionKernels *kernel = item->data.V; // Kernel used in subtraction + kernels = psArrayAdd(kernels, ARRAY_BUFFER, kernel); + } + psFree(iter); + } + psAssert((regions)->n == (kernels)->n, "Number of match regions and kernels should match"); +} + +bool saveChiSquare() { + // Record chi^2 + { + double sum = 0.0; // Sum of chi^2 + int num = 0; // Number of measurements of chi^2 + psString regex = NULL; // Regular expression + psStringAppend(®ex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL); + psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex); + psFree(regex); + psMetadataItem *item = NULL;// Item from iteration + while ((item = psMetadataGetAndIncrement(iter))) { + assert(item->type == PS_DATA_UNKNOWN); + pmSubtractionKernels *kernels = item->data.V; // Convolution kernels + sum += kernels->mean; + num++; + } + psFree(iter); + options->matchChi2->data.F32[index] = sum / (psImageCovarianceFactor(readout->covariance) * num); + } + +} + +bool renormKernel() { + // Kernel normalisation + { + double sum = 0.0; // Sum of chi^2 + int num = 0; // Number of measurements of chi^2 + psString regex = NULL; // Regular expression + psStringAppend(®ex, "^%s$", PM_SUBTRACTION_ANALYSIS_NORM); + psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex); + psFree(regex); + psMetadataItem *item = NULL;// Item from iteration + while ((item = psMetadataGetAndIncrement(iter))) { + assert(item->type == PS_TYPE_F32); + float norm = item->data.F32; // Normalisation + sum += norm; + num++; + } + psFree(iter); + float conv = sum/num; // Mean normalisation from convolution + float stars = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation from stars + float renorm = stars / conv; // Renormalisation to apply + psLogMsg("ppStack", PS_LOG_INFO, "Renormalising image %d by %f (kernel: %f, stars: %f)\n", + index, renorm, conv, stars); + psBinaryOp(readout->image, readout->image, "*", psScalarAlloc(renorm, PS_TYPE_F32)); + psBinaryOp(readout->variance, readout->variance, "*", psScalarAlloc(PS_SQR(renorm), PS_TYPE_F32)); + } + +} Index: branches/pap/psphot/src/psphotVisual.c =================================================================== --- branches/pap/psphot/src/psphotVisual.c (revision 27726) +++ branches/pap/psphot/src/psphotVisual.c (revision 28003) @@ -446,4 +446,6 @@ KapaSetFont (myKapa, "courier", 14); + section.bg = KapaColorByName ("none"); // XXX probably should be 'none' + float SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); @@ -538,5 +540,10 @@ KapaSetLimits (myKapa, &graphdata); + graphdata.padXm = NAN; + graphdata.padYm = NAN; + graphdata.padXp = 0.5; + graphdata.padYp = 0.5; KapaBox (myKapa, &graphdata); + KapaSendLabel (myKapa, "M_xx| (pixels)", KAPA_LABEL_XM); KapaSendLabel (myKapa, "M_yy| (pixels)", KAPA_LABEL_YM); @@ -547,4 +554,5 @@ graphdata.style = 2; KapaPrepPlot (myKapa, nF, &graphdata); + KapaPlotVector (myKapa, nF, xFaint->data.F32, "x"); KapaPlotVector (myKapa, nF, yFaint->data.F32, "y"); @@ -562,5 +570,5 @@ section.dy = 0.25; section.x = 0.00; - section.y = 0.80; + section.y = 0.75; section.name = psStringCopy ("MagMyy"); KapaSetSection (myKapa, §ion); @@ -574,4 +582,8 @@ KapaSetLimits (myKapa, &graphdata); + graphdata.padXm = 0.5; + graphdata.padYm = NAN; + graphdata.padXp = NAN; + graphdata.padYp = 0.5; strcpy (graphdata.labels, "0210"); KapaBox (myKapa, &graphdata); @@ -598,5 +610,5 @@ section.dx = 0.25; section.dy = 0.75; - section.x = 0.80; + section.x = 0.75; section.y = 0.00; section.name = psStringCopy ("MagMxx"); @@ -611,4 +623,8 @@ KapaSetLimits (myKapa, &graphdata); + graphdata.padXm = NAN; + graphdata.padYm = 0.5; + graphdata.padXp = 0.5; + graphdata.padYp = NAN; strcpy (graphdata.labels, "2001"); KapaBox (myKapa, &graphdata); @@ -1350,4 +1366,6 @@ assert (maskVal); + section.bg = KapaColorByName ("none"); // XXX probably should be 'none' + KapaClearPlots (myKapa); // first section : mag vs CR nSigma @@ -1586,4 +1604,6 @@ KapaSetFont (myKapa, "courier", 14); + section.bg = KapaColorByName ("none"); // XXX probably should be 'none' + // select the max psfX,Y values for the plot limits float Xmin = 1000.0, Xmax = 0.0; @@ -1754,4 +1774,8 @@ KapaSetLimits (myKapa, &graphdata); + graphdata.padXm = NAN; + graphdata.padYm = NAN; + graphdata.padXp = 0.5; + graphdata.padYp = 0.5; KapaBox (myKapa, &graphdata); KapaSendLabel (myKapa, "M_xx| (pixels)", KAPA_LABEL_XM); @@ -1822,4 +1846,8 @@ KapaSetLimits (myKapa, &graphdata); + graphdata.padXm = 0.5; + graphdata.padYm = NAN; + graphdata.padXp = NAN; + graphdata.padYp = 0.5; strcpy (graphdata.labels, "0210"); KapaBox (myKapa, &graphdata); @@ -1870,5 +1898,5 @@ section.dx = 0.25; section.dy = 0.60; - section.x = 0.80; + section.x = 0.75; section.y = 0.00; section.name = psStringCopy ("MagMxx"); @@ -1883,4 +1911,8 @@ KapaSetLimits (myKapa, &graphdata); + graphdata.padXm = NAN; + graphdata.padYm = 0.5; + graphdata.padXp = 0.5; + graphdata.padYp = NAN; strcpy (graphdata.labels, "2001"); KapaBox (myKapa, &graphdata); @@ -1930,7 +1962,7 @@ // fourth section: MagSigma section.dx = 0.75; - section.dy = 0.15; + section.dy = 0.20; section.x = 0.00; - section.y = 0.65; + section.y = 0.60; section.name = psStringCopy ("MagSigma"); KapaSetSection (myKapa, §ion); @@ -1944,4 +1976,8 @@ KapaSetLimits (myKapa, &graphdata); + graphdata.padXm = 0.5; + graphdata.padYm = NAN; + graphdata.padXp = 0.5; + graphdata.padYp = 0.5; strcpy (graphdata.labels, "0100"); KapaBox (myKapa, &graphdata); @@ -2053,4 +2089,9 @@ psFree (mDEF); psFree (sDEF); + + psFree (xLOW); + psFree (yLOW); + psFree (mLOW); + psFree (sLOW); psFree (xCR); @@ -2297,29 +2338,28 @@ if (!source) continue; if (!source->extpars) continue; - if (!source->extpars->profile) continue; - if (!source->extpars->petrosian_80) continue; - - pmSourceRadialProfile *profile = source->extpars->profile; - pmSourceExtendedFlux *petrosian = source->extpars->petrosian_80; + if (!source->extpars->petProfile) continue; + + float petrosianRadius = source->extpars->petrosianRadius; + psEllipseAxes *axes = &source->extpars->axes; overlay[Noverlay].type = KII_OVERLAY_CIRCLE; overlay[Noverlay].x = source->peak->xf; overlay[Noverlay].y = source->peak->yf; - overlay[Noverlay].dx = 2.0*petrosian->radius; - overlay[Noverlay].dy = 2.0*petrosian->radius*profile->axes.minor/profile->axes.major; - overlay[Noverlay].angle = profile->axes.theta * PS_DEG_RAD; + overlay[Noverlay].dx = 1.0*petrosianRadius; + overlay[Noverlay].dy = 1.0*petrosianRadius*axes->minor/axes->major; + overlay[Noverlay].angle = axes->theta * PS_DEG_RAD; overlay[Noverlay].text = NULL; Noverlay ++; CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); - overlay[Noverlay].type = KII_OVERLAY_CIRCLE; - overlay[Noverlay].x = source->peak->xf; - overlay[Noverlay].y = source->peak->yf; - overlay[Noverlay].dx = 4.0*petrosian->radius; - overlay[Noverlay].dy = 4.0*petrosian->radius*profile->axes.minor/profile->axes.major; - overlay[Noverlay].angle = profile->axes.theta * PS_DEG_RAD; - overlay[Noverlay].text = NULL; - Noverlay ++; - CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); + // overlay[Noverlay].type = KII_OVERLAY_CIRCLE; + // overlay[Noverlay].x = source->peak->xf; + // overlay[Noverlay].y = source->peak->yf; + // overlay[Noverlay].dx = 2.0*petrosianRadius; + // overlay[Noverlay].dy = 2.0*petrosianRadius*axes->minor/axes->major; + // overlay[Noverlay].angle = axes->theta * PS_DEG_RAD; + // overlay[Noverlay].text = NULL; + // Noverlay ++; + // CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); }