IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 25275


Ignore:
Timestamp:
Sep 4, 2009, 2:18:54 PM (17 years ago)
Author:
eugene
Message:

updates to petrosian analysis

Location:
branches/eam_branches/20090715/psphot/src
Files:
1 added
10 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/20090715/psphot/src/Makefile.am

    r25178 r25275  
    2525libpsphot_la_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
    2626
    27 # bin_PROGRAMS = psphot psphotTest psphotMomentsStudy psphotPetrosianStudy
    28 bin_PROGRAMS = psphotPetrosianStudy
     27bin_PROGRAMS = psphot psphotTest psphotMomentsStudy psphotPetrosianStudy
     28# bin_PROGRAMS = psphotPetrosianStudy
     29# bin_PROGRAMS = psphot
    2930
    3031psphot_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
    … …  
    6869
    6970psphotPetrosianStudy_SOURCES = \
    70         psphotPetrosianStudy.c \
    71         psphotPetrosianProfile.c \
    72         psphotRadialProfileByAngle.c \
    73         psphotRadiiFromProfiles.c \
    74         psphotEllipticalContour.c \
    75         psphotEllipticalProfile.c \
    76         psphotPetrosianRadialBins.c \
    77         psphotPetrosianVisual.c \
    78         psphotPetrosianStats.c
     71        psphotPetrosianStudy.c
    7972
    8073libpsphot_la_SOURCES = \
    … …  
    145138        psphotThreadTools.c            \
    146139        psphotAddNoise.c               \
     140        psphotPetrosianProfile.c       \
     141        psphotRadialProfileByAngle.c   \
     142        psphotRadiiFromProfiles.c      \
     143        psphotEllipticalContour.c      \
     144        psphotEllipticalProfile.c      \
     145        psphotPetrosianRadialBins.c    \
     146        psphotPetrosianVisual.c        \
     147        psphotPetrosianStats.c         \
     148        psphotPetrosianAnalysis.c      \
    147149        pmPetrosian.c
     150
    148151
    149152# dropped? psphotGrowthCurve.c

  • branches/eam_branches/20090715/psphot/src/pmPetrosian.h

    r25178 r25275  
    3535bool psphotPetrosianFreeVectors(pmPetrosian *petrosian);
    3636
    37 bool psphotPetrosianProfile (pmSource *source);
     37bool psphotPetrosianProfile (pmSource *source, float skynoise);
    3838bool psphotRadialProfilesByAngles (pmPetrosian *petro, pmSource *source, int Nsec, float Rmax);
    3939float psphotRadiusFromProfile (psVector *radius, psVector *flux, float fluxMin, float fluxMax);
    … …  
    4141bool psphotEllipticalProfile (pmSource *source, pmPetrosian *petrosian);
    4242bool psphotEllipticalContour (pmPetrosian *petrosian);
    43 bool psphotPetrosianRadialBins (pmSource *source, pmPetrosian *petrosian, float radiusMax);
     43bool psphotPetrosianRadialBins (pmSource *source, pmPetrosian *petrosian, float radiusMax, float skynoise);
    4444bool psphotPetrosianStats (pmPetrosian *petrosian);
    4545
    … …  
    4848
    4949bool psphotPetrosianVisualProfileByAngle (psVector *radius, psVector *flux);
    50 bool psphotPetrosianVisualProfileRadii (psVector *radius, psVector *flux, psVector *radiusBin, psVector *fluxBin, float RadiusRef);
     50bool psphotPetrosianVisualProfileRadii (psVector *radius, psVector *flux, psVector *radiusBin, psVector *fluxBin, float peakFlux, float RadiusRef);
    5151bool psphotPetrosianVisualEllipticalContour (pmPetrosian *petrosian);
    5252
    … …  
    5757                                 float petFlux, float radiusForFlux);
    5858
    59 bool pmVisualLimitsFromVectors (Graphdata *graphdata, psVector *xVec, psVector *yVec);
    60 
    6159/// @}
    6260

  • branches/eam_branches/20090715/psphot/src/psphot.h

    r25022 r25275  
    88#include <psmodules.h>
    99#include "psphotErrorCodes.h"
     10#include "pmPetrosian.h"
    1011
    1112#define PSPHOT_RECIPE "PSPHOT" // Name of the recipe to use
    … …  
    192193bool psphotVisualPlotSourceSize (psArray *sources);
    193194
     195bool psphotVisualShowPetrosian (pmSource *source, pmPetrosian *petrosian);
     196bool psphotPetrosianAnalysis (pmReadout *readout, psArray *sources, psMetadata *recipe);
     197
    194198bool psphotImageQuality (psMetadata *recipe, psArray *sources);
    195199

  • branches/eam_branches/20090715/psphot/src/psphotInternal.h

    r25105 r25275  
    1414#include <psmodules.h>
    1515#include "psphot.h"
    16 #include "pmPetrosian.h"
    1716
    1817#endif

  • branches/eam_branches/20090715/psphot/src/psphotPetrosianProfile.c

    r25178 r25275  
    77// structure to something the pmRadialProfile
    88
    9 bool psphotPetrosianProfile (pmSource *source) {
     9bool psphotPetrosianProfile (pmSource *source, float skynoise) {
    1010
    1111  // container to hold results from the radial profile analysis
    … …  
    4848  // integrate the radial profile for radial bins defined for the petrosian measurement:
    4949  // SB_i (r_i) where \alpha r_i < r < \beta r_i
    50   if (!psphotPetrosianRadialBins (source, petrosian, Rmax)) {
     50  if (!psphotPetrosianRadialBins (source, petrosian, Rmax, skynoise)) {
    5151    psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile");
    5252    return false;
    … …  
    5959  }
    6060 
     61  // XXX this will only work in the psphot context, not the psphotPetrosianStudy...
     62  // XXX add the petrosian to the pmSource structure...
     63  psphotVisualShowPetrosian (source, petrosian);
     64
    6165  psphotPetrosianFreeVectors(petrosian);
    6266

  • branches/eam_branches/20090715/psphot/src/psphotPetrosianRadialBins.c

    r25178 r25275  
    1313// track the non-overlapping radius values.
    1414
    15 bool psphotPetrosianRadialBins (pmSource *source, pmPetrosian *petrosian, float radiusMax) {
     15bool psphotPetrosianRadialBins (pmSource *source, pmPetrosian *petrosian, float radiusMax, float skynoise) {
     16
     17    // XXX for testing, let's just set this to a value
     18   
     19    float skyModelErrorSQ = PS_SQR(skynoise);
     20
     21# if (0)
     22    // if backModel or backStdev are missing, the values of sky and/or skyErr will be set to NAN
     23    pmReadout *backModel = psphotSelectBackground (config, view);
     24    pmReadout *backStdev = psphotSelectBackgroundStdev (config, view);
     25# endif
    1626
    1727    psVector *radius = petrosian->radiusElliptical;
    … …  
    104114    psVector *values = psVectorAllocEmpty (flux->n, PS_TYPE_F32);
    105115    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     116    // psStats *stats = psStatsAlloc(PS_STAT_FITTED_MEAN_V4 | PS_STAT_FITTED_STDEV_V4);
    106117
    107118    // integrate flux, radius for each of these bins.  since flux is sorted by radius,
    … …  
    123134            // binSB->data.F32[nOut] = stats->sampleMedian;
    124135            binSB->data.F32[nOut] = stats->robustMedian;
    125             binSBstdev->data.F32[nOut] = stats->robustStdev / sqrt(values->n);
     136            binSBstdev->data.F32[nOut] = sqrt(PS_SQR(stats->robustStdev) / values->n + skyModelErrorSQ);
     137            // binSB->data.F32[nOut] = stats->fittedMean;
     138            // binSBstdev->data.F32[nOut] = sqrt(PS_SQR(stats->fittedStdev) / values->n + skyModelErrorSQ);
     139
     140            // error in the SB is the stdev per bin / sqrt (number of pixels)
     141            // added in quadrature to a fraction of the local sky (not the
     142            // residual flux, but the sky from the sky model)
    126143
    127144            if (1) {
    … …  
    153170    petrosian->binSBstdev = binSBstdev;
    154171
    155     psphotPetrosianVisualProfileRadii (radius, flux, binRad, binSB, 0.0);
     172    psphotPetrosianVisualProfileRadii (radius, flux, binRad, binSB, source->peak->flux, 0.0);
    156173
    157174    psFree(radMin);

  • branches/eam_branches/20090715/psphot/src/psphotPetrosianStats.c

    r25204 r25275  
    3232    float dFsum2 = 0.0;
    3333
     34    float nSigma = 3.0;
     35    int lowestSignificantRadius = 0;
     36    float lowestSignificantRatio = 1.0;
     37
    3438    int nOut = 0;
    3539    for (int i = 0; i < binSB->n; i++) {
    … …  
    4852            areaInner += areaSum->data.F32[nOut-1];
    4953            fluxInner += fluxSum->data.F32[nOut-1];
    50             fluxInnerErr2 = PS_SQR(fluxSumErr2->data.F32[nOut-1] * areaSum->data.F32[nOut-1]);
     54            fluxInnerErr2 += fluxSumErr2->data.F32[nOut-1];
    5155        }
    5256
    … …  
    6468
    6569        float meanSBerr = sqrt(fluxInnerErr2) / areaInner;
    66         float ratioErr = ratio * sqrt(PS_SQR(binSBstdev->data.F32[i]/binSB->data.F32[i]) + PS_SQR(meanSBerr/meanSB->data.F32[nOut]));
     70        float ratioErr = fabs(ratio) * sqrt(PS_SQR(binSBstdev->data.F32[i]/binSB->data.F32[i]) + PS_SQR(meanSBerr/meanSB->data.F32[nOut]));
    6771
    6872        psVectorAppend(petRatioErr, ratioErr);
    … …  
    7478
    7579        if (1) {
    76             fprintf (stderr, "%3d : %5.2f : %5.3f  %5.3f : %5.3f %5.3f : %5.1f %5.1f\n",
    77                      i, refRadius->data.F32[nOut], binSB->data.F32[i], meanSB->data.F32[nOut], petRatio->data.F32[nOut], petRatioErr->data.F32[nOut], fluxSum->data.F32[nOut], areaSum->data.F32[nOut]);
     80            fprintf (stderr, "%3d : %5.2f : %5.3f %5.3f : %5.3f %5.3f : %5.3f %5.3f : %5.3f %5.3f : %5.1f %5.1f\n",
     81                     i, refRadius->data.F32[nOut],
     82                     binSB->data.F32[i], binSBstdev->data.F32[i],
     83                     meanSB->data.F32[nOut], meanSBerr,
     84                     petRatio->data.F32[nOut], petRatioErr->data.F32[nOut],
     85                     fluxSum->data.F32[nOut], sqrt(fluxSumErr2->data.F32[nOut]), areaSum->data.F32[nOut], areaInner);
    7886        }
    7987   
    8088        // anytime we transition below the PETROSIAN_RATIO, calculate the radius and flux
    8189        // we will keep and report the last (largest radius) value
    82         if (above && (petRatio->data.F32[nOut] < PETROSIAN_RATIO)) {
     90        if (above && (petRatio->data.F32[nOut] < PETROSIAN_RATIO) && (petRatio->data.F32[nOut] > nSigma*petRatioErr->data.F32[nOut])) {
    8391            // interpolate Rvec between i-1 and i to PETROSIAN_RATIO to get flux (Fvec) and radius (rvec)
    8492            if (i == 0) {
    … …  
    93101        }
    94102   
     103        // anytime we transition below the PETROSIAN_RATIO, calculate the radius and flux
     104        // we will keep and report the last (largest radius) value
     105        // find the last signficant measurement of the petrosian ratio
     106        if (above && (petRatio->data.F32[nOut] < lowestSignificantRatio) && (petRatio->data.F32[nOut] > nSigma*petRatioErr->data.F32[nOut])) {
     107            lowestSignificantRadius = nOut;
     108            lowestSignificantRatio = petRatio->data.F32[nOut];
     109        }
     110   
    95111        // reset on transitions up, but do not re-calculate rad_90, flux_90
    96112        if (!above && (petRatio->data.F32[nOut] >= PETROSIAN_RATIO)) {
    … …  
    101117
    102118    if (!anyPetro) {
    103         // set default radius:
     119        // interpolate Rvec between i-1 and i to PETROSIAN_RATIO to get flux (Fvec) and radius (rvec)
     120        if (lowestSignificantRadius == 0) {
     121            // assume Fmax @ R = 0.0
     122            petRadius = InterpolateValues (1.0, 0.0, petRatio->data.F32[lowestSignificantRadius], refRadius->data.F32[lowestSignificantRadius], PETROSIAN_RATIO);
     123        } else {
     124            petRadius = InterpolateValues (petRatio->data.F32[lowestSignificantRadius-1], refRadius->data.F32[lowestSignificantRadius-1], petRatio->data.F32[lowestSignificantRadius], refRadius->data.F32[lowestSignificantRadius], PETROSIAN_RATIO);
     125        }
    104126    }
    105127
    … …  
    110132        if (refRadius->data.F32[i] > apRadius) {
    111133            if (i == 0) {
    112                 psAbort ("does this case make any sense?");
     134                psWarning ("does this case make any sense? (refRadius[0] > apRadius)");
     135                continue;
    113136            } else {
    114137                petFlux = InterpolateValues (refRadius->data.F32[i-1], fluxSum->data.F32[i-1], refRadius->data.F32[i], fluxSum->data.F32[i], apRadius);
    … …  
    119142
    120143    // save petRadius, petFlux
     144    // XXX save flags (anyPetro, manyPetro)
    121145    petrosian->petrosianRadius = petRadius;
    122146    petrosian->petrosianFlux   = petFlux;

  • branches/eam_branches/20090715/psphot/src/psphotPetrosianStudy.c

    r25178 r25275  
    136136  }
    137137
    138   psphotPetrosianProfile (source);
     138  psphotPetrosianProfile (source, skynoise);
    139139
    140140  psFree (source);

  • branches/eam_branches/20090715/psphot/src/psphotPetrosianVisual.c

    r25178 r25275  
    9393}
    9494
    95 bool psphotPetrosianVisualProfileRadii (psVector *radius, psVector *flux, psVector *radiusBin, psVector *fluxBin, float RadiusRef) {
     95bool psphotPetrosianVisualProfileRadii (psVector *radius, psVector *flux, psVector *radiusBin, psVector *fluxBin, float peakFlux, float RadiusRef) {
    9696
    9797    float FluxRef = 500.0;
    … …  
    115115    KapaSetFont (kapa, "courier", 14);
    116116
    117     pmVisualLimitsFromVectors (&graphdata, radius, flux);
     117    graphdata.ymax = +1.05*peakFlux;
     118    graphdata.ymin = -0.05*peakFlux;
     119    pmVisualLimitsFromVectors (&graphdata, radius, NULL);
    118120    KapaSetLimits (kapa, &graphdata);
    119121

  • branches/eam_branches/20090715/psphot/src/psphotRadiiFromProfiles.c

    r25178 r25275  
    5656        }
    5757    }
     58    Rbin = 3;
    5859
    5960    psVector *fluxBinned = NULL;
    … …  
    117118            // XXX is there a macro in psLib that does this interpolation?
    118119            if (i == 0) {
    119               psAbort ("impossible condition f[0] < Fo");
     120                psWarning ("impossible condition f[0] < Fo");
     121                continue;
    120122            } else {
    121123                Ro = radiusBinned->data.F32[i-1] + (radiusBinned->data.F32[i] - radiusBinned->data.F32[i-1]) * (Fo - fluxBinned->data.F32[i-1]) / (fluxBinned->data.F32[i] - fluxBinned->data.F32[i-1]);
Note: See TracChangeset for help on using the changeset viewer.