Changeset 17396
- Timestamp:
- Apr 8, 2008, 8:36:06 AM (18 years ago)
- Location:
- trunk
- Files:
-
- 2 added
- 30 edited
-
psModules/src/objects/models/pmModel_QGAUSS.c (modified) (2 diffs)
-
psModules/src/objects/models/pmModel_SERSIC.c (modified) (5 diffs)
-
psModules/src/objects/pmPeaks.c (modified) (2 diffs)
-
psModules/src/objects/pmPeaks.h (modified) (2 diffs)
-
psModules/src/objects/pmSource.c (modified) (6 diffs)
-
psModules/src/objects/pmSource.h (modified) (6 diffs)
-
psModules/src/objects/pmSourceIO.c (modified) (5 diffs)
-
psModules/src/objects/pmSourceIO.h (modified) (2 diffs)
-
psModules/src/objects/pmSourceIO_PS1_DEV_1.c (modified) (12 diffs)
-
psModules/src/objects/pmSourcePhotometry.c (modified) (2 diffs)
-
psphot/doc/notes.txt (modified) (1 diff)
-
psphot/src/Makefile.am (modified) (1 diff)
-
psphot/src/psphot.h (modified) (3 diffs)
-
psphot/src/psphotAnnuli.c (modified) (2 diffs)
-
psphot/src/psphotBlendFit.c (modified) (2 diffs)
-
psphot/src/psphotEvalFLT.c (modified) (5 diffs)
-
psphot/src/psphotExtendedSourceAnalysis.c (added)
-
psphot/src/psphotExtendedSourceFits.c (added)
-
psphot/src/psphotFindDetections.c (modified) (2 diffs)
-
psphot/src/psphotFindPeaks.c (modified) (6 diffs)
-
psphot/src/psphotIsophotal.c (modified) (4 diffs)
-
psphot/src/psphotKernelFromPSF.c (modified) (1 diff)
-
psphot/src/psphotKron.c (modified) (1 diff)
-
psphot/src/psphotModelTest.c (modified) (1 diff)
-
psphot/src/psphotModelWithPSF.c (modified) (16 diffs)
-
psphot/src/psphotPSFConvModel.c (modified) (4 diffs)
-
psphot/src/psphotPetrosian.c (modified) (5 diffs)
-
psphot/src/psphotRadialProfile.c (modified) (1 diff)
-
psphot/src/psphotReadout.c (modified) (4 diffs)
-
psphot/src/psphotSourceFits.c (modified) (3 diffs)
-
psphot/src/psphotSourceSize.c (modified) (7 diffs)
-
psphot/src/psphotSourceStats.c (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/objects/models/pmModel_QGAUSS.c
r15834 r17396 237 237 if (!isfinite(shape.sxy)) return false; 238 238 239 PAR[PM_PAR_SKY] = moments->Sky; 240 PAR[PM_PAR_I0] = moments->Peak - moments->Sky; 239 // XXX turn this off here for now PAR[PM_PAR_SKY] = moments->Sky; 240 PAR[PM_PAR_SKY] = 0.0; 241 PAR[PM_PAR_I0] = moments->Peak; 241 242 PAR[PM_PAR_XPOS] = peak->x; 242 243 PAR[PM_PAR_YPOS] = peak->y; … … 457 458 458 459 status = true; 459 status &= (dP < 0.5);460 // status &= (dP < 0.5); 460 461 status &= (PAR[PM_PAR_I0] > 0); 461 462 status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5); 463 464 fprintf (stderr, "QGAUSS status pars: dP: %f, I0: %f, S/N: %f\n", 465 dP, PAR[PM_PAR_I0], (dPAR[PM_PAR_I0]/PAR[PM_PAR_I0])); 462 466 463 467 return status; -
trunk/psModules/src/objects/models/pmModel_SERSIC.c
r15834 r17396 19 19 20 20 note that a standard sersic model uses exp(-K*(z^(1/n) - 1). the additional elements (K, 21 the -1 offset) are absorbed in this model by the normalization, the exponen et, and the21 the -1 offset) are absorbed in this model by the normalization, the exponent, and the 22 22 radial scale. We fit the elements in this form, then re-normalize them on output. 23 23 *****************************************************************************/ … … 157 157 break; 158 158 case PM_PAR_SXX: 159 params_min = 0. 25;159 params_min = 0.05; 160 160 break; 161 161 case PM_PAR_SYY: 162 params_min = 0. 25;162 params_min = 0.05; 163 163 break; 164 164 case PM_PAR_SXY: … … 166 166 break; 167 167 case PM_PAR_7: 168 params_min = 0. 1;168 params_min = 0.05; 169 169 break; 170 170 default: … … 247 247 if (!isfinite(shape.sxy)) return false; 248 248 249 PAR[PM_PAR_SKY] = moments->Sky; 250 PAR[PM_PAR_I0] = moments->Peak - moments->Sky; 249 // XXX PAR[PM_PAR_SKY] = moments->Sky; 250 PAR[PM_PAR_SKY] = 0.0; 251 PAR[PM_PAR_I0] = moments->Peak; 251 252 PAR[PM_PAR_XPOS] = peak->x; 252 253 PAR[PM_PAR_YPOS] = peak->y; … … 442 443 443 444 status = true; 444 status &= (dP < 0.5);445 // status &= (dP < 0.5); 445 446 status &= (PAR[PM_PAR_I0] > 0); 446 447 status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5); 448 449 fprintf (stderr, "SERSIC status pars: dP: %f, I0: %f, S/N: %f\n", 450 dP, PAR[PM_PAR_I0], (dPAR[PM_PAR_I0]/PAR[PM_PAR_I0])); 447 451 448 452 return status; -
trunk/psModules/src/objects/pmPeaks.c
r15979 r17396 6 6 * @author EAM, IfA: significant modifications. 7 7 * 8 * @version $Revision: 1. 19$ $Name: not supported by cvs2svn $9 * @date $Date: 2008-0 1-02 20:38:28 $8 * @version $Revision: 1.20 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2008-04-08 18:35:38 $ 10 10 * 11 11 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 141 141 tmp->xf = x; 142 142 tmp->yf = y; 143 143 tmp->assigned = false; 144 144 tmp->type = type; 145 145 -
trunk/psModules/src/objects/pmPeaks.h
r15984 r17396 10 10 * @author GLG, MHPCC 11 11 * 12 * @version $Revision: 1.1 0$ $Name: not supported by cvs2svn $13 * @date $Date: 2008-0 1-02 20:42:46$12 * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ 13 * @date $Date: 2008-04-08 18:35:38 $ 14 14 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii 15 15 */ … … 58 58 float flux; ///< level in unsmoothed sci image 59 59 float SN; ///< S/N implied by detection level 60 pmPeakType type; ///< Description of peak. 60 bool assigned; ///< is peak assigned to a source? 61 pmPeakType type; ///< Description of peak. 61 62 } 62 63 pmPeak; -
trunk/psModules/src/objects/pmSource.c
r17049 r17396 6 6 * @author EAM, IfA: significant modifications. 7 7 * 8 * @version $Revision: 1.5 1$ $Name: not supported by cvs2svn $9 * @date $Date: 2008-0 3-19 00:51:09$8 * @version $Revision: 1.52 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2008-04-08 18:35:38 $ 10 10 * 11 11 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 50 50 psFree(tmp->modelPSF); 51 51 psFree(tmp->modelEXT); 52 psFree(tmp->modelConv); 52 psFree(tmp->modelFits); 53 psFree(tmp->extpars); 53 54 psFree(tmp->blends); 54 55 psTrace("psModules.objects", 5, "---- end ----\n"); … … 107 108 source->modelPSF = NULL; 108 109 source->modelEXT = NULL; 109 source->model Conv= NULL;110 source->modelFits = NULL; 110 111 source->type = PM_SOURCE_TYPE_UNKNOWN; 111 112 source->mode = PM_SOURCE_MODE_DEFAULT; … … 947 948 return model; 948 949 949 // XXX when should I return the modelConv ?? 950 // the 'best' extended model is saved in source->modelEXT (may be a pointer to one of 951 // the elements of source->modelFits) 950 952 case PM_SOURCE_TYPE_EXTENDED: 951 model = source->modelConv; 952 if (!model) { 953 model = source->modelEXT; 954 } 953 model = source->modelEXT; 955 954 if (!model && source->modelPSF) { 955 // XXX raise an error or warning here? 956 956 if (isPSF) { 957 957 *isPSF = true; … … 1015 1015 if (!strcasecmp (name, "DEFECT" )) return PM_SOURCE_MODE_DEFECT; 1016 1016 if (!strcasecmp (name, "SATURATED" )) return PM_SOURCE_MODE_SATURATED; 1017 if (!strcasecmp (name, "CRLIMIT" )) return PM_SOURCE_MODE_CRLIMIT; 1017 if (!strcasecmp (name, "CRLIMIT" )) return PM_SOURCE_MODE_CR_LIMIT; 1018 if (!strcasecmp (name, "EXTLIMIT" )) return PM_SOURCE_MODE_EXT_LIMIT; 1018 1019 if (!strcasecmp (name, "SUBTRACTED")) return PM_SOURCE_MODE_SUBTRACTED; 1019 1020 return PM_SOURCE_MODE_DEFAULT; … … 1036 1037 case PM_SOURCE_MODE_DEFECT : return psStringCopy ("DEFECT" ); 1037 1038 case PM_SOURCE_MODE_SATURATED : return psStringCopy ("SATURATED" ); 1038 case PM_SOURCE_MODE_CRLIMIT : return psStringCopy ("CRLIMIT"); 1039 case PM_SOURCE_MODE_CR_LIMIT : return psStringCopy ("CRLIMIT" ); 1040 case PM_SOURCE_MODE_EXT_LIMIT : return psStringCopy ("EXTLIMIT" ); 1039 1041 case PM_SOURCE_MODE_SUBTRACTED : return psStringCopy ("SUBTRACTED"); 1040 1042 default: -
trunk/psModules/src/objects/pmSource.h
r16819 r17396 3 3 * @author EAM, IfA; GLG, MHPCC 4 4 * 5 * @version $Revision: 1.2 3$ $Name: not supported by cvs2svn $6 * @date $Date: 2008-0 3-05 01:08:08 $5 * @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ 6 * @date $Date: 2008-04-08 18:35:38 $ 7 7 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii 8 8 */ … … 22 22 * source. 23 23 * 24 * XXX: The values given below are currently illustrative and will require25 * some modification as the source classification code is developed. (TBD)26 *27 24 */ 28 25 typedef enum { 29 PM_SOURCE_TYPE_UNKNOWN, ///< a cosmic-ray30 PM_SOURCE_TYPE_DEFECT, 31 PM_SOURCE_TYPE_SATURATED, ///< random saturated pixels32 PM_SOURCE_TYPE_STAR, ///< a good-quality star33 PM_SOURCE_TYPE_EXTENDED, ///< an extended object (eg, galaxy)26 PM_SOURCE_TYPE_UNKNOWN, ///< not yet classified 27 PM_SOURCE_TYPE_DEFECT, ///< a cosmic-ray 28 PM_SOURCE_TYPE_SATURATED, ///< random saturated pixels (eg, bleed trails) 29 PM_SOURCE_TYPE_STAR, ///< a good-quality star (subtracted model is PSF) 30 PM_SOURCE_TYPE_EXTENDED, ///< an extended object (eg, galaxy) (subtracted model is EXT) 34 31 } pmSourceType; 35 32 … … 48 45 PM_SOURCE_MODE_BADPSF = 0x0400, ///< Failed to get good estimate of object's PSF 49 46 PM_SOURCE_MODE_DEFECT = 0x0800, ///< Source is thought to be a defect 50 PM_SOURCE_MODE_SATURATED = 0x1000, ///< Source is thought to be saturation 51 PM_SOURCE_MODE_CRLIMIT = 0x2000, ///< Source has crNsigma above limit 52 PM_SOURCE_MODE_SUBTRACTED = 0x4000, ///< XXX this flag is actually only used internally (move) 47 PM_SOURCE_MODE_SATURATED = 0x1000, ///< Source is thought to be saturated pixels (bleed trail) 48 PM_SOURCE_MODE_CR_LIMIT = 0x2000, ///< Source has crNsigma above limit 49 PM_SOURCE_MODE_EXT_LIMIT = 0x4000, ///< Source has extNsigma above limit 50 PM_SOURCE_MODE_SUBTRACTED = 0x8000, ///< XXX this flag is actually only used internally (move) 53 51 } pmSourceMode; 54 52 … … 60 58 * 61 59 * XXX do I have to re-organize this (again!) to allow an arbitrary set of extended model fits?? 60 * XXX put the Mag and Err inside the pmModel? 61 * XXX keep the modelEXT or add to the psArray 62 * 62 63 * 63 64 */ … … 65 66 const int id; ///< Unique ID for object 66 67 int seq; ///< ID for output (generated on write) 67 pmPeak *peak;///< Description of peak pixel.68 pmPeak *peak; ///< Description of peak pixel. 68 69 psImage *pixels; ///< Rectangular region including object pixels. 69 70 psImage *weight; ///< Image variance. … … 72 73 psImage *modelFlux; ///< cached copy of the best model for this source 73 74 psImage *psfFlux; ///< cached copy of the psf model for this source 74 pmMoments *moments; ///< Basic moments measure for the object.75 pmMoments *moments; ///< Basic moments measured for the object. 75 76 pmModel *modelPSF; ///< PSF Model fit (parameters and type) 76 pmModel *modelEXT; ///< EXT (floating) Model fit (parameters and type).77 p mModel *modelConv; ///< PSF-Convolved Model fit (parameters and type).77 pmModel *modelEXT; ///< EXT Model fit used for subtraction (parameters and type) 78 psArray *modelFits; ///< collection of extended source models (best == modelEXT) 78 79 pmSourceType type; ///< Best identification of object. 79 pmSourceMode mode; ///< Best identification ofobject.80 psArray *blends; 81 float psfMag; ///< calculated from flux in modelP sf80 pmSourceMode mode; ///< analysis flags set for object. 81 psArray *blends; ///< collection of sources thought to be confused with object 82 float psfMag; ///< calculated from flux in modelPSF 82 83 float extMag; ///< calculated from flux in modelEXT 83 84 float errMag; ///< error in psfMag OR extMag (depending on type) 84 85 float apMag; ///< apMag corresponding to psfMag or extMag (depending on type) 85 86 float pixWeight; ///< model-weighted coverage of valid pixels 86 float psfChisq; 87 float psfChisq; ///< probability of PSF 87 88 float crNsigma; ///< Nsigma deviation from PSF to CR 88 89 float extNsigma; ///< Nsigma deviation from PSF to EXT 90 float sky, skyErr; ///< The sky and its error at the center of the object 89 91 psRegion region; ///< area on image covered by selected pixels 90 float sky, skyErr; ///< The sky and its error at the center of the object91 92 pmSourceExtendedPars *extpars; ///< extended source parameters 92 93 }; -
trunk/psModules/src/objects/pmSourceIO.c
r17010 r17396 3 3 * @author EAM, IfA 4 4 * 5 * @version $Revision: 1.5 5$ $Name: not supported by cvs2svn $6 * @date $Date: 2008-0 3-17 22:04:27$5 * @version $Revision: 1.56 $ $Name: not supported by cvs2svn $ 6 * @date $Date: 2008-04-08 18:35:38 $ 7 7 * 8 8 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 327 327 328 328 // if this is not TRUE, the output files only contain the psf measurements. 329 bool XSRC_OUTPUT = psMetadataLookupBool(&status, recipe, " SAVE.XSRC");330 bool XFIT_OUTPUT = psMetadataLookupBool(&status, recipe, " SAVE.XFIT");329 bool XSRC_OUTPUT = psMetadataLookupBool(&status, recipe, "EXTENDED_SOURCE_ANALYSIS"); 330 bool XFIT_OUTPUT = psMetadataLookupBool(&status, recipe, "EXTENDED_SOURCE_FITS"); 331 331 332 332 // define the EXTNAME values for the different data segments: … … 445 445 if (xsrcname) { 446 446 if (!strcmp (exttype, "PS1_DEV_1")) { 447 status = pmSourcesWrite_PS1_DEV_1_XSRC (file->fits, sources, xsrcname );447 status = pmSourcesWrite_PS1_DEV_1_XSRC (file->fits, sources, xsrcname, recipe); 448 448 } 449 449 } … … 457 457 psFree (headname); 458 458 psFree (dataname); 459 psFree (xsrcname); 460 psFree (xfitname); 459 461 psFree (outhead); 460 462 return false; … … 466 468 psFree (headname); 467 469 psFree (dataname); 470 psFree (xsrcname); 471 psFree (xfitname); 468 472 psFree (outhead); 469 473 break; -
trunk/psModules/src/objects/pmSourceIO.h
r16819 r17396 4 4 * @author EAM, IfA; GLG, MHPCC 5 5 * 6 * @version $Revision: 1.1 4$ $Name: not supported by cvs2svn $7 * @date $Date: 2008-0 3-05 01:08:08 $6 * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $ 7 * @date $Date: 2008-04-08 18:35:38 $ 8 8 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii 9 9 * … … 27 27 bool pmSourcesWrite_PS1_DEV_0 (psFits *fits, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname); 28 28 bool pmSourcesWrite_PS1_DEV_1 (psFits *fits, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, char *xsrcname); 29 bool pmSourcesWrite_PS1_DEV_1_XSRC (psFits *fits, psArray *sources, char *extname );29 bool pmSourcesWrite_PS1_DEV_1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe); 30 30 bool pmSourcesWrite_PS1_DEV_1_XFIT (psFits *fits, psArray *sources, char *extname); 31 31 -
trunk/psModules/src/objects/pmSourceIO_PS1_DEV_1.c
r16819 r17396 3 3 * @author EAM, IfA 4 4 * 5 * @version $Revision: 1. 9$ $Name: not supported by cvs2svn $6 * @date $Date: 2008-0 3-05 01:08:08 $5 * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ 6 * @date $Date: 2008-04-08 18:35:38 $ 7 7 * 8 8 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 251 251 } 252 252 253 bool pmSourcesWrite_PS1_DEV_1_XSRC (psFits *fits, psArray *sources, char *extname )253 bool pmSourcesWrite_PS1_DEV_1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe) 254 254 { 255 255 256 bool status; 256 257 psArray *table; 257 258 psMetadata *row; 258 int i;259 259 psF32 *PAR, *dPAR; 260 psEllipseAxes axes;261 260 psF32 xPos, yPos; 262 261 psF32 xErr, yErr; … … 273 272 table = psArrayAllocEmpty (sources->n); 274 273 274 // which extended source analyses should we perform? 275 bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); 276 bool doIsophotal = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL"); 277 bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); 278 bool doKron = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON"); 279 280 psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER"); 281 psVector *radialBinsUpper = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER"); 282 assert (radialBinsLower->n == radialBinsUpper->n); 283 275 284 // we write out all sources, regardless of quality. the source flags tell us the state 276 for (i = 0; i < sources->n; i++) {285 for (int i = 0; i < sources->n; i++) { 277 286 // skip source if it is not a ext sourc 278 287 // XXX we have two places that extended source parameters are measured: … … 283 292 pmSource *source = sources->data[i]; 284 293 285 // choose the convolved EXT model, if available, otherwise the simple one 286 pmModel *model = source->modelConv; 287 if (model == NULL) { 288 model = source->modelEXT; 289 } 290 if (model == NULL) continue; 291 292 // XXX do we need to do this? 294 // skip sources without measurements 293 295 if (source->extpars == NULL) continue; 296 297 // we require a PSF model fit (ignore the real crud) 298 pmModel *model = source->modelPSF; 299 if (model == NULL) continue; 294 300 295 301 // XXX I need to split the extended models from the extended aperture measurements … … 301 307 yErr = dPAR[PM_PAR_YPOS]; 302 308 303 // XXX for the aperture values, I probably can / should just refer to the psf position and not write any of this304 axes = pmPSF_ModelToAxes (PAR, 20.0);305 306 309 row = psMetadataAlloc (); 310 307 311 // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) 308 312 psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq); … … 311 315 psMetadataAdd (row, PS_LIST_TAIL, "X_EXT_SIG", PS_DATA_F32, "Sigma in EXT x coordinate", xErr); 312 316 psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG", PS_DATA_F32, "Sigma in EXT y coordinate", yErr); 313 psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG", PS_DATA_F32, "EXT fit instrumental magnitude", source->extMag); 314 psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->errMag); 315 316 // XXX these should be major and minor, not 'x' and 'y' 317 psMetadataAdd (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", PS_DATA_F32, "EXT width in x coordinate", axes.major); 318 psMetadataAdd (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", PS_DATA_F32, "EXT width in y coordinate", axes.minor); 319 psMetadataAdd (row, PS_LIST_TAIL, "EXT_THETA", PS_DATA_F32, "EXT orientation angle", axes.theta); 320 321 // XXX at the moment, this will be a fixed sized table, but perhaps have multiple output formats depending on what is measured? 322 323 // other values that I need to report: 324 // R, Mag Petrosian, ..... 325 if (source->extpars) { 326 327 // Petrosian measurements 317 318 // Petrosian measurements 319 // XXX insert header data: petrosian ref radius, flux ratio 320 if (doPetrosian) { 328 321 pmSourcePetrosianValues *petrosian = source->extpars->petrosian; 329 322 if (petrosian) { … … 332 325 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius", petrosian->rad); 333 326 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error", petrosian->radErr); 327 } else { 328 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG", PS_DATA_F32, "Petrosian Magnitude", NAN); 329 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR", PS_DATA_F32, "Petrosian Magnitude Error", NAN); 330 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS", PS_DATA_F32, "Petrosian Radius", NAN); 331 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error", NAN); 334 332 } 335 336 // Kron measurements 333 } 334 335 // Kron measurements 336 if (doKron) { 337 337 pmSourceKronValues *kron = source->extpars->kron; 338 338 if (kron) { … … 341 341 psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS", PS_DATA_F32, "Kron Radius", kron->rad); 342 342 psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error", kron->radErr); 343 } else { 344 psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG", PS_DATA_F32, "Kron Magnitude", NAN); 345 psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR", PS_DATA_F32, "Kron Magnitude Error", NAN); 346 psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS", PS_DATA_F32, "Kron Radius", NAN); 347 psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error", NAN); 343 348 } 344 345 // Isophot measurements 346 // XXX insert header data: isophotal level 349 } 350 351 // Isophot measurements 352 // XXX insert header data: isophotal level 353 if (doIsophotal) { 347 354 pmSourceIsophotalValues *isophot = source->extpars->isophot; 348 355 if (isophot) { … … 351 358 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS", PS_DATA_F32, "Isophot Radius", isophot->rad); 352 359 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error", isophot->radErr); 360 } else { 361 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG", PS_DATA_F32, "Isophot Magnitude", NAN); 362 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR", PS_DATA_F32, "Isophot Magnitude Error", NAN); 363 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS", PS_DATA_F32, "Isophot Radius", NAN); 364 psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error", NAN); 353 365 } 354 366 } 367 368 // Flux Annuli 369 if (doAnnuli) { 355 370 pmSourceAnnuli *annuli = source->extpars->annuli; 356 371 if (annuli) { … … 367 382 sprintf (name, "FLUX_VAR_R_%02d", j); 368 383 psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", fluxVar->data.F32[j]); 369 } 384 } 385 } else { 386 for (int j = 0; j < radialBinsLower->n; j++) { 387 char name[32]; 388 sprintf (name, "FLUX_VAL_R_%02d", j); 389 psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", NAN); 390 sprintf (name, "FLUX_ERR_R_%02d", j); 391 psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", NAN); 392 sprintf (name, "FLUX_VAR_R_%02d", j); 393 psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", NAN); 394 } 370 395 } 371 396 } 372 397 373 374 398 psArrayAdd (table, 100, row); 399 psFree (row); 375 400 } 376 401 377 402 if (table->n == 0) { 378 psFitsWriteBlank (fits, outhead, extname); 379 psFree (table); 380 return true; 403 psFitsWriteBlank (fits, outhead, extname); 404 psFree (outhead); 405 psFree (table); 406 return true; 381 407 } 382 408 383 409 psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); 384 410 if (!psFitsWriteTable (fits, outhead, table, extname)) { 385 psError(PS_ERR_IO, false, "writing ext data %s\n", extname); 386 psFree(table); 387 return false; 388 } 411 psError(PS_ERR_IO, false, "writing ext data %s\n", extname); 412 psFree (outhead); 413 psFree(table); 414 return false; 415 } 416 psFree (outhead); 389 417 psFree (table); 390 418 … … 397 425 psArray *table; 398 426 psMetadata *row; 399 int i;400 427 psF32 *PAR, *dPAR; 401 428 psEllipseAxes axes; 402 429 psF32 xPos, yPos; 403 430 psF32 xErr, yErr; 431 char name[64]; 404 432 405 433 // create a header to hold the output data … … 412 440 sources = psArraySort (sources, pmSourceSortBySN); 413 441 442 // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams 443 int nParamMax = 0; 444 for (int i = 0; i < sources->n; i++) { 445 pmSource *source = sources->data[i]; 446 if (source->modelFits == NULL) continue; 447 for (int j = 0; j < source->modelFits->n; j++) { 448 pmModel *model = source->modelFits->data[j]; 449 assert (model); 450 nParamMax = PS_MAX (nParamMax, model->params->n); 451 } 452 } 453 414 454 table = psArrayAllocEmpty (sources->n); 415 455 416 456 // we write out all sources, regardless of quality. the source flags tell us the state 417 for (i = 0; i < sources->n; i++) { 418 // skip source if it is not a ext sourc 419 // XXX we have two places that extended source parameters are measured: 420 // psphotExtendedSources, which measures the aperture-like parameters and (potentially) the psf-convolved extended source models, 421 // psphotFitEXT, which does the simple extended source model fit (not psf-convolved) 422 // should we require both? 457 for (int i = 0; i < sources->n; i++) { 423 458 424 459 pmSource *source = sources->data[i]; 425 460 426 // XXX need to have an array of model fits; loop over the array and write out all 427 // which we calculated 428 429 // choose the convolved EXT model, if available, otherwise the simple one 430 // XXX should not need to choose: write both out 431 pmModel *model = source->modelConv; 432 if (model == NULL) { 433 model = source->modelEXT; 461 // XXX if no model fits are saved, write out modelEXT? 462 if (source->modelFits == NULL) continue; 463 464 // We have multiple sources : need to flag the one used to subtract the light (the 'best' model) 465 for (int j = 0; j < source->modelFits->n; j++) { 466 467 // choose the convolved EXT model, if available, otherwise the simple one 468 pmModel *model = source->modelFits->data[j]; 469 assert (model); 470 471 PAR = model->params->data.F32; 472 dPAR = model->dparams->data.F32; 473 xPos = PAR[PM_PAR_XPOS]; 474 yPos = PAR[PM_PAR_YPOS]; 475 xErr = dPAR[PM_PAR_XPOS]; 476 yErr = dPAR[PM_PAR_YPOS]; 477 478 axes = pmPSF_ModelToAxes (PAR, 20.0); 479 480 row = psMetadataAlloc (); 481 482 // XXX we are not writing out the mode (flags) or the type (psf, ext, etc) 483 psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET", 0, "IPP detection identifier index", source->seq); 484 psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT", 0, "EXT model x coordinate", xPos); 485 psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT", 0, "EXT model y coordinate", yPos); 486 psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT_SIG", 0, "Sigma in EXT x coordinate", xErr); 487 psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT_SIG", 0, "Sigma in EXT y coordinate", yErr); 488 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG", 0, "EXT fit instrumental magnitude", model->mag); 489 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", 0, "Sigma of PSF instrumental magnitude", model->magErr); 490 491 psMetadataAddF32 (row, PS_LIST_TAIL, "NPARAMS", 0, "number of model parameters", model->params->n); 492 psMetadataAddStr (row, PS_LIST_TAIL, "MODEL_TYPE", 0, "name of model", pmModelClassGetName (model->type)); 493 494 // XXX these should be major and minor, not 'x' and 'y' 495 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", 0, "EXT width in x coordinate", axes.major); 496 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", 0, "EXT width in y coordinate", axes.minor); 497 psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA", 0, "EXT orientation angle", axes.theta); 498 499 // write out the other generic parameters 500 for (int k = 0; k < nParamMax; k++) { 501 if (k == PM_PAR_I0) continue; 502 if (k == PM_PAR_SKY) continue; 503 if (k == PM_PAR_XPOS) continue; 504 if (k == PM_PAR_YPOS) continue; 505 if (k == PM_PAR_SXX) continue; 506 if (k == PM_PAR_SXY) continue; 507 if (k == PM_PAR_SYY) continue; 508 509 snprintf (name, 64, "EXT_PAR_%02d", k); 510 511 if (k < model->params->n) { 512 psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "", model->params->data.F32[k]); 513 } else { 514 psMetadataAddF32 (row, PS_LIST_TAIL, name, PS_DATA_F32, "", NAN); 515 } 516 } 517 518 // XXX other parameters which may be set. 519 // XXX flag / value to define the model 520 // XXX write out the model type, fit status flags 521 522 psArrayAdd (table, 100, row); 523 psFree (row); 434 524 } 435 if (model == NULL) continue;436 437 PAR = model->params->data.F32;438 dPAR = model->dparams->data.F32;439 xPos = PAR[PM_PAR_XPOS];440 yPos = PAR[PM_PAR_YPOS];441 xErr = dPAR[PM_PAR_XPOS];442 yErr = dPAR[PM_PAR_YPOS];443 444 axes = pmPSF_ModelToAxes (PAR, 20.0);445 446 row = psMetadataAlloc ();447 // XXX we are not writing out the mode (flags) or the type (psf, ext, etc)448 psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET", PS_DATA_U32, "IPP detection identifier index", source->seq);449 psMetadataAdd (row, PS_LIST_TAIL, "X_EXT", PS_DATA_F32, "EXT model x coordinate", xPos);450 psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT", PS_DATA_F32, "EXT model y coordinate", yPos);451 psMetadataAdd (row, PS_LIST_TAIL, "X_EXT_SIG", PS_DATA_F32, "Sigma in EXT x coordinate", xErr);452 psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG", PS_DATA_F32, "Sigma in EXT y coordinate", yErr);453 psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG", PS_DATA_F32, "EXT fit instrumental magnitude", source->extMag);454 psMetadataAdd (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude", source->errMag);455 456 // XXX these should be major and minor, not 'x' and 'y'457 psMetadataAdd (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ", PS_DATA_F32, "EXT width in x coordinate", axes.major);458 psMetadataAdd (row, PS_LIST_TAIL, "EXT_WIDTH_MIN", PS_DATA_F32, "EXT width in y coordinate", axes.minor);459 psMetadataAdd (row, PS_LIST_TAIL, "EXT_THETA", PS_DATA_F32, "EXT orientation angle", axes.theta);460 461 // XXX other parameters which may be set.462 // XXX flag / value to define the model463 464 psArrayAdd (table, 100, row);465 psFree (row);466 525 } 467 526 468 527 if (table->n == 0) { 469 psFitsWriteBlank (fits, outhead, extname); 470 psFree (table); 471 return true; 528 psFitsWriteBlank (fits, outhead, extname); 529 psFree (outhead); 530 psFree (table); 531 return true; 472 532 } 473 533 474 534 psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname); 475 535 if (!psFitsWriteTable (fits, outhead, table, extname)) { 476 psError(PS_ERR_IO, false, "writing ext data %s\n", extname); 477 psFree(table); 478 return false; 479 } 536 psError(PS_ERR_IO, false, "writing ext data %s\n", extname); 537 psFree (outhead); 538 psFree(table); 539 return false; 540 } 541 psFree (outhead); 480 542 psFree (table); 481 482 543 return true; 483 544 } -
trunk/psModules/src/objects/pmSourcePhotometry.c
r17287 r17396 3 3 * @author EAM, IfA; GLG, MHPCC 4 4 * 5 * @version $Revision: 1.4 0$ $Name: not supported by cvs2svn $6 * @date $Date: 2008-04-0 2 22:40:36$5 * @version $Revision: 1.41 $ $Name: not supported by cvs2svn $ 6 * @date $Date: 2008-04-08 18:35:38 $ 7 7 * 8 8 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 119 119 } 120 120 121 // measure EXT model photometry (do both modelEXT and modelConv or just the one?) 122 status = pmSourcePhotometryModel (&source->extMag, source->modelEXT); 121 // if we have a collection of model fits, one of them is a pointer to modelEXT? 122 // XXX not guaranteed 123 if (source->modelFits) { 124 for (int i = 0; i < source->modelFits->n; i++) { 125 pmModel *model = source->modelFits->data[i]; 126 status = pmSourcePhotometryModel (&model->mag, model); 127 } 128 if (source->modelEXT) { 129 source->extMag = source->modelEXT->mag; 130 } 131 } else { 132 if (source->modelEXT) { 133 status = pmSourcePhotometryModel (&source->extMag, source->modelEXT); 134 } 135 } 123 136 124 137 // for PSFs, correct both apMag and psfMag to same system, consistent with infinite flux star in aperture RADIUS -
trunk/psphot/doc/notes.txt
r14655 r17396 1 2 2008.01.31 3 4 Most of the fixes listed below have been finished except for: 5 6 * careful handling of the r ~ 0 fluxes / effective positions (for eg sersic) 7 * the EXT output parameters have not been tested 8 9 In addition, some further work to be done: 10 11 * multiple image detection 12 * multiple image fitting (including projection) 13 * optimization 14 * multi-threading 15 16 * possibly try the 'active deblending' concept that Ken from WISE is using. 17 * instead of fitting a double star with the positions floating at a 18 slightly random location, search for significant additions within 19 a grid about the central peak (say ~ 1sigma / 0.2 pix?). 20 (note that this is a linear fit) 21 22 * for the blend groups, freeze the positions for the fainter sources 23 (consistent with the non-linear S/N threshold) 24 25 1 26 2 27 2007.08.17 -
trunk/psphot/src/Makefile.am
r16820 r17396 18 18 19 19 libpsphot_la_SOURCES = \ 20 psphotErrorCodes.c \ 21 pmFootprint.c \ 22 psphotVersion.c \ 23 psphotModelGroupInit.c \ 24 psphotMaskReadout.c \ 25 psphotDefineFiles.c \ 26 psphotReadout.c \ 27 psphotModelBackground.c \ 28 psphotSubtractBackground.c \ 29 psphotFindDetections.c \ 30 psphotFindPeaks.c \ 31 psphotSourceStats.c \ 32 psphotRoughClass.c \ 33 psphotBasicDeblend.c \ 34 psphotChoosePSF.c \ 35 psphotGuessModels.c \ 36 psphotFitSourcesLinear.c \ 37 psphotBlendFit.c \ 38 psphotReplaceUnfit.c \ 39 psphotApResid.c \ 40 psphotMagnitudes.c \ 41 psphotSkyReplace.c \ 42 psphotEvalPSF.c \ 43 psphotEvalFLT.c \ 44 psphotSourceFits.c \ 45 psphotRadiusChecks.c \ 46 psphotOutput.c \ 47 psphotGrowthCurve.c \ 48 psphotFakeSources.c \ 49 psphotModelWithPSF.c \ 50 psphotExtendedSources.c \ 51 psphotRadialProfile.c \ 52 psphotPetrosian.c \ 53 psphotIsophotal.c \ 54 psphotAnnuli.c \ 55 psphotKron.c \ 56 psphotKernelFromPSF.c \ 57 psphotPSFConvModel.c \ 58 psphotModelTest.c \ 59 psphotFitSet.c \ 60 psphotSourceFreePixels.c \ 61 psphotSummaryPlots.c \ 62 psphotMergeSources.c \ 63 psphotLoadPSF.c \ 64 psphotReadoutCleanup.c \ 65 psphotSourcePlots.c \ 66 psphotRadialPlot.c \ 67 psphotDeblendSatstars.c \ 68 psphotMosaicSubimage.c \ 69 psphotMakeResiduals.c \ 70 psphotSourceSize.c \ 71 psphotDiagnosticPlots.c \ 72 psphotMakeFluxScale.c \ 73 psphotCheckStarDistribution.c \ 20 psphotErrorCodes.c \ 21 pmFootprint.c \ 22 psphotVersion.c \ 23 psphotModelGroupInit.c \ 24 psphotMaskReadout.c \ 25 psphotDefineFiles.c \ 26 psphotReadout.c \ 27 psphotModelBackground.c \ 28 psphotSubtractBackground.c \ 29 psphotFindDetections.c \ 30 psphotFindPeaks.c \ 31 psphotSourceStats.c \ 32 psphotRoughClass.c \ 33 psphotBasicDeblend.c \ 34 psphotChoosePSF.c \ 35 psphotGuessModels.c \ 36 psphotFitSourcesLinear.c \ 37 psphotBlendFit.c \ 38 psphotReplaceUnfit.c \ 39 psphotApResid.c \ 40 psphotMagnitudes.c \ 41 psphotSkyReplace.c \ 42 psphotEvalPSF.c \ 43 psphotEvalFLT.c \ 44 psphotSourceFits.c \ 45 psphotRadiusChecks.c \ 46 psphotOutput.c \ 47 psphotGrowthCurve.c \ 48 psphotFakeSources.c \ 49 psphotModelWithPSF.c \ 50 psphotExtendedSourceAnalysis.c \ 51 psphotExtendedSourceFits.c \ 52 psphotRadialProfile.c \ 53 psphotPetrosian.c \ 54 psphotIsophotal.c \ 55 psphotAnnuli.c \ 56 psphotKron.c \ 57 psphotKernelFromPSF.c \ 58 psphotPSFConvModel.c \ 59 psphotModelTest.c \ 60 psphotFitSet.c \ 61 psphotSourceFreePixels.c \ 62 psphotSummaryPlots.c \ 63 psphotMergeSources.c \ 64 psphotLoadPSF.c \ 65 psphotReadoutCleanup.c \ 66 psphotSourcePlots.c \ 67 psphotRadialPlot.c \ 68 psphotDeblendSatstars.c \ 69 psphotMosaicSubimage.c \ 70 psphotMakeResiduals.c \ 71 psphotSourceSize.c \ 72 psphotDiagnosticPlots.c \ 73 psphotMakeFluxScale.c \ 74 psphotCheckStarDistribution.c \ 74 75 psphotAddNoise.c 75 76 -
trunk/psphot/src/psphot.h
r16820 r17396 45 45 bool psphotMagnitudes (psArray *sources, psMetadata *recipe, pmPSF *psf, pmReadout *background); 46 46 bool psphotSkyReplace (pmConfig *config, const pmFPAview *view); 47 bool psphotExtendedSources (pmReadout *readout, psArray *sources, psMetadata *recipe); 47 bool psphotExtendedSourceAnalysis (pmReadout *readout, psArray *sources, psMetadata *recipe); 48 bool psphotExtendedSourceFits (pmReadout *readout, psArray *sources, psMetadata *recipe); 48 49 49 50 // used by psphotFindDetections … … 87 88 bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf, psMaskType maskVal); 88 89 bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf, psMaskType maskVal); 89 pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, p sMaskType maskVal);90 pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, pmModelType modelType, psMaskType maskVal); 90 91 psArray *psphotFitDBL (pmReadout *readout, pmSource *source, psMaskType maskVal); 91 92 … … 117 118 bool psphotMakeResiduals (psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal); 118 119 119 bool psphotPSFConvModel (pmSource *source, psMetadata *recipe, psMaskType maskVal); 120 pmModel *psphotPSFConvModel (pmReadout *readout, pmSource *source, pmModelType modelType, psMaskType maskVal, int psfSize); 121 120 122 psKernel *psphotKernelFromPSF (pmSource *source, int nPix); 121 123 -
trunk/psphot/src/psphotAnnuli.c
r15562 r17396 1 # include "psphot .h"1 # include "psphotInternal.h" 2 2 3 3 bool psphotAnnuli (pmSource *source, psMetadata *recipe, psMaskType maskVal) { … … 55 55 source->extpars->annuli->fluxVar = fluxSquare; 56 56 57 psFree (pixelCount); 58 57 59 return true; 58 60 } -
trunk/psphot/src/psphotBlendFit.c
r16820 r17396 40 40 // skip non-astronomical objects (very likely defects) 41 41 if (source->mode & PM_SOURCE_MODE_BLEND) continue; 42 if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue; 42 43 if (source->type == PM_SOURCE_TYPE_DEFECT) continue; 43 44 if (source->type == PM_SOURCE_TYPE_SATURATED) continue; … … 79 80 // XXX re-consider conditions under which the source has EXT fit: 80 81 // I should try EXT if the source size measurement says it is large 81 if (psphotFitBlend (readout, source, psf, maskVal)) { 82 psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->moments->x, source->moments->y); 83 Npsf ++; 84 continue; 85 } 86 if (psphotFitBlob (readout, source, sources, psf, maskVal)) { 87 psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->moments->x, source->moments->y); 88 Next ++; 89 continue; 90 } 82 if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { 83 if (psphotFitBlob (readout, source, sources, psf, maskVal)) { 84 source->type = PM_SOURCE_TYPE_EXTENDED; 85 psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->moments->x, source->moments->y); 86 Next ++; 87 continue; 88 } 89 } else { 90 if (psphotFitBlend (readout, source, psf, maskVal)) { 91 source->type = PM_SOURCE_TYPE_STAR; 92 psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->moments->x, source->moments->y); 93 Npsf ++; 94 continue; 95 } 96 } 91 97 92 98 psTrace ("psphot", 5, "source at %7.1f, %7.1f failed", source->moments->x, source->moments->y); -
trunk/psphot/src/psphotEvalFLT.c
r15060 r17396 8 8 if (model == NULL) { 9 9 source->mode &= ~PM_SOURCE_MODE_FITTED; 10 psTrace ("psphot", 5, "no model fitted?\n"); 10 11 return false; 11 12 } … … 14 15 if (!(model->flags & PM_MODEL_STATUS_FITTED)) { 15 16 source->mode &= ~PM_SOURCE_MODE_FITTED; 17 psTrace ("psphot", 5, "no model fitted?\n"); 16 18 return false; 17 19 } 20 18 21 // did the model fit fail for one or another reason? 19 22 if (model->flags & (PM_MODEL_STATUS_BADARGS | … … 22 25 source->mode |= PM_SOURCE_MODE_FAIL; 23 26 psLogMsg ("psphot", 5, "EXT fail fit\n"); 27 psTrace ("psphot", 5, "EXT fail fit\n"); 28 29 if (model->flags & PM_MODEL_STATUS_OFFIMAGE) { 30 psTrace ("psphot", 5, "off image\n"); 31 } 32 if (model->flags & PM_MODEL_STATUS_BADARGS) { 33 psTrace ("psphot", 5, "bad args\n"); 34 } 35 if (model->flags & PM_MODEL_STATUS_NONCONVERGE) { 36 psTrace ("psphot", 5, "non converge\n"); 37 } 24 38 return false; 25 39 } … … 35 49 if (model->params->data.F32[PM_PAR_I0] <= 0.02) { 36 50 source->mode |= PM_SOURCE_MODE_FAIL; 51 psTrace ("psphot", 5, "model central intensity ~ zero\n"); 37 52 return false; 38 53 } … … 43 58 // poor-quality fit; only keep if nothing else works... 44 59 psLogMsg ("psphot", 5, "EXT poor fit\n"); 60 psTrace ("psphot", 5, "EXT poor fit\n"); 61 45 62 source->mode |= PM_SOURCE_MODE_POOR; 46 63 return false; -
trunk/psphot/src/psphotFindDetections.c
r16820 r17396 8 8 bool status; 9 9 int pass; 10 11 psTimerStart ("psphot"); 10 12 11 13 // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) … … 67 69 psphotCullPeaks(readout->image, readout->weight, recipe, detections->footprints); 68 70 detections->peaks = pmFootprintArrayToPeaks(detections->footprints); 69 psLogMsg ("psphot", PS_LOG_INFO, " peaks: %ld, total footprints: %ld\n", detections->peaks->n, detections->footprints->n);71 psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks, %ld total footprints: %f sec\n", detections->peaks->n, detections->footprints->n, psTimerMark ("psphot")); 70 72 71 73 return detections; -
trunk/psphot/src/psphotFindPeaks.c
r16820 r17396 11 11 12 12 // smooth the image and weight map 13 psTimerStart ("p sphot");13 psTimerStart ("peaks"); 14 14 15 15 // XXX if we have been supplied a PSF, we can use that to set the FWHM_X,FWHM_Y values … … 32 32 psImage *smooth_im = psImageCopy (NULL, readout->image, PS_TYPE_F32); 33 33 psImageSmoothMaskF32 (smooth_im, readout->mask, maskVal, SIGMA_SMTH, NSIGMA_SMTH); 34 psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth image: %f sec\n", psTimerMark ("p sphot"));34 psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth image: %f sec\n", psTimerMark ("peaks")); 35 35 36 36 // smooth the weight, applying the mask as we go 37 37 psImage *smooth_wt = psImageCopy (NULL, readout->weight, PS_TYPE_F32); 38 38 psImageSmoothMaskF32 (smooth_wt, readout->mask, maskVal, SIGMA_SMTH/sqrt(2), NSIGMA_SMTH); 39 psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth weight: %f sec\n", psTimerMark ("p sphot"));39 psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth weight: %f sec\n", psTimerMark ("peaks")); 40 40 41 41 psImage *mask = readout->mask; … … 61 61 } 62 62 } 63 psLogMsg ("psphot", PS_LOG_INFO, "built smoothed signficance image: %f sec\n", psTimerMark ("p sphot"));63 psLogMsg ("psphot", PS_LOG_INFO, "built smoothed signficance image: %f sec\n", psTimerMark ("peaks")); 64 64 65 65 // optionally save example images under trace … … 70 70 } 71 71 72 psTimerStart ("p sphot");72 psTimerStart ("peaks"); 73 73 // set peak threshold 74 74 … … 140 140 pmPeaksWriteText (peaks, output); 141 141 } 142 psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks: %f sec\n", peaks->n, psTimerMark ("p sphot"));142 psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks: %f sec\n", peaks->n, psTimerMark ("peaks")); 143 143 144 144 // If they asked us to return a set of pmFootprints, not a raw pmPeak list, … … 167 167 peaks = footprints; // well, you know what I mean 168 168 169 psLogMsg ("psphot", PS_LOG_INFO, "%ld footprints: %f sec\n", peaks->n, psTimerMark ("p sphot"));169 psLogMsg ("psphot", PS_LOG_INFO, "%ld footprints: %f sec\n", peaks->n, psTimerMark ("peaks")); 170 170 } 171 171 -
trunk/psphot/src/psphotIsophotal.c
r15562 r17396 1 # include "psphot.h" 1 # include "psphotInternal.h" 2 3 static float ISOPHOT_FLUX = NAN; 2 4 3 5 bool psphotIsophotal (pmSource *source, psMetadata *recipe, psMaskType maskVal) { … … 14 16 15 17 // flux at which to measure isophotal parameters 16 // XXX cache this? 17 float ISOPHOT_FLUX = psMetadataLookupF32 (&status, recipe, "ISOPHOTAL_FLUX"); 18 if (!isfinite(ISOPHOT_FLUX)) { 19 // XXX ISOPHOTAL_FLUX should be specified in mags, need the zero point to get counts/sec 20 ISOPHOT_FLUX = psMetadataLookupF32 (&status, recipe, "ISOPHOTAL_FLUX"); 21 assert (status); 22 } 18 23 19 24 // find the first bin below the flux level and the last above the level 20 // XXX can this be done faster with disection?25 // XXX can this be done faster with bisection? 21 26 // XXX do I need to worry about crazy outliers? 22 27 // XXX should i be smoothing or fitting the curve? … … 27 32 if ((firstBelow < 0) && (flux->data.F32[i] < ISOPHOT_FLUX)) firstBelow = i; 28 33 } 29 34 // if we don't go out far enough, we have a problem... 35 if (lastAbove == flux->n - 1) { 36 psTrace ("psphot", 5, "did not go out far enough to reach isophotal magnitude"); 37 // XXX raise a flag ? 38 return false; 39 } 40 if (firstBelow < 0) { 41 psTrace ("psphot", 5, "did not go out far enough to bound isophotal magnitude: error unmeasured"); 42 // XXX raise a flag ? 43 lastAbove = firstBelow; 44 return false; 45 } 46 30 47 // need to examine pixels in this vicinity 31 48 float isophotalFluxFirst = 0; … … 56 73 source->extpars->isophot->radErr = isophotalRadErr; 57 74 75 psTrace ("psphot", 5, "Isophot flux:%f +/- %f @ %f +/- %f for %f, %f\n", 76 source->extpars->isophot->mag, source->extpars->isophot->magErr, 77 source->extpars->isophot->rad, source->extpars->isophot->radErr, 78 source->peak->xf, source->peak->yf); 79 58 80 return true; 59 81 -
trunk/psphot/src/psphotKernelFromPSF.c
r14338 r17396 1 # include "psphot .h"1 # include "psphotInternal.h" 2 2 3 3 psKernel *psphotKernelFromPSF (pmSource *source, int nPix) { 4 4 5 assert (source);6 assert (source->psfFlux); // XXX build if needed?5 assert (source); 6 assert (source->psfFlux); // XXX build if needed? 7 7 8 int x0 = source->peak->xf - source->psfFlux->col0;9 int y0 = source->peak->yf - source->psfFlux->row0;8 int x0 = source->peak->xf - source->psfFlux->col0; 9 int y0 = source->peak->yf - source->psfFlux->row0; 10 10 11 // need to decide on the size: dynamically? statically?12 psKernel *psf = psKernelAlloc (-nPix, +nPix, -nPix, +nPix);11 // need to decide on the size: dynamically? statically? 12 psKernel *psf = psKernelAlloc (-nPix, +nPix, -nPix, +nPix); 13 13 14 for (int j = psf->yMin; j <= psf->yMax; j++) { 15 for (int i = psf->xMin; i <= psf->xMax; i++) { 16 psf->kernel[j][i] = source->psfFlux->data.F32[y0 + j][x0 + i]; 14 // XXX we should just re-construct a PSF at this location 15 // psModelAdd (psf->image, NULL, source->modelPSF, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM | PM_MODEL_OP_CENTER); 16 17 // if the realized PSF for this object does not cover the full kernel, give up for now 18 if (x0 + psf->xMin < 0) goto escape; 19 if (x0 + psf->xMax >= source->psfFlux->numCols) goto escape; 20 if (y0 + psf->yMin < 0) goto escape; 21 if (y0 + psf->yMax >= source->psfFlux->numRows) goto escape; 22 23 double sum = 0.0; 24 for (int j = psf->yMin; j <= psf->yMax; j++) { 25 for (int i = psf->xMin; i <= psf->xMax; i++) { 26 double value = source->psfFlux->data.F32[y0 + j][x0 + i]; 27 psf->kernel[j][i] = value; 28 sum += value; 29 } 17 30 } 18 } 19 return psf; 31 assert (sum > 0.0); 32 33 // psf must be normalized (integral = 1.0) 34 for (int i = 0; i < psf->image->numRows; i++) { 35 for (int j = 0; j < psf->image->numCols; j++) { 36 psf->image->data.F32[i][j] /= sum; 37 } 38 } 39 40 return psf; 41 42 escape: 43 psFree (psf); 44 return NULL; 20 45 } -
trunk/psphot/src/psphotKron.c
r13983 r17396 1 # include "psphot .h"1 # include "psphotInternal.h" 2 2 3 3 bool psphotKron (pmSource *source, psMetadata *recipe, psMaskType maskVal) { -
trunk/psphot/src/psphotModelTest.c
r16820 r17396 208 208 source->modelPSF = pmModelFromPSF (model, psf); 209 209 source->modelEXT = model; 210 status = psphotPSFConvModel (source, recipe, maskVal); 211 model = source->modelConv; 210 211 // what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box) 212 int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE"); 213 assert (status); 214 215 model = psphotPSFConvModel (readout, source, modelType, maskVal, psfSize); 212 216 params = model->params->data.F32; 213 217 } else { -
trunk/psphot/src/psphotModelWithPSF.c
r14347 r17396 1 # include "psphot.h" 1 # include "psphotInternal.h" 2 # define SAVE_IMAGES 0 2 3 3 4 bool psphotModelWithPSF_LMM ( … … 36 37 37 38 // allocate internal arrays (current vs Guess) 38 psImage *alpha = psImageAlloc(params->n, params->n, PS_TYPE_F32); 39 psImage *Alpha = psImageAlloc(params->n, params->n, PS_TYPE_F32); 40 psVector *beta = psVectorAlloc(params->n, PS_TYPE_F32); 41 psVector *Beta = psVectorAlloc(params->n, PS_TYPE_F32); 39 psImage *Alpha = NULL; 40 psVector *Beta = NULL; 41 42 // Alpha & Beta only contain elements to represent the unmasked parameters 43 if (!psMinLM_AllocAB (&Alpha, &Beta, params, paramMask)) { 44 psAbort ("programming error: no unmasked parameters to be fit\n"); 45 } 46 47 // allocate internal arrays (current vs Guess) 48 psImage *alpha = psImageAlloc(Alpha->numCols, Alpha->numRows, PS_TYPE_F32); 49 psVector *beta = psVectorAlloc(Beta->n, PS_TYPE_F32); 42 50 psVector *Params = psVectorAlloc(params->n, PS_TYPE_F32); 43 51 … … 79 87 // dump some useful info if trace is defined 80 88 if (psTraceGetLevel("psphot") >= 6) { 81 p_psImagePrint(psTraceGetDestination(), Alpha, "alpha guess (1)"); 82 p_psVectorPrint(psTraceGetDestination(), Beta, "beta guess (1)"); 89 p_psImagePrint(psTraceGetDestination(), Alpha, "Alpha guess (1)"); 90 p_psVectorPrint(psTraceGetDestination(), Beta, "Beta guess (1)"); 91 p_psVectorPrint(psTraceGetDestination(), beta, "beta current (1)"); 83 92 } 84 93 if (psTraceGetLevel("psphot") >= 5) { … … 116 125 beta = psVectorCopy(beta, Beta, PS_TYPE_F32); 117 126 params = psVectorCopy(params, Params, PS_TYPE_F32); 118 lambda *= 0. 1;127 lambda *= 0.25; 119 128 120 129 // save the new convolved model image … … 130 139 // construct & return the covariance matrix (if requested) 131 140 if (covar != NULL) { 132 if (!psMinLM_GuessABP( covar, Beta, Params, alpha, beta, params, paramMask, NULL, 0.0, NULL)) {141 if (!psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask, NULL, 0.0, NULL)) { 133 142 psTrace ("psphot", 5, "failure to calculate covariance matrix\n"); 134 143 } 144 // set covar values which are not masked 145 psImageInit (covar, 0.0); 146 for (int j = 0, J = 0; j < params->n; j++) { 147 if (paramMask && (paramMask->data.U8[j])) { 148 covar->data.F32[j][j] = 1.0; 149 continue; 150 } 151 for (int k = 0, K = 0; k < params->n; k++) { 152 if (paramMask && (paramMask->data.U8[k])) continue; 153 covar->data.F32[j][k] = Alpha->data.F32[J][K]; 154 K++; 155 } 156 J++; 157 } 135 158 } 136 159 … … 177 200 } 178 201 179 // generate the model and derivative images for this parameter set202 // 1 *** generate the model and derivative images for this parameter set 180 203 181 204 // storage for model derivatives … … 223 246 224 247 for (int n = 0; n < params->n; n++) { 225 226 227 248 if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; } 249 psImage *dmodel = pcm->dmodels->data[n]; 250 dmodel->data.F32[i][j] = deriv->data.F32[n]; 228 251 } 229 252 } … … 231 254 psFree(coord); 232 255 psFree(deriv); 233 234 256 235 257 // convolve model and dmodel arrays with PSF 236 258 psImageConvolveDirect (pcm->modelConv, pcm->model, psf); 237 259 for (int n = 0; n < pcm->dmodels->n; n++) { 238 239 240 241 260 if (pcm->dmodels->data[n] == NULL) continue; 261 psImage *dmodel = pcm->dmodels->data[n]; 262 psImage *dmodelConv = pcm->dmodelsConv->data[n]; 263 psImageConvolveDirect (dmodelConv, dmodel, psf); 242 264 } 243 265 244 266 // XXX TEST : SAVE IMAGES 245 # if ( 1)267 # if (SAVE_IMAGES) 246 268 psphotSaveImage (NULL, psf->image, "psf.fits"); 247 269 psphotSaveImage (NULL, pcm->model, "model.fits"); … … 252 274 # endif 253 275 276 // 2 *** accumulate alpha & beta 277 254 278 // zero alpha and beta for summing below 255 279 psImageInit (alpha, 0.0); … … 259 283 for (psS32 i = 0; i < source->pixels->numRows; i++) { 260 284 for (psS32 j = 0; j < source->pixels->numCols; j++) { 261 // XXX are we doing the right thing with the mask?285 // XXX are we doing the right thing with the mask? 262 286 // skip masked points 263 287 if (source->maskObj->data.U8[i][j]) { … … 279 303 chisq += PS_SQR(delta) * yweight; 280 304 281 if (isnan(delta)) 282 psAbort("nan in delta"); 283 if (isnan(chisq)) 284 psAbort("nan in chisq"); 285 286 for (psS32 n1 = 0; n1 < params->n; n1++) { 287 if ((paramMask != NULL) && (paramMask->data.U8[n1])) { 288 continue; 289 } 290 psImage *dmodel = pcm->dmodelsConv->data[n1]; 291 float weight = dmodel->data.F32[i][j] * yweight; 292 for (psS32 n2 = 0; n2 <= n1; n2++) { 293 if ((paramMask != NULL) && (paramMask->data.U8[n2])) { 294 continue; 295 } 296 dmodel = pcm->dmodelsConv->data[n2]; 297 alpha->data.F32[n1][n2] += weight * dmodel->data.F32[i][j]; 298 } 299 beta->data.F32[n1] += weight * delta; 305 if (isnan(delta)) psAbort("nan in delta"); 306 if (isnan(chisq)) psAbort("nan in chisq"); 307 308 // alpha & beta only contain unmasked elements 309 for (int n1 = 0, N1 = 0; n1 < params->n; n1++) { 310 if ((paramMask != NULL) && (paramMask->data.U8[n1])) continue; 311 psImage *dmodel = pcm->dmodelsConv->data[n1]; 312 float weight = dmodel->data.F32[i][j] * yweight; 313 for (int n2 = 0, N2 = 0; n2 <= n1; n2++) { 314 if ((paramMask != NULL) && (paramMask->data.U8[n2])) continue; 315 dmodel = pcm->dmodelsConv->data[n2]; 316 alpha->data.F32[N1][N2] += weight * dmodel->data.F32[i][j]; 317 N2++; 318 } 319 beta->data.F32[N1] += weight * delta; 320 N1++; 300 321 } 301 322 } … … 303 324 304 325 // calculate lower-left half of alpha 305 for (psS32 j = 1; j < params->n; j++) {326 for (psS32 j = 1; j < alpha->numCols; j++) { 306 327 for (psS32 k = 0; k < j; k++) { 307 328 alpha->data.F32[k][j] = alpha->data.F32[j][k]; 308 }309 }310 311 // fill in pivots if we apply a mask312 if (paramMask != NULL) {313 for (psS32 j = 0; j < params->n; j++) {314 if (paramMask->data.U8[j]) {315 alpha->data.F32[j][j] = 1;316 beta->data.F32[j] = 1;317 }318 329 } 319 330 } … … 345 356 pcm->dmodels = psArrayAlloc (params->n); 346 357 for (psS32 n = 0; n < params->n; n++) { 347 if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; } 348 pcm->dmodels->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32); 358 pcm->dmodels->data[n] = NULL; 359 if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; } 360 pcm->dmodels->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32); 349 361 } 350 362 … … 353 365 pcm->dmodelsConv = psArrayAlloc (params->n); 354 366 for (psS32 n = 0; n < params->n; n++) { 355 if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; } 356 pcm->dmodelsConv->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32); 367 pcm->dmodelsConv->data[n] = NULL; 368 if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; } 369 pcm->dmodelsConv->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32); 357 370 } 358 371 … … 378 391 379 392 while () { 380 381 382 383 393 GuessABP (Alpha, Beta, Params, alpha, beta, params, paramMask, checkLimits, lambda) 394 dLinear = dLinear(Beta, beta, lambda); 395 chisq = SetABX (alpha, beta, params, paramMask, x, y, dy, func) 396 convergence tests... 384 397 } 385 398 … … 388 401 ** GuessABP: 389 402 390 391 392 df_c/dp_o = d/dp_o [sum_i (kern_i * func (x_i; p_o))]393 394 df_c/dp_o = sum_i (d/dp_o [kern_i * func (x_i; p_o)])395 396 df_c/dp_o = sum_i (kern_i * d/dp_o [func (x_i; p_o)])397 398 - generate image arrays for func, dfunc/dp_j (not masked)399 - convolve each with psf400 - measure delta = f_conv - data401 - etc403 f_c = sum_i (kern_i * func (x_i; p_o)) 404 405 df_c/dp_o = d/dp_o [sum_i (kern_i * func (x_i; p_o))] 406 407 df_c/dp_o = sum_i (d/dp_o [kern_i * func (x_i; p_o)]) 408 409 df_c/dp_o = sum_i (kern_i * d/dp_o [func (x_i; p_o)]) 410 411 - generate image arrays for func, dfunc/dp_j (not masked) 412 - convolve each with psf 413 - measure delta = f_conv - data 414 - etc 402 415 */ 403 416 -
trunk/psphot/src/psphotPSFConvModel.c
r14655 r17396 1 # include "psphot.h" 1 # include "psphotInternal.h" 2 # define USE_DELTA_PSF 0 2 3 3 4 // save as static values so they may be set externally 4 5 static psF32 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = 15; 5 6 static psF32 PM_SOURCE_FIT_MODEL_TOLERANCE = 0.1; 6 // static psF32 PM_SOURCE_FIT_MODEL_WEIGHT = 1.0;7 // static bool PM_SOURCE_FIT_MODEL_PIX_WEIGHTS = true;8 7 9 8 // input source has both modelPSF and modelEXT. on successful exit, we set the 10 9 // modelConv to contain the fitted parameters, and the modelFlux to contain the 11 10 // convolved model image. 12 bool psphotPSFConvModel (pmSource *source, psMetadata *recipe, psMaskType maskVal) {11 pmModel *psphotPSFConvModel (pmReadout *readout, pmSource *source, pmModelType modelType, psMaskType maskVal, int psfSize) { 13 12 14 bool status;15 16 int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE");17 if (!status) {18 psfSize = 2;19 }20 21 13 // make sure we save a cached copy of the psf flux 22 14 pmSourceCachePSF (source, maskVal); 23 15 24 16 // convert the cached cached psf model for this source to a psKernel 25 // XXX for the moment, hard-wire the kernel to be 5x5 (2 pix radius)26 // XXX for the moment, hard-wire the kernel to be 9x9 (4 pix radius)27 17 psKernel *psf = psphotKernelFromPSF (source, psfSize); 18 if (!psf) return NULL; 28 19 29 // psf must be normalized (integral = 1.0) 30 double sum = 0.0; 31 for (int i = 0; i < psf->image->numRows; i++) { 32 for (int j = 0; j < psf->image->numCols; j++) { 33 sum += psf->image->data.F32[i][j]; 34 } 35 } 36 assert (sum > 0.0); 37 for (int i = 0; i < psf->image->numRows; i++) { 38 for (int j = 0; j < psf->image->numCols; j++) { 39 psf->image->data.F32[i][j] /= sum; 40 } 41 } 42 43 # if (0) 44 // XXX sanity check: convolve with delta function should behave like unconvolved version 45 for (int i = 0; i < psf->image->numRows; i++) { 46 for (int j = 0; j < psf->image->numCols; j++) { 47 psf->image->data.F32[i][j] = 0.0; 48 } 49 } 50 psf->image->data.F32[2][2] = 1.0; 20 # if (USE_DELTA_PSF) 21 psImageInit (psf->image, 0.0); 22 psf->image->data.F32[(int)(0.5*psf->image->numRows)][(int)(0.5*psf->image->numCols)] = 1.0; 51 23 # endif 52 24 … … 54 26 // XXX we could modify the parameter values or even the model 55 27 // here based on the observed seeing (some lookup table...) 56 pmModel *modelConv = pmModelCopy (source->modelEXT); 28 29 // use the source moments, etc to guess basic model parameters 30 pmModel *modelConv = pmSourceModelGuess (source, modelType); 31 if (!modelConv) { 32 psFree (psf); 33 return NULL; 34 } 35 36 // XXX test : modify the Io, SXX, SYY terms based on the psf SXX, SYY terms: 37 psEllipseShape psfShape; 38 psfShape.sx = source->modelPSF->params->data.F32[PM_PAR_SXX] / M_SQRT2; 39 psfShape.sxy = source->modelPSF->params->data.F32[PM_PAR_SXY]; 40 psfShape.sy = source->modelPSF->params->data.F32[PM_PAR_SYY] / M_SQRT2; 41 psEllipseAxes psfAxes = psEllipseShapeToAxes (psfShape, 20.0); 42 43 // XXX test : modify the Io, SXX, SYY terms based on the psf SXX, SYY terms: 44 psEllipseShape extShape; 45 extShape.sx = modelConv->params->data.F32[PM_PAR_SXX] / M_SQRT2; 46 extShape.sxy = modelConv->params->data.F32[PM_PAR_SXY]; 47 extShape.sy = modelConv->params->data.F32[PM_PAR_SYY] / M_SQRT2; 48 psEllipseAxes extAxes = psEllipseShapeToAxes (extShape, 20.0); 49 50 // decrease the initial guess ellipse by psf_minor axis: 51 psEllipseAxes extAxesMod; 52 extAxesMod.major = sqrt (PS_MAX (1.0, PS_SQR(extAxes.major) - PS_SQR(psfAxes.minor))); 53 extAxesMod.minor = sqrt (PS_MAX (1.0, PS_SQR(extAxes.minor) - PS_SQR(psfAxes.minor))); 54 extAxesMod.theta = extAxes.theta; 55 56 psEllipseShape extShapeMod = psEllipseAxesToShape (extAxesMod); 57 modelConv->params->data.F32[PM_PAR_SXX] = extShapeMod.sx * M_SQRT2; 58 modelConv->params->data.F32[PM_PAR_SXY] = extShapeMod.sxy; 59 modelConv->params->data.F32[PM_PAR_SYY] = extShapeMod.sy * M_SQRT2; 60 61 // increase the initial guess central intensity by 2pi r^2: 62 modelConv->params->data.F32[PM_PAR_I0] *= (1.0 + PS_SQR(psfAxes.minor) / PS_SQR(extAxesMod.minor)); 63 57 64 psVector *params = modelConv->params; 58 65 psVector *dparams = modelConv->dparams; 66 67 psphotCheckRadiusEXT (readout, source, modelConv); 59 68 60 69 // create the minimization constraints … … 91 100 psTrace ("psphot", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value); 92 101 93 // renormalize output model image 102 // renormalize output model image (generated by fitting process) 94 103 float Io = params->data.F32[PM_PAR_I0]; 95 104 for (int iy = 0; iy < source->modelFlux->numRows; iy++) { … … 115 124 onPic &= (params->data.F32[PM_PAR_YPOS] >= source->pixels->row0); 116 125 onPic &= (params->data.F32[PM_PAR_YPOS] < source->pixels->row0 + source->pixels->numRows); 117 if (!onPic) { 118 modelConv->flags |= PM_MODEL_STATUS_OFFIMAGE; 119 } 126 if (!onPic) modelConv->flags |= PM_MODEL_STATUS_OFFIMAGE; 120 127 121 source->mode |= PM_SOURCE_MODE_FITTED; 122 source->modelConv = modelConv; 128 source->mode |= PM_SOURCE_MODE_FITTED; // XXX is this needed? 123 129 130 psFree(psf); 124 131 psFree(myMin); 125 132 psFree(covar); 126 133 psFree(constraint); 127 134 128 bool retval = (onPic && fitStatus); 129 psTrace("psphot", 5, "---- %s(%d) end ----\n", __func__, retval); 130 return(retval); 135 return modelConv; 131 136 } -
trunk/psphot/src/psphotPetrosian.c
r15800 r17396 1 # include "psphot.h" 1 # include "psphotInternal.h" 2 3 static float PETROSIAN_R0 = NAN; 4 static float PETROSIAN_RF = NAN; 2 5 3 6 bool psphotPetrosian (pmSource *source, psMetadata *recipe, psMaskType maskVal) { … … 14 17 15 18 // flux at which to measure isophotal parameters 16 // XXX cache this? 17 float PETROSIAN_R0 = psMetadataLookupF32 (&status, recipe, "PETROSIAN_R0"); 18 float PETROSIAN_RF = psMetadataLookupF32 (&status, recipe, "PETROSIAN_FLUX_RATIO"); 19 if (!isfinite(PETROSIAN_R0)) { 20 PETROSIAN_R0 = psMetadataLookupF32 (&status, recipe, "PETROSIAN_R0"); 21 PETROSIAN_RF = psMetadataLookupF32 (&status, recipe, "PETROSIAN_FLUX_RATIO"); 22 assert (status); 23 } 19 24 20 25 // first find flux at R0 21 int first Below= -1;22 int last Above= -1;26 int firstAbove = -1; 27 int lastBelow = -1; 23 28 for (int i = 0; i < radius->n; i++) { 24 if (radius->data.F32[i] > PETROSIAN_R0) lastAbove = i; 25 if ((firstBelow < 0) && (flux->data.F32[i] < PETROSIAN_R0)) firstBelow = i; 29 if (radius->data.F32[i] < PETROSIAN_R0) lastBelow = i; 30 if ((firstAbove < 0) && (radius->data.F32[i] > PETROSIAN_R0)) firstAbove = i; 31 } 32 // if we don't go out far enough, we have a problem... 33 if (lastBelow == radius->n - 1) { 34 psTrace ("psphot", 5, "did not go out far enough to reach petrosian reference radius..."); 35 // XXX skip object? raise a flag ? 36 return false; 37 } 38 if (firstAbove < 0) { 39 psTrace ("psphot", 5, "did not go out far enough to bound petrosian reference radius"); 40 // XXX raise a flag ? 41 return false; 26 42 } 27 43 … … 29 45 float fluxR0 = 0.0; 30 46 int fluxRn = 0; 31 for (int i = PS_MIN(first Below, lastAbove); i <= PS_MAX(firstBelow, lastAbove); i++) {47 for (int i = PS_MIN(firstAbove, lastBelow); i <= PS_MAX(firstAbove, lastBelow); i++) { 32 48 fluxR0 += flux->data.F32[i]; 33 49 fluxRn ++; … … 42 58 // XXX do I need to worry about crazy outliers? 43 59 // XXX should i be smoothing or fitting the curve? 44 firstBelow = -1;45 lastAbove = -1;60 int firstBelow = -1; 61 int lastAbove = -1; 46 62 for (int i = 0; i < flux->n; i++) { 47 63 if (flux->data.F32[i] > fluxRP) lastAbove = i; 48 64 if ((firstBelow < 0) && (flux->data.F32[i] < fluxRP)) firstBelow = i; 65 } 66 // if we don't go out far enough, we have a problem... 67 if (lastAbove == radius->n - 1) { 68 psTrace ("psphot", 5, "did not go out far enough to reach petrosian radius..."); 69 // XXX skip object? raise a flag ? 70 return false; 71 } 72 if (firstBelow < 0) { 73 psTrace ("psphot", 5, "did not go out far enough to bound petrosian radius"); 74 // XXX raise a flag ? 75 return false; 49 76 } 50 77 … … 78 105 source->extpars->petrosian->radErr = radErr; 79 106 107 psTrace ("psphot", 5, "Petrosian flux:%f +/- %f @ %f +/- %f for %f, %f\n", 108 source->extpars->petrosian->mag, source->extpars->petrosian->magErr, 109 source->extpars->petrosian->rad, source->extpars->petrosian->radErr, 110 source->peak->xf, source->peak->yf); 111 80 112 return true; 81 113 -
trunk/psphot/src/psphotRadialProfile.c
r15819 r17396 41 41 42 42 int n = 0; 43 float Xo = source->modelEXT->params->data.F32[PM_PAR_XPOS] - source->pixels->col0; 44 float Yo = source->modelEXT->params->data.F32[PM_PAR_YPOS] - source->pixels->row0; 43 44 float Xo = 0.0; 45 float Yo = 0.0; 46 47 if (source->modelEXT) { 48 Xo = source->modelEXT->params->data.F32[PM_PAR_XPOS] - source->pixels->col0; 49 Yo = source->modelEXT->params->data.F32[PM_PAR_YPOS] - source->pixels->row0; 50 } else { 51 Xo = source->peak->xf - source->pixels->col0; 52 Yo = source->peak->yf - source->pixels->row0; 53 } 45 54 for (int iy = 0; iy < source->pixels->numRows; iy++) { 46 55 for (int ix = 0; ix < source->pixels->numCols; ix++) { -
trunk/psphot/src/psphotReadout.c
r17112 r17396 117 117 psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE); 118 118 119 // XXX test the CR/EXT measurement 120 // XXX we need to call this here and after the second pass: option for starting number? 119 // identify CRs and extended sources 121 120 psphotSourceSize (readout, sources, recipe, 0); 122 121 if (!strcasecmp (breakPt, "ENSEMBLE")) { … … 127 126 psphotBlendFit (readout, sources, recipe, psf); 128 127 129 // replace all sources (make this part of psphotFitSourcesLinear?)128 // replace all sources 130 129 psphotReplaceAll (sources, recipe); 131 130 132 131 // linear fit to include all sources 133 132 psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE); 133 134 // if we only do one pass, skip to extended source analysis 134 135 if (!strcasecmp (breakPt, "PASS1")) { 135 goto finish;136 goto pass1finish; 136 137 } 137 138 … … 145 146 psphotAddNoise (readout, sources, recipe); 146 147 148 // find fainter sources (pass 2) 147 149 detections = psphotFindDetections (detections, readout, recipe); 148 150 149 // remove noise for subtracted objects 151 // remove noise for subtracted objects (ie, return to normal noise level) 150 152 psphotSubNoise (readout, sources, recipe); 151 153 … … 172 174 psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE); 173 175 176 pass1finish: 177 174 178 // measure source size for the remaining sources 175 179 psphotSourceSize (readout, sources, recipe, 0); 176 180 177 psphotExtendedSources (readout, sources, recipe); 181 psphotExtendedSourceAnalysis (readout, sources, recipe); 182 183 psphotExtendedSourceFits (readout, sources, recipe); 178 184 179 185 finish: -
trunk/psphot/src/psphotSourceFits.c
r16820 r17396 220 220 pmSource *tmpSrc = pmSourceAlloc (); 221 221 222 pmModel *EXT = psphotFitEXT (readout, source, m askVal);222 pmModel *EXT = psphotFitEXT (readout, source, modelTypeEXT, maskVal); 223 223 okEXT = psphotEvalEXT (tmpSrc, EXT); 224 224 chiEXT = EXT->chisq / EXT->nDOF; … … 267 267 // save new model 268 268 source->modelEXT = EXT; 269 source->type = PM_SOURCE_TYPE_EXTENDED; 270 source->mode |= PM_SOURCE_MODE_EXTMODEL; 269 271 270 272 // build cached model and subtract … … 349 351 } 350 352 351 pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, p sMaskType maskVal) {353 pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, pmModelType modelType, psMaskType maskVal) { 352 354 353 355 NfitEXT ++; 354 356 355 357 // use the source moments, etc to guess basic model parameters 356 pmModel *EXT = pmSourceModelGuess (source, modelType EXT);358 pmModel *EXT = pmSourceModelGuess (source, modelType); 357 359 PS_ASSERT (EXT, NULL); 358 360 -
trunk/psphot/src/psphotSourceSize.c
r17113 r17396 1 1 # include "psphotInternal.h" 2 2 # include <gsl/gsl_sf_gamma.h> 3 4 float psphotModelContour (psImage *image, psImage *weight, psImage *mask, pmModel *model, float Ro); 3 5 4 6 // we need to call this function after sources have been fitted to the PSF model and … … 14 16 psTimerStart ("psphot"); 15 17 16 float CR_NSIGMA_LIMIT = psMetadataLookupF32 (&status, recipe, "PSPHOT.CRNSIGMA.LIMIT"); 18 float CR_NSIGMA_LIMIT = psMetadataLookupF32 (&status, recipe, "PSPHOT.CR.NSIGMA.LIMIT"); 19 assert (status); 20 21 float EXT_NSIGMA_LIMIT = psMetadataLookupF32 (&status, recipe, "PSPHOT.EXT.NSIGMA.LIMIT"); 17 22 assert (status); 18 23 … … 24 29 if (isfinite(source->crNsigma)) continue; 25 30 26 source->crNsigma = -1.0;27 source->extNsigma = -1.0;28 29 31 // source must have been subtracted 30 source->crNsigma = -3.0;31 32 if (!(source->mode & PM_SOURCE_MODE_SUBTRACTED)) continue; 32 33 33 psF32 **resid = source->pixels->data.F32;34 psF32 **resid = source->pixels->data.F32; 34 35 psF32 **weight = source->weight->data.F32; 35 psU8 **mask = source->maskObj->data.U8; 36 psU8 **mask = source->maskObj->data.U8; 37 38 // check for extendedness: measure the delta flux significance at the 1 sigma contour 39 source->extNsigma = psphotModelContour (source->pixels, source->weight, source->maskObj, source->modelPSF, 1.0); 40 41 // XXX prevent a source from being both CR and EXT? 42 if (source->extNsigma > EXT_NSIGMA_LIMIT) { 43 source->mode |= PM_SOURCE_MODE_EXT_LIMIT; 44 } 36 45 37 46 int xPeak = source->peak->xf - source->pixels->col0 + 0.5; 38 47 int yPeak = source->peak->yf - source->pixels->row0 + 0.5; 39 48 40 // skip sources which are too close to a boundary41 source->crNsigma = -4.0;49 // XXX for now, skip sources which are too close to a boundary 50 // XXX raise a flag? 42 51 if (xPeak < 1) continue; 43 52 if (xPeak > source->pixels->numCols - 2) continue; … … 46 55 47 56 // XXX for now, just skip any sources with masked pixels 48 source->crNsigma = -5.0;57 // XXX raise a flag? 49 58 bool keep = true; 50 59 for (int iy = -1; (iy <= +1) && keep; iy++) { … … 113 122 } 114 123 source->crNsigma = (nCR > 0) ? fCR / nCR : 0.0; 115 source->extNsigma = (nEXT > 0) ? fabs(fEXT) / nEXT : 0.0;116 124 // NOTE: abs needed to make the Nsigma value positive 117 125 … … 122 130 123 131 if (source->crNsigma > CR_NSIGMA_LIMIT) { 124 source->mode |= PM_SOURCE_MODE_CR LIMIT;132 source->mode |= PM_SOURCE_MODE_CR_LIMIT; 125 133 } 126 134 } … … 169 177 */ 170 178 179 180 // given the PSF ellipse parameters, navigate around the 1sigma contour, return the total 181 // deviation in sigmas. This is measure on the residual image - should we ignore negative 182 // deviations? 183 float psphotModelContour (psImage *image, psImage *weight, psImage *mask, pmModel *model, float Ro) { 184 185 psF32 *PAR = model->params->data.F32; 186 187 // Ro = (x / SXX)^2 + (y / SYY)^2 + x y SXY; 188 // y^2 (1/SYY^2) + y (x SXY) + (x / SXX)^2 - Ro = 0; 189 // y = [-(x SXY) +/- sqrt ((x SXY)^2 - 4 (1/SYY^2) ((x/SXX)^2 - Ro))] * [SYY^2 / 2]; 190 // y = [-B +/- sqrt (B^2 - 4 A C)] / [2 A]; 191 192 // min/max value of x is where T -> 0 193 // solve this for x2: 194 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); 195 if (Q < 0.0) return NAN; // ellipse is imaginary 196 197 int xMax = sqrt(Q); 198 int xMin = -1.0*xMax; 199 200 int nPts = 0; 201 float nSigma = 0.0; 202 203 for (int x = xMin; x <= xMax; x++) { 204 float A = PS_SQR (1.0 / PAR[PM_PAR_SYY]); 205 float B = x * PAR[PM_PAR_SXY]; 206 float C = PS_SQR (x / PAR[PM_PAR_SXX]) - Ro; 207 208 float T = PS_SQR(B) - 4*A*C; 209 if (T < 0.0) continue; 210 211 float yP = (-B + sqrt (T)) / (2.0 * A); 212 float yM = (-B - sqrt (T)) / (2.0 * A); 213 214 int xPix = x + PAR[PM_PAR_XPOS] - image->col0 + 0.5; 215 int yPixM = yM + PAR[PM_PAR_YPOS] - image->row0 + 0.5; 216 int yPixP = yP + PAR[PM_PAR_YPOS] - image->row0 + 0.5; 217 218 if (xPix < 0) continue; 219 if (xPix >= image->numCols) continue; 220 221 if ((yPixM >= 0) && (yPixM < image->numRows)) { 222 if (!mask || !mask->data.U8[yPixM][xPix]) { 223 float dSigma = image->data.F32[yPixM][xPix] / sqrt (weight->data.F32[yPixM][xPix]); 224 nSigma += dSigma; 225 nPts ++; 226 } 227 } 228 229 if (yPixM == yPixP) continue; 230 231 if ((yPixP >= 0) && (yPixP < image->numRows)) { 232 if (!mask || !mask->data.U8[yPixP][xPix]) { 233 float dSigma = image->data.F32[yPixP][xPix] / sqrt (weight->data.F32[yPixP][xPix]); 234 nSigma += dSigma; 235 nPts ++; 236 } 237 } 238 } 239 nSigma /= nPts; 240 return nSigma; 241 } -
trunk/psphot/src/psphotSourceStats.c
r16820 r17396 37 37 for (int i = 0; i < peaks->n; i++) { 38 38 39 pmPeak *peak = peaks->data[i]; 40 if (peak->assigned) continue; 41 39 42 // create a new source, add peak 40 43 pmSource *source = pmSourceAlloc(); 41 source->peak = (pmPeak *)psMemIncrRefCounter(peaks->data[i]);44 source->peak = psMemIncrRefCounter(peak); 42 45 43 46 // allocate image, weight, mask arrays for each peak (square of radius OUTER) 44 47 pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER); 45 48 if (!strcasecmp (breakPt, "PEAKS")) { 49 peak->assigned = true; 46 50 psArrayAdd (sources, 100, source); 47 51 psFree (source); … … 49 53 } 50 54 51 // skip faint sources 55 // skip faint sources for moments measurement 52 56 if (source->peak->SN < MIN_SN) { 57 peak->assigned = true; 53 58 psArrayAdd (sources, 100, source); 54 59 psFree (source); … … 57 62 58 63 // measure a local sky value 59 // XXX EAM : this should use ROBUST not SAMPLE median, but it is broken 64 // the local sky is now ignored; kept here for reference only 60 65 status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, mark); 61 66 if (!status) { … … 80 85 if (status) { 81 86 // add to the source array 87 peak->assigned = true; 82 88 psArrayAdd (sources, 100, source); 83 89 psFree (source); … … 93 99 if (status) { 94 100 // add to the source array 101 peak->assigned = true; 95 102 psArrayAdd (sources, 100, source); 96 103 psFree (source);
Note:
See TracChangeset
for help on using the changeset viewer.
