Index: trunk/psphot/src/Makefile.am =================================================================== --- trunk/psphot/src/Makefile.am (revision 27695) +++ trunk/psphot/src/Makefile.am (revision 27819) @@ -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: trunk/psphot/src/psphot.h =================================================================== --- trunk/psphot/src/psphot.h (revision 27695) +++ trunk/psphot/src/psphot.h (revision 27819) @@ -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: trunk/psphot/src/psphotEllipticalContour.c =================================================================== --- trunk/psphot/src/psphotEllipticalContour.c (revision 27695) +++ trunk/psphot/src/psphotEllipticalContour.c (revision 27819) @@ -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: trunk/psphot/src/psphotEllipticalProfile.c =================================================================== --- trunk/psphot/src/psphotEllipticalProfile.c (revision 27695) +++ trunk/psphot/src/psphotEllipticalProfile.c (revision 27819) @@ -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: trunk/psphot/src/psphotExtendedSourceAnalysis.c =================================================================== --- trunk/psphot/src/psphotExtendedSourceAnalysis.c (revision 27695) +++ trunk/psphot/src/psphotExtendedSourceAnalysis.c (revision 27819) @@ -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: trunk/psphot/src/psphotPetrosian.c =================================================================== --- trunk/psphot/src/psphotPetrosian.c (revision 27695) +++ trunk/psphot/src/psphotPetrosian.c (revision 27819) @@ -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: trunk/psphot/src/psphotPetrosianRadialBins.c =================================================================== --- trunk/psphot/src/psphotPetrosianRadialBins.c (revision 27695) +++ trunk/psphot/src/psphotPetrosianRadialBins.c (revision 27819) @@ -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: trunk/psphot/src/psphotPetrosianStats.c =================================================================== --- trunk/psphot/src/psphotPetrosianStats.c (revision 27695) +++ trunk/psphot/src/psphotPetrosianStats.c (revision 27819) @@ -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: trunk/psphot/src/psphotPetrosianVisual.c =================================================================== --- trunk/psphot/src/psphotPetrosianVisual.c (revision 27695) +++ trunk/psphot/src/psphotPetrosianVisual.c (revision 27819) @@ -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: trunk/psphot/src/psphotRadialBins.c =================================================================== --- trunk/psphot/src/psphotRadialBins.c (revision 27819) +++ trunk/psphot/src/psphotRadialBins.c (revision 27819) @@ -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: trunk/psphot/src/psphotRadialProfile.c =================================================================== --- trunk/psphot/src/psphotRadialProfile.c (revision 27695) +++ trunk/psphot/src/psphotRadialProfile.c (revision 27819) @@ -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: trunk/psphot/src/psphotRadialProfileByAngles.c =================================================================== --- trunk/psphot/src/psphotRadialProfileByAngles.c (revision 27695) +++ trunk/psphot/src/psphotRadialProfileByAngles.c (revision 27819) @@ -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: trunk/psphot/src/psphotRadiiFromProfiles.c =================================================================== --- trunk/psphot/src/psphotRadiiFromProfiles.c (revision 27695) +++ trunk/psphot/src/psphotRadiiFromProfiles.c (revision 27819) @@ -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: trunk/psphot/src/psphotVisual.c =================================================================== --- trunk/psphot/src/psphotVisual.c (revision 27695) +++ trunk/psphot/src/psphotVisual.c (revision 27819) @@ -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); }