Index: /branches/eam_branches/20090715/psphot/doc/notes.20090523.txt =================================================================== --- /branches/eam_branches/20090715/psphot/doc/notes.20090523.txt (revision 25031) +++ /branches/eam_branches/20090715/psphot/doc/notes.20090523.txt (revision 25032) @@ -1,2 +1,30 @@ + +20090809 : more on extended sources: + + my algorithm for getting the petrosian radii and fluxes is: + + * measure radial profile by interpolation to specific locations along the radial line + (at low radii, this gives much more accurate results; at high radii, we may need to + average a group of pixels -- say, for 15 deg separations, choose a box that is + < 5 degree in width -- this is > 1 pixel at a distance of 12 pixels. + + psphotRadialProfilesByAngle + + * find intersection of profile with isophot + + psphotRadiusFromProfile -> r50(theta) + + * generate r50x, r50y and fit to ellipse + + psphotEllipticalContour + + * generate elliptical profile + + psphotEllipticalProfile + + * convert to \alpha r_i < r < \beta r_i bins + + psphotPetrosianRadialBins + 20090725 : extended source radial profiles Index: /branches/eam_branches/20090715/psphot/src/pmPetrosian.c =================================================================== --- /branches/eam_branches/20090715/psphot/src/pmPetrosian.c (revision 25032) +++ /branches/eam_branches/20090715/psphot/src/pmPetrosian.c (revision 25032) @@ -0,0 +1,38 @@ +/* @file pmPetrosian.c + * low-level petrosian functions + * + * @author EAM, IfA + * + * @version $Revision: $ + * @date $Date: $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +# include "psphotInternal.h" + +static void pmPetrosianFree(pmPetrosian *petro) +{ + if (!petro) { + return; + } + psFree(petrosian->radii); + psFree(petrosian->fluxes); + psFree(petrosian->theta); + psFree(petrosian->isophotalRadii); +} + +pmPetrosian *pmPetrosianAlloc() +{ + pmPetrosian *petrosian = (pmPetrosian *)psAlloc(sizeof(pmPetrosian)); + psMemSetDeallocator(petrosian, (psFreeFunc) pmPetrosianFree); + + petrosian->radii = NULL; + petrosian->fluxes = NULL; + petrosian->theta = NULL; + petrosian->isophotalRadii = NULL; + + petrosian->axes = {0.0, 0.0, 0.0}; + + return(petrosian); +} + Index: /branches/eam_branches/20090715/psphot/src/pmPetrosian.h =================================================================== --- /branches/eam_branches/20090715/psphot/src/pmPetrosian.h (revision 25032) +++ /branches/eam_branches/20090715/psphot/src/pmPetrosian.h (revision 25032) @@ -0,0 +1,24 @@ +/* @file pmPetrosian.h + * + * @author EAM, IfA + * + * @version $Revision: $ + * @date $Date: $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#ifndef PM_PETROSIAN_H +#define PM_PETROSIAN_H + +typedef struct { + psArray *radii; // radii for raw radial profiles at evenly-spaced angles + psArray *fluxes; // fluxes measured at above radii + psVector *theta; // angles corresponding to above radial profiles + psVector *isophotalRadii; + psEllipseAxes axes; +} pmPetrosian; + +pmPetrosian *pmPetrosianAlloc(); + +/// @} +# endif /* PM_PETROSIAN_H */ Index: /branches/eam_branches/20090715/psphot/src/psphotEllipticalContour.c =================================================================== --- /branches/eam_branches/20090715/psphot/src/psphotEllipticalContour.c (revision 25032) +++ /branches/eam_branches/20090715/psphot/src/psphotEllipticalContour.c (revision 25032) @@ -0,0 +1,133 @@ +# include "psphotInternal.h" + +// XXX consistency : theta in radians here and in calling functions + +// model parameters +enum {PAR_PHI, PAR_EPSILON, PAR_RMIN}; + +bool psphotEllipticalContour (pmPetrosian *petrosian) { + + // use LMM to fit theta vs radius to an ellipse + psVector *theta = petrosian->theta; + psVector *radius = petrosian->isophotalRadius; + + // find Rmin and Rmax for the initial guess + float Rmin = radius->data.F32[0]; + float Rmax = radius->data.F32[0]; + + // arrays to hold the data to be fitted + // we fit x and y vs theta in separate passes. + psArray *x = psArrayAlloc(2*radius->n); + psVector *y = psVectorAlloc(2*radius->n, PS_TYPE_F32); + + for (int i = 0; i < radius->n; i++) { + + psVector *coord = NULL; + + // Rx coordinate value + coord = psVectorAlloc (2, PS_TYPE_F32); + coord->data.F32[1] = 0.0; + coord->data.F32[0] = radius->data.F32[i]*cos(theta->data.F32[i]); + x->data[n] = coord; + y->data.F32[n] = theta->data.F32[i]; + n++; + + // Ry coordinate value + coord = psVectorAlloc (2, PS_TYPE_F32); + coord->data.F32[1] = 1.0; + coord->data.F32[0] = radius->data.F32[i]*sin(theta->data.F32[i]); + x->data[n] = coord; + y->data.F32[n] = theta->data.F32[i]; + n++; + + // check the radius range + Rmin = MIN (Rmin, radius->data.F32[i]); + Rmax = MAX (Rmax, radius->data.F32[i]); + } + assert (x->n == n); + assert (y->n == n); + + psVector *params = psVectorAlloc (3, PS_TYPE_F32); + + // create the minimization constraints + psMinConstraint *constraint = psMinConstraintAlloc(); + + // XXX for now, no parameter masks, skip checkLimits + // constraint->checkLimits = psastroModelBoresiteLimits; + + params->data.F32[PAR_PHI] = 0.0; + params->data.F32[PAR_EPSILON] = Rmin / Rmax; + params->data.F32[PAR_RMIN] = Rmin; + + psMinimization *myMin = psMinimizationAlloc (25, 0.001); + psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32); + + // XXX skip the weights for now + psMinimizeLMChi2(myMin, covar, params, constraint, x, y, NULL, psphotEllipticalContourFunc); + + fprintf (stderr, "# fitted values:\n"); + fprintf (stderr, "Po: %f\n", params->data.F32[PAR_PHI]); + fprintf (stderr, "Ep: %f\n", params->data.F32[PAR_EPSILON]); + fprintf (stderr, "Rm: %f\n", params->data.F32[PAR_RMIN]); + + /// XXX rationalize? if epsilon > 1, flip major and minor axes (rotate by 90 degrees) + petrosian->axes.major = params->data.F32[PAR_RMIN] / params->data.F32[PAR_EPSILON]; + petrosian->axes.minor = params->data.F32[PAR_RMIN]; + petrosian->axes.theta = params->data.F32[PAR_PHI]; + + return true; +} + +/** + * the full chisq is built of two associated sums over coordinates: + * chisq = sum ((Rx_obs - Rx_fit(t))^2 + (Ry_obs - Ry_fit(t))^2) + * we use split this into a 2x long vector and use coord[1] to distinguish the X and Y terms: + * coord[0] = measured X or measured Y + * coord[1] = 0 or 1 + */ +psF32 psphotEllipticalContourFunc (psVector *deriv, const psVector *params, const psVector *coord) { + + psF32 *par = params->data.F32; + + float cs_alpha = cos(alpha); + float sn_alpha = sin(alpha); + + float alpha = coord->data.F32[0]; + float cs_phi = cos(alpha - par[PAR_PHI]); + float sn_phi = sin(alpha - par[PAR_PHI]); + + float r = 1.0 / sqrt(SQ(sn_phi) + SQ(par[PAR_EPSILON]*cs_phi)); + float r3 = pow(r, 3.0); + float drdE = -0.5 * r3 * SQ(cs_phi) * 2.0 * par[PAR_EPSILON]; + float drdP = -0.5 * r3 * (SQ(par[PAR_EPSILON]) - 1) * 2.0 * cs_phi * sn_phi; + + // value is X + if (coord->data.F32[1] == 0) { + + float value = par[PAR_RMIN]*cs_alpha*r; + + if (deriv) { + psF32 *dPAR = deriv->data.F32; + dpar[PAR_RMIN] = r*cs_alpha; + dpar[PAR_EPSILON] = par[PAR_RMIN]*cs_alpha*drdE; + dpar[PAR_PHI] = 4.0*par[PAR_RMIN]*cs_alpha*drdP; + } + return (value); + } + + // value is Y + if (coord->data.F32[1] == 1) { + + float value = par[PAR_RMIN]*sn_alpha*r; + + if (deriv) { + psF32 *dPAR = deriv->data.F32; + dpar[PAR_RMIN] = r*sn_alpha; + dpar[PAR_EPSILON] = par[PAR_RMIN]*sn_alpha*drdE; + dpar[PAR_PHI] = 4.0*par[PAR_RMIN]*sn_alpha*drdP; + } + return (value); + } + + psAbort ("programming error: invalid coordinate"); +} Index: /branches/eam_branches/20090715/psphot/src/psphotEllipticalProfile.c =================================================================== --- /branches/eam_branches/20090715/psphot/src/psphotEllipticalProfile.c (revision 25032) +++ /branches/eam_branches/20090715/psphot/src/psphotEllipticalProfile.c (revision 25032) @@ -0,0 +1,38 @@ +# include "psphotInternal.h" + +bool psphotEllipticalProfile (pmSource *source, pmPetrosian *petrosian) { + + petrosian->radiusElli = psArrayAllocEmpty(100); + petrosian->fluxElli = psArrayAllocEmpty(100); + + psVector *radius = petrosian->radiusElli; + psVector *flux = petrosian->fluxElli; + + // the psEllipse functions use z = 0.5(x/Sxx)^2 + 0.5(y/Syy)^2 + x y Sxy + // which are converted to z = 0.5(x/a)^2 + 0.5(y/b)^2 + // we have major and minor axes of a specific ellipse with r^2 = (x/A)^2 + (y/B)^2 + // a = A / sqrt(2) + + psEllipseAxes axes; + axes.major = petrosian->axes.major * M_SQRT1_2; + axes.minor = petrosian->axes.minor * M_SQRT1_2; + axes.theta = petrosian->axes.theta; + psEllipseShape shape = psEllipseShapeFromAxes (petrosian->axes); + + float Sxx = shape.sx; + float Sxy = shape.sxy; + float Syy = shape.sy; + + for (int iy = 0; iy < Ny; iy++) { + for (int ix = 0; ix < Nx; ix++) { + + float x = ix - source->peak->xf + source->pixels->col0; + float y = iy - source->peak->yf + source->pixels->row0; + + psVectorAppend(radius, sqrt(0.5*PS_SQR(x/Sxx) + 0.5*PS_SQR(y/Syy) + x*y*Sxy)); + psVectorAppend(flux, pixels[iy][ix]); + } + } + + return true; +} Index: /branches/eam_branches/20090715/psphot/src/psphotPetrosianProfile.c =================================================================== --- /branches/eam_branches/20090715/psphot/src/psphotPetrosianProfile.c (revision 25032) +++ /branches/eam_branches/20090715/psphot/src/psphotPetrosianProfile.c (revision 25032) @@ -0,0 +1,23 @@ +# include "psphotInternal.h" + +bool psphotPetrosianProfile (pmSource *source) { + + pmPetrosian *petrosian = pmPetrosianAlloc(); + + if (!psphotRadialProfilesByAngle (petrosian, source, Nsec, Rmax)) { + psError (PS_ERR_UNKNOWN, false, "failed to measure radial profile for petrosian"); + return false; + } + + if (!psphotRadiiFromProfiles (petrosian, fluxMin, fluxMax)) { + psError (PS_ERR_UNKNOWN, false, "failed to measure isophotal radii from profiles"); + return false; + } + + if (!psphotEllipticalContour (petrosian)) { + psError (PS_ERR_UNKNOWN, false, "failed to measure elliptical contour"); + return false; + } + + +} Index: /branches/eam_branches/20090715/psphot/src/psphotPetrosianRadialBins.c =================================================================== --- /branches/eam_branches/20090715/psphot/src/psphotPetrosianRadialBins.c (revision 25032) +++ /branches/eam_branches/20090715/psphot/src/psphotPetrosianRadialBins.c (revision 25032) @@ -0,0 +1,58 @@ +# include "psphotInternal.h" + +// convert the flux vs elliptical radius to annular bins +bool psphotPetrosianRadialBins (pmSource *source, pmPetrosian *petrosian) { + + psVector *radius = petrosian->radiusElli; + psVector *flux = petrosian->fluxElli; + + sort (radius, flux); + + // radBin stores the centers of the radial bins, + // radMin, radMax store the bounds + psVector *radBin = psVectorAllocEmpty(100, PS_TYPE_F32); + psVector *radMin = psVectorAllocEmpty(100, PS_TYPE_F32); + psVector *radMax = psVectorAllocEmpty(100, PS_TYPE_F32); + + // generate radial bin bounds + radMin->data.F32[0] = 0.0; + radMax->data.F32[0] = 1.0; + + radMin->data.F32[1] = 1.0; + radMax->data.F32[1] = 1.5; + + radMin->data.F32[2] = 1.5; + radMax->data.F32[2] = 2.0; + + for (int i = 3; i < rmax; i++) { + radMin->data.F32[i] = (i - 1); + radMax->data.F32[i] = i; + } + + // generate radial area-weighted mean radius & non-overlapping areas + float Rprev = 0.0; + for (int i = 0; i < radMin->n; i++) { + float rMin = radMin->data.F32[i]; + float rMax = radMax->data.F32[i]; + + float rMin2 = rMin*rMin; + float rMin3 = rMin2*rMin; + + float rMax2 = rMax*rMax; + float rMax3 = rMax2*rMax; + + float rBin = (rMax3 - rMin3) / (rMax2 - rMin2) / 3.0; + + radBin->data.F32[i] = rBin; + + area->data.F32[i] = + } +} + +// the area-weighted mean radius is given by: + +// integral r * 2 pi r dr / integral 2 pi r dr + +// = 2/3 pi (r_max^3 - r_min^3) / pi (r_max^2 - r_min^2) +// = 2/3 (r_max^3 - r_min^3) / (r_max^2 - r_min^2) + Index: /branches/eam_branches/20090715/psphot/src/psphotRadialProfileByAngle.c =================================================================== --- /branches/eam_branches/20090715/psphot/src/psphotRadialProfileByAngle.c (revision 25032) +++ /branches/eam_branches/20090715/psphot/src/psphotRadialProfileByAngle.c (revision 25032) @@ -0,0 +1,55 @@ +# include "psphotInternal.h" + +// Given a source at (x,y), generate a collection of radial profiles at even angular +// separations + +// XXX choices for rMax and Nsec? +// XXX where does the output go? +// XXX how does this fit in context? + +bool psphotRadialProfilesByAngles (pmPetrosianData *petro, pmSource *source, int Nsec, float Rmax) { + + float dtheta = 2.0*M_PI / Nsec; + + psFree(petro->radii); + psFree(petro->fluxes); + psFree(petro->theta); + + petro->radii = psArrayAllocEmpty(Nsec); + petro->fluxes = psArrayAllocEmpty(Nsec); + petro->theta = psVectorAllocEmpty(Nsec, PS_TYPE_F32); + + + for (int i = 0; i < Nsec; i++) { + + float theta = i*dtheta; + + psVector *radius = psVectorAllocEmpty(PS_TYPE_F32, Rmax); + psVector *flux = psVectorAllocEmpty(PS_TYPE_F32, Rmax); + + // start at Xo,Yo and find the x,y locations for r_i, theta where r_i increments by 1 pixel + // XXX at large radii ( > 10-15) use stats in a patch rather than sub-pixel interpolation + for (float r = 0; r < Rmax; r += 1.0) { + + float Xo = source->peak->xf - source->pixels->col0; + float Yo = source->peak->yf - source->pixels->row0; + + // Xo,Yo are referenced to pixels with bounds i+0.0, i+1.0 + x = r * cos (theta) + Xo + y = r * sin (theta) + Yo; + + // value is NAN if we run off the image + value = psImageInterpolatePixelBilinear(x, y, image); + if (isnan(value)) continue; + + psVectorAppend (radius, r); + psVectorAppend (flux, value); + } + + psArrayAdd (petro->radii, 100, radius); + psArrayAdd (petro->fluxes, 100, flux); + psVectorAppend (petro->theta, theta); + } + + return true; +} Index: /branches/eam_branches/20090715/psphot/src/psphotRadiiFromProfiles.c =================================================================== --- /branches/eam_branches/20090715/psphot/src/psphotRadiiFromProfiles.c (revision 25032) +++ /branches/eam_branches/20090715/psphot/src/psphotRadiiFromProfiles.c (revision 25032) @@ -0,0 +1,125 @@ +# include "psphotInternal.h" + +// Given the Petrosian data (radii, fluxes) determine the radius for each profile at the desisred isophote + +bool psphotRadiiFromProfiles (pmPetrosian *petrosian, float fluxMin, float fluxMax) { + + petrosian->isophotalRadii = psVectorAlloc(petrosian->theta->n, PS_TYPE_F32); + + for (int i = 0; i < petrosian->theta->n; i++) { + psVector *radii = petrosian->radii->data[i]; + psVector *fluxes = petrosian->fluxes->data[i]; + float radius = = psphotRadiusFromProfile (radii, fluxes, fluxMin, fluxMax); + + // warn on NAN? + petrosian->isophotalRadii->data.F32[i] = radius; + } + return true; +} + +float psphotRadiusFromProfile (psVector *radius, psVector *flux, float fluxMin, float fluxMax) { + + // 'flux' is a noisy sample of the galaxy radial profile at points 'radius' + // rebin flux into samples defined by the isophote Fo = 0.5*(fluxMax + fluxMin). the noisy + // sample is cleaned by rebinning to a well-matched radial binning + + // base selections on fluxes defined by the flux range dF + float fluxRange = fluxMax - fluxMin; + + // examine data in the two ranges Fm - Fo and Fo - Fp to define the bin size + float Fm = fluxMin + 0.25*fluxRange; + float Fp = fluxMin + 0.75*fluxRange; + float Fo = fluxMin + 0.50*fluxRange; + + // find the mean radius of the points in the flux range Fm - Fp: + float Rsum = 0; + float Rnpt = 0; + for (int i = 0; i < flux->n; i++) { + if (flux->data.F32[i] < Fm) continue; + if (flux->data.F32[i] > Fp) continue; + + Rsum += radius->data.F32[i]; + Rnpt ++; + } + + // rebin with bins 1/2 the size of the mean radius + if (Rnpt == 0) { + Rbin = 1; + } else { + Rbin = MAX(1, 0.5*(Rsum / Rnpt)); + } + + psVector *fluxBinned = NULL; + psVector *radiusBinned = NULL; + + // do not bother rebinning if the bin size is only 2 or less + if (Rbin <= 2) { + fluxBinned = psMemIncrRefCounter (flux); + radiusBinned = psMemIncrRefCounter (radius); + } else { + // storage vector for sorting + psVector *values = psVectorAllocEmpty (PS_TYPE_F32, flux->n); + + // rebinned vectors + fluxBinned = psVectorAllocEmpty (PS_TYPE_F32, flux->n); + radiusBinned = psVectorAllocEmpty (PS_TYPE_F32, flux->n); + + // sort the flux by the radius + sort2vec (radius, flux); + + int nOut = 0; + radiusBinned->data.F32[nOut] = (nOut + 0.5)*Rbin; + Rmin = radiusBinned->data.F32[nOut] - 0.5*Rbin; + Rmax = radiusBinned->data.F32[nOut] + 0.5*Rbin; + + for (int i = 0; i < flux->n; i++) { + if (radius->data.F32[i] < Rmin) { + // XXX not sure how we can hit this, if there is full coverage of radiusBinned + continue; + } + if (radius->data.F32[i] > Rmax) { + // calculate the value for the nOut bin + psStats (stats, values); + fluxBinned->data.F32[nOut] = stats->value; + nOut ++; + radiusBinned->data.F32[nOut] = (nOut + 0.5)*Rbin; + Rmin = radiusBinned->data.F32[nOut] - 0.5*Rbin; + Rmax = radiusBinned->data.F32[nOut] + 0.5*Rbin; + i--; + continue; + // XXX the sequencing here is a bit ugly -- can we unfold this? + } + psVectorAppend (values, flux->data.F32[i]); + } + fluxBinned->n = nOut; + radiusBinned->n = nOut; + psFree (values); + } + + Ro = NAN; + above = TRUE; + for (int i = 0; i < fluxBinned->n; i++) { + + // find the largest radius that matches the flux transition + if (above && (fluxBinned->data.F32[i] < Fo)) { + // XXX is there a macro in psLib that does this interpolation? + if (i == 0) { + // assume Fmax @ R = 0.0 + Ro = radiusBinned->data.F32[i] * (Fo - Fmax) / (fluxBinned->data.F32[i] - Fmax); + } else { + 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]); + } + above = FALSE; + } + + if (!above && (fluxBinned->data.F32[i] >= Fo)) { + above = TRUE; + } + } + + // XXX call a visual function to show the unbinned & binned vectors + result radius + + psFree(fluxBinned); + psFree(radiusBinned); + return Ro; +}