IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 30707


Ignore:
Timestamp:
Feb 19, 2011, 10:40:50 AM (15 years ago)
Author:
eugene
Message:

remove invalid references to isophot and kron from extended source analysis; use the measured sky stdev instead of the user guess; optionally allow stars as well as galaxies to get petrosians; fix target vs input seeing in psphotStack; allow radial apertures for single-image analysis as well as stack; calculate the radial aperture flux errors; calculate the petrosian parameter errors; turn on threads for ext source fits; save covar matrix for ext source fits; calculate petrosian fill factor

Location:
branches/eam_branches/ipp-20110213/psphot
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20110213/psphot

    • Property svn:mergeinfo changed (with no actual effect on merging)

  • branches/eam_branches/ipp-20110213/psphot/src/psphot.h

    r30624 r30707  
    430430bool psphotRadialApertures (pmConfig *config, const pmFPAview *view, const char *filerule);
    431431bool psphotRadialAperturesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
    432 bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *radMax, int entry);
     432bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, psImageMaskType maskVal, const psVector *radMax, int entry);
    433433
    434434bool psphotExtendedSourceAnalysisByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule);

  • branches/eam_branches/ipp-20110213/psphot/src/psphotBlendFit.c

    r30624 r30707  
    176176            }
    177177            psFree(job);
    178             }
     178        }
    179179    }
    180180    psFree (cellGroups);

  • branches/eam_branches/ipp-20110213/psphot/src/psphotExtendedSourceAnalysis.c

    r30624 r30707  
    11# include "psphotInternal.h"
    22
    3 // ?? these cannot happen here --> we would need to do this in psphotExtendedSourceAnalysis
    4 // XXX option to choose a consistent position
    5 // XXX option to choose a consistent elliptical contour
    6 // XXX SDSS uses the r-band petrosian radius to measure petrosian fluxes in all bands
    7 // XXX consistent choice of extendedness...
     3// measure the elliptical radial profile and use this to measure the petrosian parameters for the sources
     4// XXX this function needs to be threaded
    85
    96// for now, let's store the detections on the readout->analysis for each readout
    … …  
    4037    int Next = 0;
    4138    int Npetro = 0;
    42     int Nisophot = 0;
    4339    int Nannuli = 0;
    44     int Nkron = 0;
    4540
    4641    psTimerStart ("psphot.extended");
    … …  
    6863    assert (maskVal);
    6964
    70     // XXX require petrosian analysis for non-linear fits?
    71 
    72     // XXX temporary user-supplied systematic sky noise measurement (derive from background model)
    73     float skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE");
     65    // get the sky noise from the background analysis; if this is missing, get the user-supplied value
     66    float skynoise = psMetadataLookupF32 (&status, readout->analysis, "SKY_STDEV");
     67    if (!status) {
     68        skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE");
     69        psWarning ("failed to get sky noise level from background analysis; defaulting to user supplied value of %f\n", skynoise);
     70    }
    7471
    7572    // S/N limit to perform full non-linear fits
    … …  
    7875    // which extended source analyses should we perform?
    7976    bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
    80     bool doIsophotal    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL");
    8177    bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
    82     bool doKron         = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON");
    83 
    84 # if (0)
    85     // if backModel or backStdev are missing, the values of sky and/or skyErr will be set to NAN
    86     // XXX use this to set skynoise
    87     pmReadout *backModel = psphotSelectBackground (config, view);
    88     pmReadout *backStdev = psphotSelectBackgroundStdev (config, view);
    89 # endif
     78    bool doPetroStars   = psMetadataLookupBool (&status, recipe, "PETROSIAN_FOR_STARS");
    9079
    9180    // source analysis is done in S/N order (brightest first)
    … …  
    10392
    10493        // skip PSF-like and non-astronomical objects
    105         if (source->type == PM_SOURCE_TYPE_STAR) continue;
    10694        if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
    10795        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
    10896        if (source->mode & PM_SOURCE_MODE_DEFECT) continue;
    10997        if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;
    110         if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue;
     98
     99        // optionally allow non-extended objects to get petrosians as well
     100        if (!doPetroStars) {
     101            if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue;
     102            if (source->type == PM_SOURCE_TYPE_STAR) continue;
     103        }
    111104
    112105        // limit selection to some SN limit
    … …  
    119112        if (source->peak->x > AnalysisRegion.x1) continue;
    120113        if (source->peak->y > AnalysisRegion.y1) continue;
    121 
    122         // fprintf (stderr, "xsrc: %f, %f : %f\n", source->peak->xf, source->peak->yf, source->peak->SN);
    123114
    124115        // replace object in image
    … …  
    136127
    137128        // if we request any of these measurements, we require the radial profile
    138         if (doPetrosian || doIsophotal || doAnnuli || doKron) {
     129        if (doPetrosian || doAnnuli) {
    139130            if (!psphotRadialProfile (source, recipe, skynoise, maskVal)) {
    140131                // all measurements below require the radial profile; skip them all
    … …  
    144135                continue;
    145136            }
     137            Nannuli ++;
    146138            source->mode |= PM_SOURCE_MODE_RADIAL_FLUX;
    147139        }
    … …  
    169161    psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot.extended"), Next);
    170162    psLogMsg ("psphot", PS_LOG_INFO, "  %d petrosian\n", Npetro);
    171     psLogMsg ("psphot", PS_LOG_INFO, "  %d isophotal\n", Nisophot);
    172163    psLogMsg ("psphot", PS_LOG_INFO, "  %d annuli\n", Nannuli);
    173     psLogMsg ("psphot", PS_LOG_INFO, "  %d kron\n", Nkron);
    174164
    175165    psphotVisualShowResidualImage (readout, false);

  • branches/eam_branches/ipp-20110213/psphot/src/psphotExtendedSourceFits.c

    r30624 r30707  
    178178            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfain
    179179
    180             if (false && !psThreadJobAddPending(job)) {
     180// set this to 0 to run without threading
     181# if (1)           
     182            if (!psThreadJobAddPending(job)) {
    181183                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    182184                psFree(AnalysisRegion);
    183185                return false;
    184             } else {
    185                 if (!psphotExtendedSourceFits_Threaded(job)) {
    186                     psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    187                     psFree(AnalysisRegion);
    188                     return false;
    189                 }
    190                 psScalar *scalar = NULL;
    191                 scalar = job->args->data[7];
    192                 Next += scalar->data.S32;
    193                 scalar = job->args->data[8];
    194                 Nconvolve += scalar->data.S32;
    195                 scalar = job->args->data[9];
    196                 NconvolvePass += scalar->data.S32;
    197                 scalar = job->args->data[10];
    198                 Nplain += scalar->data.S32;
    199                 scalar = job->args->data[11];
    200                 NplainPass += scalar->data.S32;
    201                 scalar = job->args->data[12];
    202                 Nfaint += scalar->data.S32;
    203                 psFree(job);
     186            }
     187# else
     188            if (!psphotExtendedSourceFits_Threaded(job)) {
     189                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     190                psFree(AnalysisRegion);
     191                return false;
    204192            }
    205         }
     193            psScalar *scalar = NULL;
     194            scalar = job->args->data[7];
     195            Next += scalar->data.S32;
     196            scalar = job->args->data[8];
     197            Nconvolve += scalar->data.S32;
     198            scalar = job->args->data[9];
     199            NconvolvePass += scalar->data.S32;
     200            scalar = job->args->data[10];
     201            Nplain += scalar->data.S32;
     202            scalar = job->args->data[11];
     203            NplainPass += scalar->data.S32;
     204            scalar = job->args->data[12];
     205            Nfaint += scalar->data.S32;
     206            psFree(job);
     207# endif
     208        }
    206209
    207210        // wait for the threads to finish and manage results
    … …  
    268271    psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA);
    269272
     273    pthread_t tid = pthread_self();     // Thread identifier
     274
    270275    // Define source fitting parameters for extended source fits
    271276    pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc();
    272277    fitOptions->mode          = PM_SOURCE_FIT_EXT;
     278    fitOptions->saveCovariance = true;  // XXX make this a user option?
     279
    273280    // XXX for now, use the defaults for the rest:
    274281    // fitOptions->nIter         = fitIter;
    … …  
    296303        // XXX use the parameters guessed from moments
    297304        // if (source->modelEXT == NULL) continue;
     305
     306        fprintf (stderr, "fit %d,%d in thread %d\n", source->peak->x, source->peak->y, (int) tid);
    298307
    299308        // replace object in image

  • branches/eam_branches/ipp-20110213/psphot/src/psphotPetrosianRadialBins.c

    r28013 r30707  
    5151    psVector *binRad     = psVectorAllocEmpty(nMax, PS_TYPE_F32); // mean radius of radial bin
    5252    psVector *binArea    = psVectorAllocEmpty(nMax, PS_TYPE_F32); // area of radial bin (contiguous, non-overlapping)
     53    psVector *binFill    = psVectorAllocEmpty(nMax, PS_TYPE_F32); // fraction of radial bin with valid pixels
    5354
    5455    psVectorInit (binSB, 0.0);
    … …  
    144145            binSB->data.F32[nOut] = value;
    145146            binSBstdev->data.F32[nOut] = sqrt(PS_SQR(dvalue) / values->n + skyModelErrorSQ);
     147            binFill->data.F32[nOut] = values->n / binArea->data.F32[nOut];
    146148
    147149            // error in the SB is the stdev per bin / sqrt (number of pixels)
    … …  
    176178        psFree(binRad);
    177179        psFree(binArea);
     180        psFree(binFill);
    178181        psFree(radMin);
    179182        psFree(radMax);
    … …  
    202205    psFree(profile->radialBins);
    203206    psFree(profile->area);
     207    psFree(profile->binFill);
    204208
    205209    // save the vectors
    206210    profile->radialBins = binRad;
    207211    profile->area       = binArea;
     212    profile->binFill    = binFill;
    208213    profile->binSB      = binSB;
    209214    profile->binSBstdev = binSBstdev;

  • branches/eam_branches/ipp-20110213/psphot/src/psphotPetrosianStats.c

    r28013 r30707  
    66// generate the Petrosian radius and flux from the mean surface brightness (r_i)
    77
    8 float InterpolateValues (float X0, float Y0, float X1, float Y1, float X);
     8float InterpolateValues     (float X0, float Y0, float X1, float Y1, float X);
     9float InterpolateValuesErrX (float X0, float Y0, float X1, float Y1, float X, float dX0, float dX1);
     10float InterpolateValuesErrY (float X0, float Y0, float X1, float Y1, float X, float dY0, float dY1);
    911
    1012bool psphotPetrosianStats (pmSource *source) {
    … …  
    2527    psVector *binRad     = profile->radialBins;
    2628    psVector *area       = profile->area;
     29    psVector *binFill    = profile->binFill;
    2730
    2831    psVector *fluxSum     = psVectorAllocEmpty(binSB->n, PS_TYPE_F32);
    … …  
    3336    psVector *meanSB      = psVectorAllocEmpty(binSB->n, PS_TYPE_F32);
    3437    psVector *areaSum     = psVectorAllocEmpty(binSB->n, PS_TYPE_F32);
     38    psVector *apixSum     = psVectorAllocEmpty(binSB->n, PS_TYPE_F32);
    3539
    3640    float petRadius = NAN;
     41    float petRadiusErr = NAN;
    3742    float petFlux = NAN;
     43    float petFluxErr = NAN;
     44    float petArea = NAN;
     45    float petApix = NAN;
    3846
    3947    bool anyPetro = false;
    … …  
    4149    bool above = true;
    4250    float Asum = 0.0;
     51    float Psum = 0.0;
    4352    float Fsum = 0.0;
    4453    float dFsum2 = 0.0;
    … …  
    5665
    5766        float Area = area->data.F32[i];
    58         Asum += Area;
     67        Asum += Area;                   // Asum is the cumulative area interior to this bin
    5968        Fsum += binSB->data.F32[i] * Area;
    6069        dFsum2 += PS_SQR(binSBstdev->data.F32[i] * Area);
     70        Psum += Area*binFill->data.F32[i]; // Psum is the cumulative number of pixels interior to this bin
    6171
    6272        float areaInner = 0.5 * Area;
    … …  
    8797
    8898        psVectorAppend(areaSum, Asum);
     99        psVectorAppend(apixSum, Psum);
    89100        psVectorAppend(fluxSum, Fsum);
    90101        psVectorAppend(fluxSumErr2, dFsum2);
    91102        psVectorAppend(refRadius, binRad->data.F32[i]);
    92103
    93         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",
     104        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",
    94105                 i, refRadius->data.F32[nOut],
    95106                 binSB->data.F32[i], binSBstdev->data.F32[i],
    96107                 meanSB->data.F32[nOut], meanSBerr,
    97108                 petRatio->data.F32[nOut], petRatioErr->data.F32[nOut],
    98                  fluxSum->data.F32[nOut], sqrt(fluxSumErr2->data.F32[nOut]), areaSum->data.F32[nOut], areaInner);
     109                 fluxSum->data.F32[nOut], sqrt(fluxSumErr2->data.F32[nOut]), areaSum->data.F32[nOut], apixSum->data.F32[nOut], areaInner);
    99110   
    100111        // anytime we transition below the PETROSIAN_RATIO, calculate the radius and flux
    … …  
    104115            if (i == 0) {
    105116                // assume Fmax @ R = 0.0
    106                 petRadius = InterpolateValues (1.0, 0.0, petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO);
    107             } else {
    108                 petRadius = InterpolateValues (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO);
     117                petRadius    = InterpolateValues     (1.0, 0.0, petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO);
     118                petRadiusErr = InterpolateValuesErrX (1.0, 0.0, petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO, 0.0, petRatioErr->data.F32[nOut]);
     119            } else {
     120                petRadius    = InterpolateValues     (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO);
     121                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]);
    109122            }
    110123            above = false;
    … …  
    128141    }
    129142
     143    // if we failed to reach the PETROSIAN_RATIO, use the lowest significant ratio instead (flag this!)
    130144    if (!anyPetro) {
    131145        // interpolate Rvec between i-1 and i to PETROSIAN_RATIO to get flux (Fvec) and radius (rvec)
    132146        if (lowestSignificantRadius == 0) {
    133147            // assume Fmax @ R = 0.0
    134             petRadius = InterpolateValues (1.0, 0.0, petRatio->data.F32[lowestSignificantRadius], refRadius->data.F32[lowestSignificantRadius], PETROSIAN_RATIO);
     148            petRadius    = InterpolateValues     (1.0, 0.0, petRatio->data.F32[lowestSignificantRadius], refRadius->data.F32[lowestSignificantRadius], PETROSIAN_RATIO);
     149            petRadiusErr = InterpolateValuesErrX (1.0, 0.0, petRatio->data.F32[lowestSignificantRadius], refRadius->data.F32[lowestSignificantRadius], PETROSIAN_RATIO, 0.0, petRatioErr->data.F32[lowestSignificantRadius]);
     150
    135151        } else {
    136             petRadius = InterpolateValues (petRatio->data.F32[lowestSignificantRadius-1], refRadius->data.F32[lowestSignificantRadius-1], petRatio->data.F32[lowestSignificantRadius], refRadius->data.F32[lowestSignificantRadius], PETROSIAN_RATIO);
     152            int n0 = lowestSignificantRadius-1;
     153            int n1 = lowestSignificantRadius;
     154            petRadius    = InterpolateValues     (petRatio->data.F32[n0], refRadius->data.F32[n0], petRatio->data.F32[n1], refRadius->data.F32[n1], PETROSIAN_RATIO);
     155            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]);
    137156        }
    138157    }
    … …  
    147166                continue;
    148167            } else {
    149                 petFlux = InterpolateValues (refRadius->data.F32[i-1], fluxSum->data.F32[i-1], refRadius->data.F32[i], fluxSum->data.F32[i], apRadius);
     168                petFlux    = InterpolateValues     (refRadius->data.F32[i-1], fluxSum->data.F32[i-1], refRadius->data.F32[i], fluxSum->data.F32[i], apRadius);
     169                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]));
     170                petArea    = InterpolateValues     (refRadius->data.F32[i-1], areaSum->data.F32[i-1], refRadius->data.F32[i], areaSum->data.F32[i], apRadius);
     171                petApix    = InterpolateValues     (refRadius->data.F32[i-1], apixSum->data.F32[i-1], refRadius->data.F32[i], apixSum->data.F32[i], apRadius);
    150172                break;
    151173            }
    … …  
    158180    float R50 = NAN;
    159181    float R90 = NAN;
     182    float R50err = NAN;
     183    float R90err = NAN;
    160184    bool found50 = false;
    161185    bool found90 = false;
    … …  
    167191                continue;
    168192            } else {
    169                 R50 = InterpolateValues (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux50);
     193                R50    = InterpolateValues     (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux50);
     194                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]));
    170195                found50 = true;
    171196            }
    … …  
    176201                continue;
    177202            } else {
    178                 R90 = InterpolateValues (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux90);
     203                R90    = InterpolateValues     (fluxSum->data.F32[i-1], refRadius->data.F32[i-1], fluxSum->data.F32[i], refRadius->data.F32[i], flux90);
     204                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]));
    179205                found90 = true;
    180206            }
    … …  
    188214    source->extpars->petrosianR50    = R50;
    189215    source->extpars->petrosianR90    = R90;
     216    source->extpars->petrosianFill      = petApix / petArea;
    190217   
    191218    // XXX add the errors
    192     source->extpars->petrosianRadiusErr = NAN;
    193     source->extpars->petrosianFluxErr   = NAN;
    194     source->extpars->petrosianR50Err    = NAN;
    195     source->extpars->petrosianR90Err    = NAN;
     219    source->extpars->petrosianRadiusErr = petRadiusErr;
     220    source->extpars->petrosianFluxErr   = petFluxErr;
     221    source->extpars->petrosianR50Err    = R50err;
     222    source->extpars->petrosianR90Err    = R90err;
    196223
    197224    // fprintf (stderr, "source @ %f,%f\n", source->peak->xf, source->peak->yf);
    … …  
    205232    psFree(meanSB);
    206233    psFree(areaSum);
     234    psFree(apixSum);
    207235
    208236    return true;
    … …  
    210238
    211239float InterpolateValues (float X0, float Y0, float X1, float Y1, float X) {
    212     float Y = Y0 + (Y1 - Y0) * (X - X0) / (X1 - X0);
     240    float dydx = (Y1 - Y0) / (X1 - X0);
     241    float Y = Y0 + dydx * (X - X0);
    213242    return Y;
    214243}
    215244
     245float InterpolateValuesErrX (float X0, float Y0, float X1, float Y1, float X, float dX0, float dX1) {
     246
     247    float dydx = (Y1 - Y0) / (X1 - X0);
     248    float dxdx = (X  - X0) / (X1 - X0);
     249   
     250    float dY = sqrt(PS_SQR(dX1*dydx*dxdx) + PS_SQR(dX0*dydx*(dxdx - 1.0)));
     251    return dY;
     252}
     253
     254float InterpolateValuesErrY (float X0, float Y0, float X1, float Y1, float X, float dY0, float dY1) {
     255
     256    float dxdx = (X  - X0) / (X1 - X0);
     257   
     258    float dY = sqrt(PS_SQR(dY1*dxdx) + PS_SQR(dY0*(1.0 - dxdx)));
     259    return dY;
     260}
     261

  • branches/eam_branches/ipp-20110213/psphot/src/psphotRadialApertures.c

    r30624 r30707  
    3737    int Nradial = 0;
    3838
     39    // perform full non-linear fits / extended source analysis?
     40    if (!psMetadataLookupBool (&status, recipe, "RADIAL_APERTURES")) {
     41        psLogMsg ("psphot", PS_LOG_INFO, "skipping radial apertures\n");
     42        return true;
     43    }
     44
    3945    psTimerStart ("psphot.radial");
    4046
    … …  
    6268    psAssert (radMax, "annular bins (RADIAL.ANNULAR.BINS.UPPER) are not defined in the recipe");
    6369    psAssert (radMax->n, "no valid annular bins (RADIAL.ANNULAR.BINS.UPPER) are define");
     70    float outerRadius = radMax->data.F32[radMax->n - 1];
    6471
    6572    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    6673    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
    6774    assert (maskVal);
    68 
    69     // XXX temporary user-supplied systematic sky noise measurement (derive from background model)
    70     float skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE");
    7175
    7276    // S/N limit to perform full non-linear fits
    … …  
    112116            pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    113117        }
     118
     119        // we need to change the view for the radial aperture analysis, but we want to recover exactly
     120        // the original view; the following elements get destroyed by pmSourceRedefinePixels so save them:
     121        psImage *oldMaskObj   = psMemIncrRefCounter(source->maskObj);
     122        psImage *oldModelFlux = psMemIncrRefCounter(source->modelFlux);
     123        psImage *oldPSFimage  = psMemIncrRefCounter(source->psfImage);
     124        psRegion oldRegion    = source->region;
     125
    114126        Nradial ++;
    115127
    116128        // force source image to be a bit larger...
    117         float radius = source->peak->xf - source->pixels->col0;
    118         radius = PS_MAX (radius, source->peak->yf - source->pixels->row0);
    119         radius = PS_MAX (radius, source->pixels->numRows - source->peak->yf + source->pixels->row0);
    120         radius = PS_MAX (radius, source->pixels->numCols - source->peak->xf + source->pixels->col0);
    121         pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, 1.5*radius);
    122 
    123         if (!psphotRadialApertureSource (source, recipe, skynoise, maskVal, radMax, 0)) {
     129        pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, outerRadius + 2);
     130
     131        if (!psphotRadialApertureSource (source, recipe, maskVal, radMax, 0)) {
    124132            psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
    125133        } else {
    … …  
    127135        }
    128136
     137        pmSourceRedefinePixelsByRegion (source, readout, oldRegion);
     138        psFree(source->maskObj);   source->maskObj   = oldMaskObj;
     139        psFree(source->modelFlux); source->modelFlux = oldModelFlux;
     140        psFree(source->psfImage);  source->psfImage  = oldPSFimage;
     141       
    129142        // re-subtract the object, leave local sky
    130143        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    … …  
    135148}
    136149
    137 bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *aperRadii, int entry) {
     150bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, psImageMaskType maskVal, const psVector *aperRadii, int entry) {
    138151
    139152    // if we are a child source, save the results to the parent source radial aperture array
    … …  
    203216    psVector *flux    = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin
    204217    psVector *fluxErr = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin
     218    psVector *fluxStd = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin
    205219    psVector *fill    = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin
    206220
    207221    psVectorInit (flux,    0.0);
     222    psVectorInit (fluxStd, 0.0);
    208223    psVectorInit (fluxErr, 0.0);
    209224    psVectorInit (fill,    0.0);
    … …  
    215230        for (int j = 0; (*rPix2 < *aRad2) && (j < aperRadii2->n); j++, aRad2++) {
    216231            flux->data.F32[j]    += pixFlux->data.F32[i];
     232            fluxStd->data.F32[j] += PS_SQR(pixFlux->data.F32[i]);
    217233            fluxErr->data.F32[j] += pixVar->data.F32[i];
    218234            fill->data.F32[j]    += 1.0;
    219235        }
    220236    }
     237
     238    /* for each radial bin, R(i), we measure:
     239       1) the flux within that aperture: F(i) = \sum_{r_j<R_i}(F_j)
     240       2) the fractional fill factor (count of valid pixels / effective area of the aperture
     241       3) the error on the flux within that aperture
     242     */
    221243
    222244    for (int i = 0; i < flux->n; i++) {
    223245        // calculate the total flux for bin 'nOut'
    224246        float Area = M_PI*aperRadii2->data.F32[i];
    225         fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]);
     247
     248        int nPix = fill->data.F32[i];
     249        float SBmean = flux->data.F32[i] / nPix;
     250        float SBstdv = sqrt((fluxStd->data.F32[i] / nPix) - PS_SQR(SBmean));
     251
     252        flux->data.F32[i]    = SBmean * Area;
     253        fluxStd->data.F32[i] = SBstdv * Area;
     254        fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]) * Area / nPix;
     255
    226256        fill->data.F32[i] /= Area;
    227257        psTrace ("psphot", 5, "radial bins: %3d  %5.1f : %8.1f +/- %7.2f : %4.2f %6.1f\n",
    … …  
    230260   
    231261    radialAper->flux = flux;
     262    radialAper->fluxStdev = fluxStd;
    232263    radialAper->fluxErr = fluxErr;
    233264    radialAper->fill = fill;
    … …  
    245276static int nPix = 0;
    246277
    247 bool psphotRadialApertureSource_With_Sort (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *radMax, int entry) {
     278bool psphotRadialApertureSource_With_Sort (pmSource *source, psMetadata *recipe, psImageMaskType maskVal, const psVector *radMax, int entry) {
    248279
    249280    psAssert(source->radialAper->data[entry] == NULL, "why is this already defined?");

  • branches/eam_branches/ipp-20110213/psphot/src/psphotRadialAperturesByObject.c

    r30624 r30707  
    4242    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
    4343    assert (maskVal);
    44 
    45     // XXX temporary user-supplied systematic sky noise measurement (derive from background model)
    46     float skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE");
    4744
    4845    // S/N limit to perform full non-linear fits
    … …  
    147144
    148145            // force source image to be a bit larger...
    149             // float radius = source->peak->xf - source->pixels->col0;
    150             // radius = PS_MAX (radius, source->peak->yf - source->pixels->row0);
    151             // radius = PS_MAX (radius, source->pixels->numRows - source->peak->yf + source->pixels->row0);
    152             // radius = PS_MAX (radius, source->pixels->numCols - source->peak->xf + source->pixels->col0);
    153146            pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, outerRadius + 2);
    154147
    155             if (!psphotRadialApertureSource (source, recipe, skynoise, maskVal, radMax, nMatchedPSF)) {
     148            if (!psphotRadialApertureSource (source, recipe, maskVal, radMax, nMatchedPSF)) {
    156149                psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
    157150            } else {

  • branches/eam_branches/ipp-20110213/psphot/src/psphotReadout.c

    r30624 r30707  
    203203    psphotExtendedSourceAnalysis (config, view, filerule); // pass 1 (detections->allSources)
    204204    psphotExtendedSourceFits (config, view, filerule); // pass 1 (detections->allSources)
     205    psphotRadialApertures(config, view, filerule);
    205206
    206207finish:

  • branches/eam_branches/ipp-20110213/psphot/src/psphotStackMatchPSFsUtils.c

    r30624 r30707  
    363363
    364364        // we need to register the FWHM values for use downstream
    365         pmSubtractionSetFWHMs(options->inputSeeing->data.F32[index], options->targetSeeing);
     365        pmSubtractionSetFWHMs(options->targetSeeing, options->inputSeeing->data.F32[index]);
    366366
    367367        pmSubtractionParamScaleOptions(scale, scaleRef, scaleMin, scaleMax);
Note: See TracChangeset for help on using the changeset viewer.