Index: /trunk/psModules/src/objects/models/pmModel_QGAUSS.c =================================================================== --- /trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 17395) +++ /trunk/psModules/src/objects/models/pmModel_QGAUSS.c (revision 17396) @@ -237,6 +237,7 @@ if (!isfinite(shape.sxy)) return false; - PAR[PM_PAR_SKY] = moments->Sky; - PAR[PM_PAR_I0] = moments->Peak - moments->Sky; + // XXX turn this off here for now PAR[PM_PAR_SKY] = moments->Sky; + PAR[PM_PAR_SKY] = 0.0; + PAR[PM_PAR_I0] = moments->Peak; PAR[PM_PAR_XPOS] = peak->x; PAR[PM_PAR_YPOS] = peak->y; @@ -457,7 +458,10 @@ status = true; - status &= (dP < 0.5); +// status &= (dP < 0.5); status &= (PAR[PM_PAR_I0] > 0); status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5); + + fprintf (stderr, "QGAUSS status pars: dP: %f, I0: %f, S/N: %f\n", + dP, PAR[PM_PAR_I0], (dPAR[PM_PAR_I0]/PAR[PM_PAR_I0])); return status; Index: /trunk/psModules/src/objects/models/pmModel_SERSIC.c =================================================================== --- /trunk/psModules/src/objects/models/pmModel_SERSIC.c (revision 17395) +++ /trunk/psModules/src/objects/models/pmModel_SERSIC.c (revision 17396) @@ -19,5 +19,5 @@ note that a standard sersic model uses exp(-K*(z^(1/n) - 1). the additional elements (K, - the -1 offset) are absorbed in this model by the normalization, the exponenet, and the + the -1 offset) are absorbed in this model by the normalization, the exponent, and the radial scale. We fit the elements in this form, then re-normalize them on output. *****************************************************************************/ @@ -157,8 +157,8 @@ break; case PM_PAR_SXX: - params_min = 0.25; + params_min = 0.05; break; case PM_PAR_SYY: - params_min = 0.25; + params_min = 0.05; break; case PM_PAR_SXY: @@ -166,5 +166,5 @@ break; case PM_PAR_7: - params_min = 0.1; + params_min = 0.05; break; default: @@ -247,6 +247,7 @@ if (!isfinite(shape.sxy)) return false; - PAR[PM_PAR_SKY] = moments->Sky; - PAR[PM_PAR_I0] = moments->Peak - moments->Sky; + // XXX PAR[PM_PAR_SKY] = moments->Sky; + PAR[PM_PAR_SKY] = 0.0; + PAR[PM_PAR_I0] = moments->Peak; PAR[PM_PAR_XPOS] = peak->x; PAR[PM_PAR_YPOS] = peak->y; @@ -442,7 +443,10 @@ status = true; - status &= (dP < 0.5); +// status &= (dP < 0.5); status &= (PAR[PM_PAR_I0] > 0); status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5); + + fprintf (stderr, "SERSIC status pars: dP: %f, I0: %f, S/N: %f\n", + dP, PAR[PM_PAR_I0], (dPAR[PM_PAR_I0]/PAR[PM_PAR_I0])); return status; Index: /trunk/psModules/src/objects/pmPeaks.c =================================================================== --- /trunk/psModules/src/objects/pmPeaks.c (revision 17395) +++ /trunk/psModules/src/objects/pmPeaks.c (revision 17396) @@ -6,6 +6,6 @@ * @author EAM, IfA: significant modifications. * - * @version $Revision: 1.19 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-01-02 20:38:28 $ + * @version $Revision: 1.20 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-04-08 18:35:38 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -141,5 +141,5 @@ tmp->xf = x; tmp->yf = y; - + tmp->assigned = false; tmp->type = type; Index: /trunk/psModules/src/objects/pmPeaks.h =================================================================== --- /trunk/psModules/src/objects/pmPeaks.h (revision 17395) +++ /trunk/psModules/src/objects/pmPeaks.h (revision 17396) @@ -10,6 +10,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-01-02 20:42:46 $ + * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-04-08 18:35:38 $ * Copyright 2004 Maui High Performance Computing Center, University of Hawaii */ @@ -58,5 +58,6 @@ float flux; ///< level in unsmoothed sci image float SN; ///< S/N implied by detection level - pmPeakType type; ///< Description of peak. + bool assigned; ///< is peak assigned to a source? + pmPeakType type; ///< Description of peak. } pmPeak; Index: /trunk/psModules/src/objects/pmSource.c =================================================================== --- /trunk/psModules/src/objects/pmSource.c (revision 17395) +++ /trunk/psModules/src/objects/pmSource.c (revision 17396) @@ -6,6 +6,6 @@ * @author EAM, IfA: significant modifications. * - * @version $Revision: 1.51 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-03-19 00:51:09 $ + * @version $Revision: 1.52 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-04-08 18:35:38 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -50,5 +50,6 @@ psFree(tmp->modelPSF); psFree(tmp->modelEXT); - psFree(tmp->modelConv); + psFree(tmp->modelFits); + psFree(tmp->extpars); psFree(tmp->blends); psTrace("psModules.objects", 5, "---- end ----\n"); @@ -107,5 +108,5 @@ source->modelPSF = NULL; source->modelEXT = NULL; - source->modelConv = NULL; + source->modelFits = NULL; source->type = PM_SOURCE_TYPE_UNKNOWN; source->mode = PM_SOURCE_MODE_DEFAULT; @@ -947,11 +948,10 @@ return model; -// XXX when should I return the modelConv ?? + // the 'best' extended model is saved in source->modelEXT (may be a pointer to one of + // the elements of source->modelFits) case PM_SOURCE_TYPE_EXTENDED: - model = source->modelConv; - if (!model) { - model = source->modelEXT; - } + model = source->modelEXT; if (!model && source->modelPSF) { + // XXX raise an error or warning here? if (isPSF) { *isPSF = true; @@ -1015,5 +1015,6 @@ if (!strcasecmp (name, "DEFECT" )) return PM_SOURCE_MODE_DEFECT; if (!strcasecmp (name, "SATURATED" )) return PM_SOURCE_MODE_SATURATED; - if (!strcasecmp (name, "CRLIMIT" )) return PM_SOURCE_MODE_CRLIMIT; + if (!strcasecmp (name, "CRLIMIT" )) return PM_SOURCE_MODE_CR_LIMIT; + if (!strcasecmp (name, "EXTLIMIT" )) return PM_SOURCE_MODE_EXT_LIMIT; if (!strcasecmp (name, "SUBTRACTED")) return PM_SOURCE_MODE_SUBTRACTED; return PM_SOURCE_MODE_DEFAULT; @@ -1036,5 +1037,6 @@ case PM_SOURCE_MODE_DEFECT : return psStringCopy ("DEFECT" ); case PM_SOURCE_MODE_SATURATED : return psStringCopy ("SATURATED" ); - case PM_SOURCE_MODE_CRLIMIT : return psStringCopy ("CRLIMIT"); + case PM_SOURCE_MODE_CR_LIMIT : return psStringCopy ("CRLIMIT" ); + case PM_SOURCE_MODE_EXT_LIMIT : return psStringCopy ("EXTLIMIT" ); case PM_SOURCE_MODE_SUBTRACTED : return psStringCopy ("SUBTRACTED"); default: Index: /trunk/psModules/src/objects/pmSource.h =================================================================== --- /trunk/psModules/src/objects/pmSource.h (revision 17395) +++ /trunk/psModules/src/objects/pmSource.h (revision 17396) @@ -3,6 +3,6 @@ * @author EAM, IfA; GLG, MHPCC * - * @version $Revision: 1.23 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-03-05 01:08:08 $ + * @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-04-08 18:35:38 $ * Copyright 2004 Maui High Performance Computing Center, University of Hawaii */ @@ -22,14 +22,11 @@ * source. * - * XXX: The values given below are currently illustrative and will require - * some modification as the source classification code is developed. (TBD) - * */ typedef enum { - PM_SOURCE_TYPE_UNKNOWN, ///< a cosmic-ray - PM_SOURCE_TYPE_DEFECT, ///< a cosmic-ray - PM_SOURCE_TYPE_SATURATED, ///< random saturated pixels - PM_SOURCE_TYPE_STAR, ///< a good-quality star - PM_SOURCE_TYPE_EXTENDED, ///< an extended object (eg, galaxy) + PM_SOURCE_TYPE_UNKNOWN, ///< not yet classified + PM_SOURCE_TYPE_DEFECT, ///< a cosmic-ray + PM_SOURCE_TYPE_SATURATED, ///< random saturated pixels (eg, bleed trails) + PM_SOURCE_TYPE_STAR, ///< a good-quality star (subtracted model is PSF) + PM_SOURCE_TYPE_EXTENDED, ///< an extended object (eg, galaxy) (subtracted model is EXT) } pmSourceType; @@ -48,7 +45,8 @@ PM_SOURCE_MODE_BADPSF = 0x0400, ///< Failed to get good estimate of object's PSF PM_SOURCE_MODE_DEFECT = 0x0800, ///< Source is thought to be a defect - PM_SOURCE_MODE_SATURATED = 0x1000, ///< Source is thought to be saturation - PM_SOURCE_MODE_CRLIMIT = 0x2000, ///< Source has crNsigma above limit - PM_SOURCE_MODE_SUBTRACTED = 0x4000, ///< XXX this flag is actually only used internally (move) + PM_SOURCE_MODE_SATURATED = 0x1000, ///< Source is thought to be saturated pixels (bleed trail) + PM_SOURCE_MODE_CR_LIMIT = 0x2000, ///< Source has crNsigma above limit + PM_SOURCE_MODE_EXT_LIMIT = 0x4000, ///< Source has extNsigma above limit + PM_SOURCE_MODE_SUBTRACTED = 0x8000, ///< XXX this flag is actually only used internally (move) } pmSourceMode; @@ -60,4 +58,7 @@ * * XXX do I have to re-organize this (again!) to allow an arbitrary set of extended model fits?? + * XXX put the Mag and Err inside the pmModel? + * XXX keep the modelEXT or add to the psArray + * * */ @@ -65,5 +66,5 @@ const int id; ///< Unique ID for object int seq; ///< ID for output (generated on write) - pmPeak *peak; ///< Description of peak pixel. + pmPeak *peak; ///< Description of peak pixel. psImage *pixels; ///< Rectangular region including object pixels. psImage *weight; ///< Image variance. @@ -72,21 +73,21 @@ psImage *modelFlux; ///< cached copy of the best model for this source psImage *psfFlux; ///< cached copy of the psf model for this source - pmMoments *moments; ///< Basic moments measure for the object. + pmMoments *moments; ///< Basic moments measured for the object. pmModel *modelPSF; ///< PSF Model fit (parameters and type) - pmModel *modelEXT; ///< EXT (floating) Model fit (parameters and type). - pmModel *modelConv; ///< PSF-Convolved Model fit (parameters and type). + pmModel *modelEXT; ///< EXT Model fit used for subtraction (parameters and type) + psArray *modelFits; ///< collection of extended source models (best == modelEXT) pmSourceType type; ///< Best identification of object. - pmSourceMode mode; ///< Best identification of object. - psArray *blends; - float psfMag; ///< calculated from flux in modelPsf + pmSourceMode mode; ///< analysis flags set for object. + psArray *blends; ///< collection of sources thought to be confused with object + float psfMag; ///< calculated from flux in modelPSF float extMag; ///< calculated from flux in modelEXT float errMag; ///< error in psfMag OR extMag (depending on type) float apMag; ///< apMag corresponding to psfMag or extMag (depending on type) float pixWeight; ///< model-weighted coverage of valid pixels - float psfChisq; ///< probability of PSF + float psfChisq; ///< probability of PSF float crNsigma; ///< Nsigma deviation from PSF to CR float extNsigma; ///< Nsigma deviation from PSF to EXT + float sky, skyErr; ///< The sky and its error at the center of the object psRegion region; ///< area on image covered by selected pixels - float sky, skyErr; ///< The sky and its error at the center of the object pmSourceExtendedPars *extpars; ///< extended source parameters }; Index: /trunk/psModules/src/objects/pmSourceIO.c =================================================================== --- /trunk/psModules/src/objects/pmSourceIO.c (revision 17395) +++ /trunk/psModules/src/objects/pmSourceIO.c (revision 17396) @@ -3,6 +3,6 @@ * @author EAM, IfA * - * @version $Revision: 1.55 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-03-17 22:04:27 $ + * @version $Revision: 1.56 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-04-08 18:35:38 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -327,6 +327,6 @@ // if this is not TRUE, the output files only contain the psf measurements. - bool XSRC_OUTPUT = psMetadataLookupBool(&status, recipe, "SAVE.XSRC"); - bool XFIT_OUTPUT = psMetadataLookupBool(&status, recipe, "SAVE.XFIT"); + bool XSRC_OUTPUT = psMetadataLookupBool(&status, recipe, "EXTENDED_SOURCE_ANALYSIS"); + bool XFIT_OUTPUT = psMetadataLookupBool(&status, recipe, "EXTENDED_SOURCE_FITS"); // define the EXTNAME values for the different data segments: @@ -445,5 +445,5 @@ if (xsrcname) { if (!strcmp (exttype, "PS1_DEV_1")) { - status = pmSourcesWrite_PS1_DEV_1_XSRC (file->fits, sources, xsrcname); + status = pmSourcesWrite_PS1_DEV_1_XSRC (file->fits, sources, xsrcname, recipe); } } @@ -457,4 +457,6 @@ psFree (headname); psFree (dataname); + psFree (xsrcname); + psFree (xfitname); psFree (outhead); return false; @@ -466,4 +468,6 @@ psFree (headname); psFree (dataname); + psFree (xsrcname); + psFree (xfitname); psFree (outhead); break; Index: /trunk/psModules/src/objects/pmSourceIO.h =================================================================== --- /trunk/psModules/src/objects/pmSourceIO.h (revision 17395) +++ /trunk/psModules/src/objects/pmSourceIO.h (revision 17396) @@ -4,6 +4,6 @@ * @author EAM, IfA; GLG, MHPCC * - * @version $Revision: 1.14 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-03-05 01:08:08 $ + * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-04-08 18:35:38 $ * Copyright 2004 Maui High Performance Computing Center, University of Hawaii * @@ -27,5 +27,5 @@ bool pmSourcesWrite_PS1_DEV_0 (psFits *fits, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname); bool pmSourcesWrite_PS1_DEV_1 (psFits *fits, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, char *xsrcname); -bool pmSourcesWrite_PS1_DEV_1_XSRC (psFits *fits, psArray *sources, char *extname); +bool pmSourcesWrite_PS1_DEV_1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe); bool pmSourcesWrite_PS1_DEV_1_XFIT (psFits *fits, psArray *sources, char *extname); Index: /trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c =================================================================== --- /trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c (revision 17395) +++ /trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c (revision 17396) @@ -3,6 +3,6 @@ * @author EAM, IfA * - * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-03-05 01:08:08 $ + * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-04-08 18:35:38 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -251,12 +251,11 @@ } -bool pmSourcesWrite_PS1_DEV_1_XSRC (psFits *fits, psArray *sources, char *extname) +bool pmSourcesWrite_PS1_DEV_1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe) { + bool status; psArray *table; psMetadata *row; - int i; psF32 *PAR, *dPAR; - psEllipseAxes axes; psF32 xPos, yPos; psF32 xErr, yErr; @@ -273,6 +272,16 @@ table = psArrayAllocEmpty (sources->n); + // 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"); + + psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); + psVector *radialBinsUpper = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); + assert (radialBinsLower->n == radialBinsUpper->n); + // we write out all sources, regardless of quality. the source flags tell us the state - for (i = 0; i < sources->n; i++) { + for (int i = 0; i < sources->n; i++) { // skip source if it is not a ext sourc // XXX we have two places that extended source parameters are measured: @@ -283,13 +292,10 @@ pmSource *source = sources->data[i]; - // choose the convolved EXT model, if available, otherwise the simple one - pmModel *model = source->modelConv; - if (model == NULL) { - model = source->modelEXT; - } - if (model == NULL) continue; - - // XXX do we need to do this? + // skip sources without measurements if (source->extpars == NULL) continue; + + // we require a PSF model fit (ignore the real crud) + pmModel *model = source->modelPSF; + if (model == NULL) continue; // XXX I need to split the extended models from the extended aperture measurements @@ -301,8 +307,6 @@ yErr = dPAR[PM_PAR_YPOS]; - // XXX for the aperture values, I probably can / should just refer to the psf position and not write any of this - axes = pmPSF_ModelToAxes (PAR, 20.0); - row = psMetadataAlloc (); + // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); @@ -311,19 +315,8 @@ psMetadataAdd (row, PS_LIST_TAIL, "X_EXT_SIG", PS_DATA_F32, "Sigma in EXT x coordinate", xErr); psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG", PS_DATA_F32, "Sigma in EXT y coordinate", yErr); - psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG", PS_DATA_F32, "EXT fit instrumental magnitude", source->extMag); - psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->errMag); - - // XXX these should be major and minor, not 'x' and 'y' - psMetadataAdd (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", PS_DATA_F32, "EXT width in x coordinate", axes.major); - psMetadataAdd (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", PS_DATA_F32, "EXT width in y coordinate", axes.minor); - psMetadataAdd (row, PS_LIST_TAIL, "EXT_THETA", PS_DATA_F32, "EXT orientation angle", axes.theta); - - // XXX at the moment, this will be a fixed sized table, but perhaps have multiple output formats depending on what is measured? - - // other values that I need to report: - // R, Mag Petrosian, ..... - if (source->extpars) { - - // Petrosian measurements + + // Petrosian measurements + // XXX insert header data: petrosian ref radius, flux ratio + if (doPetrosian) { pmSourcePetrosianValues *petrosian = source->extpars->petrosian; if (petrosian) { @@ -332,7 +325,14 @@ psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius", petrosian->rad); psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error", petrosian->radErr); + } else { + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error", NAN); } - - // Kron measurements + } + + // Kron measurements + if (doKron) { pmSourceKronValues *kron = source->extpars->kron; if (kron) { @@ -341,8 +341,15 @@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS", PS_DATA_F32, "Kron Radius", kron->rad); psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error", kron->radErr); + } else { + psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG", PS_DATA_F32, "Kron Magnitude", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR", PS_DATA_F32, "Kron Magnitude Error", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS", PS_DATA_F32, "Kron Radius", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error", NAN); } - - // Isophot measurements - // XXX insert header data: isophotal level + } + + // Isophot measurements + // XXX insert header data: isophotal level + if (doIsophotal) { pmSourceIsophotalValues *isophot = source->extpars->isophot; if (isophot) { @@ -351,6 +358,14 @@ psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS", PS_DATA_F32, "Isophot Radius", isophot->rad); psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error", isophot->radErr); + } else { + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG", PS_DATA_F32, "Isophot Magnitude", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR", PS_DATA_F32, "Isophot Magnitude Error", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS", PS_DATA_F32, "Isophot Radius", NAN); + psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error", NAN); } - + } + + // Flux Annuli + if (doAnnuli) { pmSourceAnnuli *annuli = source->extpars->annuli; if (annuli) { @@ -367,24 +382,37 @@ sprintf (name, "FLUX_VAR_R_%02d", j); psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", fluxVar->data.F32[j]); - } + } + } else { + for (int j = 0; j < radialBinsLower->n; j++) { + char name[32]; + sprintf (name, "FLUX_VAL_R_%02d", j); + psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", NAN); + sprintf (name, "FLUX_ERR_R_%02d", j); + psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", NAN); + sprintf (name, "FLUX_VAR_R_%02d", j); + psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", NAN); + } } } - psArrayAdd (table, 100, row); - psFree (row); + psArrayAdd (table, 100, row); + psFree (row); } if (table->n == 0) { - psFitsWriteBlank (fits, outhead, extname); - psFree (table); - return true; + psFitsWriteBlank (fits, outhead, extname); + psFree (outhead); + psFree (table); + return true; } psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); if (!psFitsWriteTable (fits, outhead, table, extname)) { - psError(PS_ERR_IO, false, "writing ext data %s\n", extname); - psFree(table); - return false; - } + psError(PS_ERR_IO, false, "writing ext data %s\n", extname); + psFree (outhead); + psFree(table); + return false; + } + psFree (outhead); psFree (table); @@ -397,9 +425,9 @@ psArray *table; psMetadata *row; - int i; psF32 *PAR, *dPAR; psEllipseAxes axes; psF32 xPos, yPos; psF32 xErr, yErr; + char name[64]; // create a header to hold the output data @@ -412,72 +440,105 @@ sources = psArraySort (sources, pmSourceSortBySN); + // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams + int nParamMax = 0; + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + if (source->modelFits == NULL) continue; + for (int j = 0; j < source->modelFits->n; j++) { + pmModel *model = source->modelFits->data[j]; + assert (model); + nParamMax = PS_MAX (nParamMax, model->params->n); + } + } + table = psArrayAllocEmpty (sources->n); // we write out all sources, regardless of quality. the source flags tell us the state - for (i = 0; i < sources->n; i++) { - // skip source if it is not a ext sourc - // XXX we have two places that extended source parameters are measured: - // psphotExtendedSources, which measures the aperture-like parameters and (potentially) the psf-convolved extended source models, - // psphotFitEXT, which does the simple extended source model fit (not psf-convolved) - // should we require both? + for (int i = 0; i < sources->n; i++) { pmSource *source = sources->data[i]; - // XXX need to have an array of model fits; loop over the array and write out all - // which we calculated - - // choose the convolved EXT model, if available, otherwise the simple one - // XXX should not need to choose: write both out - pmModel *model = source->modelConv; - if (model == NULL) { - model = source->modelEXT; + // XXX if no model fits are saved, write out modelEXT? + if (source->modelFits == NULL) continue; + + // We have multiple sources : need to flag the one used to subtract the light (the 'best' model) + for (int j = 0; j < source->modelFits->n; j++) { + + // choose the convolved EXT model, if available, otherwise the simple one + pmModel *model = source->modelFits->data[j]; + assert (model); + + PAR = model->params->data.F32; + dPAR = model->dparams->data.F32; + xPos = PAR[PM_PAR_XPOS]; + yPos = PAR[PM_PAR_YPOS]; + xErr = dPAR[PM_PAR_XPOS]; + yErr = dPAR[PM_PAR_YPOS]; + + axes = pmPSF_ModelToAxes (PAR, 20.0); + + row = psMetadataAlloc (); + + // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) + psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET", 0, "IPP detection identifier index", source->seq); + psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT", 0, "EXT model x coordinate", xPos); + psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT", 0, "EXT model y coordinate", yPos); + psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT_SIG", 0, "Sigma in EXT x coordinate", xErr); + psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT_SIG", 0, "Sigma in EXT y coordinate", yErr); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG", 0, "EXT fit instrumental magnitude", model->mag); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", 0, "Sigma of PSF instrumental magnitude", model->magErr); + + psMetadataAddF32 (row, PS_LIST_TAIL, "NPARAMS", 0, "number of model parameters", model->params->n); + psMetadataAddStr (row, PS_LIST_TAIL, "MODEL_TYPE", 0, "name of model", pmModelClassGetName (model->type)); + + // XXX these should be major and minor, not 'x' and 'y' + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", 0, "EXT width in x coordinate", axes.major); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", 0, "EXT width in y coordinate", axes.minor); + psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA", 0, "EXT orientation angle", axes.theta); + + // write out the other generic parameters + for (int k = 0; k < nParamMax; k++) { + if (k == PM_PAR_I0) continue; + if (k == PM_PAR_SKY) continue; + if (k == PM_PAR_XPOS) continue; + if (k == PM_PAR_YPOS) continue; + if (k == PM_PAR_SXX) continue; + if (k == PM_PAR_SXY) continue; + if (k == PM_PAR_SYY) continue; + + snprintf (name, 64, "EXT_PAR_%02d", k); + + if (k < model->params->n) { + psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "", model->params->data.F32[k]); + } else { + psMetadataAddF32 (row, PS_LIST_TAIL, name, PS_DATA_F32, "", NAN); + } + } + + // XXX other parameters which may be set. + // XXX flag / value to define the model + // XXX write out the model type, fit status flags + + psArrayAdd (table, 100, row); + psFree (row); } - if (model == NULL) continue; - - PAR = model->params->data.F32; - dPAR = model->dparams->data.F32; - xPos = PAR[PM_PAR_XPOS]; - yPos = PAR[PM_PAR_YPOS]; - xErr = dPAR[PM_PAR_XPOS]; - yErr = dPAR[PM_PAR_YPOS]; - - axes = pmPSF_ModelToAxes (PAR, 20.0); - - row = psMetadataAlloc (); - // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) - psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); - psMetadataAdd (row, PS_LIST_TAIL, "X_EXT", PS_DATA_F32, "EXT model x coordinate", xPos); - psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT", PS_DATA_F32, "EXT model y coordinate", yPos); - psMetadataAdd (row, PS_LIST_TAIL, "X_EXT_SIG", PS_DATA_F32, "Sigma in EXT x coordinate", xErr); - psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG", PS_DATA_F32, "Sigma in EXT y coordinate", yErr); - psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG", PS_DATA_F32, "EXT fit instrumental magnitude", source->extMag); - psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->errMag); - - // XXX these should be major and minor, not 'x' and 'y' - psMetadataAdd (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", PS_DATA_F32, "EXT width in x coordinate", axes.major); - psMetadataAdd (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", PS_DATA_F32, "EXT width in y coordinate", axes.minor); - psMetadataAdd (row, PS_LIST_TAIL, "EXT_THETA", PS_DATA_F32, "EXT orientation angle", axes.theta); - - // XXX other parameters which may be set. - // XXX flag / value to define the model - - psArrayAdd (table, 100, row); - psFree (row); } if (table->n == 0) { - psFitsWriteBlank (fits, outhead, extname); - psFree (table); - return true; + psFitsWriteBlank (fits, outhead, extname); + psFree (outhead); + psFree (table); + return true; } psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); if (!psFitsWriteTable (fits, outhead, table, extname)) { - psError(PS_ERR_IO, false, "writing ext data %s\n", extname); - psFree(table); - return false; - } + psError(PS_ERR_IO, false, "writing ext data %s\n", extname); + psFree (outhead); + psFree(table); + return false; + } + psFree (outhead); psFree (table); - return true; } Index: /trunk/psModules/src/objects/pmSourcePhotometry.c =================================================================== --- /trunk/psModules/src/objects/pmSourcePhotometry.c (revision 17395) +++ /trunk/psModules/src/objects/pmSourcePhotometry.c (revision 17396) @@ -3,6 +3,6 @@ * @author EAM, IfA; GLG, MHPCC * - * @version $Revision: 1.40 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-04-02 22:40:36 $ + * @version $Revision: 1.41 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-04-08 18:35:38 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -119,6 +119,19 @@ } - // measure EXT model photometry (do both modelEXT and modelConv or just the one?) - status = pmSourcePhotometryModel (&source->extMag, source->modelEXT); + // if we have a collection of model fits, one of them is a pointer to modelEXT? + // XXX not guaranteed + if (source->modelFits) { + for (int i = 0; i < source->modelFits->n; i++) { + pmModel *model = source->modelFits->data[i]; + status = pmSourcePhotometryModel (&model->mag, model); + } + if (source->modelEXT) { + source->extMag = source->modelEXT->mag; + } + } else { + if (source->modelEXT) { + status = pmSourcePhotometryModel (&source->extMag, source->modelEXT); + } + } // for PSFs, correct both apMag and psfMag to same system, consistent with infinite flux star in aperture RADIUS Index: /trunk/psphot/doc/notes.txt =================================================================== --- /trunk/psphot/doc/notes.txt (revision 17395) +++ /trunk/psphot/doc/notes.txt (revision 17396) @@ -1,2 +1,27 @@ + +2008.01.31 + + Most of the fixes listed below have been finished except for: + + * careful handling of the r ~ 0 fluxes / effective positions (for eg sersic) + * the EXT output parameters have not been tested + + In addition, some further work to be done: + + * multiple image detection + * multiple image fitting (including projection) + * optimization + * multi-threading + + * possibly try the 'active deblending' concept that Ken from WISE is using. + * instead of fitting a double star with the positions floating at a + slightly random location, search for significant additions within + a grid about the central peak (say ~ 1sigma / 0.2 pix?). + (note that this is a linear fit) + + * for the blend groups, freeze the positions for the fainter sources + (consistent with the non-linear S/N threshold) + + 2007.08.17 Index: /trunk/psphot/src/Makefile.am =================================================================== --- /trunk/psphot/src/Makefile.am (revision 17395) +++ /trunk/psphot/src/Makefile.am (revision 17396) @@ -18,58 +18,59 @@ libpsphot_la_SOURCES = \ - psphotErrorCodes.c \ - pmFootprint.c \ - psphotVersion.c \ - psphotModelGroupInit.c \ - psphotMaskReadout.c \ - psphotDefineFiles.c \ - psphotReadout.c \ - psphotModelBackground.c \ - psphotSubtractBackground.c \ - psphotFindDetections.c \ - psphotFindPeaks.c \ - psphotSourceStats.c \ - psphotRoughClass.c \ - psphotBasicDeblend.c \ - psphotChoosePSF.c \ - psphotGuessModels.c \ - psphotFitSourcesLinear.c \ - psphotBlendFit.c \ - psphotReplaceUnfit.c \ - psphotApResid.c \ - psphotMagnitudes.c \ - psphotSkyReplace.c \ - psphotEvalPSF.c \ - psphotEvalFLT.c \ - psphotSourceFits.c \ - psphotRadiusChecks.c \ - psphotOutput.c \ - psphotGrowthCurve.c \ - psphotFakeSources.c \ - psphotModelWithPSF.c \ - psphotExtendedSources.c \ - psphotRadialProfile.c \ - psphotPetrosian.c \ - psphotIsophotal.c \ - psphotAnnuli.c \ - psphotKron.c \ - psphotKernelFromPSF.c \ - psphotPSFConvModel.c \ - psphotModelTest.c \ - psphotFitSet.c \ - psphotSourceFreePixels.c \ - psphotSummaryPlots.c \ - psphotMergeSources.c \ - psphotLoadPSF.c \ - psphotReadoutCleanup.c \ - psphotSourcePlots.c \ - psphotRadialPlot.c \ - psphotDeblendSatstars.c \ - psphotMosaicSubimage.c \ - psphotMakeResiduals.c \ - psphotSourceSize.c \ - psphotDiagnosticPlots.c \ - psphotMakeFluxScale.c \ - psphotCheckStarDistribution.c \ + psphotErrorCodes.c \ + pmFootprint.c \ + psphotVersion.c \ + psphotModelGroupInit.c \ + psphotMaskReadout.c \ + psphotDefineFiles.c \ + psphotReadout.c \ + psphotModelBackground.c \ + psphotSubtractBackground.c \ + psphotFindDetections.c \ + psphotFindPeaks.c \ + psphotSourceStats.c \ + psphotRoughClass.c \ + psphotBasicDeblend.c \ + psphotChoosePSF.c \ + psphotGuessModels.c \ + psphotFitSourcesLinear.c \ + psphotBlendFit.c \ + psphotReplaceUnfit.c \ + psphotApResid.c \ + psphotMagnitudes.c \ + psphotSkyReplace.c \ + psphotEvalPSF.c \ + psphotEvalFLT.c \ + psphotSourceFits.c \ + psphotRadiusChecks.c \ + psphotOutput.c \ + psphotGrowthCurve.c \ + psphotFakeSources.c \ + psphotModelWithPSF.c \ + psphotExtendedSourceAnalysis.c \ + psphotExtendedSourceFits.c \ + psphotRadialProfile.c \ + psphotPetrosian.c \ + psphotIsophotal.c \ + psphotAnnuli.c \ + psphotKron.c \ + psphotKernelFromPSF.c \ + psphotPSFConvModel.c \ + psphotModelTest.c \ + psphotFitSet.c \ + psphotSourceFreePixels.c \ + psphotSummaryPlots.c \ + psphotMergeSources.c \ + psphotLoadPSF.c \ + psphotReadoutCleanup.c \ + psphotSourcePlots.c \ + psphotRadialPlot.c \ + psphotDeblendSatstars.c \ + psphotMosaicSubimage.c \ + psphotMakeResiduals.c \ + psphotSourceSize.c \ + psphotDiagnosticPlots.c \ + psphotMakeFluxScale.c \ + psphotCheckStarDistribution.c \ psphotAddNoise.c Index: /trunk/psphot/src/psphot.h =================================================================== --- /trunk/psphot/src/psphot.h (revision 17395) +++ /trunk/psphot/src/psphot.h (revision 17396) @@ -45,5 +45,6 @@ bool psphotMagnitudes (psArray *sources, psMetadata *recipe, pmPSF *psf, pmReadout *background); bool psphotSkyReplace (pmConfig *config, const pmFPAview *view); -bool psphotExtendedSources (pmReadout *readout, psArray *sources, psMetadata *recipe); +bool psphotExtendedSourceAnalysis (pmReadout *readout, psArray *sources, psMetadata *recipe); +bool psphotExtendedSourceFits (pmReadout *readout, psArray *sources, psMetadata *recipe); // used by psphotFindDetections @@ -87,5 +88,5 @@ bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf, psMaskType maskVal); bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf, psMaskType maskVal); -pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, psMaskType maskVal); +pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, pmModelType modelType, psMaskType maskVal); psArray *psphotFitDBL (pmReadout *readout, pmSource *source, psMaskType maskVal); @@ -117,5 +118,6 @@ bool psphotMakeResiduals (psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal); -bool psphotPSFConvModel (pmSource *source, psMetadata *recipe, psMaskType maskVal); +pmModel *psphotPSFConvModel (pmReadout *readout, pmSource *source, pmModelType modelType, psMaskType maskVal, int psfSize); + psKernel *psphotKernelFromPSF (pmSource *source, int nPix); Index: /trunk/psphot/src/psphotAnnuli.c =================================================================== --- /trunk/psphot/src/psphotAnnuli.c (revision 17395) +++ /trunk/psphot/src/psphotAnnuli.c (revision 17396) @@ -1,3 +1,3 @@ -# include "psphot.h" +# include "psphotInternal.h" bool psphotAnnuli (pmSource *source, psMetadata *recipe, psMaskType maskVal) { @@ -55,4 +55,6 @@ source->extpars->annuli->fluxVar = fluxSquare; + psFree (pixelCount); + return true; } Index: /trunk/psphot/src/psphotBlendFit.c =================================================================== --- /trunk/psphot/src/psphotBlendFit.c (revision 17395) +++ /trunk/psphot/src/psphotBlendFit.c (revision 17396) @@ -40,4 +40,5 @@ // skip non-astronomical objects (very likely defects) if (source->mode & PM_SOURCE_MODE_BLEND) continue; + if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue; if (source->type == PM_SOURCE_TYPE_DEFECT) continue; if (source->type == PM_SOURCE_TYPE_SATURATED) continue; @@ -79,14 +80,19 @@ // XXX re-consider conditions under which the source has EXT fit: // I should try EXT if the source size measurement says it is large - if (psphotFitBlend (readout, source, psf, maskVal)) { - psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->moments->x, source->moments->y); - Npsf ++; - continue; - } - if (psphotFitBlob (readout, source, sources, psf, maskVal)) { - psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->moments->x, source->moments->y); - Next ++; - continue; - } + if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { + if (psphotFitBlob (readout, source, sources, psf, maskVal)) { + source->type = PM_SOURCE_TYPE_EXTENDED; + psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->moments->x, source->moments->y); + Next ++; + continue; + } + } else { + if (psphotFitBlend (readout, source, psf, maskVal)) { + source->type = PM_SOURCE_TYPE_STAR; + psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->moments->x, source->moments->y); + Npsf ++; + continue; + } + } psTrace ("psphot", 5, "source at %7.1f, %7.1f failed", source->moments->x, source->moments->y); Index: /trunk/psphot/src/psphotEvalFLT.c =================================================================== --- /trunk/psphot/src/psphotEvalFLT.c (revision 17395) +++ /trunk/psphot/src/psphotEvalFLT.c (revision 17396) @@ -8,4 +8,5 @@ if (model == NULL) { source->mode &= ~PM_SOURCE_MODE_FITTED; + psTrace ("psphot", 5, "no model fitted?\n"); return false; } @@ -14,6 +15,8 @@ if (!(model->flags & PM_MODEL_STATUS_FITTED)) { source->mode &= ~PM_SOURCE_MODE_FITTED; + psTrace ("psphot", 5, "no model fitted?\n"); return false; } + // did the model fit fail for one or another reason? if (model->flags & (PM_MODEL_STATUS_BADARGS | @@ -22,4 +25,15 @@ source->mode |= PM_SOURCE_MODE_FAIL; psLogMsg ("psphot", 5, "EXT fail fit\n"); + psTrace ("psphot", 5, "EXT fail fit\n"); + + if (model->flags & PM_MODEL_STATUS_OFFIMAGE) { + psTrace ("psphot", 5, "off image\n"); + } + if (model->flags & PM_MODEL_STATUS_BADARGS) { + psTrace ("psphot", 5, "bad args\n"); + } + if (model->flags & PM_MODEL_STATUS_NONCONVERGE) { + psTrace ("psphot", 5, "non converge\n"); + } return false; } @@ -35,4 +49,5 @@ if (model->params->data.F32[PM_PAR_I0] <= 0.02) { source->mode |= PM_SOURCE_MODE_FAIL; + psTrace ("psphot", 5, "model central intensity ~ zero\n"); return false; } @@ -43,4 +58,6 @@ // poor-quality fit; only keep if nothing else works... psLogMsg ("psphot", 5, "EXT poor fit\n"); + psTrace ("psphot", 5, "EXT poor fit\n"); + source->mode |= PM_SOURCE_MODE_POOR; return false; Index: /trunk/psphot/src/psphotExtendedSourceAnalysis.c =================================================================== --- /trunk/psphot/src/psphotExtendedSourceAnalysis.c (revision 17396) +++ /trunk/psphot/src/psphotExtendedSourceAnalysis.c (revision 17396) @@ -0,0 +1,131 @@ +# include "psphotInternal.h" + +// aperture-like measurements for extended sources +bool psphotExtendedSourceAnalysis (pmReadout *readout, psArray *sources, psMetadata *recipe) { + + bool status; + int Next = 0; + int Npetro = 0; + int Nisophot = 0; + int Nannuli = 0; + int Nkron = 0; + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (maskVal); + + // perform full non-linear fits / extended source analysis? + if (!psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANALYSIS")) { + psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source measurements\n"); + return true; + } + + // 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"); + + // S/N limit to perform full non-linear fits + float SN_LIM = psMetadataLookupF32 (&status, recipe, "EXTENDED_SOURCE_SN_LIM"); + + // 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"); + + // source analysis is done in S/N order (brightest first) + sources = psArraySort (sources, pmSourceSortBySN); + + // XXX some init functions for the extended source recipe options? + + // choose the sources of interest + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + + // skip PSF-like and non-astronomical objects + if (source->type == PM_SOURCE_TYPE_STAR) continue; + if (source->type == PM_SOURCE_TYPE_DEFECT) continue; + if (source->type == PM_SOURCE_TYPE_SATURATED) continue; + + // limit selection to some SN limit + assert (source->peak); // how can a source not have a peak? + if (source->peak->SN < SN_LIM) continue; + + // limit selection by analysis region + if (source->peak->x < AnalysisRegion.x0) continue; + if (source->peak->y < AnalysisRegion.y0) continue; + if (source->peak->x > AnalysisRegion.x1) continue; + if (source->peak->y > AnalysisRegion.y1) continue; + + // replace object in image + if (source->mode & PM_SOURCE_MODE_SUBTRACTED) { + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + } + Next ++; + + // if we request any of these measurements, we require the radial profile + if (doPetrosian || doIsophotal || doAnnuli || doKron) { + if (!psphotRadialProfile (source, recipe, maskVal)) { + // all measurements below require the radial profile; skip them all + // re-subtract the object, leave local sky + psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->x, source->moments->y); + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + source->mode |= PM_SOURCE_MODE_SUBTRACTED; + continue; + } + } + + // Isophotal Mags + if (doIsophotal) { + if (!psphotIsophotal (source, recipe, maskVal)) { + psTrace ("psphot", 5, "failed to measure isophotal mags for source at %7.1f, %7.1f", source->moments->x, source->moments->y); + } else { + psTrace ("psphot", 5, "measured isophotal mags for source at %7.1f, %7.1f", source->moments->x, source->moments->y); + Nisophot ++; + } + } + + // 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->x, source->moments->y); + } else { + psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->x, source->moments->y); + Npetro ++; + } + } + + // Kron Mags + if (doKron) { + if (!psphotKron (source, recipe, maskVal)) { + psTrace ("psphot", 5, "failed to measure kron mags for source at %7.1f, %7.1f", source->moments->x, source->moments->y); + } else { + psTrace ("psphot", 5, "measure kron mags for source at %7.1f, %7.1f", source->moments->x, source->moments->y); + Nkron ++; + } + } + + // 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->x, source->moments->y); + Nannuli ++; + } + + // re-subtract the object, leave local sky + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + source->mode |= PM_SOURCE_MODE_SUBTRACTED; + } + + psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot"), Next); + psLogMsg ("psphot", PS_LOG_INFO, " %d petrosian\n", Npetro); + psLogMsg ("psphot", PS_LOG_INFO, " %d isophotal\n", Nisophot); + psLogMsg ("psphot", PS_LOG_INFO, " %d annuli\n", Nannuli); + psLogMsg ("psphot", PS_LOG_INFO, " %d kron\n", Nkron); + return true; +} Index: /trunk/psphot/src/psphotExtendedSourceFits.c =================================================================== --- /trunk/psphot/src/psphotExtendedSourceFits.c (revision 17396) +++ /trunk/psphot/src/psphotExtendedSourceFits.c (revision 17396) @@ -0,0 +1,270 @@ +# include "psphotInternal.h" + +// non-linear model fitting for extended sources +bool psphotExtendedSourceFits (pmReadout *readout, psArray *sources, psMetadata *recipe) { + + bool status; + int Next = 0; + int Nconvolve = 0; + int NconvolvePass = 0; + int Nplain = 0; + int NplainPass = 0; + bool savePics = false; + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (maskVal); + + // perform full extended source non-linear fits? + if (!psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_FITS")) { + psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source measurements\n"); + return true; + } + + // select the collection of desired models + psMetadata *models = psMetadataLookupMetadata (&status, recipe, "EXTENDED_SOURCE_MODELS"); + if (!status) { + psWarning ("extended source model fits requested but model model is missing (EXTENDED_SOURCE_MODELS)\n"); + return true; + } + if (models->list->n == 0) { + psWarning ("extended source model fits requested but no models are specified\n"); + return true; + } + + // validate the model entries + psMetadataIterator *iter = psMetadataIteratorAlloc (models, PS_LIST_HEAD, NULL); + psMetadataItem *item = NULL; + while ((item = psMetadataGetAndIncrement (iter)) != NULL) { + + if (item->type != PS_DATA_METADATA) { + psAbort ("Invalid type for EXTENDED_SOURCE_MODEL entry %s, not a metadata folder", item->name); + // XXX we could cull the bad entries or build a validated model folder + } + + psMetadata *model = (psMetadata *) item->data.md; + + // check on the model type + char *modelName = psMetadataLookupStr (&status, model, "MODEL"); + int modelType = pmModelClassGetType (modelName); + if (modelType < 0) { + psAbort ("Unknown model class for EXTENDED_SOURCE_MODEL entry %s: %s", item->name, modelName); + } + psMetadataAddS32 (model, PS_LIST_TAIL, "MODEL_TYPE", PS_META_REPLACE, "", modelType); + + // check on the SNLIM, set a float value + char *SNword = psMetadataLookupStr (&status, model, "SNLIM"); + if (!status) { + psAbort("SNLIM not defined for extended source model %s\n", item->name); + } + float SNlim = atof (SNword); + psMetadataAddF32 (model, PS_LIST_TAIL, "SNLIM_VALUE", PS_META_REPLACE, "", SNlim); + + // check on the PSF-Convolution status + char *convolvedWord = psMetadataLookupStr (&status, model, "PSF_CONVOLVED"); + if (!status || (strcasecmp (convolvedWord, "true") && strcasecmp (convolvedWord, "false"))) { + psAbort ("PSF_CONVOLVED entry invalid or missing for EXTENDED_SOURCE_MODEL entry %s", item->name); + } + bool convolved = !strcasecmp (convolvedWord, "true"); + psMetadataAddBool (model, PS_LIST_TAIL, "PSF_CONVOLVED_VALUE", PS_META_REPLACE, "", convolved); + } + psFree (iter); + + // 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"); + + // what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box) + int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE"); + assert (status); + + // source analysis is done in S/N order (brightest first) + sources = psArraySort (sources, pmSourceSortBySN); + + // choose the sources of interest + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + + // skip PSF-like and non-astronomical objects + if (source->type == PM_SOURCE_TYPE_STAR) continue; + if (source->type == PM_SOURCE_TYPE_DEFECT) continue; + if (source->type == PM_SOURCE_TYPE_SATURATED) continue; + + // XXX this should use peak? + if (source->peak->x < AnalysisRegion.x0) continue; + if (source->peak->y < AnalysisRegion.y0) continue; + if (source->peak->x > AnalysisRegion.x1) continue; + if (source->peak->y > AnalysisRegion.y1) continue; + + // if model is NULL, we don't have a starting guess + // XXX use the parameters guessed from moments + // if (source->modelEXT == NULL) continue; + + // replace object in image + if (source->mode & PM_SOURCE_MODE_SUBTRACTED) { + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + } + Next ++; + + // save the modelFlux here in case we need to subtract it (for failure) + psImage *modelFluxStart = psMemIncrRefCounter (source->modelFlux); + + if (savePics) { + psphotSaveImage (NULL, readout->image, "image.xp.fits"); + } + + // array to store the pointers to the model flux images while the models are being fitted + psArray *modelFluxes = psArrayAllocEmpty (4); + + // allocate the array to store the model fits + if (source->modelFits == NULL) { + source->modelFits = psArrayAllocEmpty (4); + } + + // loop here over the models chosen for each source (exclude by S/N) + psMetadataIterator *iter = psMetadataIteratorAlloc (models, PS_LIST_HEAD, NULL); + psMetadataItem *item = NULL; + while ((item = psMetadataGetAndIncrement (iter)) != NULL) { + + // XXX this should have been forced above + assert (item->type == PS_DATA_METADATA); + psMetadata *model = (psMetadata *) item->data.md; + + // check the SNlim and skip model if source is too faint + float SNlim = psMetadataLookupF32 (&status, model, "SNLIM_VALUE"); + assert (status); + + // limit selection to some SN limit + assert (source->peak); // how can a source not have a peak? + if (source->peak->SN < SNlim) continue; + + // check on the model type + pmModelType modelType = psMetadataLookupS32 (&status, model, "MODEL_TYPE"); + assert (status); + + // check on the PSF-Convolution status + bool convolved = psMetadataLookupBool (&status, model, "PSF_CONVOLVED_VALUE"); + assert (status); + + // XXX psTraceSetLevel ("psLib.math.psMinimizeLMChi2", 6); + // XXX psTraceSetLevel ("psphot.psphotModelWithPSF_LMM", 6); + + // fit the model as convolved or not + pmModel *modelFit = NULL; + if (convolved) { + modelFit = psphotPSFConvModel (readout, source, modelType, maskVal, psfSize); + if (!modelFit) { + psTrace ("psphot", 5, "failed to fit psf-conv model for object at %f, %f", source->moments->x, source->moments->y); + continue; + } + psTrace ("psphot", 4, "fit psf-conv model for %f, %f : %s chisq = %f\n", source->moments->x, source->moments->y, pmModelClassGetName (modelFit->type), modelFit->chisq); + Nconvolve ++; + if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) { + NconvolvePass ++; + } + } else { + psFree (source->modelFlux); + source->modelFlux = NULL; + modelFit = psphotFitEXT (readout, source, modelType, maskVal); + if (!modelFit) { + psTrace ("psphot", 5, "failed to fit plain model for object at %f, %f", source->moments->x, source->moments->y); + continue; + } + pmSourceCacheModel (source, maskVal); + psTrace ("psphot", 4, "fit plain model for %f, %f : %s chisq = %f\n", source->moments->x, source->moments->y, pmModelClassGetName (modelFit->type), modelFit->chisq); + Nplain ++; + if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) { + NplainPass ++; + } + } + + // save each of the model flux images and store the best + psArrayAdd (modelFluxes, 4, source->modelFlux); + + // test for fit quality / result + psArrayAdd (source->modelFits, 4, modelFit); + + psFree (modelFit); + } + psFree (iter); + + // evaluate the relative quality of the models, choose one + float minChisq = NAN; + int minModel = -1; + for (int i = 0; i < source->modelFits->n; i++) { + pmModel *model = source->modelFits->data[i]; + + if (!(model->flags & PM_MODEL_STATUS_FITTED)) continue; + + if (model->flags & (PM_MODEL_STATUS_BADARGS)) continue; + if (model->flags & (PM_MODEL_STATUS_NONCONVERGE)) continue; + if (model->flags & (PM_MODEL_STATUS_OFFIMAGE)) continue; + + if ((minModel < 0) || (model->chisq < minChisq)) { + minChisq = model->chisq; + minModel = i; + } + } + + if (minModel == -1) { + // re-subtract the object, leave local sky + psTrace ("psphot", 5, "failed to fit extended source model to object at %f, %f", source->moments->x, source->moments->y); + + // replace original model, subtract it + psFree (source->modelFlux); + source->modelFlux = modelFluxStart; + + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + source->mode |= PM_SOURCE_MODE_SUBTRACTED; + + psFree (modelFluxes); + + if (savePics) { + psphotSaveImage (NULL, readout->image, "image.xp.fits"); + char key[10]; + fprintf (stdout, "continue? "); + fgets (key, 8, stdin); + if (key[0] == 'n') { + savePics = false; + } + } + continue; + } + + // save the best extended model in modelEXT + psFree (source->modelEXT); + source->modelEXT = psMemIncrRefCounter (source->modelFits->data[minModel]); + + // save the modelFlux for the best model + psFree (source->modelFlux); + source->modelFlux = psMemIncrRefCounter (modelFluxes->data[minModel]); + + // subtract the best fit from the object, leave local sky + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + source->mode |= PM_SOURCE_MODE_SUBTRACTED; + + // the initial model flux is no longer needed + psFree (modelFluxStart); + psFree (modelFluxes); + + psTrace ("psphot", 4, "best ext model for %f, %f : %s chisq = %f\n", source->moments->x, source->moments->y, pmModelClassGetName (source->modelEXT->type), source->modelEXT->chisq); + psTrace ("psphot", 5, "extended source model for source at %7.1f, %7.1f", source->moments->x, source->moments->y); + + if (savePics) { + psphotSaveImage (NULL, readout->image, "image.xm.fits"); + char key[10]; + fprintf (stdout, "continue? "); + fgets (key, 8, stdin); + if (key[0] == 'n') { + savePics = false; + } + } + } + + psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot"), Next); + psLogMsg ("psphot", PS_LOG_INFO, " %d convolved models (%d passed)\n", Nconvolve, NconvolvePass); + psLogMsg ("psphot", PS_LOG_INFO, " %d plain models (%d passed)\n", Nplain, NplainPass); + return true; +} Index: /trunk/psphot/src/psphotFindDetections.c =================================================================== --- /trunk/psphot/src/psphotFindDetections.c (revision 17395) +++ /trunk/psphot/src/psphotFindDetections.c (revision 17396) @@ -8,4 +8,6 @@ bool status; int pass; + + psTimerStart ("psphot"); // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) @@ -67,5 +69,5 @@ psphotCullPeaks(readout->image, readout->weight, recipe, detections->footprints); detections->peaks = pmFootprintArrayToPeaks(detections->footprints); - psLogMsg ("psphot", PS_LOG_INFO, "peaks: %ld, total footprints: %ld\n", detections->peaks->n, detections->footprints->n); + psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks, %ld total footprints: %f sec\n", detections->peaks->n, detections->footprints->n, psTimerMark ("psphot")); return detections; Index: /trunk/psphot/src/psphotFindPeaks.c =================================================================== --- /trunk/psphot/src/psphotFindPeaks.c (revision 17395) +++ /trunk/psphot/src/psphotFindPeaks.c (revision 17396) @@ -11,5 +11,5 @@ // smooth the image and weight map - psTimerStart ("psphot"); + psTimerStart ("peaks"); // XXX if we have been supplied a PSF, we can use that to set the FWHM_X,FWHM_Y values @@ -32,10 +32,10 @@ psImage *smooth_im = psImageCopy (NULL, readout->image, PS_TYPE_F32); psImageSmoothMaskF32 (smooth_im, readout->mask, maskVal, SIGMA_SMTH, NSIGMA_SMTH); - psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth image: %f sec\n", psTimerMark ("psphot")); + psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth image: %f sec\n", psTimerMark ("peaks")); // smooth the weight, applying the mask as we go psImage *smooth_wt = psImageCopy (NULL, readout->weight, PS_TYPE_F32); psImageSmoothMaskF32 (smooth_wt, readout->mask, maskVal, SIGMA_SMTH/sqrt(2), NSIGMA_SMTH); - psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth weight: %f sec\n", psTimerMark ("psphot")); + psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth weight: %f sec\n", psTimerMark ("peaks")); psImage *mask = readout->mask; @@ -61,5 +61,5 @@ } } - psLogMsg ("psphot", PS_LOG_INFO, "built smoothed signficance image: %f sec\n", psTimerMark ("psphot")); + psLogMsg ("psphot", PS_LOG_INFO, "built smoothed signficance image: %f sec\n", psTimerMark ("peaks")); // optionally save example images under trace @@ -70,5 +70,5 @@ } - psTimerStart ("psphot"); + psTimerStart ("peaks"); // set peak threshold @@ -140,5 +140,5 @@ pmPeaksWriteText (peaks, output); } - psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks: %f sec\n", peaks->n, psTimerMark ("psphot")); + psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks: %f sec\n", peaks->n, psTimerMark ("peaks")); // If they asked us to return a set of pmFootprints, not a raw pmPeak list, @@ -167,5 +167,5 @@ peaks = footprints; // well, you know what I mean - psLogMsg ("psphot", PS_LOG_INFO, "%ld footprints: %f sec\n", peaks->n, psTimerMark ("psphot")); + psLogMsg ("psphot", PS_LOG_INFO, "%ld footprints: %f sec\n", peaks->n, psTimerMark ("peaks")); } Index: /trunk/psphot/src/psphotIsophotal.c =================================================================== --- /trunk/psphot/src/psphotIsophotal.c (revision 17395) +++ /trunk/psphot/src/psphotIsophotal.c (revision 17396) @@ -1,3 +1,5 @@ -# include "psphot.h" +# include "psphotInternal.h" + +static float ISOPHOT_FLUX = NAN; bool psphotIsophotal (pmSource *source, psMetadata *recipe, psMaskType maskVal) { @@ -14,9 +16,12 @@ // flux at which to measure isophotal parameters - // XXX cache this? - float ISOPHOT_FLUX = psMetadataLookupF32 (&status, recipe, "ISOPHOTAL_FLUX"); + if (!isfinite(ISOPHOT_FLUX)) { + // XXX ISOPHOTAL_FLUX should be specified in mags, need the zero point to get counts/sec + ISOPHOT_FLUX = psMetadataLookupF32 (&status, recipe, "ISOPHOTAL_FLUX"); + assert (status); + } // find the first bin below the flux level and the last above the level - // XXX can this be done faster with disection? + // 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? @@ -27,5 +32,17 @@ if ((firstBelow < 0) && (flux->data.F32[i] < ISOPHOT_FLUX)) firstBelow = i; } - + // if we don't go out far enough, we have a problem... + if (lastAbove == flux->n - 1) { + psTrace ("psphot", 5, "did not go out far enough to reach isophotal magnitude"); + // XXX raise a flag ? + return false; + } + if (firstBelow < 0) { + psTrace ("psphot", 5, "did not go out far enough to bound isophotal magnitude: error unmeasured"); + // XXX raise a flag ? + lastAbove = firstBelow; + return false; + } + // need to examine pixels in this vicinity float isophotalFluxFirst = 0; @@ -56,4 +73,9 @@ source->extpars->isophot->radErr = isophotalRadErr; + psTrace ("psphot", 5, "Isophot flux:%f +/- %f @ %f +/- %f for %f, %f\n", + source->extpars->isophot->mag, source->extpars->isophot->magErr, + source->extpars->isophot->rad, source->extpars->isophot->radErr, + source->peak->xf, source->peak->yf); + return true; Index: /trunk/psphot/src/psphotKernelFromPSF.c =================================================================== --- /trunk/psphot/src/psphotKernelFromPSF.c (revision 17395) +++ /trunk/psphot/src/psphotKernelFromPSF.c (revision 17396) @@ -1,20 +1,45 @@ -# include "psphot.h" +# include "psphotInternal.h" psKernel *psphotKernelFromPSF (pmSource *source, int nPix) { - assert (source); - assert (source->psfFlux); // XXX build if needed? + assert (source); + assert (source->psfFlux); // XXX build if needed? - int x0 = source->peak->xf - source->psfFlux->col0; - int y0 = source->peak->yf - source->psfFlux->row0; + int x0 = source->peak->xf - source->psfFlux->col0; + int y0 = source->peak->yf - source->psfFlux->row0; - // need to decide on the size: dynamically? statically? - psKernel *psf = psKernelAlloc (-nPix, +nPix, -nPix, +nPix); + // need to decide on the size: dynamically? statically? + psKernel *psf = psKernelAlloc (-nPix, +nPix, -nPix, +nPix); - for (int j = psf->yMin; j <= psf->yMax; j++) { - for (int i = psf->xMin; i <= psf->xMax; i++) { - psf->kernel[j][i] = source->psfFlux->data.F32[y0 + j][x0 + i]; + // XXX we should just re-construct a PSF at this location + // psModelAdd (psf->image, NULL, source->modelPSF, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM | PM_MODEL_OP_CENTER); + + // if the realized PSF for this object does not cover the full kernel, give up for now + if (x0 + psf->xMin < 0) goto escape; + if (x0 + psf->xMax >= source->psfFlux->numCols) goto escape; + if (y0 + psf->yMin < 0) goto escape; + if (y0 + psf->yMax >= source->psfFlux->numRows) goto escape; + + double sum = 0.0; + for (int j = psf->yMin; j <= psf->yMax; j++) { + for (int i = psf->xMin; i <= psf->xMax; i++) { + double value = source->psfFlux->data.F32[y0 + j][x0 + i]; + psf->kernel[j][i] = value; + sum += value; + } } - } - return psf; + assert (sum > 0.0); + + // psf must be normalized (integral = 1.0) + for (int i = 0; i < psf->image->numRows; i++) { + for (int j = 0; j < psf->image->numCols; j++) { + psf->image->data.F32[i][j] /= sum; + } + } + + return psf; + +escape: + psFree (psf); + return NULL; } Index: /trunk/psphot/src/psphotKron.c =================================================================== --- /trunk/psphot/src/psphotKron.c (revision 17395) +++ /trunk/psphot/src/psphotKron.c (revision 17396) @@ -1,3 +1,3 @@ -# include "psphot.h" +# include "psphotInternal.h" bool psphotKron (pmSource *source, psMetadata *recipe, psMaskType maskVal) { Index: /trunk/psphot/src/psphotModelTest.c =================================================================== --- /trunk/psphot/src/psphotModelTest.c (revision 17395) +++ /trunk/psphot/src/psphotModelTest.c (revision 17396) @@ -208,6 +208,10 @@ source->modelPSF = pmModelFromPSF (model, psf); source->modelEXT = model; - status = psphotPSFConvModel (source, recipe, maskVal); - model = source->modelConv; + + // what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box) + int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE"); + assert (status); + + model = psphotPSFConvModel (readout, source, modelType, maskVal, psfSize); params = model->params->data.F32; } else { Index: /trunk/psphot/src/psphotModelWithPSF.c =================================================================== --- /trunk/psphot/src/psphotModelWithPSF.c (revision 17395) +++ /trunk/psphot/src/psphotModelWithPSF.c (revision 17396) @@ -1,3 +1,4 @@ -# include "psphot.h" +# include "psphotInternal.h" +# define SAVE_IMAGES 0 bool psphotModelWithPSF_LMM ( @@ -36,8 +37,15 @@ // allocate internal arrays (current vs Guess) - psImage *alpha = psImageAlloc(params->n, params->n, PS_TYPE_F32); - psImage *Alpha = psImageAlloc(params->n, params->n, PS_TYPE_F32); - psVector *beta = psVectorAlloc(params->n, PS_TYPE_F32); - psVector *Beta = psVectorAlloc(params->n, PS_TYPE_F32); + psImage *Alpha = NULL; + psVector *Beta = NULL; + + // Alpha & Beta only contain elements to represent the unmasked parameters + if (!psMinLM_AllocAB (&Alpha, &Beta, params, paramMask)) { + psAbort ("programming error: no unmasked parameters to be fit\n"); + } + + // allocate internal arrays (current vs Guess) + psImage *alpha = psImageAlloc(Alpha->numCols, Alpha->numRows, PS_TYPE_F32); + psVector *beta = psVectorAlloc(Beta->n, PS_TYPE_F32); psVector *Params = psVectorAlloc(params->n, PS_TYPE_F32); @@ -79,6 +87,7 @@ // dump some useful info if trace is defined if (psTraceGetLevel("psphot") >= 6) { - p_psImagePrint(psTraceGetDestination(), Alpha, "alpha guess (1)"); - p_psVectorPrint(psTraceGetDestination(), Beta, "beta guess (1)"); + p_psImagePrint(psTraceGetDestination(), Alpha, "Alpha guess (1)"); + p_psVectorPrint(psTraceGetDestination(), Beta, "Beta guess (1)"); + p_psVectorPrint(psTraceGetDestination(), beta, "beta current (1)"); } if (psTraceGetLevel("psphot") >= 5) { @@ -116,5 +125,5 @@ beta = psVectorCopy(beta, Beta, PS_TYPE_F32); params = psVectorCopy(params, Params, PS_TYPE_F32); - lambda *= 0.1; + lambda *= 0.25; // save the new convolved model image @@ -130,7 +139,21 @@ // construct & return the covariance matrix (if requested) if (covar != NULL) { - if (!psMinLM_GuessABP(covar, Beta, Params, alpha, beta, params, paramMask, NULL, 0.0, NULL)) { + if (!psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask, NULL, 0.0, NULL)) { psTrace ("psphot", 5, "failure to calculate covariance matrix\n"); } + // set covar values which are not masked + psImageInit (covar, 0.0); + for (int j = 0, J = 0; j < params->n; j++) { + if (paramMask && (paramMask->data.U8[j])) { + covar->data.F32[j][j] = 1.0; + continue; + } + for (int k = 0, K = 0; k < params->n; k++) { + if (paramMask && (paramMask->data.U8[k])) continue; + covar->data.F32[j][k] = Alpha->data.F32[J][K]; + K++; + } + J++; + } } @@ -177,5 +200,5 @@ } - // generate the model and derivative images for this parameter set + // 1 *** generate the model and derivative images for this parameter set // storage for model derivatives @@ -223,7 +246,7 @@ for (int n = 0; n < params->n; n++) { - if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; } - psImage *dmodel = pcm->dmodels->data[n]; - dmodel->data.F32[i][j] = deriv->data.F32[n]; + if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; } + psImage *dmodel = pcm->dmodels->data[n]; + dmodel->data.F32[i][j] = deriv->data.F32[n]; } } @@ -231,17 +254,16 @@ psFree(coord); psFree(deriv); - // convolve model and dmodel arrays with PSF psImageConvolveDirect (pcm->modelConv, pcm->model, psf); for (int n = 0; n < pcm->dmodels->n; n++) { - if (pcm->dmodels->data[n] == NULL) continue; - psImage *dmodel = pcm->dmodels->data[n]; - psImage *dmodelConv = pcm->dmodelsConv->data[n]; - psImageConvolveDirect (dmodelConv, dmodel, psf); + if (pcm->dmodels->data[n] == NULL) continue; + psImage *dmodel = pcm->dmodels->data[n]; + psImage *dmodelConv = pcm->dmodelsConv->data[n]; + psImageConvolveDirect (dmodelConv, dmodel, psf); } // XXX TEST : SAVE IMAGES -# if (1) +# if (SAVE_IMAGES) psphotSaveImage (NULL, psf->image, "psf.fits"); psphotSaveImage (NULL, pcm->model, "model.fits"); @@ -252,4 +274,6 @@ # endif + // 2 *** accumulate alpha & beta + // zero alpha and beta for summing below psImageInit (alpha, 0.0); @@ -259,5 +283,5 @@ for (psS32 i = 0; i < source->pixels->numRows; i++) { for (psS32 j = 0; j < source->pixels->numCols; j++) { - // XXX are we doing the right thing with the mask? + // XXX are we doing the right thing with the mask? // skip masked points if (source->maskObj->data.U8[i][j]) { @@ -279,23 +303,20 @@ chisq += PS_SQR(delta) * yweight; - if (isnan(delta)) - psAbort("nan in delta"); - if (isnan(chisq)) - psAbort("nan in chisq"); - - for (psS32 n1 = 0; n1 < params->n; n1++) { - if ((paramMask != NULL) && (paramMask->data.U8[n1])) { - continue; - } - psImage *dmodel = pcm->dmodelsConv->data[n1]; - float weight = dmodel->data.F32[i][j] * yweight; - for (psS32 n2 = 0; n2 <= n1; n2++) { - if ((paramMask != NULL) && (paramMask->data.U8[n2])) { - continue; - } - dmodel = pcm->dmodelsConv->data[n2]; - alpha->data.F32[n1][n2] += weight * dmodel->data.F32[i][j]; - } - beta->data.F32[n1] += weight * delta; + if (isnan(delta)) psAbort("nan in delta"); + if (isnan(chisq)) psAbort("nan in chisq"); + + // alpha & beta only contain unmasked elements + for (int n1 = 0, N1 = 0; n1 < params->n; n1++) { + if ((paramMask != NULL) && (paramMask->data.U8[n1])) continue; + psImage *dmodel = pcm->dmodelsConv->data[n1]; + float weight = dmodel->data.F32[i][j] * yweight; + for (int n2 = 0, N2 = 0; n2 <= n1; n2++) { + if ((paramMask != NULL) && (paramMask->data.U8[n2])) continue; + dmodel = pcm->dmodelsConv->data[n2]; + alpha->data.F32[N1][N2] += weight * dmodel->data.F32[i][j]; + N2++; + } + beta->data.F32[N1] += weight * delta; + N1++; } } @@ -303,17 +324,7 @@ // calculate lower-left half of alpha - for (psS32 j = 1; j < params->n; j++) { + for (psS32 j = 1; j < alpha->numCols; j++) { for (psS32 k = 0; k < j; k++) { alpha->data.F32[k][j] = alpha->data.F32[j][k]; - } - } - - // fill in pivots if we apply a mask - if (paramMask != NULL) { - for (psS32 j = 0; j < params->n; j++) { - if (paramMask->data.U8[j]) { - alpha->data.F32[j][j] = 1; - beta->data.F32[j] = 1; - } } } @@ -345,6 +356,7 @@ pcm->dmodels = psArrayAlloc (params->n); for (psS32 n = 0; n < params->n; n++) { - if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; } - pcm->dmodels->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32); + pcm->dmodels->data[n] = NULL; + if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; } + pcm->dmodels->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32); } @@ -353,6 +365,7 @@ pcm->dmodelsConv = psArrayAlloc (params->n); for (psS32 n = 0; n < params->n; n++) { - if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; } - pcm->dmodelsConv->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32); + pcm->dmodelsConv->data[n] = NULL; + if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; } + pcm->dmodelsConv->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32); } @@ -378,8 +391,8 @@ while () { - GuessABP (Alpha, Beta, Params, alpha, beta, params, paramMask, checkLimits, lambda) - dLinear = dLinear(Beta, beta, lambda); - chisq = SetABX (alpha, beta, params, paramMask, x, y, dy, func) - convergence tests... + GuessABP (Alpha, Beta, Params, alpha, beta, params, paramMask, checkLimits, lambda) + dLinear = dLinear(Beta, beta, lambda); + chisq = SetABX (alpha, beta, params, paramMask, x, y, dy, func) + convergence tests... } @@ -388,16 +401,16 @@ ** GuessABP: - f_c = sum_i (kern_i * func (x_i; p_o)) - -df_c/dp_o = d/dp_o [sum_i (kern_i * func (x_i; p_o))] - -df_c/dp_o = sum_i (d/dp_o [kern_i * func (x_i; p_o)]) - -df_c/dp_o = sum_i (kern_i * d/dp_o [func (x_i; p_o)]) - -- generate image arrays for func, dfunc/dp_j (not masked) -- convolve each with psf -- measure delta = f_conv - data -- etc + f_c = sum_i (kern_i * func (x_i; p_o)) + + df_c/dp_o = d/dp_o [sum_i (kern_i * func (x_i; p_o))] + + df_c/dp_o = sum_i (d/dp_o [kern_i * func (x_i; p_o)]) + + df_c/dp_o = sum_i (kern_i * d/dp_o [func (x_i; p_o)]) + + - generate image arrays for func, dfunc/dp_j (not masked) + - convolve each with psf + - measure delta = f_conv - data + - etc */ Index: /trunk/psphot/src/psphotPSFConvModel.c =================================================================== --- /trunk/psphot/src/psphotPSFConvModel.c (revision 17395) +++ /trunk/psphot/src/psphotPSFConvModel.c (revision 17396) @@ -1,52 +1,24 @@ -# include "psphot.h" +# include "psphotInternal.h" +# define USE_DELTA_PSF 0 // save as static values so they may be set externally static psF32 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = 15; static psF32 PM_SOURCE_FIT_MODEL_TOLERANCE = 0.1; -// static psF32 PM_SOURCE_FIT_MODEL_WEIGHT = 1.0; -// static bool PM_SOURCE_FIT_MODEL_PIX_WEIGHTS = true; // input source has both modelPSF and modelEXT. on successful exit, we set the // modelConv to contain the fitted parameters, and the modelFlux to contain the // convolved model image. -bool psphotPSFConvModel (pmSource *source, psMetadata *recipe, psMaskType maskVal) { +pmModel *psphotPSFConvModel (pmReadout *readout, pmSource *source, pmModelType modelType, psMaskType maskVal, int psfSize) { - bool status; - - int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE"); - if (!status) { - psfSize = 2; - } - // make sure we save a cached copy of the psf flux pmSourceCachePSF (source, maskVal); // convert the cached cached psf model for this source to a psKernel - // XXX for the moment, hard-wire the kernel to be 5x5 (2 pix radius) - // XXX for the moment, hard-wire the kernel to be 9x9 (4 pix radius) psKernel *psf = psphotKernelFromPSF (source, psfSize); + if (!psf) return NULL; - // psf must be normalized (integral = 1.0) - double sum = 0.0; - for (int i = 0; i < psf->image->numRows; i++) { - for (int j = 0; j < psf->image->numCols; j++) { - sum += psf->image->data.F32[i][j]; - } - } - assert (sum > 0.0); - for (int i = 0; i < psf->image->numRows; i++) { - for (int j = 0; j < psf->image->numCols; j++) { - psf->image->data.F32[i][j] /= sum; - } - } - -# if (0) - // XXX sanity check: convolve with delta function should behave like unconvolved version - for (int i = 0; i < psf->image->numRows; i++) { - for (int j = 0; j < psf->image->numCols; j++) { - psf->image->data.F32[i][j] = 0.0; - } - } - psf->image->data.F32[2][2] = 1.0; +# if (USE_DELTA_PSF) + psImageInit (psf->image, 0.0); + psf->image->data.F32[(int)(0.5*psf->image->numRows)][(int)(0.5*psf->image->numCols)] = 1.0; # endif @@ -54,7 +26,44 @@ // XXX we could modify the parameter values or even the model // here based on the observed seeing (some lookup table...) - pmModel *modelConv = pmModelCopy (source->modelEXT); + + // use the source moments, etc to guess basic model parameters + pmModel *modelConv = pmSourceModelGuess (source, modelType); + if (!modelConv) { + psFree (psf); + return NULL; + } + + // XXX test : modify the Io, SXX, SYY terms based on the psf SXX, SYY terms: + psEllipseShape psfShape; + psfShape.sx = source->modelPSF->params->data.F32[PM_PAR_SXX] / M_SQRT2; + psfShape.sxy = source->modelPSF->params->data.F32[PM_PAR_SXY]; + psfShape.sy = source->modelPSF->params->data.F32[PM_PAR_SYY] / M_SQRT2; + psEllipseAxes psfAxes = psEllipseShapeToAxes (psfShape, 20.0); + + // XXX test : modify the Io, SXX, SYY terms based on the psf SXX, SYY terms: + psEllipseShape extShape; + extShape.sx = modelConv->params->data.F32[PM_PAR_SXX] / M_SQRT2; + extShape.sxy = modelConv->params->data.F32[PM_PAR_SXY]; + extShape.sy = modelConv->params->data.F32[PM_PAR_SYY] / M_SQRT2; + psEllipseAxes extAxes = psEllipseShapeToAxes (extShape, 20.0); + + // decrease the initial guess ellipse by psf_minor axis: + psEllipseAxes extAxesMod; + extAxesMod.major = sqrt (PS_MAX (1.0, PS_SQR(extAxes.major) - PS_SQR(psfAxes.minor))); + extAxesMod.minor = sqrt (PS_MAX (1.0, PS_SQR(extAxes.minor) - PS_SQR(psfAxes.minor))); + extAxesMod.theta = extAxes.theta; + + psEllipseShape extShapeMod = psEllipseAxesToShape (extAxesMod); + modelConv->params->data.F32[PM_PAR_SXX] = extShapeMod.sx * M_SQRT2; + modelConv->params->data.F32[PM_PAR_SXY] = extShapeMod.sxy; + modelConv->params->data.F32[PM_PAR_SYY] = extShapeMod.sy * M_SQRT2; + + // increase the initial guess central intensity by 2pi r^2: + modelConv->params->data.F32[PM_PAR_I0] *= (1.0 + PS_SQR(psfAxes.minor) / PS_SQR(extAxesMod.minor)); + psVector *params = modelConv->params; psVector *dparams = modelConv->dparams; + + psphotCheckRadiusEXT (readout, source, modelConv); // create the minimization constraints @@ -91,5 +100,5 @@ psTrace ("psphot", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value); - // renormalize output model image + // renormalize output model image (generated by fitting process) float Io = params->data.F32[PM_PAR_I0]; for (int iy = 0; iy < source->modelFlux->numRows; iy++) { @@ -115,17 +124,13 @@ onPic &= (params->data.F32[PM_PAR_YPOS] >= source->pixels->row0); onPic &= (params->data.F32[PM_PAR_YPOS] < source->pixels->row0 + source->pixels->numRows); - if (!onPic) { - modelConv->flags |= PM_MODEL_STATUS_OFFIMAGE; - } + if (!onPic) modelConv->flags |= PM_MODEL_STATUS_OFFIMAGE; - source->mode |= PM_SOURCE_MODE_FITTED; - source->modelConv = modelConv; + source->mode |= PM_SOURCE_MODE_FITTED; // XXX is this needed? + psFree(psf); psFree(myMin); psFree(covar); psFree(constraint); - bool retval = (onPic && fitStatus); - psTrace("psphot", 5, "---- %s(%d) end ----\n", __func__, retval); - return(retval); + return modelConv; } Index: /trunk/psphot/src/psphotPetrosian.c =================================================================== --- /trunk/psphot/src/psphotPetrosian.c (revision 17395) +++ /trunk/psphot/src/psphotPetrosian.c (revision 17396) @@ -1,3 +1,6 @@ -# include "psphot.h" +# include "psphotInternal.h" + +static float PETROSIAN_R0 = NAN; +static float PETROSIAN_RF = NAN; bool psphotPetrosian (pmSource *source, psMetadata *recipe, psMaskType maskVal) { @@ -14,14 +17,27 @@ // flux at which to measure isophotal parameters - // XXX cache this? - float PETROSIAN_R0 = psMetadataLookupF32 (&status, recipe, "PETROSIAN_R0"); - float PETROSIAN_RF = psMetadataLookupF32 (&status, recipe, "PETROSIAN_FLUX_RATIO"); + if (!isfinite(PETROSIAN_R0)) { + PETROSIAN_R0 = psMetadataLookupF32 (&status, recipe, "PETROSIAN_R0"); + PETROSIAN_RF = psMetadataLookupF32 (&status, recipe, "PETROSIAN_FLUX_RATIO"); + assert (status); + } // first find flux at R0 - int firstBelow = -1; - int lastAbove = -1; + int firstAbove = -1; + int lastBelow = -1; for (int i = 0; i < radius->n; i++) { - if (radius->data.F32[i] > PETROSIAN_R0) lastAbove = i; - if ((firstBelow < 0) && (flux->data.F32[i] < PETROSIAN_R0)) firstBelow = 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; } @@ -29,5 +45,5 @@ float fluxR0 = 0.0; int fluxRn = 0; - for (int i = PS_MIN(firstBelow, lastAbove); i <= PS_MAX(firstBelow, lastAbove); i++) { + for (int i = PS_MIN(firstAbove, lastBelow); i <= PS_MAX(firstAbove, lastBelow); i++) { fluxR0 += flux->data.F32[i]; fluxRn ++; @@ -42,9 +58,20 @@ // XXX do I need to worry about crazy outliers? // XXX should i be smoothing or fitting the curve? - firstBelow = -1; - lastAbove = -1; + 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; } @@ -78,4 +105,9 @@ 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; Index: /trunk/psphot/src/psphotRadialProfile.c =================================================================== --- /trunk/psphot/src/psphotRadialProfile.c (revision 17395) +++ /trunk/psphot/src/psphotRadialProfile.c (revision 17396) @@ -41,6 +41,15 @@ int n = 0; - float Xo = source->modelEXT->params->data.F32[PM_PAR_XPOS] - source->pixels->col0; - float Yo = source->modelEXT->params->data.F32[PM_PAR_YPOS] - source->pixels->row0; + + 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; + } for (int iy = 0; iy < source->pixels->numRows; iy++) { for (int ix = 0; ix < source->pixels->numCols; ix++) { Index: /trunk/psphot/src/psphotReadout.c =================================================================== --- /trunk/psphot/src/psphotReadout.c (revision 17395) +++ /trunk/psphot/src/psphotReadout.c (revision 17396) @@ -117,6 +117,5 @@ psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE); - // XXX test the CR/EXT measurement - // XXX we need to call this here and after the second pass: option for starting number? + // identify CRs and extended sources psphotSourceSize (readout, sources, recipe, 0); if (!strcasecmp (breakPt, "ENSEMBLE")) { @@ -127,11 +126,13 @@ psphotBlendFit (readout, sources, recipe, psf); - // replace all sources (make this part of psphotFitSourcesLinear?) + // replace all sources psphotReplaceAll (sources, recipe); // linear fit to include all sources psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE); + + // if we only do one pass, skip to extended source analysis if (!strcasecmp (breakPt, "PASS1")) { - goto finish; + goto pass1finish; } @@ -145,7 +146,8 @@ psphotAddNoise (readout, sources, recipe); + // find fainter sources (pass 2) detections = psphotFindDetections (detections, readout, recipe); - // remove noise for subtracted objects + // remove noise for subtracted objects (ie, return to normal noise level) psphotSubNoise (readout, sources, recipe); @@ -172,8 +174,12 @@ psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE); +pass1finish: + // measure source size for the remaining sources psphotSourceSize (readout, sources, recipe, 0); - psphotExtendedSources (readout, sources, recipe); + psphotExtendedSourceAnalysis (readout, sources, recipe); + + psphotExtendedSourceFits (readout, sources, recipe); finish: Index: /trunk/psphot/src/psphotSourceFits.c =================================================================== --- /trunk/psphot/src/psphotSourceFits.c (revision 17395) +++ /trunk/psphot/src/psphotSourceFits.c (revision 17396) @@ -220,5 +220,5 @@ pmSource *tmpSrc = pmSourceAlloc (); - pmModel *EXT = psphotFitEXT (readout, source, maskVal); + pmModel *EXT = psphotFitEXT (readout, source, modelTypeEXT, maskVal); okEXT = psphotEvalEXT (tmpSrc, EXT); chiEXT = EXT->chisq / EXT->nDOF; @@ -267,4 +267,6 @@ // save new model source->modelEXT = EXT; + source->type = PM_SOURCE_TYPE_EXTENDED; + source->mode |= PM_SOURCE_MODE_EXTMODEL; // build cached model and subtract @@ -349,10 +351,10 @@ } -pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, psMaskType maskVal) { +pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, pmModelType modelType, psMaskType maskVal) { NfitEXT ++; // use the source moments, etc to guess basic model parameters - pmModel *EXT = pmSourceModelGuess (source, modelTypeEXT); + pmModel *EXT = pmSourceModelGuess (source, modelType); PS_ASSERT (EXT, NULL); Index: /trunk/psphot/src/psphotSourceSize.c =================================================================== --- /trunk/psphot/src/psphotSourceSize.c (revision 17395) +++ /trunk/psphot/src/psphotSourceSize.c (revision 17396) @@ -1,4 +1,6 @@ # include "psphotInternal.h" # include + +float psphotModelContour (psImage *image, psImage *weight, psImage *mask, pmModel *model, float Ro); // we need to call this function after sources have been fitted to the PSF model and @@ -14,5 +16,8 @@ psTimerStart ("psphot"); - float CR_NSIGMA_LIMIT = psMetadataLookupF32 (&status, recipe, "PSPHOT.CRNSIGMA.LIMIT"); + float CR_NSIGMA_LIMIT = psMetadataLookupF32 (&status, recipe, "PSPHOT.CR.NSIGMA.LIMIT"); + assert (status); + + float EXT_NSIGMA_LIMIT = psMetadataLookupF32 (&status, recipe, "PSPHOT.EXT.NSIGMA.LIMIT"); assert (status); @@ -24,20 +29,24 @@ if (isfinite(source->crNsigma)) continue; - source->crNsigma = -1.0; - source->extNsigma = -1.0; - // source must have been subtracted - source->crNsigma = -3.0; if (!(source->mode & PM_SOURCE_MODE_SUBTRACTED)) continue; - psF32 **resid = source->pixels->data.F32; + psF32 **resid = source->pixels->data.F32; psF32 **weight = source->weight->data.F32; - psU8 **mask = source->maskObj->data.U8; + psU8 **mask = source->maskObj->data.U8; + + // check for extendedness: measure the delta flux significance at the 1 sigma contour + source->extNsigma = psphotModelContour (source->pixels, source->weight, source->maskObj, source->modelPSF, 1.0); + + // XXX prevent a source from being both CR and EXT? + if (source->extNsigma > EXT_NSIGMA_LIMIT) { + source->mode |= PM_SOURCE_MODE_EXT_LIMIT; + } int xPeak = source->peak->xf - source->pixels->col0 + 0.5; int yPeak = source->peak->yf - source->pixels->row0 + 0.5; - // skip sources which are too close to a boundary - source->crNsigma = -4.0; + // XXX for now, skip sources which are too close to a boundary + // XXX raise a flag? if (xPeak < 1) continue; if (xPeak > source->pixels->numCols - 2) continue; @@ -46,5 +55,5 @@ // XXX for now, just skip any sources with masked pixels - source->crNsigma = -5.0; + // XXX raise a flag? bool keep = true; for (int iy = -1; (iy <= +1) && keep; iy++) { @@ -113,5 +122,4 @@ } source->crNsigma = (nCR > 0) ? fCR / nCR : 0.0; - source->extNsigma = (nEXT > 0) ? fabs(fEXT) / nEXT : 0.0; // NOTE: abs needed to make the Nsigma value positive @@ -122,5 +130,5 @@ if (source->crNsigma > CR_NSIGMA_LIMIT) { - source->mode |= PM_SOURCE_MODE_CRLIMIT; + source->mode |= PM_SOURCE_MODE_CR_LIMIT; } } @@ -169,2 +177,65 @@ */ + +// given the PSF ellipse parameters, navigate around the 1sigma contour, return the total +// deviation in sigmas. This is measure on the residual image - should we ignore negative +// deviations? +float psphotModelContour (psImage *image, psImage *weight, psImage *mask, pmModel *model, float Ro) { + + psF32 *PAR = model->params->data.F32; + + // Ro = (x / SXX)^2 + (y / SYY)^2 + x y SXY; + // y^2 (1/SYY^2) + y (x SXY) + (x / SXX)^2 - Ro = 0; + // y = [-(x SXY) +/- sqrt ((x SXY)^2 - 4 (1/SYY^2) ((x/SXX)^2 - Ro))] * [SYY^2 / 2]; + // y = [-B +/- sqrt (B^2 - 4 A C)] / [2 A]; + + // min/max value of x is where T -> 0 + // solve this for x2: + float Q = Ro * PS_SQR(PAR[PM_PAR_SXX]) / (1.0 - PS_SQR(PAR[PM_PAR_SXX]*PAR[PM_PAR_SYY]*PAR[PM_PAR_SXY]) / 4.0); + if (Q < 0.0) return NAN; // ellipse is imaginary + + int xMax = sqrt(Q); + int xMin = -1.0*xMax; + + int nPts = 0; + float nSigma = 0.0; + + for (int x = xMin; x <= xMax; x++) { + float A = PS_SQR (1.0 / PAR[PM_PAR_SYY]); + float B = x * PAR[PM_PAR_SXY]; + float C = PS_SQR (x / PAR[PM_PAR_SXX]) - Ro; + + float T = PS_SQR(B) - 4*A*C; + if (T < 0.0) continue; + + float yP = (-B + sqrt (T)) / (2.0 * A); + float yM = (-B - sqrt (T)) / (2.0 * A); + + int xPix = x + PAR[PM_PAR_XPOS] - image->col0 + 0.5; + int yPixM = yM + PAR[PM_PAR_YPOS] - image->row0 + 0.5; + int yPixP = yP + PAR[PM_PAR_YPOS] - image->row0 + 0.5; + + if (xPix < 0) continue; + if (xPix >= image->numCols) continue; + + if ((yPixM >= 0) && (yPixM < image->numRows)) { + if (!mask || !mask->data.U8[yPixM][xPix]) { + float dSigma = image->data.F32[yPixM][xPix] / sqrt (weight->data.F32[yPixM][xPix]); + nSigma += dSigma; + nPts ++; + } + } + + if (yPixM == yPixP) continue; + + if ((yPixP >= 0) && (yPixP < image->numRows)) { + if (!mask || !mask->data.U8[yPixP][xPix]) { + float dSigma = image->data.F32[yPixP][xPix] / sqrt (weight->data.F32[yPixP][xPix]); + nSigma += dSigma; + nPts ++; + } + } + } + nSigma /= nPts; + return nSigma; +} Index: /trunk/psphot/src/psphotSourceStats.c =================================================================== --- /trunk/psphot/src/psphotSourceStats.c (revision 17395) +++ /trunk/psphot/src/psphotSourceStats.c (revision 17396) @@ -37,11 +37,15 @@ for (int i = 0; i < peaks->n; i++) { + pmPeak *peak = peaks->data[i]; + if (peak->assigned) continue; + // create a new source, add peak pmSource *source = pmSourceAlloc(); - source->peak = (pmPeak *)psMemIncrRefCounter(peaks->data[i]); + source->peak = psMemIncrRefCounter(peak); // allocate image, weight, mask arrays for each peak (square of radius OUTER) pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER); if (!strcasecmp (breakPt, "PEAKS")) { + peak->assigned = true; psArrayAdd (sources, 100, source); psFree (source); @@ -49,6 +53,7 @@ } - // skip faint sources + // skip faint sources for moments measurement if (source->peak->SN < MIN_SN) { + peak->assigned = true; psArrayAdd (sources, 100, source); psFree (source); @@ -57,5 +62,5 @@ // measure a local sky value - // XXX EAM : this should use ROBUST not SAMPLE median, but it is broken + // the local sky is now ignored; kept here for reference only status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, mark); if (!status) { @@ -80,4 +85,5 @@ if (status) { // add to the source array + peak->assigned = true; psArrayAdd (sources, 100, source); psFree (source); @@ -93,4 +99,5 @@ if (status) { // add to the source array + peak->assigned = true; psArrayAdd (sources, 100, source); psFree (source);