Index: /branches/eam_branches/ipp-20110213/psphot/src/psphot.h =================================================================== --- /branches/eam_branches/ipp-20110213/psphot/src/psphot.h (revision 30706) +++ /branches/eam_branches/ipp-20110213/psphot/src/psphot.h (revision 30707) @@ -430,5 +430,5 @@ bool psphotRadialApertures (pmConfig *config, const pmFPAview *view, const char *filerule); bool psphotRadialAperturesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe); -bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *radMax, int entry); +bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, psImageMaskType maskVal, const psVector *radMax, int entry); bool psphotExtendedSourceAnalysisByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule); Index: /branches/eam_branches/ipp-20110213/psphot/src/psphotBlendFit.c =================================================================== --- /branches/eam_branches/ipp-20110213/psphot/src/psphotBlendFit.c (revision 30706) +++ /branches/eam_branches/ipp-20110213/psphot/src/psphotBlendFit.c (revision 30707) @@ -176,5 +176,5 @@ } psFree(job); - } + } } psFree (cellGroups); Index: /branches/eam_branches/ipp-20110213/psphot/src/psphotExtendedSourceAnalysis.c =================================================================== --- /branches/eam_branches/ipp-20110213/psphot/src/psphotExtendedSourceAnalysis.c (revision 30706) +++ /branches/eam_branches/ipp-20110213/psphot/src/psphotExtendedSourceAnalysis.c (revision 30707) @@ -1,9 +1,6 @@ # include "psphotInternal.h" -// ?? these cannot happen here --> we would need to do this in psphotExtendedSourceAnalysis -// XXX option to choose a consistent position -// XXX option to choose a consistent elliptical contour -// XXX SDSS uses the r-band petrosian radius to measure petrosian fluxes in all bands -// XXX consistent choice of extendedness... +// measure the elliptical radial profile and use this to measure the petrosian parameters for the sources +// XXX this function needs to be threaded // for now, let's store the detections on the readout->analysis for each readout @@ -40,7 +37,5 @@ int Next = 0; int Npetro = 0; - int Nisophot = 0; int Nannuli = 0; - int Nkron = 0; psTimerStart ("psphot.extended"); @@ -68,8 +63,10 @@ assert (maskVal); - // XXX require petrosian analysis for non-linear fits? - - // XXX temporary user-supplied systematic sky noise measurement (derive from background model) - float skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE"); + // get the sky noise from the background analysis; if this is missing, get the user-supplied value + float skynoise = psMetadataLookupF32 (&status, readout->analysis, "SKY_STDEV"); + if (!status) { + skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE"); + psWarning ("failed to get sky noise level from background analysis; defaulting to user supplied value of %f\n", skynoise); + } // S/N limit to perform full non-linear fits @@ -78,14 +75,6 @@ // which extended source analyses should we perform? bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); - bool doIsophotal = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL"); bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); - bool doKron = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON"); - -# if (0) - // if backModel or backStdev are missing, the values of sky and/or skyErr will be set to NAN - // XXX use this to set skynoise - pmReadout *backModel = psphotSelectBackground (config, view); - pmReadout *backStdev = psphotSelectBackgroundStdev (config, view); -# endif + bool doPetroStars = psMetadataLookupBool (&status, recipe, "PETROSIAN_FOR_STARS"); // source analysis is done in S/N order (brightest first) @@ -103,10 +92,14 @@ // skip PSF-like and non-astronomical objects - if (source->type == PM_SOURCE_TYPE_STAR) continue; if (source->type == PM_SOURCE_TYPE_DEFECT) continue; if (source->type == PM_SOURCE_TYPE_SATURATED) continue; if (source->mode & PM_SOURCE_MODE_DEFECT) continue; if (source->mode & PM_SOURCE_MODE_SATSTAR) continue; - if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue; + + // optionally allow non-extended objects to get petrosians as well + if (!doPetroStars) { + if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue; + if (source->type == PM_SOURCE_TYPE_STAR) continue; + } // limit selection to some SN limit @@ -119,6 +112,4 @@ if (source->peak->x > AnalysisRegion.x1) continue; if (source->peak->y > AnalysisRegion.y1) continue; - - // fprintf (stderr, "xsrc: %f, %f : %f\n", source->peak->xf, source->peak->yf, source->peak->SN); // replace object in image @@ -136,5 +127,5 @@ // if we request any of these measurements, we require the radial profile - if (doPetrosian || doIsophotal || doAnnuli || doKron) { + if (doPetrosian || doAnnuli) { if (!psphotRadialProfile (source, recipe, skynoise, maskVal)) { // all measurements below require the radial profile; skip them all @@ -144,4 +135,5 @@ continue; } + Nannuli ++; source->mode |= PM_SOURCE_MODE_RADIAL_FLUX; } @@ -169,7 +161,5 @@ psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next); psLogMsg ("psphot", PS_LOG_INFO, " %d petrosian\n", Npetro); - psLogMsg ("psphot", PS_LOG_INFO, " %d isophotal\n", Nisophot); psLogMsg ("psphot", PS_LOG_INFO, " %d annuli\n", Nannuli); - psLogMsg ("psphot", PS_LOG_INFO, " %d kron\n", Nkron); psphotVisualShowResidualImage (readout, false); Index: /branches/eam_branches/ipp-20110213/psphot/src/psphotExtendedSourceFits.c =================================================================== --- /branches/eam_branches/ipp-20110213/psphot/src/psphotExtendedSourceFits.c (revision 30706) +++ /branches/eam_branches/ipp-20110213/psphot/src/psphotExtendedSourceFits.c (revision 30707) @@ -178,30 +178,33 @@ PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfain - if (false && !psThreadJobAddPending(job)) { +// set this to 0 to run without threading +# if (1) + if (!psThreadJobAddPending(job)) { psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); psFree(AnalysisRegion); return false; - } else { - if (!psphotExtendedSourceFits_Threaded(job)) { - psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - psFree(AnalysisRegion); - return false; - } - psScalar *scalar = NULL; - scalar = job->args->data[7]; - Next += scalar->data.S32; - scalar = job->args->data[8]; - Nconvolve += scalar->data.S32; - scalar = job->args->data[9]; - NconvolvePass += scalar->data.S32; - scalar = job->args->data[10]; - Nplain += scalar->data.S32; - scalar = job->args->data[11]; - NplainPass += scalar->data.S32; - scalar = job->args->data[12]; - Nfaint += scalar->data.S32; - psFree(job); + } +# else + if (!psphotExtendedSourceFits_Threaded(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + psFree(AnalysisRegion); + return false; } - } + psScalar *scalar = NULL; + scalar = job->args->data[7]; + Next += scalar->data.S32; + scalar = job->args->data[8]; + Nconvolve += scalar->data.S32; + scalar = job->args->data[9]; + NconvolvePass += scalar->data.S32; + scalar = job->args->data[10]; + Nplain += scalar->data.S32; + scalar = job->args->data[11]; + NplainPass += scalar->data.S32; + scalar = job->args->data[12]; + Nfaint += scalar->data.S32; + psFree(job); +# endif + } // wait for the threads to finish and manage results @@ -268,7 +271,11 @@ psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA); + pthread_t tid = pthread_self(); // Thread identifier + // Define source fitting parameters for extended source fits pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc(); fitOptions->mode = PM_SOURCE_FIT_EXT; + fitOptions->saveCovariance = true; // XXX make this a user option? + // XXX for now, use the defaults for the rest: // fitOptions->nIter = fitIter; @@ -296,4 +303,6 @@ // XXX use the parameters guessed from moments // if (source->modelEXT == NULL) continue; + + fprintf (stderr, "fit %d,%d in thread %d\n", source->peak->x, source->peak->y, (int) tid); // replace object in image Index: /branches/eam_branches/ipp-20110213/psphot/src/psphotPetrosianRadialBins.c =================================================================== --- /branches/eam_branches/ipp-20110213/psphot/src/psphotPetrosianRadialBins.c (revision 30706) +++ /branches/eam_branches/ipp-20110213/psphot/src/psphotPetrosianRadialBins.c (revision 30707) @@ -51,4 +51,5 @@ psVector *binRad = psVectorAllocEmpty(nMax, PS_TYPE_F32); // mean radius of radial bin psVector *binArea = psVectorAllocEmpty(nMax, PS_TYPE_F32); // area of radial bin (contiguous, non-overlapping) + psVector *binFill = psVectorAllocEmpty(nMax, PS_TYPE_F32); // fraction of radial bin with valid pixels psVectorInit (binSB, 0.0); @@ -144,4 +145,5 @@ binSB->data.F32[nOut] = value; binSBstdev->data.F32[nOut] = sqrt(PS_SQR(dvalue) / values->n + skyModelErrorSQ); + binFill->data.F32[nOut] = values->n / binArea->data.F32[nOut]; // error in the SB is the stdev per bin / sqrt (number of pixels) @@ -176,4 +178,5 @@ psFree(binRad); psFree(binArea); + psFree(binFill); psFree(radMin); psFree(radMax); @@ -202,8 +205,10 @@ psFree(profile->radialBins); psFree(profile->area); + psFree(profile->binFill); // save the vectors profile->radialBins = binRad; profile->area = binArea; + profile->binFill = binFill; profile->binSB = binSB; profile->binSBstdev = binSBstdev; Index: /branches/eam_branches/ipp-20110213/psphot/src/psphotPetrosianStats.c =================================================================== --- /branches/eam_branches/ipp-20110213/psphot/src/psphotPetrosianStats.c (revision 30706) +++ /branches/eam_branches/ipp-20110213/psphot/src/psphotPetrosianStats.c (revision 30707) @@ -6,5 +6,7 @@ // generate the Petrosian radius and flux from the mean surface brightness (r_i) -float InterpolateValues (float X0, float Y0, float X1, float Y1, float X); +float InterpolateValues (float X0, float Y0, float X1, float Y1, float X); +float InterpolateValuesErrX (float X0, float Y0, float X1, float Y1, float X, float dX0, float dX1); +float InterpolateValuesErrY (float X0, float Y0, float X1, float Y1, float X, float dY0, float dY1); bool psphotPetrosianStats (pmSource *source) { @@ -25,4 +27,5 @@ psVector *binRad = profile->radialBins; psVector *area = profile->area; + psVector *binFill = profile->binFill; psVector *fluxSum = psVectorAllocEmpty(binSB->n, PS_TYPE_F32); @@ -33,7 +36,12 @@ psVector *meanSB = psVectorAllocEmpty(binSB->n, PS_TYPE_F32); psVector *areaSum = psVectorAllocEmpty(binSB->n, PS_TYPE_F32); + psVector *apixSum = psVectorAllocEmpty(binSB->n, PS_TYPE_F32); float petRadius = NAN; + float petRadiusErr = NAN; float petFlux = NAN; + float petFluxErr = NAN; + float petArea = NAN; + float petApix = NAN; bool anyPetro = false; @@ -41,4 +49,5 @@ bool above = true; float Asum = 0.0; + float Psum = 0.0; float Fsum = 0.0; float dFsum2 = 0.0; @@ -56,7 +65,8 @@ float Area = area->data.F32[i]; - Asum += Area; + Asum += Area; // Asum is the cumulative area interior to this bin Fsum += binSB->data.F32[i] * Area; dFsum2 += PS_SQR(binSBstdev->data.F32[i] * Area); + Psum += Area*binFill->data.F32[i]; // Psum is the cumulative number of pixels interior to this bin float areaInner = 0.5 * Area; @@ -87,14 +97,15 @@ psVectorAppend(areaSum, Asum); + psVectorAppend(apixSum, Psum); psVectorAppend(fluxSum, Fsum); psVectorAppend(fluxSumErr2, dFsum2); psVectorAppend(refRadius, binRad->data.F32[i]); - psTrace ("psphot", 4, "%3d : %5.2f : %5.3f %5.3f : %5.3f %5.3f : %5.3f %5.3f : %5.3f %5.3f : %5.1f %5.1f\n", + psTrace ("psphot", 4, "%3d : %5.2f : %5.3f %5.3f : %5.3f %5.3f : %5.3f %5.3f : %5.3f %5.3f : %5.1f %5.1f %5.1f\n", i, refRadius->data.F32[nOut], binSB->data.F32[i], binSBstdev->data.F32[i], meanSB->data.F32[nOut], meanSBerr, petRatio->data.F32[nOut], petRatioErr->data.F32[nOut], - fluxSum->data.F32[nOut], sqrt(fluxSumErr2->data.F32[nOut]), areaSum->data.F32[nOut], areaInner); + fluxSum->data.F32[nOut], sqrt(fluxSumErr2->data.F32[nOut]), areaSum->data.F32[nOut], apixSum->data.F32[nOut], areaInner); // anytime we transition below the PETROSIAN_RATIO, calculate the radius and flux @@ -104,7 +115,9 @@ if (i == 0) { // assume Fmax @ R = 0.0 - petRadius = InterpolateValues (1.0, 0.0, petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO); - } else { - petRadius = InterpolateValues (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO); + petRadius = InterpolateValues (1.0, 0.0, petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO); + petRadiusErr = InterpolateValuesErrX (1.0, 0.0, petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO, 0.0, petRatioErr->data.F32[nOut]); + } else { + petRadius = InterpolateValues (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO); + petRadiusErr = InterpolateValuesErrX (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO, petRatioErr->data.F32[nOut-1], petRatioErr->data.F32[nOut]); } above = false; @@ -128,11 +141,17 @@ } + // if we failed to reach the PETROSIAN_RATIO, use the lowest significant ratio instead (flag this!) if (!anyPetro) { // interpolate Rvec between i-1 and i to PETROSIAN_RATIO to get flux (Fvec) and radius (rvec) if (lowestSignificantRadius == 0) { // assume Fmax @ R = 0.0 - petRadius = InterpolateValues (1.0, 0.0, petRatio->data.F32[lowestSignificantRadius], refRadius->data.F32[lowestSignificantRadius], PETROSIAN_RATIO); + petRadius = InterpolateValues (1.0, 0.0, petRatio->data.F32[lowestSignificantRadius], refRadius->data.F32[lowestSignificantRadius], PETROSIAN_RATIO); + petRadiusErr = InterpolateValuesErrX (1.0, 0.0, petRatio->data.F32[lowestSignificantRadius], refRadius->data.F32[lowestSignificantRadius], PETROSIAN_RATIO, 0.0, petRatioErr->data.F32[lowestSignificantRadius]); + } else { - petRadius = InterpolateValues (petRatio->data.F32[lowestSignificantRadius-1], refRadius->data.F32[lowestSignificantRadius-1], petRatio->data.F32[lowestSignificantRadius], refRadius->data.F32[lowestSignificantRadius], PETROSIAN_RATIO); + int n0 = lowestSignificantRadius-1; + int n1 = lowestSignificantRadius; + petRadius = InterpolateValues (petRatio->data.F32[n0], refRadius->data.F32[n0], petRatio->data.F32[n1], refRadius->data.F32[n1], PETROSIAN_RATIO); + petRadiusErr = InterpolateValuesErrX (petRatio->data.F32[n0], refRadius->data.F32[n0], petRatio->data.F32[n1], refRadius->data.F32[n1], PETROSIAN_RATIO, petRatioErr->data.F32[n0], petRatioErr->data.F32[n1]); } } @@ -147,5 +166,8 @@ continue; } else { - petFlux = InterpolateValues (refRadius->data.F32[i-1], fluxSum->data.F32[i-1], refRadius->data.F32[i], fluxSum->data.F32[i], apRadius); + petFlux = InterpolateValues (refRadius->data.F32[i-1], fluxSum->data.F32[i-1], refRadius->data.F32[i], fluxSum->data.F32[i], apRadius); + petFluxErr = InterpolateValuesErrY (refRadius->data.F32[i-1], fluxSum->data.F32[i-1], refRadius->data.F32[i], fluxSum->data.F32[i], apRadius, sqrt(fluxSumErr2->data.F32[i-1]), sqrt(fluxSumErr2->data.F32[i])); + petArea = InterpolateValues (refRadius->data.F32[i-1], areaSum->data.F32[i-1], refRadius->data.F32[i], areaSum->data.F32[i], apRadius); + petApix = InterpolateValues (refRadius->data.F32[i-1], apixSum->data.F32[i-1], refRadius->data.F32[i], apixSum->data.F32[i], apRadius); break; } @@ -158,4 +180,6 @@ float R50 = NAN; float R90 = NAN; + float R50err = NAN; + float R90err = NAN; bool found50 = false; bool found90 = false; @@ -167,5 +191,6 @@ continue; } else { - R50 = InterpolateValues (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux50); + R50 = InterpolateValues (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux50); + R50err = InterpolateValuesErrX (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux50, sqrt(fluxSumErr2->data.F32[i-1]), sqrt(fluxSumErr2->data.F32[i])); found50 = true; } @@ -176,5 +201,6 @@ continue; } else { - R90 = InterpolateValues (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux90); + R90 = InterpolateValues (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux90); + R90err = InterpolateValuesErrX (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux90, sqrt(fluxSumErr2->data.F32[i-1]), sqrt(fluxSumErr2->data.F32[i])); found90 = true; } @@ -188,10 +214,11 @@ source->extpars->petrosianR50 = R50; source->extpars->petrosianR90 = R90; + source->extpars->petrosianFill = petApix / petArea; // XXX add the errors - source->extpars->petrosianRadiusErr = NAN; - source->extpars->petrosianFluxErr = NAN; - source->extpars->petrosianR50Err = NAN; - source->extpars->petrosianR90Err = NAN; + source->extpars->petrosianRadiusErr = petRadiusErr; + source->extpars->petrosianFluxErr = petFluxErr; + source->extpars->petrosianR50Err = R50err; + source->extpars->petrosianR90Err = R90err; // fprintf (stderr, "source @ %f,%f\n", source->peak->xf, source->peak->yf); @@ -205,4 +232,5 @@ psFree(meanSB); psFree(areaSum); + psFree(apixSum); return true; @@ -210,6 +238,24 @@ float InterpolateValues (float X0, float Y0, float X1, float Y1, float X) { - float Y = Y0 + (Y1 - Y0) * (X - X0) / (X1 - X0); + float dydx = (Y1 - Y0) / (X1 - X0); + float Y = Y0 + dydx * (X - X0); return Y; } +float InterpolateValuesErrX (float X0, float Y0, float X1, float Y1, float X, float dX0, float dX1) { + + float dydx = (Y1 - Y0) / (X1 - X0); + float dxdx = (X - X0) / (X1 - X0); + + float dY = sqrt(PS_SQR(dX1*dydx*dxdx) + PS_SQR(dX0*dydx*(dxdx - 1.0))); + return dY; +} + +float InterpolateValuesErrY (float X0, float Y0, float X1, float Y1, float X, float dY0, float dY1) { + + float dxdx = (X - X0) / (X1 - X0); + + float dY = sqrt(PS_SQR(dY1*dxdx) + PS_SQR(dY0*(1.0 - dxdx))); + return dY; +} + Index: /branches/eam_branches/ipp-20110213/psphot/src/psphotRadialApertures.c =================================================================== --- /branches/eam_branches/ipp-20110213/psphot/src/psphotRadialApertures.c (revision 30706) +++ /branches/eam_branches/ipp-20110213/psphot/src/psphotRadialApertures.c (revision 30707) @@ -37,4 +37,10 @@ int Nradial = 0; + // perform full non-linear fits / extended source analysis? + if (!psMetadataLookupBool (&status, recipe, "RADIAL_APERTURES")) { + psLogMsg ("psphot", PS_LOG_INFO, "skipping radial apertures\n"); + return true; + } + psTimerStart ("psphot.radial"); @@ -62,11 +68,9 @@ psAssert (radMax, "annular bins (RADIAL.ANNULAR.BINS.UPPER) are not defined in the recipe"); psAssert (radMax->n, "no valid annular bins (RADIAL.ANNULAR.BINS.UPPER) are define"); + float outerRadius = radMax->data.F32[radMax->n - 1]; // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels assert (maskVal); - - // XXX temporary user-supplied systematic sky noise measurement (derive from background model) - float skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE"); // S/N limit to perform full non-linear fits @@ -112,14 +116,18 @@ pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); } + + // we need to change the view for the radial aperture analysis, but we want to recover exactly + // the original view; the following elements get destroyed by pmSourceRedefinePixels so save them: + psImage *oldMaskObj = psMemIncrRefCounter(source->maskObj); + psImage *oldModelFlux = psMemIncrRefCounter(source->modelFlux); + psImage *oldPSFimage = psMemIncrRefCounter(source->psfImage); + psRegion oldRegion = source->region; + Nradial ++; // force source image to be a bit larger... - float radius = source->peak->xf - source->pixels->col0; - radius = PS_MAX (radius, source->peak->yf - source->pixels->row0); - radius = PS_MAX (radius, source->pixels->numRows - source->peak->yf + source->pixels->row0); - radius = PS_MAX (radius, source->pixels->numCols - source->peak->xf + source->pixels->col0); - pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, 1.5*radius); - - if (!psphotRadialApertureSource (source, recipe, skynoise, maskVal, radMax, 0)) { + pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, outerRadius + 2); + + if (!psphotRadialApertureSource (source, recipe, maskVal, radMax, 0)) { psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); } else { @@ -127,4 +135,9 @@ } + pmSourceRedefinePixelsByRegion (source, readout, oldRegion); + psFree(source->maskObj); source->maskObj = oldMaskObj; + psFree(source->modelFlux); source->modelFlux = oldModelFlux; + psFree(source->psfImage); source->psfImage = oldPSFimage; + // re-subtract the object, leave local sky pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); @@ -135,5 +148,5 @@ } -bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *aperRadii, int entry) { +bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, psImageMaskType maskVal, const psVector *aperRadii, int entry) { // if we are a child source, save the results to the parent source radial aperture array @@ -203,7 +216,9 @@ psVector *flux = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin psVector *fluxErr = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin + psVector *fluxStd = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin psVector *fill = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin psVectorInit (flux, 0.0); + psVectorInit (fluxStd, 0.0); psVectorInit (fluxErr, 0.0); psVectorInit (fill, 0.0); @@ -215,13 +230,28 @@ for (int j = 0; (*rPix2 < *aRad2) && (j < aperRadii2->n); j++, aRad2++) { flux->data.F32[j] += pixFlux->data.F32[i]; + fluxStd->data.F32[j] += PS_SQR(pixFlux->data.F32[i]); fluxErr->data.F32[j] += pixVar->data.F32[i]; fill->data.F32[j] += 1.0; } } + + /* for each radial bin, R(i), we measure: + 1) the flux within that aperture: F(i) = \sum_{r_jn; i++) { // calculate the total flux for bin 'nOut' float Area = M_PI*aperRadii2->data.F32[i]; - fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]); + + int nPix = fill->data.F32[i]; + float SBmean = flux->data.F32[i] / nPix; + float SBstdv = sqrt((fluxStd->data.F32[i] / nPix) - PS_SQR(SBmean)); + + flux->data.F32[i] = SBmean * Area; + fluxStd->data.F32[i] = SBstdv * Area; + fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]) * Area / nPix; + fill->data.F32[i] /= Area; psTrace ("psphot", 5, "radial bins: %3d %5.1f : %8.1f +/- %7.2f : %4.2f %6.1f\n", @@ -230,4 +260,5 @@ radialAper->flux = flux; + radialAper->fluxStdev = fluxStd; radialAper->fluxErr = fluxErr; radialAper->fill = fill; @@ -245,5 +276,5 @@ static int nPix = 0; -bool psphotRadialApertureSource_With_Sort (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *radMax, int entry) { +bool psphotRadialApertureSource_With_Sort (pmSource *source, psMetadata *recipe, psImageMaskType maskVal, const psVector *radMax, int entry) { psAssert(source->radialAper->data[entry] == NULL, "why is this already defined?"); Index: /branches/eam_branches/ipp-20110213/psphot/src/psphotRadialAperturesByObject.c =================================================================== --- /branches/eam_branches/ipp-20110213/psphot/src/psphotRadialAperturesByObject.c (revision 30706) +++ /branches/eam_branches/ipp-20110213/psphot/src/psphotRadialAperturesByObject.c (revision 30707) @@ -42,7 +42,4 @@ psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels assert (maskVal); - - // XXX temporary user-supplied systematic sky noise measurement (derive from background model) - float skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE"); // S/N limit to perform full non-linear fits @@ -147,11 +144,7 @@ // force source image to be a bit larger... - // float radius = source->peak->xf - source->pixels->col0; - // radius = PS_MAX (radius, source->peak->yf - source->pixels->row0); - // radius = PS_MAX (radius, source->pixels->numRows - source->peak->yf + source->pixels->row0); - // radius = PS_MAX (radius, source->pixels->numCols - source->peak->xf + source->pixels->col0); pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, outerRadius + 2); - if (!psphotRadialApertureSource (source, recipe, skynoise, maskVal, radMax, nMatchedPSF)) { + if (!psphotRadialApertureSource (source, recipe, maskVal, radMax, nMatchedPSF)) { psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); } else { Index: /branches/eam_branches/ipp-20110213/psphot/src/psphotReadout.c =================================================================== --- /branches/eam_branches/ipp-20110213/psphot/src/psphotReadout.c (revision 30706) +++ /branches/eam_branches/ipp-20110213/psphot/src/psphotReadout.c (revision 30707) @@ -203,4 +203,5 @@ psphotExtendedSourceAnalysis (config, view, filerule); // pass 1 (detections->allSources) psphotExtendedSourceFits (config, view, filerule); // pass 1 (detections->allSources) + psphotRadialApertures(config, view, filerule); finish: Index: /branches/eam_branches/ipp-20110213/psphot/src/psphotStackMatchPSFsUtils.c =================================================================== --- /branches/eam_branches/ipp-20110213/psphot/src/psphotStackMatchPSFsUtils.c (revision 30706) +++ /branches/eam_branches/ipp-20110213/psphot/src/psphotStackMatchPSFsUtils.c (revision 30707) @@ -363,5 +363,5 @@ // we need to register the FWHM values for use downstream - pmSubtractionSetFWHMs(options->inputSeeing->data.F32[index], options->targetSeeing); + pmSubtractionSetFWHMs(options->targetSeeing, options->inputSeeing->data.F32[index]); pmSubtractionParamScaleOptions(scale, scaleRef, scaleMin, scaleMax);