IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 25105


Ignore:
Timestamp:
Aug 18, 2009, 6:24:53 AM (17 years ago)
Author:
eugene
Message:

updates for petrosian analysis study

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

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/20090715/psphot/src

    • Property svn:ignore

      •  

        old new  
        1818psphotVersionDefinitions.h
        1919psphotMomentsStudy
         20psphotPetrosianStudy

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

    r24586 r25105  
    2525libpsphot_la_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
    2626
    27 bin_PROGRAMS = psphot psphotTest psphotMomentsStudy
     27bin_PROGRAMS = psphot psphotTest psphotMomentsStudy psphotPetrosianStudy
    2828
    2929psphot_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
    … …  
    3838psphotMomentsStudy_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
    3939psphotMomentsStudy_LDADD = libpsphot.la
     40
     41psphotPetrosianStudy_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
     42psphotPetrosianStudy_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
     43psphotPetrosianStudy_LDADD = libpsphot.la
    4044
    4145psphot_SOURCES = \
    … …  
    6165psphotMomentsStudy_SOURCES = \
    6266        psphotMomentsStudy.c
     67
     68psphotPetrosianStudy_SOURCES = \
     69        psphotPetrosianStudy.c \
     70        psphotPetrosianProfile.c \
     71        psphotRadialProfileByAngle.c \
     72        psphotRadiiFromProfiles.c \
     73        psphotEllipticalContour.c \
     74        psphotEllipticalProfile.c \
     75        psphotPetrosianRadialBins.c \
     76        psphotPetrosianVisual.c \
     77        psphotPetrosianStats.c
    6378
    6479libpsphot_la_SOURCES = \
    … …  
    128143        psphotCheckStarDistribution.c  \
    129144        psphotThreadTools.c            \
    130         psphotAddNoise.c
     145        psphotAddNoise.c               \
     146        pmPetrosian.c
    131147
    132148# dropped? psphotGrowthCurve.c
    … …  
    134150include_HEADERS = \
    135151        psphot.h \
     152        pmPetrosian.h \
    136153        psphotErrorCodes.h
    137154

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

    r25032 r25105  
    1111# include "psphotInternal.h"
    1212
    13 static void pmPetrosianFree(pmPetrosian *petro)
     13static void pmPetrosianFree(pmPetrosian *petrosian)
    1414{
    15    if (!petro) {
     15    if (!petrosian) {
    1616        return;
    17    }
    18    psFree(petrosian->radii);
    19    psFree(petrosian->fluxes);
    20    psFree(petrosian->theta);
    21    psFree(petrosian->isophotalRadii);
     17    }
     18    psFree(petrosian->radii);
     19    psFree(petrosian->fluxes);
     20    psFree(petrosian->theta);
     21    psFree(petrosian->isophotalRadii);
     22
     23    psFree(petrosian->radiusElliptical);
     24    psFree(petrosian->fluxElliptical);
     25
     26    psFree(petrosian->binnedFlux);
     27    psFree(petrosian->radialBins);
     28    psFree(petrosian->area);
    2229}
    2330
    … …  
    3239    petrosian->isophotalRadii = NULL;
    3340
    34     petrosian->axes = {0.0, 0.0, 0.0};
     41    petrosian->radiusElliptical = NULL;
     42    petrosian->fluxElliptical = NULL;
    3543
    36     return(petrosian);
     44    petrosian->radialBins = NULL;
     45    petrosian->area = NULL;
     46    petrosian->binnedFlux = NULL;
     47
     48    petrosian->petrosianRadius = NAN;
     49    petrosian->petrosianFlux = NAN;
     50
     51    // petrosian->axes = {0.0, 0.0, 0.0};
     52
     53    return petrosian;
    3754}
    3855
     56bool psphotPetrosianFreeVectors(pmPetrosian *petrosian) {
     57
     58    psFree(petrosian->radii);
     59    psFree(petrosian->fluxes);
     60    psFree(petrosian->theta);
     61    psFree(petrosian->isophotalRadii);
     62
     63    psFree(petrosian->radiusElliptical);
     64    psFree(petrosian->fluxElliptical);
     65
     66    psFree(petrosian->binnedFlux);
     67    psFree(petrosian->radialBins);
     68    psFree(petrosian->area);
     69
     70    petrosian->radii = NULL;
     71    petrosian->fluxes = NULL;
     72    petrosian->theta = NULL;
     73    petrosian->isophotalRadii = NULL;
     74
     75    petrosian->radiusElliptical = NULL;
     76    petrosian->fluxElliptical = NULL;
     77
     78    petrosian->radialBins = NULL;
     79    petrosian->area = NULL;
     80    petrosian->binnedFlux = NULL;
     81   
     82    return true;
     83}
     84
     85# define COMPARE_INDEX(A,B) (index->data.F32[A] < index->data.F32[B])
     86# define SWAP_INDEX(TYPE,A,B) { \
     87  float tmp; \
     88  if (A != B) { \
     89    tmp = index->data.F32[A]; \
     90    index->data.F32[A] = index->data.F32[B]; \
     91    index->data.F32[B] = tmp; \
     92    tmp = extra->data.F32[A]; \
     93    extra->data.F32[A] = extra->data.F32[B]; \
     94    extra->data.F32[B] = tmp; \
     95  } \
     96}
     97
     98bool psphotPetrosianSortPair (psVector *index, psVector *extra) {
     99
     100    // sort the vector set by the radius
     101    PSSORT (index->n, COMPARE_INDEX, SWAP_INDEX, NONE);
     102    return true;
     103}

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

    r25032 r25105  
    1212
    1313typedef struct {
    14     psArray *radii;                     // radii for raw radial profiles at evenly-spaced angles
    15     psArray *fluxes;                    // fluxes measured at above radii
     14    psArray  *radii;                    // radii for raw radial profiles at evenly-spaced angles
     15    psArray  *fluxes;                   // fluxes measured at above radii
    1616    psVector *theta;                    // angles corresponding to above radial profiles
    17     psVector *isophotalRadii;
    18     psEllipseAxes axes;
     17    psVector *isophotalRadii;           // isophotal radius for the above angles
     18
     19    psVector *radiusElliptical;         // normalized radial coordinates for all relevant pixels
     20    psVector *fluxElliptical;           // flux for the above radial coordinates
     21
     22    psVector *binnedFlux;               // mean surface brightness within radial bins
     23    psVector *radialBins;               // radii corresponding to above binnedBlux
     24    psVector *area;                     // differential area of the non-overlapping radial bins
     25
     26    psEllipseAxes axes;                 // shape of elliptical contour
     27
     28    float petrosianRadius;
     29    float petrosianFlux;
     30
    1931} pmPetrosian;
    2032
    2133pmPetrosian *pmPetrosianAlloc();
     34bool psphotPetrosianFreeVectors(pmPetrosian *petrosian);
     35
     36bool psphotPetrosianProfile (pmSource *source);
     37bool psphotRadialProfilesByAngles (pmPetrosian *petro, pmSource *source, int Nsec, float Rmax);
     38float psphotRadiusFromProfile (psVector *radius, psVector *flux, float fluxMin, float fluxMax);
     39bool psphotRadiiFromProfiles (pmPetrosian *petrosian, float fluxMin, float fluxMax);
     40bool psphotEllipticalProfile (pmSource *source, pmPetrosian *petrosian);
     41bool psphotEllipticalContour (pmPetrosian *petrosian);
     42bool psphotPetrosianRadialBins (pmSource *source, pmPetrosian *petrosian, float radiusMax);
     43bool psphotPetrosianStats (pmPetrosian *petrosian);
     44
     45bool psphotPetrosianSortPair (psVector *index, psVector *extra);
     46
     47
     48bool psphotPetrosianVisualProfileByAngle (psVector *radius, psVector *flux);
     49bool psphotPetrosianVisualProfileRadii (psVector *radius, psVector *flux, psVector *radiusBin, psVector *fluxBin, float RadiusRef);
     50bool psphotPetrosianVisualEllipticalContour (pmPetrosian *petrosian);
     51
     52bool psphotPetrosianVisualStats (psVector *radBin, psVector *fluxBin,
     53                                 psVector *refRadius, psVector *meanSB,
     54                                 psVector *petRatio, psVector *fluxSum,
     55                                 float petRadius, float petFlux);
     56
     57bool pmVisualLimitsFromVectors (Graphdata *graphdata, psVector *xVec, psVector *yVec);
    2258
    2359/// @}
     60
     61
    2462# endif /* PM_PETROSIAN_H */

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

    r25032 r25105  
    55// model parameters
    66enum {PAR_PHI, PAR_EPSILON, PAR_RMIN};
     7psF32 psphotEllipticalContourFunc (psVector *deriv, const psVector *params, const psVector *coord);
    78
    89bool psphotEllipticalContour (pmPetrosian *petrosian) {
    … …  
    1011    // use LMM to fit theta vs radius to an ellipse
    1112    psVector *theta = petrosian->theta;
    12     psVector *radius = petrosian->isophotalRadius;
     13    psVector *radius = petrosian->isophotalRadii;
    1314
    1415    // find Rmin and Rmax for the initial guess
    … …  
    2021    psArray *x = psArrayAlloc(2*radius->n);
    2122    psVector *y = psVectorAlloc(2*radius->n, PS_TYPE_F32);
     23    psVector *yErr = psVectorAlloc(2*radius->n, PS_TYPE_F32);
    2224
     25    int n = 0;
    2326    for (int i = 0; i < radius->n; i++) {
    2427
    … …  
    2831        coord = psVectorAlloc (2, PS_TYPE_F32);
    2932        coord->data.F32[1] = 0.0;
    30         coord->data.F32[0] = radius->data.F32[i]*cos(theta->data.F32[i]);
     33        coord->data.F32[0] = theta->data.F32[i];
    3134        x->data[n] = coord;
    32         y->data.F32[n] = theta->data.F32[i];
     35        y->data.F32[n] = radius->data.F32[i]*cos(theta->data.F32[i]);
     36        yErr->data.F32[n] = 1000.0;
    3337        n++;
    3438
    … …  
    3640        coord = psVectorAlloc (2, PS_TYPE_F32);
    3741        coord->data.F32[1] = 1.0;
    38         coord->data.F32[0] = radius->data.F32[i]*sin(theta->data.F32[i]);
     42        coord->data.F32[0] = theta->data.F32[i];
    3943        x->data[n] = coord;
    40         y->data.F32[n] = theta->data.F32[i];
     44        y->data.F32[n] = radius->data.F32[i]*sin(theta->data.F32[i]);
     45        yErr->data.F32[n] = 1000.0;
    4146        n++;
    4247
    … …  
    4954
    5055    psVector *params = psVectorAlloc (3, PS_TYPE_F32);
     56   
     57    // psTraceSetLevel ("psLib.math.psMinimizeLMChi2", 7);
    5158   
    5259    // create the minimization constraints
    … …  
    6471   
    6572    // XXX skip the weights for now
    66     psMinimizeLMChi2(myMin, covar, params, constraint, x, y, NULL, psphotEllipticalContourFunc);
     73    psMinimizeLMChi2(myMin, covar, params, constraint, x, y, yErr, psphotEllipticalContourFunc);
    6774
    6875    fprintf (stderr, "# fitted values:\n");
    69     fprintf (stderr, "Po:  %f\n", params->data.F32[PAR_PHI]);
     76    fprintf (stderr, "Po:  %f\n", params->data.F32[PAR_PHI]*PS_DEG_RAD);
    7077    fprintf (stderr, "Ep:  %f\n", params->data.F32[PAR_EPSILON]);
    7178    fprintf (stderr, "Rm:  %f\n", params->data.F32[PAR_RMIN]);
    … …  
    7582    petrosian->axes.minor = params->data.F32[PAR_RMIN];
    7683    petrosian->axes.theta = params->data.F32[PAR_PHI];
     84
     85    // show the results
     86    psphotPetrosianVisualEllipticalContour (petrosian);
     87
     88    psFree (x);
     89    psFree (y);
     90    psFree (yErr);
    7791
    7892    return true;
    … …  
    88102psF32 psphotEllipticalContourFunc (psVector *deriv, const psVector *params, const psVector *coord) {
    89103
     104    static int pass = 0;
     105
    90106    psF32 *par = params->data.F32;
     107
     108    float alpha = coord->data.F32[0];
    91109
    92110    float cs_alpha = cos(alpha);
    93111    float sn_alpha = sin(alpha);
    94112
    95     float alpha = coord->data.F32[0];
    96113    float cs_phi = cos(alpha - par[PAR_PHI]);
    97114    float sn_phi = sin(alpha - par[PAR_PHI]);
    … …  
    103120
    104121    // value is X
    105     if (coord->data.F32[1] == 0) {
     122    // if (coord->data.F32[1] == 0) {
     123    if (pass == 0) {
     124        pass = 1;
    106125
    107126        float value = par[PAR_RMIN]*cs_alpha*r;
    … …  
    109128        if (deriv) {
    110129            psF32 *dPAR = deriv->data.F32;
    111             dpar[PAR_RMIN]    = r*cs_alpha;
    112             dpar[PAR_EPSILON] = par[PAR_RMIN]*cs_alpha*drdE;
    113             dpar[PAR_PHI]     = 4.0*par[PAR_RMIN]*cs_alpha*drdP;
     130            dPAR[PAR_RMIN]    = r*cs_alpha;
     131            dPAR[PAR_EPSILON] = par[PAR_RMIN]*cs_alpha*drdE;
     132            dPAR[PAR_PHI]     = 4.0*par[PAR_RMIN]*cs_alpha*drdP;
    114133        }
    115134        return (value);
    … …  
    117136
    118137    // value is Y
    119     if (coord->data.F32[1] == 1) {
     138    // if (coord->data.F32[1] == 1) {
     139    if (pass == 1) {
     140        pass = 0;
    120141
    121142        float value = par[PAR_RMIN]*sn_alpha*r;
    … …  
    123144        if (deriv) {
    124145            psF32 *dPAR = deriv->data.F32;
    125             dpar[PAR_RMIN]    = r*sn_alpha;
    126             dpar[PAR_EPSILON] = par[PAR_RMIN]*sn_alpha*drdE;
    127             dpar[PAR_PHI]     = 4.0*par[PAR_RMIN]*sn_alpha*drdP;
     146            dPAR[PAR_RMIN]    = r*sn_alpha;
     147            dPAR[PAR_EPSILON] = par[PAR_RMIN]*sn_alpha*drdE;
     148            dPAR[PAR_PHI]     = 4.0*par[PAR_RMIN]*sn_alpha*drdP;
    128149        }
    129150        return (value);

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

    r25032 r25105  
    33bool psphotEllipticalProfile (pmSource *source, pmPetrosian *petrosian) {
    44
    5     petrosian->radiusElli = psArrayAllocEmpty(100);
    6     petrosian->fluxElli = psArrayAllocEmpty(100);
     5    petrosian->radiusElliptical = psVectorAllocEmpty(100, PS_TYPE_F32);
     6    petrosian->fluxElliptical = psVectorAllocEmpty(100, PS_TYPE_F32);
    77
    8     psVector *radius = petrosian->radiusElli;
    9     psVector *flux = petrosian->fluxElli;
     8    psVector *radius = petrosian->radiusElliptical;
     9    psVector *flux = petrosian->fluxElliptical;
    1010
    1111    // the psEllipse functions use z = 0.5(x/Sxx)^2 + 0.5(y/Syy)^2 + x y Sxy
    … …  
    1414    // a = A / sqrt(2)
    1515
     16    // we have the shape parameters of the elliptical contour at the reference isophote.
     17    // use the axis ratio (major/minor) to rescale the radial profile so that 1 pixel
     18    // along the major axis is 1 pixel, and a smaller amount on the minor axis
     19
    1620    psEllipseAxes axes;
    17     axes.major = petrosian->axes.major * M_SQRT1_2;
    18     axes.minor = petrosian->axes.minor * M_SQRT1_2;
     21    axes.major = M_SQRT1_2;
     22    axes.minor = M_SQRT1_2 * (petrosian->axes.minor / petrosian->axes.major);
     23
     24    // axes.major = 1.0;
     25    // axes.minor = petrosian->axes.minor / petrosian->axes.major;
     26
    1927    axes.theta = petrosian->axes.theta;
    20     psEllipseShape shape = psEllipseShapeFromAxes (petrosian->axes);
     28    psEllipseShape shape = psEllipseAxesToShape (axes);
    2129
    2230    float Sxx = shape.sx;
    … …  
    2432    float Syy = shape.sy;
    2533
    26     for (int iy = 0; iy < Ny; iy++) {
    27         for (int ix = 0; ix < Nx; ix++) {
     34    psVector *radiusRaw = psVectorAllocEmpty(100, PS_TYPE_F32);
     35    psVector *fluxRaw = psVectorAllocEmpty(100, PS_TYPE_F32);
     36
     37    for (int iy = 0; iy < source->pixels->numRows; iy++) {
     38        for (int ix = 0; ix < source->pixels->numCols; ix++) {
    2839
    2940            float x = ix - source->peak->xf + source->pixels->col0;
    3041            float y = iy - source->peak->yf + source->pixels->row0;
    3142
    32             psVectorAppend(radius, sqrt(0.5*PS_SQR(x/Sxx) + 0.5*PS_SQR(y/Syy) + x*y*Sxy));
    33             psVectorAppend(flux, pixels[iy][ix]);
     43            float r2 = 0.5*PS_SQR(x/Sxx) + 0.5*PS_SQR(y/Syy) + x*y*Sxy;
     44            // float r2 = PS_SQR(x/Sxx) + PS_SQR(y/Syy) + x*y*Sxy;
     45
     46            psVectorAppend(radius, sqrt(r2));
     47            psVectorAppend(flux, source->pixels->data.F32[iy][ix]);
     48
     49            float Rraw = hypot(x, y);
     50            psVectorAppend(radiusRaw, Rraw);
     51            psVectorAppend(fluxRaw, source->pixels->data.F32[iy][ix]);
    3452        }
    3553    }
    3654
     55    // psVector *radiusRaw = psVectorAllocEmpty(100, PS_TYPE_F32);
     56    // psVector *fluxRaw = psVectorAllocEmpty(100, PS_TYPE_F32);
     57    // for (int i = 0; i < petrosian->radii->n; i++) {
     58    //   psVector *r = petrosian->radii->data[i];
     59    //   psVector *f = petrosian->fluxes->data[i];
     60    //   for (int j = 0; j < r->n; j++) {
     61    //  psVectorAppend(radiusRaw, r->data.F32[j]);
     62    //  psVectorAppend(fluxRaw, f->data.F32[j]);
     63    //   }
     64    // }
     65
     66    // psphotPetrosianVisualProfileRadii (radius, flux, radiusRaw, fluxRaw, 0.0);
     67    psphotPetrosianVisualProfileByAngle (radius, flux);
    3768    return true;
    3869}

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

    r12805 r25105  
    1414#include <psmodules.h>
    1515#include "psphot.h"
     16#include "pmPetrosian.h"
    1617
    1718#endif

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

    r25032 r25105  
    11# include "psphotInternal.h"
     2
     3// generate the Petrosian radius and flux using elliptical contours
     4
     5// XXX much of this function is focused on generating the clean contours, which can be used by
     6// any number of aperture-like measurements.  probably will want to rename the pmPetrosian
     7// structure to something the pmRadialProfile
    28
    39bool psphotPetrosianProfile (pmSource *source) {
    410
     11  // container to hold results from the radial profile analysis
    512  pmPetrosian *petrosian = pmPetrosianAlloc();
    613
    7   if (!psphotRadialProfilesByAngle (petrosian, source, Nsec, Rmax)) {
     14  int Nsec = 24;
     15  float Rmax = 64;
     16  float fluxMin = 0.0;
     17  float fluxMax = 1000.0;
     18
     19  // generate a series of radial profiles at Nsec evenly spaced angles.  the profile flux
     20  // is measured by interpolation for small radii; for large radii, the pixels in a box
     21  // are averaged to increase the S/N (XXX not yet done)
     22  if (!psphotRadialProfilesByAngles (petrosian, source, Nsec, Rmax)) {
    823    psError (PS_ERR_UNKNOWN, false, "failed to measure radial profile for petrosian");
    924    return false;
    1025  }
    1126
     27  // use the radial profiles to determine the radius of a given isophote.  this isophote
     28  // is used to determine the elliptical shape of the object, so it has a relatively high
     29  // value (nominally 50% of the peak)
    1230  if (!psphotRadiiFromProfiles (petrosian, fluxMin, fluxMax)) {
    1331    psError (PS_ERR_UNKNOWN, false, "failed to measure isophotal radii from profiles");
    … …  
    1533  }
    1634
     35  // convert the isophotal radius vs angle measurements to an elliptical contour
    1736  if (!psphotEllipticalContour (petrosian)) {
    1837    psError (PS_ERR_UNKNOWN, false, "failed to measure elliptical contour");
    … …  
    2039  }
    2140 
     41  // generate a single, normalized radial profile following the elliptical contours.
     42  // the radius is normalized by the axis ratio so that on the major axis, 1 pixel = 1 pixel
     43  if (!psphotEllipticalProfile (source, petrosian)) {
     44    psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile");
     45    return false;
     46  }
     47 
     48  // integrate the radial profile for radial bins defined for the petrosian measurement:
     49  // SB_i (r_i) where \alpha r_i < r < \beta r_i
     50  if (!psphotPetrosianRadialBins (source, petrosian, Rmax)) {
     51    psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile");
     52    return false;
     53  }
     54 
     55  // use the SB_i from above to calculate the petrosian radius and the flux within that radius
     56  if (!psphotPetrosianStats (petrosian)) {
     57    psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile");
     58    return false;
     59  }
     60 
     61  psphotPetrosianFreeVectors(petrosian);
    2262
     63  fprintf (stderr, "petrosian radius: %f, flux: %f, axis ratio: %f, angle: %f\n",
     64           petrosian->petrosianRadius, petrosian->petrosianFlux, petrosian->axes.minor/petrosian->axes.major, PS_DEG_RAD*petrosian->axes.theta);
     65
     66  return true;
    2367}

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

    r25032 r25105  
    22
    33// convert the flux vs elliptical radius to annular bins
    4 bool psphotPetrosianRadialBins (pmSource *source, pmPetrosian *petrosian) {
    54
    6     psVector *radius = petrosian->radiusElli;
    7     psVector *flux = petrosian->fluxElli;
     5// we are guaranteed to be limited by either the seeing (1 - few pixels) or by the pixels
     6// themselves.  this function does not attempt to measure the radial profiles accurately
     7// for radii that are smaller than a minimum (currently 1.0 pixels). 
    88
    9     sort (radius, flux);
     9// for small radii, we are measuring the mean surface brightness in non-overlapping radial
     10// bins.  for large radii (r > 2 pixels), we are measuring the surface brightness for a
     11// radius range \alpha r_i < i < \beta r_i, but performing this measurement for radii more
     12// finely spaced than r_{i+1} = r_i * \beta / \alpha.  for the integration, we need to
     13// track the non-overlapping radius values.
     14
     15# define PETROSIAN_ALPHA 0.8
     16# define PETROSIAN_BETA 1.25
     17
     18bool psphotPetrosianRadialBins (pmSource *source, pmPetrosian *petrosian, float radiusMax) {
     19
     20    psVector *radius = petrosian->radiusElliptical;
     21    psVector *flux = petrosian->fluxElliptical;
     22
     23    // sort incoming vectors by radius
     24    psphotPetrosianSortPair (radius, flux);
     25
     26    int nMax = radiusMax;
    1027
    1128    // radBin stores the centers of the radial bins,
    1229    // radMin, radMax store the bounds
    13     psVector *radBin = psVectorAllocEmpty(100, PS_TYPE_F32);
    14     psVector *radMin = psVectorAllocEmpty(100, PS_TYPE_F32);
    15     psVector *radMax = psVectorAllocEmpty(100, PS_TYPE_F32);
     30    psVector *radMin  = psVectorAllocEmpty(nMax, PS_TYPE_F32);
     31    psVector *radMax  = psVectorAllocEmpty(nMax, PS_TYPE_F32);
     32    psVector *radAlp  = psVectorAllocEmpty(nMax, PS_TYPE_F32);
     33    psVector *radBet  = psVectorAllocEmpty(nMax, PS_TYPE_F32);
     34
     35    psVector *fluxBin = psVectorAllocEmpty(nMax, PS_TYPE_F32);
     36    psVector *radBin  = psVectorAllocEmpty(nMax, PS_TYPE_F32);
     37    psVector *area    = psVectorAllocEmpty(nMax, PS_TYPE_F32);
     38
     39    psVectorInit (fluxBin, 0.0);
     40    psVectorInit (radBin, 0.0);
    1641
    1742    // generate radial bin bounds
    1843    radMin->data.F32[0] = 0.0;
    1944    radMax->data.F32[0] = 1.0;
     45    radAlp->data.F32[0] = 0.0;
     46    radBet->data.F32[0] = 1.0;
    2047   
    2148    radMin->data.F32[1] = 1.0;
    2249    radMax->data.F32[1] = 1.5;
     50    radAlp->data.F32[1] = 1.0;
     51    radBet->data.F32[1] = 1.5;
    2352   
    2453    radMin->data.F32[2] = 1.5;
    2554    radMax->data.F32[2] = 2.0;
     55    radAlp->data.F32[2] = 1.5;
     56    radBet->data.F32[2] = 2.0;
    2657   
    27     for (int i = 3; i < rmax; i++) {
    28         radMin->data.F32[i] = (i - 1);
    29         radMax->data.F32[i] = i;
     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++;
    3063    }
     64    radMin->n = radMax->n = radAlp->n = radBet->n = nPts;
    3165
    3266    // generate radial area-weighted mean radius & non-overlapping areas
    33     float Rprev = 0.0;
    3467    for (int i = 0; i < radMin->n; i++) {
    3568        float rMin = radMin->data.F32[i];
    … …  
    4275        float rMax3 = rMax2*rMax;
    4376
    44         float rBin = (rMax3 - rMin3) / (rMax2 - rMin2) / 3.0;
     77        float rBin = 2.0 * (rMax3 - rMin3) / (rMax2 - rMin2) / 3.0;
    4578       
     79        // XXX calculate area-weighted radius rather than asserting?
    4680        radBin->data.F32[i] = rBin;
     81        area->data.F32[i] = M_PI * (rMax2 - rMin2);
    4782
    48         area->data.F32[i] =
     83        if (i > 2) {
     84            radAlp->data.F32[i] = rBin*PETROSIAN_ALPHA;
     85            radBet->data.F32[i] = rBin*PETROSIAN_BETA;
     86        }
    4987    }
     88
     89    // storage vector for stats
     90    psVector *values = psVectorAllocEmpty (flux->n, PS_TYPE_F32);
     91    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
     92
     93    // integrate flux, radius for each of these bins.  since flux is sorted by radius,
     94    // we can do this fairly quickly
     95
     96    bool done = false;
     97    int nOut = 0;
     98    float Rmin = radAlp->data.F32[nOut];
     99    float Rmax = radBet->data.F32[nOut];
     100    float Rnxt = radAlp->data.F32[nOut+1];  // minimum radius for next range
     101    int iNext = 0;
     102    for (int i = 0; !done && (i < radius->n); i++) {
     103        if (radius->data.F32[i] < Rnxt) {
     104          iNext = i;
     105        }
     106        if (radius->data.F32[i] > Rmax) {
     107            // calculate the value for the nOut bin
     108            psVectorStats (stats, values, NULL, NULL, 0);
     109            fluxBin->data.F32[nOut] = stats->sampleMedian;
     110            nOut ++;
     111            if (nOut >= nMax) break;
     112            Rmin = radAlp->data.F32[nOut];
     113            Rmax = radBet->data.F32[nOut];
     114            Rnxt = (nOut < nMax - 1) ? radAlp->data.F32[nOut+1] : Rmax;  // minimum radius for next range
     115            values->n = 0;
     116            psStatsInit(stats);
     117            i = iNext;
     118        }
     119        if (radius->data.F32[i] < Rmin) {
     120            continue;
     121        }
     122        psVectorAppend (values, flux->data.F32[i]);
     123    }
     124    fluxBin->n = radBin->n = area->n = nOut;
     125    // XXX I think this misses the last radial bin -- do we care?
     126
     127    // save the vectors
     128    petrosian->radialBins = radBin;
     129    petrosian->area = area;
     130    petrosian->binnedFlux = fluxBin;
     131
     132    psphotPetrosianVisualProfileRadii (radius, flux, radBin, fluxBin, 0.0);
     133
     134    psFree(radMin);
     135    psFree(radMax);
     136    psFree(radAlp);
     137    psFree(radBet);
     138    psFree(values);
     139    psFree(stats);
     140
     141    return true;
    50142}
    51143

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

    r25032 r25105  
    88// XXX how does this fit in context?
    99
    10 bool psphotRadialProfilesByAngles (pmPetrosianData *petro, pmSource *source, int Nsec, float Rmax) {
     10bool psphotRadialProfilesByAngles (pmPetrosian *petrosian, pmSource *source, int Nsec, float Rmax) {
    1111
    1212    float dtheta = 2.0*M_PI / Nsec;
    1313
    14     psFree(petro->radii);
    15     psFree(petro->fluxes);
    16     psFree(petro->theta);
     14    psFree(petrosian->radii);
     15    psFree(petrosian->fluxes);
     16    psFree(petrosian->theta);
    1717
    18     petro->radii = psArrayAllocEmpty(Nsec);
    19     petro->fluxes = psArrayAllocEmpty(Nsec);
    20     petro->theta = psVectorAllocEmpty(Nsec, PS_TYPE_F32);
     18    petrosian->radii = psArrayAllocEmpty(Nsec);
     19    petrosian->fluxes = psArrayAllocEmpty(Nsec);
     20    petrosian->theta = psVectorAllocEmpty(Nsec, PS_TYPE_F32);
    2121
    2222
    … …  
    2525        float theta = i*dtheta;
    2626
    27         psVector *radius = psVectorAllocEmpty(PS_TYPE_F32, Rmax);
    28         psVector *flux   = psVectorAllocEmpty(PS_TYPE_F32, Rmax);
     27        psVector *radius = psVectorAllocEmpty(Rmax, PS_TYPE_F32);
     28        psVector *flux   = psVectorAllocEmpty(Rmax, PS_TYPE_F32);
    2929
    3030        // start at Xo,Yo and find the x,y locations for r_i, theta where r_i increments by 1 pixel
    … …  
    3232        for (float r = 0; r < Rmax; r += 1.0) {
    3333
    34             float Xo = source->peak->xf - source->pixels->col0;
    35             float Yo = source->peak->yf - source->pixels->row0;
     34            float Xo = source->peak->xf + 0.5;
     35            float Yo = source->peak->yf + 0.5;
    3636
    3737            // Xo,Yo are referenced to pixels with bounds i+0.0, i+1.0
    38             x = r * cos (theta) + Xo
    39             y = r * sin (theta) + Yo;
     38            float x = r * cos (theta) + Xo;
     39            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;
    4045
    4146            // value is NAN if we run off the image
    42             value = psImageInterpolatePixelBilinear(x, y, image);
     47            float value = psImageInterpolatePixelBilinear(x, y, source->pixels);
    4348            if (isnan(value)) continue;
    4449
    … …  
    4752        }
    4853
    49         psArrayAdd (petro->radii, 100, radius);
    50         psArrayAdd (petro->fluxes, 100, flux);
    51         psVectorAppend (petro->theta, theta);
     54        psArrayAdd (petrosian->radii, 100, radius);
     55        psArrayAdd (petrosian->fluxes, 100, flux);
     56        psVectorAppend (petrosian->theta, theta);
     57
     58        // psphotPetrosianVisualProfileByAngle (radius, flux);
     59
     60        psFree(radius);
     61        psFree(flux);
    5262    }
    5363

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

    r25032 r25105  
    1010      psVector *radii = petrosian->radii->data[i];
    1111      psVector *fluxes = petrosian->fluxes->data[i];
    12       float radius =  = psphotRadiusFromProfile (radii, fluxes, fluxMin, fluxMax);
     12      float radius = psphotRadiusFromProfile (radii, fluxes, fluxMin, fluxMax);
     13
     14      // psphotPetrosianVisualProfileByAngle (radii, fluxes, radius);
    1315
    1416      // warn on NAN?
    … …  
    4446
    4547    // rebin with bins 1/2 the size of the mean radius
     48    int Rbin = 0;
    4649    if (Rnpt == 0) {
    4750        Rbin = 1;
    … …  
    5861        radiusBinned = psMemIncrRefCounter (radius);
    5962    } else {
    60         // storage vector for sorting
    61         psVector *values = psVectorAllocEmpty (PS_TYPE_F32, flux->n);
     63        // storage vector for stats
     64        psVector *values = psVectorAllocEmpty (flux->n, PS_TYPE_F32);
     65        psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
    6266 
    6367        // rebinned vectors
    64         fluxBinned = psVectorAllocEmpty (PS_TYPE_F32, flux->n);
    65         radiusBinned = psVectorAllocEmpty (PS_TYPE_F32, flux->n);
     68        fluxBinned = psVectorAllocEmpty (flux->n, PS_TYPE_F32);
     69        radiusBinned = psVectorAllocEmpty (flux->n, PS_TYPE_F32);
    6670
    6771        // sort the flux by the radius
    68         sort2vec (radius, flux);
     72        psphotPetrosianSortPair (radius, flux);
    6973
    7074        int nOut = 0;
    7175        radiusBinned->data.F32[nOut] = (nOut + 0.5)*Rbin;
    72         Rmin = radiusBinned->data.F32[nOut] - 0.5*Rbin;
    73         Rmax = radiusBinned->data.F32[nOut] + 0.5*Rbin;
     76        float Rmin = radiusBinned->data.F32[nOut] - 0.5*Rbin;
     77        float Rmax = radiusBinned->data.F32[nOut] + 0.5*Rbin;
    7478
    7579        for (int i = 0; i < flux->n; i++) {
    … …  
    8084            if (radius->data.F32[i] > Rmax) {
    8185                // calculate the value for the nOut bin
    82                 psStats (stats, values);
    83                 fluxBinned->data.F32[nOut] = stats->value;
     86                // XXX need to fix this as well psStats (stats, values);
     87                psVectorStats (stats, values, NULL, NULL, 0);
     88                fluxBinned->data.F32[nOut] = stats->sampleMedian;
    8489                nOut ++;
    8590                radiusBinned->data.F32[nOut] = (nOut + 0.5)*Rbin;
    8691                Rmin = radiusBinned->data.F32[nOut] - 0.5*Rbin;
    8792                Rmax = radiusBinned->data.F32[nOut] + 0.5*Rbin;
    88                 i--;
    89                 continue;
    90                 // XXX the sequencing here is a bit ugly -- can we unfold this?
     93                values->n = 0;
     94                psStatsInit(stats);
    9195            }
    9296            psVectorAppend (values, flux->data.F32[i]);
    … …  
    9599        radiusBinned->n = nOut;
    96100        psFree (values);
     101        psFree(stats);
    97102    }
    98103
    99     Ro = NAN;
    100     above = TRUE;
     104    float Ro = NAN;
     105    bool above = true;
    101106    for (int i = 0; i < fluxBinned->n; i++) {
    102107
    … …  
    105110            // XXX is there a macro in psLib that does this interpolation?
    106111            if (i == 0) {
    107                 // assume Fmax @ R = 0.0
    108                 Ro = radiusBinned->data.F32[i] * (Fo - Fmax) / (fluxBinned->data.F32[i] - Fmax);
     112              psAbort ("impossible condition f[0] < Fo");
    109113            } else {
    110114                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]);
    … …  
    118122    }
    119123
    120     // XXX call a visual function to show the unbinned & binned vectors + result radius
     124    // show the results
     125    // psphotPetrosianVisualProfileRadii (radius, flux, radiusBinned, fluxBinned, Ro);
    121126
    122127    psFree(fluxBinned);
Note: See TracChangeset for help on using the changeset viewer.