IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 25178


Ignore:
Timestamp:
Aug 24, 2009, 8:40:34 AM (17 years ago)
Author:
eugene
Message:

various petrosian / radial profile analysis improvements

Location:
branches/eam_branches/20090715/psphot/src
Files:
11 edited

Legend:

Unmodified
Added
Removed

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

    r25105 r25178  
    1717# Force recompilation of psphotVersion.c, since it gets the version information
    1818psphotVersion.c: psphotVersionDefinitions.h
    19 psphotVersionDefinitions.h: psphotVersionDefinitions.h.in FORCE
    20         -$(RM) psphotVersionDefinitions.h
    21         $(SED) -e "s|@PSPHOT_VERSION@|\"$(PSPHOT_VERSION)\"|" -e "s|@PSPHOT_BRANCH@|\"$(PSPHOT_BRANCH)\"|" -e "s|@PSPHOT_SOURCE@|\"$(PSPHOT_SOURCE)\"|" psphotVersionDefinitions.h.in > psphotVersionDefinitions.h
    22 FORCE: ;
     19# psphotVersionDefinitions.h: psphotVersionDefinitions.h.in FORCE
     20#       -$(RM) psphotVersionDefinitions.h
     21#       $(SED) -e "s|@PSPHOT_VERSION@|\"$(PSPHOT_VERSION)\"|" -e "s|@PSPHOT_BRANCH@|\"$(PSPHOT_BRANCH)\"|" -e "s|@PSPHOT_SOURCE@|\"$(PSPHOT_SOURCE)\"|" psphotVersionDefinitions.h.in > psphotVersionDefinitions.h
     22# FORCE: ;
    2323
    2424libpsphot_la_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
    2525libpsphot_la_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
    2626
    27 bin_PROGRAMS = psphot psphotTest psphotMomentsStudy psphotPetrosianStudy
     27# bin_PROGRAMS = psphot psphotTest psphotMomentsStudy psphotPetrosianStudy
     28bin_PROGRAMS = psphotPetrosianStudy
    2829
    2930psphot_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)

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

    r25105 r25178  
    2424    psFree(petrosian->fluxElliptical);
    2525
    26     psFree(petrosian->binnedFlux);
     26    psFree(petrosian->binSB);
     27    psFree(petrosian->binSBstdev);
    2728    psFree(petrosian->radialBins);
    2829    psFree(petrosian->area);
    … …  
    4445    petrosian->radialBins = NULL;
    4546    petrosian->area = NULL;
    46     petrosian->binnedFlux = NULL;
     47    petrosian->binSB = NULL;
     48    petrosian->binSBstdev = NULL;
    4749
    4850    petrosian->petrosianRadius = NAN;
    … …  
    6466    psFree(petrosian->fluxElliptical);
    6567
    66     psFree(petrosian->binnedFlux);
     68    psFree(petrosian->binSB);
     69    psFree(petrosian->binSBstdev);
    6770    psFree(petrosian->radialBins);
    6871    psFree(petrosian->area);
    … …  
    7881    petrosian->radialBins = NULL;
    7982    petrosian->area = NULL;
    80     petrosian->binnedFlux = NULL;
     83    petrosian->binSB = NULL;
     84    petrosian->binSBstdev = NULL;
    8185   
    8286    return true;

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

    r25105 r25178  
    2020    psVector *fluxElliptical;           // flux for the above radial coordinates
    2121
    22     psVector *binnedFlux;               // mean surface brightness within radial bins
     22    psVector *binSB;                    // mean surface brightness within radial bins
     23    psVector *binSBstdev;               // scatter of mean surface brightness within radial bins
    2324    psVector *radialBins;               // radii corresponding to above binnedBlux
    2425    psVector *area;                     // differential area of the non-overlapping radial bins
    … …  
    5253bool psphotPetrosianVisualStats (psVector *radBin, psVector *fluxBin,
    5354                                 psVector *refRadius, psVector *meanSB,
    54                                  psVector *petRatio, psVector *fluxSum,
    55                                  float petRadius, float petFlux);
     55                                 psVector *petRatio, psVector *petRatioErr, psVector *fluxSum,
     56                                 float petRadius, float ratioForRadius,
     57                                 float petFlux, float radiusForFlux);
    5658
    5759bool pmVisualLimitsFromVectors (Graphdata *graphdata, psVector *xVec, psVector *yVec);

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

    r25105 r25178  
    7373    psMinimizeLMChi2(myMin, covar, params, constraint, x, y, yErr, psphotEllipticalContourFunc);
    7474
     75    /// XXX rationalize? if epsilon > 1, flip major and minor axes (rotate by 90 degrees)
     76    if (params->data.F32[PAR_EPSILON] < 1.0) {
     77        petrosian->axes.major = params->data.F32[PAR_RMIN] / params->data.F32[PAR_EPSILON];
     78        petrosian->axes.minor = params->data.F32[PAR_RMIN];
     79        petrosian->axes.theta = params->data.F32[PAR_PHI];
     80    } else {
     81        petrosian->axes.major = params->data.F32[PAR_RMIN];
     82        petrosian->axes.minor = params->data.F32[PAR_RMIN] / params->data.F32[PAR_EPSILON];
     83        petrosian->axes.theta = params->data.F32[PAR_PHI] + 0.5*M_PI;
     84    }
     85
    7586    fprintf (stderr, "# fitted values:\n");
    7687    fprintf (stderr, "Po:  %f\n", params->data.F32[PAR_PHI]*PS_DEG_RAD);
    7788    fprintf (stderr, "Ep:  %f\n", params->data.F32[PAR_EPSILON]);
    7889    fprintf (stderr, "Rm:  %f\n", params->data.F32[PAR_RMIN]);
    79 
    80     /// XXX rationalize? if epsilon > 1, flip major and minor axes (rotate by 90 degrees)
    81     petrosian->axes.major = params->data.F32[PAR_RMIN] / params->data.F32[PAR_EPSILON];
    82     petrosian->axes.minor = params->data.F32[PAR_RMIN];
    83     petrosian->axes.theta = params->data.F32[PAR_PHI];
    8490
    8591    // show the results

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

    r25128 r25178  
    6161  psphotPetrosianFreeVectors(petrosian);
    6262
    63   fprintf (stderr, "petrosian radius: %f, flux: %f, axis ratio: %f, angle: %f\n",
     63  fprintf (stdout, "\n petrosian radius: %f\n flux: %f\n axis ratio: %f\n angle: %f\n",
    6464           petrosian->petrosianRadius, petrosian->petrosianFlux, petrosian->axes.minor/petrosian->axes.major, PS_DEG_RAD*petrosian->axes.theta);
    6565

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

    r25105 r25178  
    1212// finely spaced than r_{i+1} = r_i * \beta / \alpha.  for the integration, we need to
    1313// track the non-overlapping radius values.
    14 
    15 # define PETROSIAN_ALPHA 0.8
    16 # define PETROSIAN_BETA 1.25
    1714
    1815bool psphotPetrosianRadialBins (pmSource *source, pmPetrosian *petrosian, float radiusMax) {
    … …  
    3330    psVector *radBet  = psVectorAllocEmpty(nMax, PS_TYPE_F32);
    3431
    35     psVector *fluxBin = psVectorAllocEmpty(nMax, PS_TYPE_F32);
    36     psVector *radBin  = psVectorAllocEmpty(nMax, PS_TYPE_F32);
    37     psVector *area    = psVectorAllocEmpty(nMax, PS_TYPE_F32);
     32    psVector *binSB      = psVectorAllocEmpty(nMax, PS_TYPE_F32); // surface brightness of radial bin
     33    psVector *binSBstdev = psVectorAllocEmpty(nMax, PS_TYPE_F32); // surface brightness of radial bin
     34    psVector *binRad     = psVectorAllocEmpty(nMax, PS_TYPE_F32); // mean radius of radial bin
     35    psVector *binArea    = psVectorAllocEmpty(nMax, PS_TYPE_F32); // area of radial bin (contiguous, non-overlapping)
    3836
    39     psVectorInit (fluxBin, 0.0);
    40     psVectorInit (radBin, 0.0);
     37    psVectorInit (binSB, 0.0);
     38    psVectorInit (binSBstdev, 0.0);
     39    psVectorInit (binRad, 0.0);
    4140
    4241    // generate radial bin bounds
    … …  
    5655    radBet->data.F32[2] = 2.0;
    5756   
    58     int nPts = 3;
    59     for (int i = 3; i < radiusMax; i++) {
    60         radMin->data.F32[nPts] = (i - 1);
    61         radMax->data.F32[nPts] = i;
    62         nPts++;
     57# define PETROSIAN_ALPHA 0.8
     58# define PETROSIAN_BETA 1.25
     59# define POWER_LAW_SPACING true
     60   
     61    // power-law spacing with overlapping boundaries at the geometric mid-points
     62    float rBeta = sqrt(PETROSIAN_BETA);
     63    for (int i = 3; radBet->data.F32[i-1] < radiusMax; i++) {
     64        if (POWER_LAW_SPACING) {
     65            radMin->data.F32[i] = radMax->data.F32[i-1];
     66            radMax->data.F32[i] = radMin->data.F32[i] * PETROSIAN_BETA;
     67            radAlp->data.F32[i] = radMin->data.F32[i] / rBeta;
     68            radBet->data.F32[i] = radMax->data.F32[i] * rBeta;
     69        } else {
     70            radMin->data.F32[i] = radMax->data.F32[i-1];
     71            radMax->data.F32[i] = radMin->data.F32[i] + 1;
     72            float rMid = 0.5*(radMin->data.F32[i] + radMax->data.F32[i]);
     73            radAlp->data.F32[i] = rMid * PETROSIAN_ALPHA;
     74            radBet->data.F32[i] = rMid * PETROSIAN_BETA;
     75        }
     76        radMin->n = radMax->n = radAlp->n = radBet->n = i + 1;
    6377    }
    64     radMin->n = radMax->n = radAlp->n = radBet->n = nPts;
    6578
    6679    // generate radial area-weighted mean radius & non-overlapping areas
    … …  
    7891       
    7992        // XXX calculate area-weighted radius rather than asserting?
    80         radBin->data.F32[i] = rBin;
    81         area->data.F32[i] = M_PI * (rMax2 - rMin2);
     93        binRad->data.F32[i] = rBin;
     94        binArea->data.F32[i] = M_PI * (rMax2 - rMin2);
    8295
    83         if (i > 2) {
    84             radAlp->data.F32[i] = rBin*PETROSIAN_ALPHA;
    85             radBet->data.F32[i] = rBin*PETROSIAN_BETA;
     96        if (0) {
     97            fprintf (stderr, "%3d  %5.1f %5.1f : %5.1f : %5.1f %5.1f\n",
     98                 i, radAlp->data.F32[i], radMin->data.F32[i], binRad->data.F32[i],
     99                 radMax->data.F32[i], radBet->data.F32[i]);
    86100        }
    87101    }
    … …  
    89103    // storage vector for stats
    90104    psVector *values = psVectorAllocEmpty (flux->n, PS_TYPE_F32);
    91     psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
     105    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    92106
    93107    // integrate flux, radius for each of these bins.  since flux is sorted by radius,
    … …  
    107121            // calculate the value for the nOut bin
    108122            psVectorStats (stats, values, NULL, NULL, 0);
    109             fluxBin->data.F32[nOut] = stats->sampleMedian;
     123            // binSB->data.F32[nOut] = stats->sampleMedian;
     124            binSB->data.F32[nOut] = stats->robustMedian;
     125            binSBstdev->data.F32[nOut] = stats->robustStdev / sqrt(values->n);
     126
     127            if (1) {
     128                fprintf (stderr, "%3d  %5.1f %5.1f : %5.1f  %5.2f\n",
     129                         nOut, radAlp->data.F32[nOut], radBet->data.F32[nOut], binSB->data.F32[nOut], binSBstdev->data.F32[nOut]);
     130            }
     131
    110132            nOut ++;
    111             if (nOut >= nMax) break;
     133            if (nOut >= radAlp->n) break;
    112134            Rmin = radAlp->data.F32[nOut];
    113135            Rmax = radBet->data.F32[nOut];
    … …  
    122144        psVectorAppend (values, flux->data.F32[i]);
    123145    }
    124     fluxBin->n = radBin->n = area->n = nOut;
     146    binSB->n = binSBstdev->n = binRad->n = binArea->n = nOut;
    125147    // XXX I think this misses the last radial bin -- do we care?
    126148
    127149    // save the vectors
    128     petrosian->radialBins = radBin;
    129     petrosian->area = area;
    130     petrosian->binnedFlux = fluxBin;
     150    petrosian->radialBins = binRad;
     151    petrosian->area       = binArea;
     152    petrosian->binSB      = binSB;
     153    petrosian->binSBstdev = binSBstdev;
    131154
    132     psphotPetrosianVisualProfileRadii (radius, flux, radBin, fluxBin, 0.0);
     155    psphotPetrosianVisualProfileRadii (radius, flux, binRad, binSB, 0.0);
    133156
    134157    psFree(radMin);

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

    r25105 r25178  
    11# include "psphotInternal.h"
    22
    3 # define PETROSIAN_RATIO 0.05
     3# define PETROSIAN_RATIO 0.2
     4# define PETROSIAN_RADII 2.0
    45
    56// generate the Petrosian radius and flux from the mean surface brightness (r_i)
    67
     8float InterpolateValues (float X0, float Y0, float X1, float Y1, float X);
     9
    710bool psphotPetrosianStats (pmPetrosian *petrosian) {
    811
    9   float petRadius, petFlux;
     12    float petRadius, petFlux;
    1013
    11   psVector *fluxBin = petrosian->binnedFlux;
    12   psVector *radBin  = petrosian->radialBins;
    13   psVector *area    = petrosian->area;
     14    psVector *binSB      = petrosian->binSB;
     15    psVector *binSBstdev = petrosian->binSBstdev;
     16    psVector *binRad     = petrosian->radialBins;
     17    psVector *area       = petrosian->area;
    1418
    15   psVector *fluxSum  = psVectorAllocEmpty(fluxBin->n, PS_TYPE_F32);
    16   psVector *refRadius= psVectorAllocEmpty(fluxBin->n, PS_TYPE_F32);
    17   psVector *petRatio = psVectorAllocEmpty(fluxBin->n, PS_TYPE_F32);
    18   psVector *meanSB   = psVectorAllocEmpty(fluxBin->n, PS_TYPE_F32);
    19   psVector *areaSum  = psVectorAllocEmpty(fluxBin->n, PS_TYPE_F32);
     19    psVector *fluxSum     = psVectorAllocEmpty(binSB->n, PS_TYPE_F32);
     20    psVector *refRadius   = psVectorAllocEmpty(binSB->n, PS_TYPE_F32);
     21    psVector *petRatio    = psVectorAllocEmpty(binSB->n, PS_TYPE_F32);
     22    psVector *petRatioErr = psVectorAllocEmpty(binSB->n, PS_TYPE_F32);
     23    psVector *meanSB      = psVectorAllocEmpty(binSB->n, PS_TYPE_F32);
     24    psVector *areaSum     = psVectorAllocEmpty(binSB->n, PS_TYPE_F32);
    2025
    21   bool above = true;
    22   float Fsum = 0.0;
    23   float Asum = 0.0;
    24   float Area = area->data.F32[0];
     26    bool anyPetro = false;
     27    bool manyPetro = false;
     28    bool above = true;
     29    float Fsum = 0.0;
     30    float Asum = 0.0;
    2531
    26   int nOut = 0;
    27   for (int i = 0; i < fluxBin->n; i++) {
    28     // for nan bins, we keep the area for use with the next valid bin
    29     if (!isfinite(fluxBin->data.F32[i])) {
    30       Area += area->data.F32[i];
    31       continue;
    32     }
    33     Fsum += fluxBin->data.F32[i] * Area;
    34     Asum += Area;
    35     if (i+1 < fluxBin->n) {
    36       Area = area->data.F32[i+1];
     32    int nOut = 0;
     33    for (int i = 0; i < binSB->n; i++) {
     34        // skip nan bins (do not contribute to flux or area)
     35        if (!isfinite(binSB->data.F32[i])) continue;
     36
     37        float Area = area->data.F32[i];
     38        Fsum += binSB->data.F32[i] * Area;
     39        Asum += Area;
     40
     41        float fluxInner = 0.5 * Area * binSB->data.F32[i];
     42        float areaInner = 0.5 * Area;
     43        if (nOut > 0) {
     44            fluxInner += fluxSum->data.F32[nOut-1];
     45            areaInner += areaSum->data.F32[nOut-1];
     46        }
     47
     48        psVectorAppend(meanSB, (fluxInner / areaInner));
     49        psVectorAppend(petRatio, binSB->data.F32[i] / meanSB->data.F32[nOut]);
     50        psVectorAppend(petRatioErr, binSBstdev->data.F32[i] / meanSB->data.F32[nOut]);
     51        psVectorAppend(fluxSum, Fsum);
     52        psVectorAppend(areaSum, Asum);
     53        psVectorAppend(refRadius, binRad->data.F32[i]);
     54
     55        if (1) {
     56            fprintf (stderr, "%3d : %5.2f : %5.3f  %5.3f : %5.3f %5.3f : %5.1f %5.1f\n",
     57                     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]);
     58        }
     59   
     60        // anytime we transition below the PETROSIAN_RATIO, calculate the radius and flux
     61        // we will keep and report the last (largest radius) value
     62        if (above && (petRatio->data.F32[nOut] < PETROSIAN_RATIO)) {
     63            // interpolate Rvec between i-1 and i to PETROSIAN_RATIO to get flux (Fvec) and radius (rvec)
     64            if (i == 0) {
     65                // assume Fmax @ R = 0.0
     66                petRadius = InterpolateValues (1.0, 0.0, petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO);
     67            } else {
     68                petRadius = InterpolateValues (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO);
     69            }
     70            above = false;
     71            if (anyPetro) manyPetro = true;
     72            anyPetro = true;
     73        }
     74   
     75        // reset on transitions up, but do not re-calculate rad_90, flux_90
     76        if (!above && (petRatio->data.F32[nOut] >= PETROSIAN_RATIO)) {
     77            above = true;
     78        }
     79        nOut ++;
    3780    }
    3881
    39     psVectorAppend(petRatio, Asum * fluxBin->data.F32[i] / Fsum);
    40     psVectorAppend(fluxSum, Fsum);
    41     psVectorAppend(meanSB, Fsum / Asum);
    42     psVectorAppend(areaSum, Asum);
    43     psVectorAppend(refRadius, radBin->data.F32[i]);
     82    if (!anyPetro) {
     83        // set default radius:
     84    }
    4485
    45     // anytime we transition below the PETROSIAN_RATIO, calculate the radius and flux
    46     // we will keep and report the last (largest radius) value
    47     if (above && (petRatio->data.F32[nOut] < PETROSIAN_RATIO)) {
    48       // interpolate Rvec between i-1 and i to PETROSIAN_RATIO to get flux (Fvec) and radius (rvec)
    49       if (i == 0) {
    50         // assume Fmax @ R = 0.0
    51         petRadius = radBin->data.F32[nOut]  * (PETROSIAN_RATIO - 1.0) / (petRatio->data.F32[nOut] - 1.0);
    52         petFlux   = fluxSum->data.F32[nOut] * (PETROSIAN_RATIO - 1.0) / (petRatio->data.F32[nOut] - 1.0);
    53       } else {
    54         petRadius = radBin->data.F32[nOut-1]  + (radBin->data.F32[nOut]  - radBin->data.F32[nOut-1])  * (PETROSIAN_RATIO - petRatio->data.F32[nOut-1]) / (petRatio->data.F32[nOut] - petRatio->data.F32[nOut-1]);
    55         petFlux   = fluxSum->data.F32[nOut-1] + (fluxSum->data.F32[nOut] - fluxSum->data.F32[nOut-1]) * (PETROSIAN_RATIO - petRatio->data.F32[nOut-1]) / (petRatio->data.F32[nOut] - petRatio->data.F32[nOut-1]);
    56       }
    57       above = false;
     86    // now measure the flux within PETROSIAN_RADII * petRadius
     87    float apRadius = PETROSIAN_RADII * petRadius;
     88    for (int i = 0; i < refRadius->n; i++) {
     89        // XXX use bisection to do this faster:
     90        if (refRadius->data.F32[i] > apRadius) {
     91            if (i == 0) {
     92                psAbort ("does this case make any sense?");
     93            } else {
     94                petFlux = InterpolateValues (refRadius->data.F32[i-1], fluxSum->data.F32[i-1], refRadius->data.F32[i], fluxSum->data.F32[i], apRadius);
     95                break;
     96            }
     97        }
    5898    }
    59    
    60     // reset on transitions up, but do not re-calculate rad_90, flux_90
    61     if (!above && (petRatio->data.F32[nOut] >= PETROSIAN_RATIO)) {
    62       above = true;
    63     }
    64     nOut ++;
    65   }
    6699
    67   // save petRadius, petFlux
    68   petrosian->petrosianRadius = petRadius;
    69   petrosian->petrosianFlux   = petFlux;
     100    // save petRadius, petFlux
     101    petrosian->petrosianRadius = petRadius;
     102    petrosian->petrosianFlux   = petFlux;
    70103
    71   psphotPetrosianVisualStats (radBin, fluxBin, refRadius, meanSB, petRatio, fluxSum, petRadius, petFlux);
     104    psphotPetrosianVisualStats (binRad, binSB, refRadius, meanSB, petRatio, petRatioErr, fluxSum, petRadius, PETROSIAN_RATIO, petFlux, apRadius);
    72105
    73   psFree(fluxSum);
    74   psFree(petRatio);
    75   psFree(meanSB);
    76   psFree(areaSum);
     106    psFree(fluxSum);
     107    psFree(petRatio);
     108    psFree(meanSB);
     109    psFree(areaSum);
    77110
    78   return true;
     111    return true;
    79112}
     113
     114float InterpolateValues (float X0, float Y0, float X1, float Y1, float X) {
     115    float Y = Y0 + (Y1 - Y0) * (X - X0) / (X1 - X0);
     116    return Y;
     117}
     118

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

    r25128 r25178  
    133133      source = pmSourceAlloc();
    134134      source->peak = psphotLocalPeak(readout, Xo, Yo);
    135       pmSourceDefinePixels (source, readout, Xo, Yo, 64);
     135      pmSourceDefinePixels (source, readout, Xo, Yo, 128);
    136136  }
    137137

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

    r25128 r25178  
    161161bool psphotPetrosianVisualStats (psVector *radBin, psVector *fluxBin,
    162162                                 psVector *refRadius, psVector *meanSB,
    163                                  psVector *petRatio, psVector *fluxSum,
    164                                  float petRadius, float petFlux)
     163                                 psVector *petRatio, psVector *petRatioErr,
     164                                 psVector *fluxSum,
     165                                 float petRadius, float ratioForRadius,
     166                                 float petFlux, float radiusForFlux)
    165167{
    166168    Graphdata graphdata;
    … …  
    229231
    230232    graphdata.color = KapaColorByName ("black");
    231     pmVisualLimitsFromVectors (&graphdata, radBin, petRatio);
     233    graphdata.ymax = +1.05;
     234    graphdata.ymin = -0.05;
     235    pmVisualLimitsFromVectors (&graphdata, radBin, NULL);
    232236    KapaSetLimits (kapa, &graphdata);
    233237
    … …  
    240244    graphdata.ptype = 0;
    241245    graphdata.size = 1.0;
     246    graphdata.etype = 0x01;
    242247    KapaPrepPlot (kapa, refRadius->n, &graphdata);
    243248    KapaPlotVector (kapa, refRadius->n, refRadius->data.F32, "x");
    244249    KapaPlotVector (kapa, refRadius->n, petRatio->data.F32, "y");
    245 
    246     float ratio_90 = 0.1;
     250    KapaPlotVector (kapa, refRadius->n, petRatioErr->data.F32, "dym");
     251    KapaPlotVector (kapa, refRadius->n, petRatioErr->data.F32, "dyp");
     252    graphdata.etype = 0;
     253
    247254    graphdata.color = KapaColorByName ("red");
    248255    graphdata.style = 2;
    … …  
    251258    KapaPrepPlot   (kapa, 1, &graphdata);
    252259    KapaPlotVector (kapa, 1, &petRadius, "x");
    253     KapaPlotVector (kapa, 1, &ratio_90, "y");
     260    KapaPlotVector (kapa, 1, &ratioForRadius, "y");
    254261
    255262    // *** section 3: radius vs integrated flux
    … …  
    283290    graphdata.size = 2.0;
    284291    KapaPrepPlot   (kapa, 1, &graphdata);
    285     KapaPlotVector (kapa, 1, &petRadius, "x");
     292    KapaPlotVector (kapa, 1, &radiusForFlux, "x");
    286293    KapaPlotVector (kapa, 1, &petFlux, "y");
    287294

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

    r25128 r25178  
    44// separations
    55
    6 // XXX choices for rMax and Nsec?
     6float psphotMeanSectorValue (psImage *image, int x, int y, float dL, float dW, float theta);
     7psVector *psphotBoxValues (psImage *image, float x0, float y0, float dL, float dW, float theta);
     8psVector *psphotLineValues (psImage *image, double x1, double y1, double x2, double y2, int dW);
     9psVector *psphotLineValuesBresen (psImage *image, int X1, int Y1, int X2, int Y2, int dW, int swapcoords);
    710
    811bool psphotRadialProfilesByAngles (pmPetrosian *petrosian, pmSource *source, int Nsec, float Rmax) {
    … …  
    2831        psVector *flux   = psVectorAllocEmpty(Rmax, PS_TYPE_F32);
    2932
    30         // start at Xo,Yo and find the x,y locations for r_i, theta where r_i increments by 1 pixel
    31         // XXX at large radii ( > 10-15) use stats in a patch rather than sub-pixel interpolation
    32         for (float r = 0; r < Rmax; r += 1.0) {
     33        // Start at Xo,Yo and find the x,y locations for r_i, theta where r_i initially
     34        // increments by 1 pixel.  At large radii (r*dtheta > 2) use stats in a box rather than
     35        // sub-pixel interpolation
     36
     37        int dR = 1.0;
     38        for (float r = 0; r < Rmax; r += dR) {
    3339
    3440            float Xo = source->peak->xf + 0.5;
    … …  
    3844            float x = r * cos (theta) + Xo;
    3945            float y = r * sin (theta) + Yo;
    40 
    41             if (x < 0) continue;
    42             if (y < 0) continue;
    43             if (x >= source->pixels->parent->numCols) continue;
    44             if (y >= source->pixels->parent->numRows) continue;
    45 
    46             // value is NAN if we run off the image
    47             float value = psImageInterpolatePixelBilinear(x, y, source->pixels);
    48 
    49             // keep the nan values so all vectors are matched
    50             // if (isnan(value)) continue;
    51 
     46            dR = 2*(int)(0.5*r*sin(dtheta)) + 1;
     47
     48            if (x < 0) goto badvalue;
     49            if (y < 0) goto badvalue;
     50            if (x >= source->pixels->parent->numCols) goto badvalue;
     51            if (y >= source->pixels->parent->numRows) goto badvalue;
     52
     53            float value = NAN;
     54            if (dR < 2) {
     55                // value is NAN if we run off the image
     56                value = psImageInterpolatePixelBilinear(x, y, source->pixels);
     57            } else {
     58                value = psphotMeanSectorValue(source->pixels, x, y, dR, dR, theta);
     59            }
     60
     61            // keep the all values (even NAN) so all vectors are matched in length
    5262            psVectorAppend (radius, r);
    5363            psVectorAppend (flux, value);
     64            continue;
     65           
     66        badvalue:
     67            psVectorAppend (radius, r);
     68            psVectorAppend (flux, NAN);
    5469        }
    5570
    … …  
    103118    return true;
    104119}
     120
     121float psphotMeanSectorValue (psImage *image, int x, int y, float dL, float dW, float theta) {
     122
     123    psVector *values = psphotBoxValues (image, x, y, dL, dW, theta);
     124   
     125    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
     126    psVectorStats (stats, values, NULL, NULL, 0);
     127
     128    float value = stats->sampleMedian;
     129
     130    psFree (stats);
     131    psFree (values);
     132   
     133    return value;
     134}
     135
     136psVector *psphotBoxValues (psImage *image, float x0, float y0, float dL, float dW, float theta) {
     137
     138    // extract pixels from a series of lines (from -0.5*dW to +0.5*dW) of length dL,
     139    // centered on x0, y0 in parent pixel coordinates (not pixel indicies)
     140
     141    float xs = x0 - image->col0 - 0.5*dL*cos(theta);
     142    float ys = y0 - image->row0 - 0.5*dL*sin(theta);
     143
     144    float xe = xs + 0.5*dL*cos(theta);
     145    float ye = ys + 0.5*dL*sin(theta);
     146
     147    psVector *values = psphotLineValues (image, xs, ys, xe, ye, (int) dW);
     148    return values;
     149}
     150
     151/**
     152 * identify the quadrant and draw the correct line
     153 */
     154psVector *psphotLineValues (psImage *image, double x1, double y1, double x2, double y2, int dW) {
     155
     156  int FlipDirect, FlipCoords;
     157  int X1, Y1, X2, Y2, dX, dY;
     158
     159  /* rather than draw the line from float positions, we find the closest
     160     integer end-points and draw the line between those pixels */
     161
     162  X1 = ROUND(x1);
     163  Y1 = ROUND(y1);
     164  X2 = ROUND(x2);
     165  Y2 = ROUND(y2);
     166
     167  dX = X2 - X1;
     168  dY = Y2 - Y1;
     169
     170  FlipCoords = (abs(dX) < abs(dY));
     171  FlipDirect = FlipCoords ? (y1 > y2) : (x1 > x2);
     172
     173  psVector *values = NULL;
     174  if (!FlipDirect && !FlipCoords) values = psphotLineValuesBresen (image, X1, Y1, X2, Y2, dW, FALSE);
     175  if ( FlipDirect && !FlipCoords) values = psphotLineValuesBresen (image, X2, Y2, X1, Y1, dW, FALSE);
     176  if (!FlipDirect &&  FlipCoords) values = psphotLineValuesBresen (image, Y1, X1, Y2, X2, dW, TRUE);
     177  if ( FlipDirect &&  FlipCoords) values = psphotLineValuesBresen (image, Y2, X2, Y1, X1, dW, TRUE);
     178
     179  return values;
     180}
     181
     182/**
     183 * use the Bresenham line drawing technique
     184 * integer-only Bresenham line-draw version which is fast
     185 */
     186psVector *psphotLineValuesBresen (psImage *image, int X1, int Y1, int X2, int Y2, int dW, int swapcoords) {
     187
     188    int X, Y, dX, dY;
     189    int e, e2;
     190
     191    psVector *values = psVectorAllocEmpty(100, PS_TYPE_F32);
     192
     193    dX = X2 - X1;
     194    dY = Y2 - Y1;
     195
     196    Y = Y1;
     197    e = 0;
     198    for (X = X1; X <= X2; X++) {
     199        if (X > 0) {
     200            if (swapcoords) {
     201                if (X >= image->numRows) continue;
     202                for (int y = Y - dW; y <= Y + dW; y++) {
     203                    if (y < 0) continue;
     204                    if (y >= image->numCols) continue;
     205                    psVectorAppend(values, image->data.F32[X][y]);
     206                }
     207            } else {
     208                if (X >= image->numCols) continue;
     209                for (int y = Y - dW; y <= Y + dW; y++) {
     210                    if (y < 0) continue;
     211                    if (y >= image->numRows) continue;
     212                    psVectorAppend(values, image->data.F32[y][X]);
     213                }
     214            }
     215        }
     216        e += dY;
     217        e2 = 2 * e;
     218        if (e2 > dX) {
     219            Y++;
     220            e -= dX;
     221        }
     222        if (e2 < -dX) {
     223            Y--;
     224            e += dX;
     225        }
     226    }
     227    return values;
     228}
     229

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

    r25128 r25178  
    3131    // examine data in the two ranges Fm - Fo and Fo - Fp to define the bin size
    3232    // XXX reconsider the fractional isophote value
    33     float Fm = fluxMin + 0.25*fluxRange;
    34     float Fp = fluxMin + 0.75*fluxRange;
    35     float Fo = fluxMin + 0.50*fluxRange;
     33    float Fm = fluxMin + 0.10*fluxRange;
     34    float Fo = fluxMin + 0.25*fluxRange;
     35    float Fp = fluxMin + 0.50*fluxRange;
    3636    int Rbin = 1;
    3737     
Note: See TracChangeset for help on using the changeset viewer.