Index: branches/eam_branches/20090715/ippconfig/mosaic2/camera.config =================================================================== --- branches/eam_branches/20090715/ippconfig/mosaic2/camera.config (revision 25451) +++ branches/eam_branches/20090715/ippconfig/mosaic2/camera.config (revision 25452) @@ -90,4 +90,5 @@ CMF.XSRC STR {CHIP.NAME}.xsrc # use .PSF and .EXT? CMF.XFIT STR {CHIP.NAME}.xfit # use .PSF and .EXT? + CMF.DETEFF STR {CHIP.NAME}.deteff PSF.HEAD STR {CHIP.NAME}.hdr Index: branches/eam_branches/20090715/ippconfig/mosaic2/psphot.config =================================================================== --- branches/eam_branches/20090715/ippconfig/mosaic2/psphot.config (revision 25451) +++ branches/eam_branches/20090715/ippconfig/mosaic2/psphot.config (revision 25452) @@ -24,9 +24,13 @@ USE_FOOTPRINTS BOOL TRUE # use new pmFootprint peak packaging +PSF.TREND.MODE STR MAP +PSF.TREND.NX S32 1 +PSF.TREND.NY S32 1 + PEAKS_NSIGMA_LIMIT F32 15.0 # peak significance threshold MOMENTS_SN_MIN F32 10.0 # min S/N to measure moments PSF_SN_LIM F32 20.0 # minimum S/N for stars used for PSF model -FULL_FIT_SN_LIM F32 20.0 -EXT_MIN_SN F32 20.0 # fit galaxies above this S/N limit +FULL_FIT_SN_LIM F32 10.0 +EXT_MIN_SN F32 10.0 # fit galaxies above this S/N limit AP_MIN_SN F32 10.0 @@ -42,2 +46,6 @@ PSPHOT.EXT.NSIGMA.LIMIT F32 3.0 # sources with extNsigma greater that this get tagged as likely extended sources PSPHOT.EXT.NSIGMA.MOMENTS F32 2.0 # sources with extNsigma greater that this get tagged as likely extended sources + +EXTENDED_SOURCE_SN_LIM F32 3.0 +EXTENDED_SOURCE_ANALYSIS BOOL TRUE # perform any of the aperture-like measurements? +EXTENDED_SOURCE_PETROSIAN BOOL TRUE Index: branches/eam_branches/20090715/psModules/src/objects/pmPSFtry.c =================================================================== --- branches/eam_branches/20090715/psModules/src/objects/pmPSFtry.c (revision 25451) +++ branches/eam_branches/20090715/psModules/src/objects/pmPSFtry.c (revision 25452) @@ -443,4 +443,6 @@ float dSysBright = psVectorSystematicError (bright, brightErr, 0.1); fprintf (stderr, "bright systematic error: %f\n", dSysBright); + psFree(bright); + psFree(brightErr); // XXX test dump of fitted model (dump when tracing?) @@ -1178,4 +1180,9 @@ // free local allocations + psFree (mask); + psFree (chisq); + psFree (stats); + psFree (index); + return (sqrt(S2guess)); } Index: branches/eam_branches/20090715/psModules/src/objects/pmPetrosian.c =================================================================== --- branches/eam_branches/20090715/psModules/src/objects/pmPetrosian.c (revision 25452) +++ branches/eam_branches/20090715/psModules/src/objects/pmPetrosian.c (revision 25452) @@ -0,0 +1,107 @@ +/* @file pmPetrosian.c + * low-level petrosian functions + * + * @author EAM, IfA + * + * @version $Revision: $ + * @date $Date: $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +# include "pmPetrosian.h" + +static void pmPetrosianFree(pmPetrosian *petrosian) +{ + if (!petrosian) { + return; + } + psFree(petrosian->radii); + psFree(petrosian->fluxes); + psFree(petrosian->theta); + psFree(petrosian->isophotalRadii); + + psFree(petrosian->radiusElliptical); + psFree(petrosian->fluxElliptical); + + psFree(petrosian->binSB); + psFree(petrosian->binSBstdev); + psFree(petrosian->radialBins); + psFree(petrosian->area); +} + +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->radiusElliptical = NULL; + petrosian->fluxElliptical = NULL; + + petrosian->radialBins = NULL; + petrosian->area = NULL; + petrosian->binSB = NULL; + petrosian->binSBstdev = NULL; + + petrosian->petrosianRadius = NAN; + petrosian->petrosianFlux = NAN; + + // petrosian->axes = {0.0, 0.0, 0.0}; + + return petrosian; +} + +bool pmPetrosianFreeVectors(pmPetrosian *petrosian) { + + psFree(petrosian->radii); + psFree(petrosian->fluxes); + psFree(petrosian->theta); + psFree(petrosian->isophotalRadii); + + psFree(petrosian->radiusElliptical); + psFree(petrosian->fluxElliptical); + + psFree(petrosian->binSB); + psFree(petrosian->binSBstdev); + psFree(petrosian->radialBins); + psFree(petrosian->area); + + petrosian->radii = NULL; + petrosian->fluxes = NULL; + petrosian->theta = NULL; + petrosian->isophotalRadii = NULL; + + petrosian->radiusElliptical = NULL; + petrosian->fluxElliptical = NULL; + + petrosian->radialBins = NULL; + petrosian->area = NULL; + petrosian->binSB = NULL; + petrosian->binSBstdev = NULL; + + return true; +} + +# define COMPARE_INDEX(A,B) (index->data.F32[A] < index->data.F32[B]) +# define SWAP_INDEX(TYPE,A,B) { \ + float tmp; \ + if (A != B) { \ + tmp = index->data.F32[A]; \ + index->data.F32[A] = index->data.F32[B]; \ + index->data.F32[B] = tmp; \ + tmp = extra->data.F32[A]; \ + extra->data.F32[A] = extra->data.F32[B]; \ + extra->data.F32[B] = tmp; \ + } \ +} + +bool pmPetrosianSortPair (psVector *index, psVector *extra) { + + // sort the vector set by the radius + PSSORT (index->n, COMPARE_INDEX, SWAP_INDEX, NONE); + return true; +} Index: branches/eam_branches/20090715/psModules/src/objects/pmPetrosian.h =================================================================== --- branches/eam_branches/20090715/psModules/src/objects/pmPetrosian.h (revision 25452) +++ branches/eam_branches/20090715/psModules/src/objects/pmPetrosian.h (revision 25452) @@ -0,0 +1,44 @@ +/* @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 + +/// @addtogroup Objects Object Detection / Analysis Functions +/// @{ + +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; // isophotal radius for the above angles + + psVector *radiusElliptical; // normalized radial coordinates for all relevant pixels + psVector *fluxElliptical; // flux for the above radial coordinates + + psVector *binSB; // mean surface brightness within radial bins + psVector *binSBstdev; // scatter of mean surface brightness within radial bins + psVector *radialBins; // radii corresponding to above binnedBlux + psVector *area; // differential area of the non-overlapping radial bins + + psEllipseAxes axes; // shape of elliptical contour + + float petrosianRadius; + float petrosianFlux; + +} pmPetrosian; + +pmPetrosian *pmPetrosianAlloc(); + +bool pmPetrosianFreeVectors(pmPetrosian *petrosian); +bool pmPetrosianSortPair (psVector *index, psVector *extra); + +/// @} + +# endif /* PM_PETROSIAN_H */ Index: branches/eam_branches/20090715/psModules/src/objects/pmSourceExtendedPars.c =================================================================== --- branches/eam_branches/20090715/psModules/src/objects/pmSourceExtendedPars.c (revision 25451) +++ branches/eam_branches/20090715/psModules/src/objects/pmSourceExtendedPars.c (revision 25452) @@ -17,30 +17,140 @@ #endif -#include -#include -#include +// #include +// #include +// #include #include -#include "pmHDU.h" -#include "pmFPA.h" -#include "pmFPAMaskWeight.h" -#include "pmSpan.h" -#include "pmFootprint.h" -#include "pmPeaks.h" -#include "pmMoments.h" -#include "pmResiduals.h" -#include "pmGrowthCurve.h" -#include "pmTrend2D.h" -#include "pmPSF.h" -#include "pmModel.h" -#include "pmSource.h" - +// #include "pmHDU.h" +// #include "pmFPA.h" +// #include "pmFPAMaskWeight.h" +// #include "pmSpan.h" +// #include "pmFootprint.h" +// #include "pmPeaks.h" +// #include "pmMoments.h" +// #include "pmResiduals.h" +// #include "pmGrowthCurve.h" +// #include "pmTrend2D.h" +// #include "pmPSF.h" +// #include "pmModel.h" +// #include "pmSource.h" +#include "pmSourceExtendedPars.h" + +// *** pmSourceRadialProfile describes the radial profile of a source in elliptical contours, and +// intermediate data used to measure the profile +static void pmSourceRadialProfileFree(pmSourceRadialProfile *profile) +{ + if (!profile) { + return; + } + psFree(profile->radii); + psFree(profile->fluxes); + psFree(profile->theta); + psFree(profile->isophotalRadii); + + psFree(profile->radiusElliptical); + psFree(profile->fluxElliptical); + + psFree(profile->binSB); + psFree(profile->binSBstdev); + psFree(profile->binSBerror); + + psFree(profile->radialBins); + psFree(profile->area); +} + +pmSourceRadialProfile *pmSourceRadialProfileAlloc() +{ + pmSourceRadialProfile *profile = (pmSourceRadialProfile *)psAlloc(sizeof(pmSourceRadialProfile)); + psMemSetDeallocator(profile, (psFreeFunc) pmSourceRadialProfileFree); + + profile->radii = NULL; + profile->fluxes = NULL; + profile->theta = NULL; + profile->isophotalRadii = NULL; + + profile->radiusElliptical = NULL; + profile->fluxElliptical = NULL; + + profile->binSB = NULL; + profile->binSBstdev = NULL; + profile->binSBerror = NULL; + + profile->radialBins = NULL; + profile->area = NULL; + + return profile; +} + +bool psMemCheckSourceRadialProfile(psPtr ptr) +{ + PS_ASSERT_PTR(ptr, false); + return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceRadialProfileFree); +} + + +// *** pmSourceRadialProfileFreeVectors frees the intermediate data values +bool pmSourceRadialProfileFreeVectors(pmSourceRadialProfile *profile) { + + psFree(profile->radii); + psFree(profile->fluxes); + psFree(profile->theta); + psFree(profile->isophotalRadii); + + psFree(profile->radiusElliptical); + psFree(profile->fluxElliptical); + + // psFree(profile->binSB); + // psFree(profile->binSBstdev); + // psFree(profile->binSBerror); + + // psFree(profile->radialBins); + psFree(profile->area); + + profile->radii = NULL; + profile->fluxes = NULL; + profile->theta = NULL; + profile->isophotalRadii = NULL; + + profile->radiusElliptical = NULL; + profile->fluxElliptical = NULL; + + // profile->binSB = NULL; + // profile->binSBstdev = NULL; + // profile->binSBerror = NULL; + + // profile->radialBins = NULL; + profile->area = NULL; + + return true; +} + +// *** pmSourceRadialProfileSortPair is a utility function for sorting a pair of vectors +# define COMPARE_INDEX(A,B) (index->data.F32[A] < index->data.F32[B]) +# define SWAP_INDEX(TYPE,A,B) { \ + float tmp; \ + if (A != B) { \ + tmp = index->data.F32[A]; \ + index->data.F32[A] = index->data.F32[B]; \ + index->data.F32[B] = tmp; \ + tmp = extra->data.F32[A]; \ + extra->data.F32[A] = extra->data.F32[B]; \ + extra->data.F32[B] = tmp; \ + } \ +} + +bool pmSourceRadialProfileSortPair (psVector *index, psVector *extra) { + + // sort the vector set by the radius + PSSORT (index->n, COMPARE_INDEX, SWAP_INDEX, NONE); + return true; +} + +// *** pmSourceExtendedPars describes the possible collection of extended flux measurements for a source static void pmSourceExtendedParsFree (pmSourceExtendedPars *pars) { if (!pars) return; psFree(pars->profile); - psFree(pars->annuli); - psFree(pars->isophot); - psFree(pars->petrosian); - psFree(pars->kron); + psFree(pars->petrosian_50); + psFree(pars->petrosian_80); return; } @@ -51,8 +161,6 @@ pars->profile = NULL; - pars->annuli = NULL; - pars->isophot = NULL; - pars->petrosian = NULL; - pars->kron = NULL; + pars->petrosian_50 = NULL; + pars->petrosian_80 = NULL; return pars; @@ -66,135 +174,27 @@ -static void pmSourceRadialProfileFree (pmSourceRadialProfile *profile) { - if (!profile) return; - - psFree(profile->radius); - psFree(profile->flux); - psFree(profile->variance); +// *** pmSourceExtendedFlux describes the flux within an elliptical aperture of some kind +static void pmSourceExtendedFluxFree (pmSourceExtendedFlux *flux) { + if (!flux) return; return; } -pmSourceRadialProfile *pmSourceRadialProfileAlloc (void) { - - pmSourceRadialProfile *profile = (pmSourceRadialProfile *) psAlloc(sizeof(pmSourceRadialProfile)); - psMemSetDeallocator(profile, (psFreeFunc) pmSourceRadialProfileFree); - - profile->radius = NULL; - profile->flux = NULL; - profile->variance = NULL; - - return profile; -} - -bool psMemCheckSourceRadialProfile(psPtr ptr) +pmSourceExtendedFlux *pmSourceExtendedFluxAlloc (void) { + + pmSourceExtendedFlux *flux = (pmSourceExtendedFlux *) psAlloc(sizeof(pmSourceExtendedFlux)); + psMemSetDeallocator(flux, (psFreeFunc) pmSourceExtendedFluxFree); + + flux->flux = 0.0; + flux->fluxErr = 0.0; + flux->radius = 0.0; + flux->radiusErr = 0.0; + + return flux; +} + + +bool psMemCheckSourceExtendedFlux(psPtr ptr) { PS_ASSERT_PTR(ptr, false); - return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceRadialProfileFree); -} - - -static void pmSourceIsophotalValuesFree (pmSourceIsophotalValues *isophot) { - if (!isophot) return; - return; -} - -pmSourceIsophotalValues *pmSourceIsophotalValuesAlloc (void) { - - pmSourceIsophotalValues *isophot = (pmSourceIsophotalValues *) psAlloc(sizeof(pmSourceIsophotalValues)); - psMemSetDeallocator(isophot, (psFreeFunc) pmSourceIsophotalValuesFree); - - isophot->mag = 0.0; - isophot->magErr = 0.0; - isophot->rad = 0.0; - isophot->radErr = 0.0; - - return isophot; -} - - -bool psMemCheckSourceIsophotalValues(psPtr ptr) -{ - PS_ASSERT_PTR(ptr, false); - return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceIsophotalValuesFree); -} - - -static void pmSourcePetrosianValuesFree (pmSourcePetrosianValues *petrosian) { - if (!petrosian) return; - return; -} - -pmSourcePetrosianValues *pmSourcePetrosianValuesAlloc (void) { - - pmSourcePetrosianValues *petrosian = (pmSourcePetrosianValues *) psAlloc(sizeof(pmSourcePetrosianValues)); - psMemSetDeallocator(petrosian, (psFreeFunc) pmSourcePetrosianValuesFree); - - petrosian->mag = 0.0; - petrosian->magErr = 0.0; - petrosian->rad = 0.0; - petrosian->radErr = 0.0; - - return petrosian; -} - - -bool psMemCheckSourcePetrosianValues(psPtr ptr) -{ - PS_ASSERT_PTR(ptr, false); - return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourcePetrosianValuesFree); -} - -static void pmSourceKronValuesFree (pmSourceKronValues *kron) { - if (!kron) return; - return; -} - -pmSourceKronValues *pmSourceKronValuesAlloc (void) { - - pmSourceKronValues *kron = (pmSourceKronValues *) psAlloc(sizeof(pmSourceKronValues)); - psMemSetDeallocator(kron, (psFreeFunc) pmSourceKronValuesFree); - - kron->mag = 0.0; - kron->magErr = 0.0; - kron->rad = 0.0; - kron->radErr = 0.0; - - return kron; -} - - -bool psMemCheckSourceKronValues(psPtr ptr) -{ - PS_ASSERT_PTR(ptr, false); - return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceKronValuesFree); -} - - -static void pmSourceAnnuliFree (pmSourceAnnuli *annuli) { - if (!annuli) return; - - psFree (annuli->flux); - psFree (annuli->fluxErr); - psFree (annuli->fluxVar); - - return; -} - -pmSourceAnnuli *pmSourceAnnuliAlloc (void) { - - pmSourceAnnuli *annuli = (pmSourceAnnuli *) psAlloc(sizeof(pmSourceAnnuli)); - psMemSetDeallocator(annuli, (psFreeFunc) pmSourceAnnuliFree); - - annuli->flux = NULL; - annuli->fluxErr = NULL; - annuli->fluxVar = NULL; - - return annuli; -} - - -bool psMemCheckSourceAnnuli(psPtr ptr) -{ - PS_ASSERT_PTR(ptr, false); - return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceAnnuliFree); -} + return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmSourceExtendedFluxFree); +} Index: branches/eam_branches/20090715/psModules/src/objects/pmSourceExtendedPars.h =================================================================== --- branches/eam_branches/20090715/psModules/src/objects/pmSourceExtendedPars.h (revision 25451) +++ branches/eam_branches/20090715/psModules/src/objects/pmSourceExtendedPars.h (revision 25452) @@ -15,56 +15,53 @@ typedef struct { - psVector *radius; - psVector *flux; - psVector *variance; + 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; // isophotal radius for the above angles + + psVector *radiusElliptical; // normalized radial coordinates for all relevant pixels + psVector *fluxElliptical; // flux for the above radial coordinates + + psVector *binSB; // mean surface brightness within radial bins + psVector *binSBstdev; // scatter of mean surface brightness within radial bins + psVector *binSBerror; // formal error on mean surface brightness within radial bins + + psVector *radialBins; // radii corresponding to above binnedBlux + psVector *area; // differential area of the non-overlapping radial bins + + psEllipseAxes axes; // shape of elliptical contour } pmSourceRadialProfile; typedef struct { - psVector *flux; - psVector *fluxErr; - psVector *fluxVar; -} pmSourceAnnuli; - -typedef struct { - float mag; - float magErr; - float rad; - float radErr; -} pmSourceIsophotalValues; - -typedef struct { - float mag; - float magErr; - float rad; - float radErr; -} pmSourcePetrosianValues; - -typedef struct { - float mag; - float magErr; - float rad; - float radErr; -} pmSourceKronValues; + float flux; + float fluxErr; + float radius; + float radiusErr; +} pmSourceExtendedFlux; typedef struct { pmSourceRadialProfile *profile; - pmSourceAnnuli *annuli; - pmSourceIsophotalValues *isophot; - pmSourcePetrosianValues *petrosian; - pmSourceKronValues *kron; + pmSourceExtendedFlux *petrosian_50; + pmSourceExtendedFlux *petrosian_80; } pmSourceExtendedPars; -pmSourceExtendedPars *pmSourceExtendedParsAlloc(void); +// *** pmSourceRadialProfile describes the radial profile of a source in elliptical contours, and +// intermediate data used to measure the profile +pmSourceRadialProfile *pmSourceRadialProfileAlloc(); +bool psMemCheckSourceRadialProfile(psPtr ptr); + +// *** pmSourceRadialProfileFreeVectors frees the intermediate data values +bool pmSourceRadialProfileFreeVectors(pmSourceRadialProfile *profile); + +// *** pmSourceExtendedPars describes the possible collection of extended flux measurements for a source +pmSourceExtendedPars *pmSourceExtendedParsAlloc (void); bool psMemCheckSourceExtendedPars(psPtr ptr); -pmSourceRadialProfile *pmSourceRadialProfileAlloc(void); -bool psMemCheckSourceRadialProfile(psPtr ptr); -pmSourceIsophotalValues *pmSourceIsophotalValuesAlloc(void); -bool psMemCheckSourceIsophotalValues(psPtr ptr); -pmSourcePetrosianValues *pmSourcePetrosianValuesAlloc(void); -bool psMemCheckSourcePetrosianValues(psPtr ptr); -pmSourceKronValues *pmSourceKronValuesAlloc(void); -bool psMemCheckSourceKronValues(psPtr ptr); -pmSourceAnnuli *pmSourceAnnuliAlloc(void); -bool psMemCheckSourceAnnuli(psPtr ptr); + +// *** pmSourceExtendedFlux describes the flux within an elliptical aperture of some kind +pmSourceExtendedFlux *pmSourceExtendedFluxAlloc(void); +bool psMemCheckSourceExtendedFlux(psPtr ptr); + +// *** pmSourceRadialProfileSortPair is a utility function for sorting a pair of vectors +bool pmSourceRadialProfileSortPair(psVector *index, psVector *extra); /// @} Index: branches/eam_branches/20090715/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c =================================================================== --- branches/eam_branches/20090715/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c (revision 25451) +++ branches/eam_branches/20090715/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c (revision 25452) @@ -353,8 +353,8 @@ // which extended source analyses should we perform? - bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); - bool doIsophotal = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL"); - bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); - bool doKron = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON"); + // bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); + // bool doIsophotal = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL"); + // bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); + // bool doKron = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON"); psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); @@ -396,4 +396,5 @@ psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG", PS_DATA_F32, "Sigma in EXT y coordinate", yErr); +# if (0) // Petrosian measurements // XXX insert header data: petrosian ref radius, flux ratio @@ -476,4 +477,5 @@ } +# endif psArrayAdd (table, 100, row); psFree (row); Index: branches/eam_branches/20090715/psModules/src/objects/pmSourceIO_PS1_CAL_0.c =================================================================== --- branches/eam_branches/20090715/psModules/src/objects/pmSourceIO_PS1_CAL_0.c (revision 25451) +++ branches/eam_branches/20090715/psModules/src/objects/pmSourceIO_PS1_CAL_0.c (revision 25452) @@ -349,8 +349,8 @@ // which extended source analyses should we perform? - bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); - bool doIsophotal = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL"); - bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); - bool doKron = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON"); + // bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); + // bool doIsophotal = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL"); + // bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); + // bool doKron = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON"); psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); @@ -410,4 +410,5 @@ psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG", PS_DATA_F32, "Sigma in EXT y coordinate", yErr); +# if (0) // Petrosian measurements // XXX insert header data: petrosian ref radius, flux ratio @@ -501,4 +502,5 @@ } } +# endif psArrayAdd (table, 100, row); Index: branches/eam_branches/20090715/psModules/src/objects/pmSourceIO_PS1_DEV_1.c =================================================================== --- branches/eam_branches/20090715/psModules/src/objects/pmSourceIO_PS1_DEV_1.c (revision 25451) +++ branches/eam_branches/20090715/psModules/src/objects/pmSourceIO_PS1_DEV_1.c (revision 25452) @@ -300,8 +300,8 @@ // which extended source analyses should we perform? - bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); - bool doIsophotal = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL"); - bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); - bool doKron = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON"); + // bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); + // bool doIsophotal = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL"); + // bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); + // bool doKron = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON"); psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); @@ -343,4 +343,6 @@ psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG", PS_DATA_F32, "Sigma in EXT y coordinate", yErr); +// XXX disable these outputs until we clean up the names +# if (0) // Petrosian measurements // XXX insert header data: petrosian ref radius, flux ratio @@ -422,4 +424,5 @@ } } +# endif psArrayAdd (table, 100, row); Index: branches/eam_branches/20090715/psphot/src/Makefile.am =================================================================== --- branches/eam_branches/20090715/psphot/src/Makefile.am (revision 25451) +++ branches/eam_branches/20090715/psphot/src/Makefile.am (revision 25452) @@ -25,5 +25,5 @@ libpsphot_la_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) -bin_PROGRAMS = psphot psphotTest psphotMomentsStudy psphotPetrosianStudy +bin_PROGRAMS = psphot psphotTest psphotMomentsStudy # bin_PROGRAMS = psphotPetrosianStudy # bin_PROGRAMS = psphot @@ -41,7 +41,7 @@ psphotMomentsStudy_LDADD = libpsphot.la -psphotPetrosianStudy_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) -psphotPetrosianStudy_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) -psphotPetrosianStudy_LDADD = libpsphot.la +# psphotPetrosianStudy_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) +# psphotPetrosianStudy_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) +# psphotPetrosianStudy_LDADD = libpsphot.la psphot_SOURCES = \ @@ -68,6 +68,6 @@ psphotMomentsStudy.c -psphotPetrosianStudy_SOURCES = \ - psphotPetrosianStudy.c +# psphotPetrosianStudy_SOURCES = \ +# psphotPetrosianStudy.c libpsphot_la_SOURCES = \ @@ -113,9 +113,4 @@ psphotExtendedSourceAnalysis.c \ psphotExtendedSourceFits.c \ - psphotRadialProfile.c \ - psphotPetrosian.c \ - psphotIsophotal.c \ - psphotAnnuli.c \ - psphotKron.c \ psphotKernelFromPSF.c \ psphotPSFConvModel.c \ @@ -138,19 +133,28 @@ psphotThreadTools.c \ psphotAddNoise.c \ - psphotPetrosianProfile.c \ - psphotRadialProfileByAngle.c \ + psphotRadialProfile.c \ + psphotRadialProfileByAngles.c \ psphotRadiiFromProfiles.c \ psphotEllipticalContour.c \ psphotEllipticalProfile.c \ + psphotPetrosian.c \ psphotPetrosianRadialBins.c \ - psphotPetrosianVisual.c \ psphotPetrosianStats.c \ - psphotPetrosianAnalysis.c \ - pmPetrosian.c \ psphotEfficiency.c + +# re-instate these +# psphotIsophotal.c \ +# psphotAnnuli.c \ +# psphotKron.c \ +# psphotPetrosianVisual.c \ +# + +# test versions +# psphotPetrosianProfile.c \ +# psphotPetrosianAnalysis.c \ +# include_HEADERS = \ psphot.h \ - pmPetrosian.h \ psphotErrorCodes.h Index: branches/eam_branches/20090715/psphot/src/pmPetrosian.c =================================================================== --- branches/eam_branches/20090715/psphot/src/pmPetrosian.c (revision 25451) +++ (revision ) @@ -1,107 +1,0 @@ -/* @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 *petrosian) -{ - if (!petrosian) { - return; - } - psFree(petrosian->radii); - psFree(petrosian->fluxes); - psFree(petrosian->theta); - psFree(petrosian->isophotalRadii); - - psFree(petrosian->radiusElliptical); - psFree(petrosian->fluxElliptical); - - psFree(petrosian->binSB); - psFree(petrosian->binSBstdev); - psFree(petrosian->radialBins); - psFree(petrosian->area); -} - -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->radiusElliptical = NULL; - petrosian->fluxElliptical = NULL; - - petrosian->radialBins = NULL; - petrosian->area = NULL; - petrosian->binSB = NULL; - petrosian->binSBstdev = NULL; - - petrosian->petrosianRadius = NAN; - petrosian->petrosianFlux = NAN; - - // petrosian->axes = {0.0, 0.0, 0.0}; - - return petrosian; -} - -bool psphotPetrosianFreeVectors(pmPetrosian *petrosian) { - - psFree(petrosian->radii); - psFree(petrosian->fluxes); - psFree(petrosian->theta); - psFree(petrosian->isophotalRadii); - - psFree(petrosian->radiusElliptical); - psFree(petrosian->fluxElliptical); - - psFree(petrosian->binSB); - psFree(petrosian->binSBstdev); - psFree(petrosian->radialBins); - psFree(petrosian->area); - - petrosian->radii = NULL; - petrosian->fluxes = NULL; - petrosian->theta = NULL; - petrosian->isophotalRadii = NULL; - - petrosian->radiusElliptical = NULL; - petrosian->fluxElliptical = NULL; - - petrosian->radialBins = NULL; - petrosian->area = NULL; - petrosian->binSB = NULL; - petrosian->binSBstdev = NULL; - - return true; -} - -# define COMPARE_INDEX(A,B) (index->data.F32[A] < index->data.F32[B]) -# define SWAP_INDEX(TYPE,A,B) { \ - float tmp; \ - if (A != B) { \ - tmp = index->data.F32[A]; \ - index->data.F32[A] = index->data.F32[B]; \ - index->data.F32[B] = tmp; \ - tmp = extra->data.F32[A]; \ - extra->data.F32[A] = extra->data.F32[B]; \ - extra->data.F32[B] = tmp; \ - } \ -} - -bool psphotPetrosianSortPair (psVector *index, psVector *extra) { - - // sort the vector set by the radius - PSSORT (index->n, COMPARE_INDEX, SWAP_INDEX, NONE); - return true; -} Index: branches/eam_branches/20090715/psphot/src/pmPetrosian.h =================================================================== --- branches/eam_branches/20090715/psphot/src/pmPetrosian.h (revision 25451) +++ (revision ) @@ -1,62 +1,0 @@ -/* @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; // isophotal radius for the above angles - - psVector *radiusElliptical; // normalized radial coordinates for all relevant pixels - psVector *fluxElliptical; // flux for the above radial coordinates - - psVector *binSB; // mean surface brightness within radial bins - psVector *binSBstdev; // scatter of mean surface brightness within radial bins - psVector *radialBins; // radii corresponding to above binnedBlux - psVector *area; // differential area of the non-overlapping radial bins - - psEllipseAxes axes; // shape of elliptical contour - - float petrosianRadius; - float petrosianFlux; - -} pmPetrosian; - -pmPetrosian *pmPetrosianAlloc(); -bool psphotPetrosianFreeVectors(pmPetrosian *petrosian); - -bool psphotPetrosianProfile (pmReadout *readout, pmSource *source, float skynoise); -bool psphotRadialProfilesByAngles (pmSource *source, pmPetrosian *petro, int Nsec, float Rmax); -float psphotRadiusFromProfile (pmSource *source, psVector *radius, psVector *flux, float fluxMin, float fluxMax); -bool psphotRadiiFromProfiles (pmSource *source, pmPetrosian *petrosian, float fluxMin, float fluxMax); -bool psphotEllipticalProfile (pmSource *source, pmPetrosian *petrosian); -bool psphotEllipticalContour (pmSource *source, pmPetrosian *petrosian); -bool psphotPetrosianRadialBins (pmSource *source, pmPetrosian *petrosian, float radiusMax, float skynoise); -bool psphotPetrosianStats (pmSource *source, pmPetrosian *petrosian); - -bool psphotPetrosianSortPair (psVector *index, psVector *extra); - - -bool psphotPetrosianVisualProfileByAngle (psVector *radius, psVector *flux); -bool psphotPetrosianVisualProfileRadii (psVector *radius, psVector *flux, psVector *radiusBin, psVector *fluxBin, float peakFlux, float RadiusRef); -bool psphotPetrosianVisualEllipticalContour (pmPetrosian *petrosian); - -bool psphotPetrosianVisualStats (psVector *radBin, psVector *fluxBin, - psVector *refRadius, psVector *meanSB, - psVector *petRatio, psVector *petRatioErr, psVector *fluxSum, - float petRadius, float ratioForRadius, - float petFlux, float radiusForFlux); - -/// @} - - -# endif /* PM_PETROSIAN_H */ Index: branches/eam_branches/20090715/psphot/src/psphot.h =================================================================== --- branches/eam_branches/20090715/psphot/src/psphot.h (revision 25451) +++ branches/eam_branches/20090715/psphot/src/psphot.h (revision 25452) @@ -8,5 +8,4 @@ #include #include "psphotErrorCodes.h" -#include "pmPetrosian.h" #define PSPHOT_RECIPE "PSPHOT" // Name of the recipe to use @@ -169,9 +168,11 @@ psKernel *psphotKernelFromPSF (pmSource *source, int nPix); -bool psphotRadialProfile (pmSource *source, psMetadata *recipe, psImageMaskType maskVal); -bool psphotPetrosian (pmSource *source, psMetadata *recipe, psImageMaskType maskVal); -bool psphotIsophotal (pmSource *source, psMetadata *recipe, psImageMaskType maskVal); -bool psphotAnnuli (pmSource *source, psMetadata *recipe, psImageMaskType maskVal); -bool psphotKron (pmSource *source, psMetadata *recipe, psImageMaskType maskVal); +// functions related to extended source analysis +bool psphotRadialProfile (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal); +bool psphotRadialProfilesByAngles (pmSource *source, int Nsec, float Rmax); +float psphotRadiusFromProfile (pmSource *source, psVector *radius, psVector *flux, float fluxMin, float fluxMax); +bool psphotRadiiFromProfiles (pmSource *source, float fluxMin, float fluxMax); +bool psphotEllipticalProfile (pmSource *source); +bool psphotEllipticalContour (pmSource *source); // psphotVisual functions @@ -194,7 +195,27 @@ bool psphotVisualPlotApResid (psArray *sources); bool psphotVisualPlotSourceSize (psMetadata *recipe, psArray *sources); - -bool psphotVisualShowPetrosian (pmSource *source, pmPetrosian *petrosian); -bool psphotPetrosianAnalysis (pmReadout *readout, psArray *sources, psMetadata *recipe); +bool psphotVisualShowPetrosians (psArray *sources); + +// bool psphotPetrosianAnalysis (pmReadout *readout, psArray *sources, psMetadata *recipe); +// bool psphotPetrosianProfile (pmReadout *readout, pmSource *source, float skynoise); + +bool psphotPetrosian (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal); +bool psphotPetrosianRadialBins (pmSource *source, float radiusMax, float skynoise); +bool psphotPetrosianStats (pmSource *source); + +// XXX old versions, currently disabled +// bool psphotIsophotal (pmSource *source, psMetadata *recipe, psImageMaskType maskVal); +// bool psphotAnnuli (pmSource *source, psMetadata *recipe, psImageMaskType maskVal); +// bool psphotKron (pmSource *source, psMetadata *recipe, psImageMaskType maskVal); + +// XXX visualization functions related to radial profiles (disabled) +// bool psphotPetrosianVisualProfileByAngle (psVector *radius, psVector *flux); +// bool psphotPetrosianVisualProfileRadii (psVector *radius, psVector *flux, psVector *radiusBin, psVector *fluxBin, float peakFlux, float RadiusRef); +// bool psphotPetrosianVisualEllipticalContour (pmPetrosian *petrosian); +// bool psphotPetrosianVisualStats (psVector *radBin, psVector *fluxBin, +// psVector *refRadius, psVector *meanSB, +// psVector *petRatio, psVector *petRatioErr, psVector *fluxSum, +// float petRadius, float ratioForRadius, +// float petFlux, float radiusForFlux); bool psphotImageQuality (psMetadata *recipe, psArray *sources); Index: branches/eam_branches/20090715/psphot/src/psphotEllipticalContour.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotEllipticalContour.c (revision 25451) +++ branches/eam_branches/20090715/psphot/src/psphotEllipticalContour.c (revision 25452) @@ -5,9 +5,11 @@ psF32 psphotEllipticalContourFunc (psVector *deriv, const psVector *params, const psVector *coord); -bool psphotEllipticalContour (pmSource *source, pmPetrosian *petrosian) { +bool psphotEllipticalContour (pmSource *source) { + + pmSourceRadialProfile *profile = source->extpars->profile; // use LMM to fit theta vs radius to an ellipse - psVector *theta = petrosian->theta; - psVector *radius = petrosian->isophotalRadii; + psVector *theta = profile->theta; + psVector *radius = profile->isophotalRadii; // find Rmin and Rmax for the initial guess @@ -83,17 +85,17 @@ /// XXX rationalize? if epsilon > 1, flip major and minor axes (rotate by 90 degrees) if (params->data.F32[PAR_EPSILON] < 1.0) { - 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]; + profile->axes.major = params->data.F32[PAR_RMIN] / params->data.F32[PAR_EPSILON]; + profile->axes.minor = params->data.F32[PAR_RMIN]; + profile->axes.theta = params->data.F32[PAR_PHI]; } else { - petrosian->axes.major = params->data.F32[PAR_RMIN]; - petrosian->axes.minor = params->data.F32[PAR_RMIN] / params->data.F32[PAR_EPSILON]; - petrosian->axes.theta = params->data.F32[PAR_PHI] + 0.5*M_PI; + profile->axes.major = params->data.F32[PAR_RMIN]; + profile->axes.minor = params->data.F32[PAR_RMIN] / params->data.F32[PAR_EPSILON]; + profile->axes.theta = params->data.F32[PAR_PHI] + 0.5*M_PI; } psTrace ("psphot", 4, "# fitted values:\n"); - psTrace ("psphot", 4, "Phi: %f\n", petrosian->axes.theta*PS_DEG_RAD); - psTrace ("psphot", 4, "Rmaj: %f\n", petrosian->axes.major); - psTrace ("psphot", 4, "Rmin: %f\n", petrosian->axes.minor); + psTrace ("psphot", 4, "Phi: %f\n", profile->axes.theta*PS_DEG_RAD); + psTrace ("psphot", 4, "Rmaj: %f\n", profile->axes.major); + psTrace ("psphot", 4, "Rmin: %f\n", profile->axes.minor); // show the results Index: branches/eam_branches/20090715/psphot/src/psphotEllipticalProfile.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotEllipticalProfile.c (revision 25451) +++ branches/eam_branches/20090715/psphot/src/psphotEllipticalProfile.c (revision 25452) @@ -1,11 +1,13 @@ # include "psphotInternal.h" -bool psphotEllipticalProfile (pmSource *source, pmPetrosian *petrosian) { +bool psphotEllipticalProfile (pmSource *source) { - petrosian->radiusElliptical = psVectorAllocEmpty(100, PS_TYPE_F32); - petrosian->fluxElliptical = psVectorAllocEmpty(100, PS_TYPE_F32); + pmSourceRadialProfile *profile = source->extpars->profile; - psVector *radius = petrosian->radiusElliptical; - psVector *flux = petrosian->fluxElliptical; + profile->radiusElliptical = psVectorAllocEmpty(100, PS_TYPE_F32); + profile->fluxElliptical = psVectorAllocEmpty(100, PS_TYPE_F32); + + psVector *radius = profile->radiusElliptical; + psVector *flux = profile->fluxElliptical; // the psEllipse functions use z = 0.5(x/Sxx)^2 + 0.5(y/Syy)^2 + x y Sxy @@ -20,10 +22,10 @@ psEllipseAxes axes; axes.major = M_SQRT1_2; - axes.minor = M_SQRT1_2 * (petrosian->axes.minor / petrosian->axes.major); + axes.minor = M_SQRT1_2 * (profile->axes.minor / profile->axes.major); // axes.major = 1.0; - // axes.minor = petrosian->axes.minor / petrosian->axes.major; + // axes.minor = profile->axes.minor / profile->axes.major; - axes.theta = petrosian->axes.theta; + axes.theta = profile->axes.theta; psEllipseShape shape = psEllipseAxesToShape (axes); @@ -56,7 +58,7 @@ // psVector *radiusRaw = psVectorAllocEmpty(100, PS_TYPE_F32); // psVector *fluxRaw = psVectorAllocEmpty(100, PS_TYPE_F32); - // for (int i = 0; i < petrosian->radii->n; i++) { - // psVector *r = petrosian->radii->data[i]; - // psVector *f = petrosian->fluxes->data[i]; + // for (int i = 0; i < profile->radii->n; i++) { + // psVector *r = profile->radii->data[i]; + // psVector *f = profile->fluxes->data[i]; // for (int j = 0; j < r->n; j++) { // psVectorAppend(radiusRaw, r->data.F32[j]); Index: branches/eam_branches/20090715/psphot/src/psphotExtendedSourceAnalysis.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotExtendedSourceAnalysis.c (revision 25451) +++ branches/eam_branches/20090715/psphot/src/psphotExtendedSourceAnalysis.c (revision 25452) @@ -21,8 +21,15 @@ } - // option to limit analysis to a specific region - char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION"); - psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region)); - if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined"); + // XXX require petrosian analysis for non-linear fits? + + // XXX temporary user-supplied systematic sky noise measurement (derive from background model) + float skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE"); + +# if (0) + // if backModel or backStdev are missing, the values of sky and/or skyErr will be set to NAN + // XXX use this to set skynoise + pmReadout *backModel = psphotSelectBackground (config, view); + pmReadout *backStdev = psphotSelectBackgroundStdev (config, view); +# endif // S/N limit to perform full non-linear fits @@ -38,5 +45,8 @@ sources = psArraySort (sources, pmSourceSortBySN); - // XXX some init functions for the extended source recipe options? + // option to limit analysis to a specific region + char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION"); + psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region)); + if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined"); // choose the sources of interest @@ -49,4 +59,7 @@ if (source->type == PM_SOURCE_TYPE_DEFECT) continue; if (source->type == PM_SOURCE_TYPE_SATURATED) continue; + if (source->mode & PM_SOURCE_MODE_DEFECT) continue; + if (source->mode & PM_SOURCE_MODE_SATSTAR) continue; + if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue; // limit selection to some SN limit @@ -68,5 +81,5 @@ // if we request any of these measurements, we require the radial profile if (doPetrosian || doIsophotal || doAnnuli || doKron) { - if (!psphotRadialProfile (source, recipe, maskVal)) { + if (!psphotRadialProfile (source, recipe, skynoise, maskVal)) { // all measurements below require the radial profile; skip them all // re-subtract the object, leave local sky @@ -79,4 +92,16 @@ } + // Petrosian Mags + if (doPetrosian) { + if (!psphotPetrosian (source, recipe, skynoise, maskVal)) { + psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); + } else { + psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); + Npetro ++; + source->mode |= PM_SOURCE_MODE_EXTENDED_STATS; + } + } + +# if (0) // Isophotal Mags if (doIsophotal) { @@ -89,16 +114,4 @@ } } - - // Petrosian Mags - if (doPetrosian) { - if (!psphotPetrosian (source, recipe, maskVal)) { - psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); - } else { - psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); - Npetro ++; - source->mode |= PM_SOURCE_MODE_EXTENDED_STATS; - } - } - // Kron Mags if (doKron) { @@ -111,19 +124,13 @@ } } - - // Radial Annuli - if (doAnnuli) { - if (!psphotAnnuli (source, recipe, maskVal)) { - psError(PSPHOT_ERR_UNKNOWN, false, "failure in Annuli analysis"); - return false; - } - psTrace ("psphot", 5, "measured annuli for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); - Nannuli ++; - source->mode |= PM_SOURCE_MODE_EXTENDED_STATS; - } +# endif // re-subtract the object, leave local sky pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED; + + if (source->extpars) { + pmSourceRadialProfileFreeVectors(source->extpars->profile); + } } @@ -133,4 +140,11 @@ psLogMsg ("psphot", PS_LOG_INFO, " %d annuli\n", Nannuli); psLogMsg ("psphot", PS_LOG_INFO, " %d kron\n", Nkron); + + psphotVisualShowResidualImage (readout); + + if (doPetrosian) { + psphotVisualShowPetrosians (sources); + } + return true; } Index: branches/eam_branches/20090715/psphot/src/psphotFitSourcesLinear.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotFitSourcesLinear.c (revision 25451) +++ branches/eam_branches/20090715/psphot/src/psphotFitSourcesLinear.c (revision 25452) @@ -76,5 +76,5 @@ if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) continue; } else { - if (source->mode & PM_SOURCE_MODE_BLEND) continue; + // if (source->mode & PM_SOURCE_MODE_BLEND) continue; } Index: branches/eam_branches/20090715/psphot/src/psphotPetrosian.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotPetrosian.c (revision 25451) +++ branches/eam_branches/20090715/psphot/src/psphotPetrosian.c (revision 25452) @@ -1,109 +1,32 @@ # include "psphotInternal.h" -bool psphotPetrosian (pmSource *source, psMetadata *recipe, psImageMaskType maskVal) { +bool psphotPetrosian (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal) { - bool status; + // XXX these need to go into recipe values + float Rmax = 200; - assert (source->extpars); - assert (source->extpars->profile); - assert (source->extpars->profile->radius); - assert (source->extpars->profile->flux); + psAssert (source->extpars, "need to run psphotRadialProfile first"); + psAssert (source->extpars->profile, "need to run psphotRadialProfile first"); - psVector *radius = source->extpars->profile->radius; - psVector *flux = source->extpars->profile->flux; + // integrate the radial profile for radial bins defined for the petrosian measurement: + // SB_i (r_i) where \alpha r_i < r < \beta r_i + if (!psphotPetrosianRadialBins (source, Rmax, skynoise)) { + psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile"); + return false; + } + + // use the SB_i from above to calculate the petrosian radius and the flux within that radius + if (!psphotPetrosianStats (source)) { + psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile"); + return false; + } + + psTrace ("psphot", 3, "source at %f,%f: petrosian radius: %f, flux: %f, axis ratio: %f, angle: %f", + source->peak->xf, source->peak->yf, + source->extpars->petrosian_80->radius, + source->extpars->petrosian_80->flux, + source->extpars->profile->axes.minor/source->extpars->profile->axes.major, + source->extpars->profile->axes.theta*PS_DEG_RAD); - // flux at which to measure isophotal parameters - float PETROSIAN_R0 = psMetadataLookupF32 (&status, recipe, "PETROSIAN_R0"); - float PETROSIAN_RF = psMetadataLookupF32 (&status, recipe, "PETROSIAN_FLUX_RATIO"); - assert (status); - - // first find flux at R0 - int firstAbove = -1; - int lastBelow = -1; - for (int i = 0; i < radius->n; i++) { - if (radius->data.F32[i] < PETROSIAN_R0) lastBelow = i; - if ((firstAbove < 0) && (radius->data.F32[i] > PETROSIAN_R0)) firstAbove = i; - } - // if we don't go out far enough, we have a problem... - if (lastBelow == radius->n - 1) { - psTrace ("psphot", 5, "did not go out far enough to reach petrosian reference radius..."); - // XXX skip object? raise a flag ? - return false; - } - if (firstAbove < 0) { - psTrace ("psphot", 5, "did not go out far enough to bound petrosian reference radius"); - // XXX raise a flag ? - return false; - } - - // average flux in this range - float fluxR0 = 0.0; - int fluxRn = 0; - for (int i = PS_MIN(firstAbove, lastBelow); i <= PS_MAX(firstAbove, lastBelow); i++) { - fluxR0 += flux->data.F32[i]; - fluxRn ++; - } - fluxR0 /= (float)(fluxRn); - - // target flux for petrosian radius - float fluxRP = fluxR0 * PETROSIAN_RF; - - // find the first bin below the flux level and the last above the level - // XXX can this be done faster with bisection? - // XXX do I need to worry about crazy outliers? - // XXX should i be smoothing or fitting the curve? - int firstBelow = -1; - int lastAbove = -1; - for (int i = 0; i < flux->n; i++) { - if (flux->data.F32[i] > fluxRP) lastAbove = i; - if ((firstBelow < 0) && (flux->data.F32[i] < fluxRP)) firstBelow = i; - } - // if we don't go out far enough, we have a problem... - if (lastAbove == radius->n - 1) { - psTrace ("psphot", 5, "did not go out far enough to reach petrosian radius..."); - // XXX skip object? raise a flag ? - return false; - } - if (firstBelow < 0) { - psTrace ("psphot", 5, "did not go out far enough to bound petrosian radius"); - // XXX raise a flag ? - return false; - } - - // need to examine pixels in this vicinity - float fluxFirst = 0; - float fluxLast = 0; - for (int i = 0; i <= PS_MAX(firstBelow, lastAbove); i++) { - if (i <= firstBelow) { - fluxFirst += flux->data.F32[i]; - } - if (i <= lastAbove) { - fluxLast += flux->data.F32[i]; - } - } - float fluxRPSum = 0.5*(fluxLast + fluxFirst); - float fluxRPSumErr = 0.5*fabs(fluxLast - fluxFirst); - // XXX need to use the weight appropriately here... - - float rad = 0.5*(radius->data.F32[firstBelow] + radius->data.F32[lastAbove]); - float radErr = 0.5*fabs(radius->data.F32[firstBelow] - radius->data.F32[lastAbove]); - - if (!source->extpars->petrosian) { - source->extpars->petrosian = pmSourcePetrosianValuesAlloc (); - } - - // these are uncalibrated: instrumental mags and pixel units - source->extpars->petrosian->mag = -2.5*log10(fluxRPSum); - source->extpars->petrosian->magErr = fluxRPSumErr / fluxRPSum; - - source->extpars->petrosian->rad = rad; - source->extpars->petrosian->radErr = radErr; - - psTrace ("psphot", 5, "Petrosian flux:%f +/- %f @ %f +/- %f for %f, %f\n", - source->extpars->petrosian->mag, source->extpars->petrosian->magErr, - source->extpars->petrosian->rad, source->extpars->petrosian->radErr, - source->peak->xf, source->peak->yf); - - return true; - + return true; } Index: branches/eam_branches/20090715/psphot/src/psphotPetrosianAnalysis.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotPetrosianAnalysis.c (revision 25451) +++ branches/eam_branches/20090715/psphot/src/psphotPetrosianAnalysis.c (revision 25452) @@ -66,13 +66,4 @@ psphotVisualShowResidualImage (readout); - - // pause and wait for user input: - // continue, save (provide name), ?? - char key[10]; - fprintf (stdout, "[c]ontinue? "); - if (!fgets(key, 8, stdin)) { - psWarning("Unable to read option"); - } - return true; } Index: branches/eam_branches/20090715/psphot/src/psphotPetrosianProfile.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotPetrosianProfile.c (revision 25451) +++ branches/eam_branches/20090715/psphot/src/psphotPetrosianProfile.c (revision 25452) @@ -12,4 +12,5 @@ pmPetrosian *petrosian = pmPetrosianAlloc(); + // XXX these need to go into recipe values int Nsec = 24; float Rmax = 200; Index: branches/eam_branches/20090715/psphot/src/psphotPetrosianRadialBins.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotPetrosianRadialBins.c (revision 25451) +++ branches/eam_branches/20090715/psphot/src/psphotPetrosianRadialBins.c (revision 25452) @@ -13,15 +13,16 @@ // track the non-overlapping radius values. -bool psphotPetrosianRadialBins (pmSource *source, pmPetrosian *petrosian, float radiusMax, float skynoise) { +// XXX move the resulting elements from profile to extpars->petrosian? +bool psphotPetrosianRadialBins (pmSource *source, float radiusMax, float skynoise) { - // XXX for testing, let's just set this to a value - + pmSourceRadialProfile *profile = source->extpars->profile; + float skyModelErrorSQ = PS_SQR(skynoise); - psVector *radius = petrosian->radiusElliptical; - psVector *flux = petrosian->fluxElliptical; + psVector *radius = profile->radiusElliptical; + psVector *flux = profile->fluxElliptical; // sort incoming vectors by radius - psphotPetrosianSortPair (radius, flux); + pmSourceRadialProfileSortPair (radius, flux); int nMax = radiusMax; @@ -163,8 +164,8 @@ // save the vectors - petrosian->radialBins = binRad; - petrosian->area = binArea; - petrosian->binSB = binSB; - petrosian->binSBstdev = binSBstdev; + profile->radialBins = binRad; + profile->area = binArea; + profile->binSB = binSB; + profile->binSBstdev = binSBstdev; // psphotPetrosianVisualProfileRadii (radius, flux, binRad, binSB, source->peak->flux, 0.0); Index: branches/eam_branches/20090715/psphot/src/psphotPetrosianStats.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotPetrosianStats.c (revision 25451) +++ branches/eam_branches/20090715/psphot/src/psphotPetrosianStats.c (revision 25452) @@ -8,12 +8,14 @@ float InterpolateValues (float X0, float Y0, float X1, float Y1, float X); -bool psphotPetrosianStats (pmSource *source, pmPetrosian *petrosian) { +bool psphotPetrosianStats (pmSource *source) { + + pmSourceRadialProfile *profile = source->extpars->profile; float petRadius, petFlux; - psVector *binSB = petrosian->binSB; - psVector *binSBstdev = petrosian->binSBstdev; - psVector *binRad = petrosian->radialBins; - psVector *area = petrosian->area; + psVector *binSB = profile->binSB; + psVector *binSBstdev = profile->binSBstdev; + psVector *binRad = profile->radialBins; + psVector *area = profile->area; psVector *fluxSum = psVectorAllocEmpty(binSB->n, PS_TYPE_F32); @@ -139,10 +141,14 @@ } - // save petRadius, petFlux + if (!source->extpars->petrosian_80) { + source->extpars->petrosian_80 = pmSourceExtendedFluxAlloc (); + } + pmSourceExtendedFlux *petrosian = source->extpars->petrosian_80; + // XXX save flags (anyPetro, manyPetro) - petrosian->petrosianRadius = petRadius; - petrosian->petrosianFlux = petFlux; + petrosian->radius = petRadius; + petrosian->flux = petFlux; - psphotPetrosianVisualStats (binRad, binSB, refRadius, meanSB, petRatio, petRatioErr, fluxSum, petRadius, PETROSIAN_RATIO, petFlux, apRadius); + // psphotPetrosianVisualStats (binRad, binSB, refRadius, meanSB, petRatio, petRatioErr, fluxSum, petRadius, PETROSIAN_RATIO, petFlux, apRadius); psFree(fluxSum); Index: branches/eam_branches/20090715/psphot/src/psphotRadialProfile.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotRadialProfile.c (revision 25451) +++ branches/eam_branches/20090715/psphot/src/psphotRadialProfile.c (revision 25452) @@ -1,21 +1,5 @@ # include "psphotInternal.h" -# define COMPARE_RADIUS(A,B) (radius->data.F32[A] < radius->data.F32[B]) -# define SWAP_RADIUS(TYPE,A,B) { \ - float tmp; \ - if (A != B) { \ - tmp = radius->data.F32[A]; \ - radius->data.F32[A] = radius->data.F32[B]; \ - radius->data.F32[B] = tmp; \ - tmp = flux->data.F32[A]; \ - flux->data.F32[A] = flux->data.F32[B]; \ - flux->data.F32[B] = tmp; \ - tmp = variance->data.F32[A]; \ - variance->data.F32[A] = variance->data.F32[B]; \ - variance->data.F32[B] = tmp; \ - } \ -} - -bool psphotRadialProfile (pmSource *source, psMetadata *recipe, psImageMaskType maskVal) { +bool psphotRadialProfile (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal) { // allocate pmSourceExtendedParameters, if not already defined @@ -28,44 +12,39 @@ } - int nPts = source->pixels->numRows * source->pixels->numCols; - source->extpars->profile->radius = psVectorAllocEmpty (nPts, PS_TYPE_F32); - source->extpars->profile->flux = psVectorAllocEmpty (nPts, PS_TYPE_F32); - source->extpars->profile->variance = psVectorAllocEmpty (nPts, PS_TYPE_F32); + // XXX these need to go into recipe values + int Nsec = 24; + float Rmax = 200; + float fluxMin = 0.0; + float fluxMax = source->peak->flux; - psVector *radius = source->extpars->profile->radius; - psVector *flux = source->extpars->profile->flux; - psVector *variance = source->extpars->profile->variance; + // generate a series of radial profiles at Nsec evenly spaced angles. the profile flux + // is measured by interpolation for small radii; for large radii, the pixels in a box + // are averaged to increase the S/N + if (!psphotRadialProfilesByAngles (source, Nsec, Rmax)) { + psError (PS_ERR_UNKNOWN, false, "failed to measure radial profile for petrosian"); + return false; + } - // XXX use the extended source model here for Xo, Yo? - // XXX define a radius scaled to the elliptical contour? + // use the radial profiles to determine the radius of a given isophote. this isophote + // is used to determine the elliptical shape of the object, so it has a relatively high + // value (nominally 50% of the peak) + if (!psphotRadiiFromProfiles (source, fluxMin, fluxMax)) { + psError (PS_ERR_UNKNOWN, false, "failed to measure isophotal radii from profiles"); + return false; + } - int n = 0; - - float Xo = 0.0; - float Yo = 0.0; - - if (source->modelEXT) { - Xo = source->modelEXT->params->data.F32[PM_PAR_XPOS] - source->pixels->col0; - Yo = source->modelEXT->params->data.F32[PM_PAR_YPOS] - source->pixels->row0; - } else { - Xo = source->peak->xf - source->pixels->col0; - Yo = source->peak->yf - source->pixels->row0; + // convert the isophotal radius vs angle measurements to an elliptical contour + if (!psphotEllipticalContour (source)) { + psLogMsg ("psphot", 3, "failed to measure elliptical contour"); + return false; } - for (int iy = 0; iy < source->pixels->numRows; iy++) { - for (int ix = 0; ix < source->pixels->numCols; ix++) { - if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) continue; - radius->data.F32[n] = hypot (ix - Xo, iy - Yo) ; - flux->data.F32[n] = source->pixels->data.F32[iy][ix]; - variance->data.F32[n] = source->variance->data.F32[iy][ix]; - n++; - } + + // generate a single, normalized radial profile following the elliptical contours. + // the radius is normalized by the axis ratio so that on the major axis, 1 pixel = 1 pixel + if (!psphotEllipticalProfile (source)) { + psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile"); + return false; } - radius->n = n; - variance->n = n; - flux->n = n; - - // sort the vector set by the radius - PSSORT (radius->n, COMPARE_RADIUS, SWAP_RADIUS, NONE); - + return true; } Index: branches/eam_branches/20090715/psphot/src/psphotRadialProfileByAngle.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotRadialProfileByAngle.c (revision 25451) +++ (revision ) @@ -1,238 +1,0 @@ -# include "psphotInternal.h" - -// Given a source at (x,y), generate a collection of radial profiles at even angular separations - -// These functions are used to calculate the stats in a rectangle at arbitrary orientation. -// XXX Move these elsewhere (psLib?) -float psphotMeanSectorValue (psImage *image, float x, float y, float dL, float dW, float theta); -psVector *psphotBoxValues (psImage *image, float x0, float y0, float dL, float dW, float theta); -psVector *psphotLineValues (psImage *image, double x1, double y1, double x2, double y2, int dW); -psVector *psphotLineValuesBresen (psImage *image, int X1, int Y1, int X2, int Y2, int dW, int swapcoords); - -bool psphotRadialProfilesByAngles (pmSource *source, pmPetrosian *petrosian, int Nsec, float Rmax) { - - // we want to have an even number of sectors so we can do 180 deg symmetrizing - Nsec = (Nsec % 2) ? Nsec + 1 : Nsec; - float dtheta = 2.0*M_PI / Nsec; - - psFree(petrosian->radii); - psFree(petrosian->fluxes); - psFree(petrosian->theta); - - petrosian->radii = psArrayAllocEmpty(Nsec); - petrosian->fluxes = psArrayAllocEmpty(Nsec); - petrosian->theta = psVectorAllocEmpty(Nsec, PS_TYPE_F32); - - - for (int i = 0; i < Nsec; i++) { - - float theta = i*dtheta; - - psVector *radius = psVectorAllocEmpty(Rmax, PS_TYPE_F32); - psVector *flux = psVectorAllocEmpty(Rmax, PS_TYPE_F32); - - // Start at Xo,Yo and find the x,y locations for r_i, theta where r_i initially - // increments by 1 pixel. At large radii (r*dtheta > 2) use stats in a box rather than - // sub-pixel interpolation - - int dR = 1.0; - for (float r = 0; r < Rmax; r += dR) { - - float Xo = source->peak->xf; - float Yo = source->peak->yf; - - // Xo,Yo are referenced to pixels with bounds i+0.0, i+1.0 - float x = r * cos (theta) + Xo; - float y = r * sin (theta) + Yo; - dR = 2*(int)(0.5*r*sin(dtheta)) + 1; - - if (x < 0) goto badvalue; - if (y < 0) goto badvalue; - if (x >= source->pixels->parent->numCols) goto badvalue; - if (y >= source->pixels->parent->numRows) goto badvalue; - - float value = NAN; - if (dR < 2) { - // value is NAN if we run off the image - // 0.5 PIX: this function takes pixel coords; source peak is in pixel coords - value = psImageInterpolatePixelBilinear(x, y, source->pixels); - } else { - // 0.5 PIX: this function takes pixel coords; source peak is in pixel coords - value = psphotMeanSectorValue(source->pixels, x, y, dR, dR, theta); - } - - // keep the all values (even NAN) so all vectors are matched in length - psVectorAppend (radius, r); - psVectorAppend (flux, value); - continue; - - badvalue: - psVectorAppend (radius, r); - psVectorAppend (flux, NAN); - } - - psArrayAdd (petrosian->radii, 100, radius); - psArrayAdd (petrosian->fluxes, 100, flux); - psVectorAppend (petrosian->theta, theta); - - // psphotPetrosianVisualProfileByAngle (radius, flux); - - psFree(radius); - psFree(flux); - } - - for (int i = 0; i < Nsec / 2; i++) { - - psVector *r1 = petrosian->radii->data[i]; - psVector *r2 = petrosian->radii->data[i+Nsec/2]; - - psVector *f1 = petrosian->fluxes->data[i]; - psVector *f2 = petrosian->fluxes->data[i+Nsec/2]; - - psAssert (r1->n == r2->n, "mis-matched vectors"); - psAssert (f1->n == f2->n, "mis-matched vectors"); - - // we have a pair of vectors i, i+Nsec/2; replace them with the finite minimum of the pair - for (int j = 0; j < r1->n; j++) { - - float flux; - - if (!isfinite(f1->data.F32[j]) && !isfinite(f2->data.F32[j])) { - flux = NAN; - goto setflux; - } - - if (!isfinite(f1->data.F32[j])) { - flux = f2->data.F32[j]; - goto setflux; - } - if (!isfinite(f2->data.F32[j])) { - flux = f1->data.F32[j]; - goto setflux; - } - - flux = PS_MIN(f1->data.F32[j], f2->data.F32[j]); - - setflux: - f1->data.F32[j] = flux; - f2->data.F32[j] = flux; - } - } - return true; -} - -float psphotMeanSectorValue (psImage *image, float x, float y, float dL, float dW, float theta) { - - psVector *values = psphotBoxValues (image, x, y, dL, dW, theta); - if (!values) goto escape; - if (!values->n) goto escape; - - psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); - psVectorStats (stats, values, NULL, NULL, 0); - - float value = stats->sampleMedian; - - psFree (stats); - psFree (values); - - return value; - -escape: - psFree(values); - return NAN; -} - -psVector *psphotBoxValues (psImage *image, float x0, float y0, float dL, float dW, float theta) { - - // extract pixels from a series of lines (from -0.5*dW to +0.5*dW) of length dL, - // centered on x0, y0 in parent pixel coordinates (not pixel indicies) - - float xs = x0 - image->col0 - 0.5*dL*cos(theta); - float ys = y0 - image->row0 - 0.5*dL*sin(theta); - - float xe = xs + 0.5*dL*cos(theta); - float ye = ys + 0.5*dL*sin(theta); - - psVector *values = psphotLineValues (image, xs, ys, xe, ye, (int) dW); - return values; -} - -/** - * identify the quadrant and draw the correct line - */ -psVector *psphotLineValues (psImage *image, double x1, double y1, double x2, double y2, int dW) { - - int FlipDirect, FlipCoords; - int X1, Y1, X2, Y2, dX, dY; - - /* rather than draw the line from float positions, we find the closest - integer end-points and draw the line between those pixels */ - - X1 = ROUND(x1); - Y1 = ROUND(y1); - X2 = ROUND(x2); - Y2 = ROUND(y2); - - dX = X2 - X1; - dY = Y2 - Y1; - - FlipCoords = (abs(dX) < abs(dY)); - FlipDirect = FlipCoords ? (y1 > y2) : (x1 > x2); - - psVector *values = NULL; - if (!FlipDirect && !FlipCoords) values = psphotLineValuesBresen (image, X1, Y1, X2, Y2, dW, FALSE); - if ( FlipDirect && !FlipCoords) values = psphotLineValuesBresen (image, X2, Y2, X1, Y1, dW, FALSE); - if (!FlipDirect && FlipCoords) values = psphotLineValuesBresen (image, Y1, X1, Y2, X2, dW, TRUE); - if ( FlipDirect && FlipCoords) values = psphotLineValuesBresen (image, Y2, X2, Y1, X1, dW, TRUE); - - return values; -} - -/** - * use the Bresenham line drawing technique - * integer-only Bresenham line-draw version which is fast - */ -psVector *psphotLineValuesBresen (psImage *image, int X1, int Y1, int X2, int Y2, int dW, int swapcoords) { - - int X, Y, dX, dY; - int e, e2; - - psVector *values = psVectorAllocEmpty(100, PS_TYPE_F32); - - dX = X2 - X1; - dY = Y2 - Y1; - - Y = Y1; - e = 0; - for (X = X1; X <= X2; X++) { - if (X > 0) { - if (swapcoords) { - if (X >= image->numRows) continue; - for (int y = Y - dW; y <= Y + dW; y++) { - if (y < 0) continue; - if (y >= image->numCols) continue; - psVectorAppend(values, image->data.F32[X][y]); - } - } else { - if (X >= image->numCols) continue; - for (int y = Y - dW; y <= Y + dW; y++) { - if (y < 0) continue; - if (y >= image->numRows) continue; - psVectorAppend(values, image->data.F32[y][X]); - } - } - } - e += dY; - e2 = 2 * e; - if (e2 > dX) { - Y++; - e -= dX; - } - if (e2 < -dX) { - Y--; - e += dX; - } - } - return values; -} - Index: branches/eam_branches/20090715/psphot/src/psphotRadialProfileByAngles.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotRadialProfileByAngles.c (revision 25452) +++ branches/eam_branches/20090715/psphot/src/psphotRadialProfileByAngles.c (revision 25452) @@ -0,0 +1,239 @@ +# include "psphotInternal.h" + +// Given a source at (x,y), generate a collection of radial profiles at even angular separations + +// These functions are used to calculate the stats in a rectangle at arbitrary orientation. +// XXX Move these elsewhere (psLib?) +float psphotMeanSectorValue (psImage *image, float x, float y, float dL, float dW, float theta); +psVector *psphotBoxValues (psImage *image, float x0, float y0, float dL, float dW, float theta); +psVector *psphotLineValues (psImage *image, double x1, double y1, double x2, double y2, int dW); +psVector *psphotLineValuesBresen (psImage *image, int X1, int Y1, int X2, int Y2, int dW, int swapcoords); + +bool psphotRadialProfilesByAngles (pmSource *source, int Nsec, float Rmax) { + + // we want to have an even number of sectors so we can do 180 deg symmetrizing + Nsec = (Nsec % 2) ? Nsec + 1 : Nsec; + float dtheta = 2.0*M_PI / Nsec; + + pmSourceRadialProfile *profile = source->extpars->profile; + psFree(profile->radii); + psFree(profile->fluxes); + psFree(profile->theta); + + profile->radii = psArrayAllocEmpty(Nsec); + profile->fluxes = psArrayAllocEmpty(Nsec); + profile->theta = psVectorAllocEmpty(Nsec, PS_TYPE_F32); + + + for (int i = 0; i < Nsec; i++) { + + float theta = i*dtheta; + + psVector *radius = psVectorAllocEmpty(Rmax, PS_TYPE_F32); + psVector *flux = psVectorAllocEmpty(Rmax, PS_TYPE_F32); + + // Start at Xo,Yo and find the x,y locations for r_i, theta where r_i initially + // increments by 1 pixel. At large radii (r*dtheta > 2) use stats in a box rather than + // sub-pixel interpolation + + int dR = 1.0; + for (float r = 0; r < Rmax; r += dR) { + + float Xo = source->peak->xf; + float Yo = source->peak->yf; + + // Xo,Yo are referenced to pixels with bounds i+0.0, i+1.0 + float x = r * cos (theta) + Xo; + float y = r * sin (theta) + Yo; + dR = 2*(int)(0.5*r*sin(dtheta)) + 1; + + if (x < 0) goto badvalue; + if (y < 0) goto badvalue; + if (x >= source->pixels->parent->numCols) goto badvalue; + if (y >= source->pixels->parent->numRows) goto badvalue; + + float value = NAN; + if (dR < 2) { + // value is NAN if we run off the image + // 0.5 PIX: this function takes pixel coords; source peak is in pixel coords + value = psImageInterpolatePixelBilinear(x, y, source->pixels); + } else { + // 0.5 PIX: this function takes pixel coords; source peak is in pixel coords + value = psphotMeanSectorValue(source->pixels, x, y, dR, dR, theta); + } + + // keep the all values (even NAN) so all vectors are matched in length + psVectorAppend (radius, r); + psVectorAppend (flux, value); + continue; + + badvalue: + psVectorAppend (radius, r); + psVectorAppend (flux, NAN); + } + + psArrayAdd (profile->radii, 100, radius); + psArrayAdd (profile->fluxes, 100, flux); + psVectorAppend (profile->theta, theta); + + // psphotPetrosianVisualProfileByAngle (radius, flux); + + psFree(radius); + psFree(flux); + } + + for (int i = 0; i < Nsec / 2; i++) { + + psVector *r1 = profile->radii->data[i]; + psVector *r2 = profile->radii->data[i+Nsec/2]; + + psVector *f1 = profile->fluxes->data[i]; + psVector *f2 = profile->fluxes->data[i+Nsec/2]; + + psAssert (r1->n == r2->n, "mis-matched vectors"); + psAssert (f1->n == f2->n, "mis-matched vectors"); + + // we have a pair of vectors i, i+Nsec/2; replace them with the finite minimum of the pair + for (int j = 0; j < r1->n; j++) { + + float flux; + + if (!isfinite(f1->data.F32[j]) && !isfinite(f2->data.F32[j])) { + flux = NAN; + goto setflux; + } + + if (!isfinite(f1->data.F32[j])) { + flux = f2->data.F32[j]; + goto setflux; + } + if (!isfinite(f2->data.F32[j])) { + flux = f1->data.F32[j]; + goto setflux; + } + + flux = PS_MIN(f1->data.F32[j], f2->data.F32[j]); + + setflux: + f1->data.F32[j] = flux; + f2->data.F32[j] = flux; + } + } + return true; +} + +float psphotMeanSectorValue (psImage *image, float x, float y, float dL, float dW, float theta) { + + psVector *values = psphotBoxValues (image, x, y, dL, dW, theta); + if (!values) goto escape; + if (!values->n) goto escape; + + psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); + psVectorStats (stats, values, NULL, NULL, 0); + + float value = stats->sampleMedian; + + psFree (stats); + psFree (values); + + return value; + +escape: + psFree(values); + return NAN; +} + +psVector *psphotBoxValues (psImage *image, float x0, float y0, float dL, float dW, float theta) { + + // extract pixels from a series of lines (from -0.5*dW to +0.5*dW) of length dL, + // centered on x0, y0 in parent pixel coordinates (not pixel indicies) + + float xs = x0 - image->col0 - 0.5*dL*cos(theta); + float ys = y0 - image->row0 - 0.5*dL*sin(theta); + + float xe = xs + 0.5*dL*cos(theta); + float ye = ys + 0.5*dL*sin(theta); + + psVector *values = psphotLineValues (image, xs, ys, xe, ye, (int) dW); + return values; +} + +/** + * identify the quadrant and draw the correct line + */ +psVector *psphotLineValues (psImage *image, double x1, double y1, double x2, double y2, int dW) { + + int FlipDirect, FlipCoords; + int X1, Y1, X2, Y2, dX, dY; + + /* rather than draw the line from float positions, we find the closest + integer end-points and draw the line between those pixels */ + + X1 = ROUND(x1); + Y1 = ROUND(y1); + X2 = ROUND(x2); + Y2 = ROUND(y2); + + dX = X2 - X1; + dY = Y2 - Y1; + + FlipCoords = (abs(dX) < abs(dY)); + FlipDirect = FlipCoords ? (y1 > y2) : (x1 > x2); + + psVector *values = NULL; + if (!FlipDirect && !FlipCoords) values = psphotLineValuesBresen (image, X1, Y1, X2, Y2, dW, FALSE); + if ( FlipDirect && !FlipCoords) values = psphotLineValuesBresen (image, X2, Y2, X1, Y1, dW, FALSE); + if (!FlipDirect && FlipCoords) values = psphotLineValuesBresen (image, Y1, X1, Y2, X2, dW, TRUE); + if ( FlipDirect && FlipCoords) values = psphotLineValuesBresen (image, Y2, X2, Y1, X1, dW, TRUE); + + return values; +} + +/** + * use the Bresenham line drawing technique + * integer-only Bresenham line-draw version which is fast + */ +psVector *psphotLineValuesBresen (psImage *image, int X1, int Y1, int X2, int Y2, int dW, int swapcoords) { + + int X, Y, dX, dY; + int e, e2; + + psVector *values = psVectorAllocEmpty(100, PS_TYPE_F32); + + dX = X2 - X1; + dY = Y2 - Y1; + + Y = Y1; + e = 0; + for (X = X1; X <= X2; X++) { + if (X > 0) { + if (swapcoords) { + if (X >= image->numRows) continue; + for (int y = Y - dW; y <= Y + dW; y++) { + if (y < 0) continue; + if (y >= image->numCols) continue; + psVectorAppend(values, image->data.F32[X][y]); + } + } else { + if (X >= image->numCols) continue; + for (int y = Y - dW; y <= Y + dW; y++) { + if (y < 0) continue; + if (y >= image->numRows) continue; + psVectorAppend(values, image->data.F32[y][X]); + } + } + } + e += dY; + e2 = 2 * e; + if (e2 > dX) { + Y++; + e -= dX; + } + if (e2 < -dX) { + Y--; + e += dX; + } + } + return values; +} + Index: branches/eam_branches/20090715/psphot/src/psphotRadiiFromProfiles.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotRadiiFromProfiles.c (revision 25451) +++ branches/eam_branches/20090715/psphot/src/psphotRadiiFromProfiles.c (revision 25452) @@ -1,21 +1,24 @@ # include "psphotInternal.h" -// Given the Petrosian data (radii, fluxes) determine the radius for each profile at the desisred isophote +// Given the Radial Profiles (radii, fluxes) determine the radius for each profile at the desired isophote -bool psphotRadiiFromProfiles (pmSource *source, pmPetrosian *petrosian, float fluxMin, float fluxMax) { +bool psphotRadiiFromProfiles (pmSource *source, float fluxMin, float fluxMax) { - petrosian->isophotalRadii = psVectorAlloc(petrosian->theta->n, PS_TYPE_F32); + pmSourceRadialProfile *profile = source->extpars->profile; - for (int i = 0; i < petrosian->theta->n; i++) { - psVector *radii = petrosian->radii->data[i]; - psVector *fluxes = petrosian->fluxes->data[i]; - float radius = psphotRadiusFromProfile (source, radii, fluxes, fluxMin, fluxMax); + psFree(profile->isophotalRadii); + profile->isophotalRadii = psVectorAlloc(profile->theta->n, PS_TYPE_F32); - // psphotPetrosianVisualProfileByAngle (radii, fluxes, radius); + for (int i = 0; i < profile->theta->n; i++) { + psVector *radii = profile->radii->data[i]; + psVector *fluxes = profile->fluxes->data[i]; + float radius = psphotRadiusFromProfile (source, radii, fluxes, fluxMin, fluxMax); - // warn on NAN? - petrosian->isophotalRadii->data.F32[i] = radius; - } - return true; + // psphotPetrosianVisualProfileByAngle (radii, fluxes, radius); + + // warn on NAN? + profile->isophotalRadii->data.F32[i] = radius; + } + return true; } @@ -79,5 +82,5 @@ // sort the flux by the radius - psphotPetrosianSortPair (radius, flux); + pmSourceRadialProfileSortPair (radius, flux); int nOut = 0; @@ -128,5 +131,5 @@ // XXX is there a macro in psLib that does this interpolation? if (i == 0) { - psLogMsg ("psphot", 3, "bogus radial profile for ..., skipping"); + psTrace ("psphot", 4, "bogus radial profile for source at %f, %f, skipping", source->peak->xf, source->peak->yf); psFree (fluxBinned); psFree (radiusBinned); Index: branches/eam_branches/20090715/psphot/src/psphotReadout.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotReadout.c (revision 25451) +++ branches/eam_branches/20090715/psphot/src/psphotReadout.c (revision 25452) @@ -212,5 +212,6 @@ finish: - psphotPetrosianAnalysis (readout, sources, recipe); + // XXX drop this : test dev version + // psphotPetrosianAnalysis (readout, sources, recipe); // plot positive sources Index: branches/eam_branches/20090715/psphot/src/psphotSourceSize.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotSourceSize.c (revision 25451) +++ branches/eam_branches/20090715/psphot/src/psphotSourceSize.c (revision 25452) @@ -334,4 +334,7 @@ bool isPSF = (fabs(nSigma) < options->nSigmaApResid) && (fabs(Mxx - psfClump->X) < options->nSigmaMoments*psfClump->dX) && (fabs(Myy - psfClump->Y) < options->nSigmaMoments*psfClump->dY); if (isPSF) { + if (Mxx > 3.0) { + fprintf (stderr, "!"); + } Npsf ++; continue; @@ -354,5 +357,5 @@ // XXX allow the Mxx, Myy to be less than psfClump->X,Y (by some nSigma)? - bool isEXT = (nSigma > options->nSigmaApResid) && (Mxx > psfClump->X) && (Myy > psfClump->Y); + bool isEXT = (nSigma > options->nSigmaApResid) || ((Mxx > psfClump->X) && (Myy > psfClump->Y)); if (isEXT) { source->mode |= PM_SOURCE_MODE_EXT_LIMIT; @@ -361,4 +364,5 @@ } + fprintf (stderr, "miss %f,%f : %f %f -- %f\n", source->peak->xf, source->peak->yf, Mxx, Myy, nSigma); Nmiss ++; } Index: branches/eam_branches/20090715/psphot/src/psphotVisual.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotVisual.c (revision 25451) +++ branches/eam_branches/20090715/psphot/src/psphotVisual.c (revision 25452) @@ -2094,9 +2094,8 @@ } -bool psphotVisualShowPetrosian (pmSource *source, pmPetrosian *petrosian) { - - KiiOverlay overlay[2]; - - if (source == NULL) return true; +bool psphotVisualShowPetrosians (psArray *sources) { + + int Noverlay, NOVERLAY; + KiiOverlay *overlay; // if (!pmVisualIsVisual()) return true; @@ -2105,26 +2104,50 @@ if (kapa == -1) return false; - if (kapa == -1) { - fprintf (stderr, "kapa not opened, skipping\n"); - return false; - } - - overlay[0].type = KII_OVERLAY_CIRCLE; - overlay[0].x = source->peak->xf; - overlay[0].y = source->peak->yf; - overlay[0].dx = 2.0*petrosian->petrosianRadius; - overlay[0].dy = 2.0*petrosian->petrosianRadius*petrosian->axes.minor/petrosian->axes.major; - overlay[0].angle = petrosian->axes.theta * PS_DEG_RAD; - overlay[0].text = NULL; - - overlay[1].type = KII_OVERLAY_CIRCLE; - overlay[1].x = source->peak->xf; - overlay[1].y = source->peak->yf; - overlay[1].dx = 4.0*petrosian->petrosianRadius; - overlay[1].dy = 4.0*petrosian->petrosianRadius*petrosian->axes.minor/petrosian->axes.major; - overlay[1].angle = petrosian->axes.theta * PS_DEG_RAD; - overlay[1].text = NULL; - - KiiLoadOverlay (kapa, overlay, 2, "red"); + Noverlay = 0; + NOVERLAY = 100; + ALLOCATE (overlay, KiiOverlay, NOVERLAY); + + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + + if (!source) continue; + if (!source->extpars) continue; + if (!source->extpars->profile) continue; + if (!source->extpars->petrosian_80) continue; + + pmSourceRadialProfile *profile = source->extpars->profile; + pmSourceExtendedFlux *petrosian = source->extpars->petrosian_80; + + overlay[Noverlay].type = KII_OVERLAY_CIRCLE; + overlay[Noverlay].x = source->peak->xf; + overlay[Noverlay].y = source->peak->yf; + overlay[Noverlay].dx = 2.0*petrosian->radius; + overlay[Noverlay].dy = 2.0*petrosian->radius*profile->axes.minor/profile->axes.major; + overlay[Noverlay].angle = profile->axes.theta * PS_DEG_RAD; + overlay[Noverlay].text = NULL; + Noverlay ++; + CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); + + overlay[Noverlay].type = KII_OVERLAY_CIRCLE; + overlay[Noverlay].x = source->peak->xf; + overlay[Noverlay].y = source->peak->yf; + overlay[Noverlay].dx = 4.0*petrosian->radius; + overlay[Noverlay].dy = 4.0*petrosian->radius*profile->axes.minor/profile->axes.major; + overlay[Noverlay].angle = profile->axes.theta * PS_DEG_RAD; + overlay[Noverlay].text = NULL; + Noverlay ++; + CHECK_REALLOCATE (overlay, KiiOverlay, NOVERLAY, Noverlay, 100); + } + + KiiLoadOverlay (kapa, overlay, Noverlay, "red"); + FREE (overlay); + + // pause and wait for user input: + // continue, save (provide name), ?? + char key[10]; + fprintf (stdout, "[c]ontinue? "); + if (!fgets(key, 8, stdin)) { + psWarning("Unable to read option"); + } return true; @@ -2152,37 +2175,37 @@ # if (0) - // *** make a histogram of the source counts in the x and y directions - psHistogram *nX = psHistogramAlloc (graphdata.xmin, graphdata.xmax, 50.0); - psHistogram *nY = psHistogramAlloc (graphdata.ymin, graphdata.ymax, 50.0); - psVectorHistogram (nX, xFaint, NULL, NULL, 0); - psVectorHistogram (nY, yFaint, NULL, NULL, 0); - psVector *dX = psVectorAlloc (nX->nums->n, PS_TYPE_F32); - psVector *vX = psVectorAlloc (nX->nums->n, PS_TYPE_F32); - psVector *dY = psVectorAlloc (nY->nums->n, PS_TYPE_F32); - psVector *vY = psVectorAlloc (nY->nums->n, PS_TYPE_F32); - for (int i = 0; i < nX->nums->n; i++) { - dX->data.F32[i] = nX->nums->data.S32[i]; - vX->data.F32[i] = 0.5*(nX->bounds->data.F32[i] + nX->bounds->data.F32[i+1]); - } - for (int i = 0; i < nY->nums->n; i++) { - dY->data.F32[i] = nY->nums->data.S32[i]; - vY->data.F32[i] = 0.5*(nY->bounds->data.F32[i] + nY->bounds->data.F32[i+1]); - } - - graphdata.color = KapaColorByName ("black"); - graphdata.ptype = 0; - graphdata.size = 0.0; - graphdata.style = 0; - KapaPrepPlot (myKapa, dX->n, &graphdata); - KapaPlotVector (myKapa, dX->n, dX->data.F32, "x"); - KapaPlotVector (myKapa, vX->n, vX->data.F32, "y"); - - psFree (nX); - psFree (dX); - psFree (vX); - - psFree (nY); - psFree (dY); - psFree (vY); +// *** make a histogram of the source counts in the x and y directions +psHistogram *nX = psHistogramAlloc (graphdata.xmin, graphdata.xmax, 50.0); +psHistogram *nY = psHistogramAlloc (graphdata.ymin, graphdata.ymax, 50.0); +psVectorHistogram (nX, xFaint, NULL, NULL, 0); +psVectorHistogram (nY, yFaint, NULL, NULL, 0); +psVector *dX = psVectorAlloc (nX->nums->n, PS_TYPE_F32); +psVector *vX = psVectorAlloc (nX->nums->n, PS_TYPE_F32); +psVector *dY = psVectorAlloc (nY->nums->n, PS_TYPE_F32); +psVector *vY = psVectorAlloc (nY->nums->n, PS_TYPE_F32); +for (int i = 0; i < nX->nums->n; i++) { + dX->data.F32[i] = nX->nums->data.S32[i]; + vX->data.F32[i] = 0.5*(nX->bounds->data.F32[i] + nX->bounds->data.F32[i+1]); +} +for (int i = 0; i < nY->nums->n; i++) { + dY->data.F32[i] = nY->nums->data.S32[i]; + vY->data.F32[i] = 0.5*(nY->bounds->data.F32[i] + nY->bounds->data.F32[i+1]); +} + +graphdata.color = KapaColorByName ("black"); +graphdata.ptype = 0; +graphdata.size = 0.0; +graphdata.style = 0; +KapaPrepPlot (myKapa, dX->n, &graphdata); +KapaPlotVector (myKapa, dX->n, dX->data.F32, "x"); +KapaPlotVector (myKapa, vX->n, vX->data.F32, "y"); + +psFree (nX); +psFree (dX); +psFree (vX); + +psFree (nY); +psFree (dY); +psFree (vY);