Changeset 15000 for trunk/psModules/src/objects
- Timestamp:
- Sep 24, 2007, 11:27:58 AM (19 years ago)
- Location:
- trunk/psModules/src/objects
- Files:
-
- 1 added
- 14 edited
-
models/pmModel_GAUSS.c (modified) (2 diffs)
-
models/pmModel_PGAUSS.c (modified) (2 diffs)
-
models/pmModel_QGAUSS.c (modified) (2 diffs)
-
models/pmModel_RGAUSS.c (modified) (2 diffs)
-
models/pmModel_SERSIC.c (modified) (2 diffs)
-
pmPSF.c (modified) (11 diffs)
-
pmPSF.h (modified) (6 diffs)
-
pmPSF_IO.c (modified) (12 diffs)
-
pmPSF_IO.h (modified) (2 diffs)
-
pmPSFtoMetadata.c (added)
-
pmPSFtry.c (modified) (20 diffs)
-
pmPSFtry.h (modified) (4 diffs)
-
pmSourcePhotometry.c (modified) (5 diffs)
-
pmTrend2D.c (modified) (10 diffs)
-
pmTrend2D.h (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/objects/models/pmModel_GAUSS.c
r14652 r15000 274 274 out[i] = in[i]; 275 275 } else { 276 p sPolynomial2D *poly= psf->params->data[i];277 out[i] = p sPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]);276 pmTrend2D *trend = psf->params->data[i]; 277 out[i] = pmTrend2DEval(trend, in[PM_PAR_XPOS], in[PM_PAR_YPOS]); 278 278 } 279 279 } … … 331 331 if (i == PM_PAR_XPOS) continue; 332 332 if (i == PM_PAR_YPOS) continue; 333 psPolynomial2D *poly = psf->params->data[i]; 334 assert (poly); 335 PAR[i] = psPolynomial2DEval(poly, Xo, Yo); 333 pmTrend2D *trend = psf->params->data[i]; 334 PAR[i] = pmTrend2DEval(trend, Xo, Yo); 336 335 } 337 336 -
trunk/psModules/src/objects/models/pmModel_PGAUSS.c
r14652 r15000 328 328 out[i] = in[i]; 329 329 } else { 330 p sPolynomial2D *poly= psf->params->data[i];331 out[i] = p sPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]);330 pmTrend2D *trend = psf->params->data[i]; 331 out[i] = pmTrend2DEval(trend, in[PM_PAR_XPOS], in[PM_PAR_YPOS]); 332 332 } 333 333 } … … 384 384 if (i == PM_PAR_XPOS) continue; 385 385 if (i == PM_PAR_YPOS) continue; 386 psPolynomial2D *poly = psf->params->data[i]; 387 assert (poly); 388 PAR[i] = psPolynomial2DEval(poly, Xo, Yo); 386 pmTrend2D *trend = psf->params->data[i]; 387 PAR[i] = pmTrend2DEval(trend, Xo, Yo); 389 388 } 390 389 -
trunk/psModules/src/objects/models/pmModel_QGAUSS.c
r14652 r15000 358 358 out[i] = in[i]; 359 359 } else { 360 p sPolynomial2D *poly= psf->params->data[i];361 out[i] = p sPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]);360 pmTrend2D *trend = psf->params->data[i]; 361 out[i] = pmTrend2DEval(trend, in[PM_PAR_XPOS], in[PM_PAR_YPOS]); 362 362 } 363 363 } … … 412 412 if (i == PM_PAR_XPOS) continue; 413 413 if (i == PM_PAR_YPOS) continue; 414 psPolynomial2D *poly = psf->params->data[i]; 415 assert (poly); 416 PAR[i] = psPolynomial2DEval(poly, Xo, Yo); 414 pmTrend2D *trend = psf->params->data[i]; 415 PAR[i] = pmTrend2DEval(trend, Xo, Yo); 417 416 } 418 417 -
trunk/psModules/src/objects/models/pmModel_RGAUSS.c
r14652 r15000 351 351 out[i] = in[i]; 352 352 } else { 353 p sPolynomial2D *poly= psf->params->data[i];354 out[i] = p sPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]);353 pmTrend2D *trend = psf->params->data[i]; 354 out[i] = pmTrend2DEval(trend, in[PM_PAR_XPOS], in[PM_PAR_YPOS]); 355 355 } 356 356 } … … 405 405 if (i == PM_PAR_XPOS) continue; 406 406 if (i == PM_PAR_YPOS) continue; 407 psPolynomial2D *poly = psf->params->data[i]; 408 assert (poly); 409 PAR[i] = psPolynomial2DEval(poly, Xo, Yo); 407 pmTrend2D *trend = psf->params->data[i]; 408 PAR[i] = pmTrend2DEval(trend, Xo, Yo); 410 409 } 411 410 -
trunk/psModules/src/objects/models/pmModel_SERSIC.c
r14669 r15000 343 343 out[i] = in[i]; 344 344 } else { 345 p sPolynomial2D *poly= psf->params->data[i];346 out[i] = p sPolynomial2DEval(poly, in[PM_PAR_XPOS], in[PM_PAR_YPOS]);345 pmTrend2D *trend = psf->params->data[i]; 346 out[i] = pmTrend2DEval(trend, in[PM_PAR_XPOS], in[PM_PAR_YPOS]); 347 347 } 348 348 } … … 397 397 if (i == PM_PAR_XPOS) continue; 398 398 if (i == PM_PAR_YPOS) continue; 399 psPolynomial2D *poly = psf->params->data[i]; 400 assert (poly); 401 PAR[i] = psPolynomial2DEval(poly, Xo, Yo); 399 pmTrend2D *trend = psf->params->data[i]; 400 PAR[i] = pmTrend2DEval(trend, Xo, Yo); 402 401 } 403 402 -
trunk/psModules/src/objects/pmPSF.c
r14961 r15000 6 6 * @author EAM, IfA 7 7 * 8 * @version $Revision: 1.2 8$ $Name: not supported by cvs2svn $9 * @date $Date: 2007-09-2 1 02:46:25$8 * @version $Revision: 1.29 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2007-09-24 21:27:58 $ 10 10 * 11 11 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 41 41 42 42 /*****************************************************************************/ 43 /* DEFINE STATEMENTS */44 /*****************************************************************************/45 46 /*****************************************************************************/47 /* TYPE DEFINITIONS */48 /*****************************************************************************/49 50 /*****************************************************************************/51 /* GLOBAL VARIABLES */52 /*****************************************************************************/53 54 /*****************************************************************************/55 /* FILE STATIC VARIABLES */56 /*****************************************************************************/57 58 /*****************************************************************************/59 /* FUNCTION IMPLEMENTATION - LOCAL */60 /*****************************************************************************/61 62 /*****************************************************************************/63 43 /* FUNCTION IMPLEMENTATION - PUBLIC */ 64 44 /*****************************************************************************/ 45 46 static void pmPSFOptionsFree (pmPSFOptions *options) { 47 48 if (!options) return; 49 50 psFree (options->stats); 51 return; 52 } 53 54 pmPSFOptions *pmPSFOptionsAlloc () { 55 56 pmPSFOptions *options = (pmPSFOptions *) psAlloc(sizeof(pmPSFOptions)); 57 psMemSetDeallocator(options, (psFreeFunc) pmPSFOptionsFree); 58 59 options->type = 0; 60 61 options->stats = NULL; 62 63 options->psfTrendMode = PM_TREND_NONE; 64 options->psfTrendNx = 0; 65 options->psfTrendNy = 0; 66 options->psfFieldNx = 0; 67 options->psfFieldNy = 0; 68 options->psfFieldXo = 0; 69 options->psfFieldYo = 0; 70 71 options->poissonErrorsPhotLMM = true; 72 options->poissonErrorsPhotLin = false; 73 options->poissonErrorsParams = true; 74 75 return options; 76 } 65 77 66 78 /***************************************************************************** … … 74 86 75 87 psFree (psf->ChiTrend); 88 psFree (psf->psfTrendStats); 76 89 psFree (psf->ApTrend); 77 90 psFree (psf->FluxScale); … … 94 107 Object Normalization 95 108 *****************************************************************************/ 96 pmPSF *pmPSFAlloc (pm ModelType type, bool poissonErrors, psPolynomial2D *psfTrendMask)109 pmPSF *pmPSFAlloc (pmPSFOptions *options) 97 110 { 98 111 int Nparams; 99 112 100 113 pmPSF *psf = (pmPSF *) psAlloc(sizeof(pmPSF)); 101 102 psf->type = type; 114 psMemSetDeallocator(psf, (psFreeFunc) pmPSFFree); 115 116 psf->type = options->type; 103 117 psf->chisq = 0.0; 104 118 psf->ApResid = 0.0; … … 108 122 psf->nPSFstars = 0; 109 123 psf->nApResid = 0; 110 psf->poissonErrors = poissonErrors; 124 125 psf->poissonErrorsPhotLMM = options->poissonErrorsPhotLMM; 126 psf->poissonErrorsPhotLin = options->poissonErrorsPhotLin; 127 psf->poissonErrorsParams = options->poissonErrorsParams; 111 128 112 129 // the ApTrend components are (x, y). It may be represented with a polynomial or with a … … 121 138 psf->FluxScale = NULL; 122 139 123 if (psf->poissonErrors ) {140 if (psf->poissonErrorsPhotLMM) { 124 141 psf->ChiTrend = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1); 125 142 } else { … … 133 150 psf->residuals = NULL; 134 151 135 Nparams = pmModelClassParameterCount ( type);152 Nparams = pmModelClassParameterCount (options->type); 136 153 if (!Nparams) { 137 154 psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType"); … … 139 156 } 140 157 psf->params = psArrayAlloc(Nparams); 158 159 // save the trend stats on the psf for use in pmPSFFromPSFtry 160 psf->psfTrendStats = psMemIncrRefCounter (options->stats); 141 161 142 162 // the order of the PSF parameter (X,Y) fits is determined by the psfTrendMask polynomial … … 147 167 // PM_PAR_XPOS, etc) 148 168 149 if (psfTrendMask) { 169 psImageBinning *binning = psImageBinningAlloc(); 170 binning->nXruff = options->psfTrendNx; 171 binning->nYruff = options->psfTrendNy; 172 binning->nXfine = options->psfFieldNx; 173 binning->nYfine = options->psfFieldNy; 174 psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER); 175 psImageBinningSetSkipByOffset (binning, options->psfFieldXo, options->psfFieldYo); 176 177 psf->fieldNx = options->psfFieldNx; 178 psf->fieldNy = options->psfFieldNy; 179 psf->fieldXo = options->psfFieldXo; 180 psf->fieldYo = options->psfFieldYo; 181 182 // define the parameter trends 183 if (options->psfTrendMode != PM_TREND_NONE) { 150 184 for (int i = 0; i < psf->params->n; i++) { 151 if (i == PM_PAR_SKY) 152 continue; 153 if (i == PM_PAR_I0) 154 continue; 155 if (i == PM_PAR_XPOS) 156 continue; 157 if (i == PM_PAR_YPOS) 158 continue; 159 160 psPolynomial2D *param = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, psfTrendMask->nX, psfTrendMask->nY); 161 for (int nx = 0; nx < param->nX + 1; nx++) { 162 for (int ny = 0; ny < param->nY + 1; ny++) { 163 param->mask[nx][ny] = psfTrendMask->mask[nx][ny]; 164 } 165 } 166 psf->params->data[i] = param; 185 if (i == PM_PAR_SKY) continue; 186 if (i == PM_PAR_I0) continue; 187 if (i == PM_PAR_XPOS) continue; 188 if (i == PM_PAR_YPOS) continue; 189 190 psf->params->data[i] = pmTrend2DNoImageAlloc (options->psfTrendMode, binning, options->stats); 167 191 } 168 192 } 169 170 psMemSetDeallocator(psf, (psFreeFunc) pmPSFFree); 171 return(psf); 193 psFree (binning); 194 return psf; 172 195 } 173 196 … … 300 323 va_start(ap, sxy); 301 324 302 pmModelType type = pmModelClassGetType (typeName); 303 psPolynomial2D *psfTrend = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 0, 0); 304 pmPSF *psf = pmPSFAlloc (type, true, psfTrend); 325 pmPSFOptions *options = pmPSFOptionsAlloc (); 326 options->type = pmModelClassGetType (typeName); 327 options->psfTrendMode = PM_TREND_POLY_ORD; 328 options->psfTrendNx = 0; 329 options->psfTrendNy = 0; 330 331 pmPSF *psf = pmPSFAlloc (options); 305 332 306 333 psVector *par = psVectorAlloc (psf->params->n, PS_TYPE_F32); … … 311 338 psEllipsePol pol = pmPSF_ModelToFit(par->data.F32); 312 339 340 pmTrend2D *trend = NULL; 341 313 342 // set the psf shape parameters 314 psPolynomial2D *poly = NULL; 315 poly = psf->params->data[PM_PAR_E0]; 316 poly->coeff[0][0] = pol.e0; 317 318 poly = psf->params->data[PM_PAR_E1]; 319 poly->coeff[0][0] = pol.e1; 320 321 poly = psf->params->data[PM_PAR_E2]; 322 poly->coeff[0][0] = pol.e2; 343 trend = psf->params->data[PM_PAR_E0]; 344 trend->poly->coeff[0][0] = pol.e0; 345 346 trend = psf->params->data[PM_PAR_E1]; 347 trend->poly->coeff[0][0] = pol.e1; 348 349 trend = psf->params->data[PM_PAR_E2]; 350 trend->poly->coeff[0][0] = pol.e2; 323 351 324 352 for (int i = PM_PAR_SXY + 1; i < psf->params->n; i++) { 325 poly= psf->params->data[i];326 poly->coeff[0][0] = (psF32)va_arg(ap, psF64);353 trend = psf->params->data[i]; 354 trend->poly->coeff[0][0] = (psF32)va_arg(ap, psF64); 327 355 } 328 356 va_end(ap); 329 357 330 358 psFree (par); 331 psFree ( psfTrend);359 psFree (options); 332 360 return psf; 333 361 } -
trunk/psModules/src/objects/pmPSF.h
r14936 r15000 6 6 * @author EAM, IfA 7 7 * 8 * @version $Revision: 1.1 6$ $Name: not supported by cvs2svn $9 * @date $Date: 2007-09-2 1 00:05:35$8 * @version $Revision: 1.17 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2007-09-24 21:27:58 $ 10 10 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii 11 11 */ … … 19 19 // type of model carried by the pmModel structure 20 20 typedef int pmModelType; 21 22 typedef enum {23 PM_PSF_APTREND_ERROR = 0,24 PM_PSF_APTREND_NONE,25 PM_PSF_APTREND_CONSTANT,26 PM_PSF_APTREND_SKYBIAS,27 PM_PSF_APTREND_SKYSAT,28 PM_PSF_APTREND_XY_LIN,29 PM_PSF_APTREND_XY_QUAD,30 PM_PSF_APTREND_SKY_XY_LIN,31 PM_PSF_APTREND_SKYSAT_XY_LIN,32 PM_PSF_APTREND_ALL33 } pmPSFApTrendOptions;34 21 35 22 /** pmPSF data structure … … 49 36 pmModelType type; ///< PSF Model in use 50 37 psArray *params; ///< Model parameters (psPolynomial2D) 38 psStats *psfTrendStats; ///< psf parameter trend clipping stats 51 39 psPolynomial1D *ChiTrend; ///< Chisq vs flux fit (correction for systematic errors) 52 40 pmTrend2D *ApTrend; ///< ApResid vs (x,y) … … 59 47 int nPSFstars; ///< number of stars used to measure PSF 60 48 int nApResid; ///< number of stars used to measure ApResid 61 bool poissonErrors; 49 int fieldNx; 50 int fieldNy; 51 int fieldXo; 52 int fieldYo; 53 bool poissonErrorsPhotLMM; ///< use poission errors for non-linear model fitting 54 bool poissonErrorsPhotLin; ///< use poission errors for linear model fitting 55 bool poissonErrorsParams; ///< use poission errors for model parameter fitting 62 56 pmGrowthCurve *growth; ///< apMag vs Radius 63 57 pmResiduals *residuals; ///< normalized residual image (no spatial variation) 64 58 } 65 59 pmPSF; 60 61 typedef struct { 62 pmModelType type; 63 psStats *stats; // psfTrend clipping stats 64 pmTrend2DMode psfTrendMode; 65 int psfTrendNx; 66 int psfTrendNy; 67 int psfFieldNx; 68 int psfFieldNy; 69 int psfFieldXo; 70 int psfFieldYo; 71 bool poissonErrorsPhotLMM; ///< use poission errors for non-linear model fitting 72 bool poissonErrorsPhotLin; ///< use poission errors for linear model fitting 73 bool poissonErrorsParams; ///< use poission errors for model parameter fitting 74 float radius; 75 } pmPSFOptions; 66 76 67 77 # define PM_PAR_E0 PM_PAR_SXX … … 74 84 * 75 85 */ 76 pmPSF *pmPSFAlloc( 77 pmModelType type, // type of model for PSF 78 bool poissonErrors, 79 psPolynomial2D *psfTrendMask 80 ); 81 82 bool pmPSFMaskApTrend (psPolynomial4D *trend, pmPSFApTrendOptions option); 83 84 pmPSFApTrendOptions pmPSFApTrendOptionFromName (char *name); 86 pmPSF *pmPSFAlloc (pmPSFOptions *options); 87 pmPSFOptions *pmPSFOptionsAlloc (); 85 88 86 89 double pmPSF_SXYfromModel (psF32 *modelPar); … … 91 94 92 95 bool pmPSF_AxesToModel (psF32 *modelPar, psEllipseAxes axes); 96 bool pmPSF_FitToModel (psF32 *fittedPar, float minMinorAxis); 97 98 psEllipsePol pmPSF_ModelToFit (psF32 *modelPar); 93 99 psEllipseAxes pmPSF_ModelToAxes (psF32 *modelPar, double maxAR); 94 bool pmPSF_FitToModel (psF32 *fittedPar, float minMinorAxis);95 psEllipsePol pmPSF_ModelToFit (psF32 *modelPar);96 100 97 101 /// @} -
trunk/psModules/src/objects/pmPSF_IO.c
r14936 r15000 6 6 * @author EAM, IfA 7 7 * 8 * @version $Revision: 1.2 2$ $Name: not supported by cvs2svn $9 * @date $Date: 2007-09-2 1 00:05:57$8 * @version $Revision: 1.23 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2007-09-24 21:27:58 $ 10 10 * 11 11 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 148 148 149 149 // for a pmPSF supplied on the analysis metadata, we write out 150 // if needed: 151 // - a PHU blank header 150 152 // - image header : FITS Image NAXIS = 0 151 // - psf table (+header) : FITS Table 152 // - psf resid (+header) : FITS Image 153 // if needed, we also write out a PHU blank header 153 // if (trendMode == MAP) 154 // - psf resid (+header) : FITS Image 155 // else 156 // - psf table (+header) : FITS Table 154 157 bool pmPSFmodelWrite (psMetadata *analysis, const pmFPAview *view, pmFPAfile *file, const pmConfig *config) 155 158 { … … 159 162 160 163 if (!analysis) return false; 164 165 // select the current recipe 166 psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, "PSPHOT"); 167 if (!recipe) { 168 psError(PS_ERR_UNKNOWN, false, "missing recipe %s", "PSPHOT"); 169 return false; 170 } 161 171 162 172 // write a PHU? (only if input image is MEF) … … 215 225 216 226 // write out the IMAGE header segment 217 // this header block is new, write it to disk227 // if this header block is new, write it to disk 218 228 if (hdu->header != file->header) { 219 229 // add EXTNAME, EXTHEAD, EXTTYPE to header … … 250 260 psMetadataAddStr (header, PS_LIST_TAIL, "PSF_NAME", 0, "PSF model name", modelName); 251 261 252 psMetadataAddBool (header, PS_LIST_TAIL, "POISSON", 0, "Use Poisson errors in fits?", psf->poissonErrors); 262 psMetadataAddBool (header, PS_LIST_TAIL, "ERR_LMM", 0, "Use Poisson errors in fits?", psf->poissonErrorsPhotLMM); 263 psMetadataAddBool (header, PS_LIST_TAIL, "ERR_LIN", 0, "Use Poisson errors in fits?", psf->poissonErrorsPhotLin); 264 psMetadataAddBool (header, PS_LIST_TAIL, "ERR_PAR", 0, "Use Poisson errors in fits?", psf->poissonErrorsParams); 253 265 254 266 int nPar = pmModelClassParameterCount (psf->type) ; 255 267 psMetadataAdd (header, PS_LIST_TAIL, "PSF_NPAR", PS_DATA_S32, "PSF model parameter count", nPar); 268 269 psMetadataAddS32 (header, PS_LIST_TAIL, "IMAXIS1", 0, "Image X Size", psf->fieldNx); 270 psMetadataAddS32 (header, PS_LIST_TAIL, "IMAXIS2", 0, "Image Y Size", psf->fieldNy); 271 psMetadataAddS32 (header, PS_LIST_TAIL, "IMREF1", 0, "Image X Ref", psf->fieldXo); 272 psMetadataAddS32 (header, PS_LIST_TAIL, "IMREF2", 0, "Image Y Ref", psf->fieldYo); 273 274 // extract PSF Clump info 275 pmPSFClump psfClump; 276 277 psfClump.X = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.X"); assert (status); 278 psfClump.Y = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.Y"); assert (status); 279 psfClump.dX = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.DX"); assert (status); 280 psfClump.dY = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.DY"); assert (status); 281 282 psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLX", 0, "psf clump center", psfClump.X); 283 psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLY", 0, "psf clump center", psfClump.Y); 284 psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLDX", 0, "psf clump size", psfClump.dX); 285 psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLDY", 0, "psf clump size", psfClump.dY); 256 286 257 287 // save the dimensions of each parameter 258 288 for (int i = 0; i < nPar; i++) { 259 289 char name[9]; 260 psPolynomial2D *poly = psf->params->data[i]; 261 if (poly == NULL) continue; 290 int nX, nY; 291 292 pmTrend2D *trend = psf->params->data[i]; 293 if (trend == NULL) continue; 294 295 if (trend->mode == PM_TREND_MAP) { 296 nX = trend->map->map->numCols; 297 nY = trend->map->map->numRows; 298 } else { 299 nX = trend->poly->nX; 300 nY = trend->poly->nY; 301 } 262 302 snprintf (name, 9, "PAR%02d_NX", i); 263 psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", poly->nX);303 psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", nX); 264 304 snprintf (name, 9, "PAR%02d_NY", i); 265 psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", poly->nY); 305 psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", nY); 306 snprintf (name, 9, "PAR%02d_MD", i); 307 char *modeName = pmTrend2DModeToString (trend->mode); 308 psMetadataAddStr (header, PS_LIST_TAIL, name, 0, "", modeName); 309 psFree (modeName); 266 310 } 267 311 268 312 // other required information describing the PSF 269 psMetadataAddF32 270 psMetadataAddF32 271 psMetadataAddF32 272 psMetadataAddS32 313 psMetadataAddF32 (header, PS_LIST_TAIL, "AP_RESID", 0, "aperture residual", psf->ApResid); 314 psMetadataAddF32 (header, PS_LIST_TAIL, "AP_ERROR", 0, "aperture residual scatter", psf->dApResid); 315 psMetadataAddF32 (header, PS_LIST_TAIL, "CHISQ", 0, "chi-square for fit", psf->chisq); 316 psMetadataAddS32 (header, PS_LIST_TAIL, "NSTARS", 0, "number of stars used to measure PSF", psf->nPSFstars); 273 317 274 318 // XXX can we drop this now? … … 278 322 psArray *psfTable = psArrayAllocEmpty (100); 279 323 for (int i = 0; i < nPar; i++) { 280 psPolynomial2D *poly = psf->params->data[i]; 281 if (poly == NULL) continue; // skip unset parameters (eg, XPOS) 282 for (int ix = 0; ix <= poly->nX; ix++) { 283 for (int iy = 0; iy <= poly->nY; iy++) { 284 285 psMetadata *row = psMetadataAlloc (); 286 psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i); 287 psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER", 0, "", ix); 288 psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER", 0, "", iy); 289 psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE", 0, "", poly->coeff[ix][iy]); 290 psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR", 0, "", poly->coeffErr[ix][iy]); 291 psMetadataAddU8 (row, PS_LIST_TAIL, "MASK", 0, "", poly->mask[ix][iy]); 292 293 psArrayAdd (psfTable, 100, row); 294 psFree (row); 324 pmTrend2D *trend = psf->params->data[i]; 325 if (trend == NULL) continue; // skip unset parameters (eg, XPOS) 326 327 if (trend->mode == PM_TREND_MAP) { 328 // write the image components into a table: this is needed because they may each be a different size 329 psImageMap *map = trend->map; 330 for (int ix = 0; ix < map->map->numCols; ix++) { 331 for (int iy = 0; iy < map->map->numRows; iy++) { 332 psMetadata *row = psMetadataAlloc (); 333 psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i); 334 psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER", 0, "", ix); 335 psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER", 0, "", iy); 336 psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE", 0, "", map->map->data.F32[iy][ix]); 337 psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR", 0, "", map->error->data.F32[iy][ix]); 338 psMetadataAddU8 (row, PS_LIST_TAIL, "MASK", 0, "", 0); // no cells are masked 339 340 psArrayAdd (psfTable, 100, row); 341 psFree (row); 342 } 343 } 344 } else { 345 // write the polynomial components into a table 346 psPolynomial2D *poly = trend->poly; 347 for (int ix = 0; ix <= poly->nX; ix++) { 348 for (int iy = 0; iy <= poly->nY; iy++) { 349 psMetadata *row = psMetadataAlloc (); 350 psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i); 351 psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER", 0, "", ix); 352 psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER", 0, "", iy); 353 psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE", 0, "", poly->coeff[ix][iy]); 354 psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR", 0, "", poly->coeffErr[ix][iy]); 355 psMetadataAddU8 (row, PS_LIST_TAIL, "MASK", 0, "", poly->mask[ix][iy]); 356 357 psArrayAdd (psfTable, 100, row); 358 psFree (row); 359 } 295 360 } 296 361 } … … 356 421 // XXX save the growth curve 357 422 // XXX save ApTrend (as image?) 423 // XXX write the ApTrend with the same API as will be used for the PSF parameters above 358 424 359 425 # if (0) … … 474 540 psTrace ("psModules.objects", 5, "read psf model for %s\n", file->filename); 475 541 542 // select the current recipe 543 psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, "PSPHOT"); 544 if (!recipe) { 545 psError(PS_ERR_UNKNOWN, false, "missing recipe %s", "PSPHOT"); 546 return false; 547 } 548 476 549 // Menu of EXTNAME rules 477 550 psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES"); … … 506 579 if (!header) psAbort("cannot read table header"); 507 580 581 pmPSFOptions *options = pmPSFOptionsAlloc(); 582 508 583 // load the PSF model parameters from the FITS table 509 584 char *modelName = psMetadataLookupStr (&status, header, "PSF_NAME"); 510 pmModelTypetype = pmModelClassGetType (modelName);511 if ( type == -1) {585 options->type = pmModelClassGetType (modelName); 586 if (options->type == -1) { 512 587 psError(PS_ERR_UNKNOWN, true, "invalid model name %s in psf file %s", modelName, file->filename); 513 588 return false; 514 589 } 515 590 516 bool poissonErrors = psMetadataLookupBool (&status, header, "POISSON"); 591 // psf clump data 592 pmPSFClump psfClump; 593 594 psfClump.X = psMetadataLookupF32 (&status, header, "PSF_CLX" ); assert(status); 595 psfClump.Y = psMetadataLookupF32 (&status, header, "PSF_CLY" ); assert(status); 596 psfClump.dX = psMetadataLookupF32 (&status, header, "PSF_CLDX"); assert(status); 597 psfClump.dY = psMetadataLookupF32 (&status, header, "PSF_CLDY"); assert(status); 598 599 psMetadataAddF32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.X" , 0, "psf clump center", psfClump.X); 600 psMetadataAddF32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.Y" , 0, "psf clump center", psfClump.Y); 601 psMetadataAddF32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.DX", 0, "psf clump size", psfClump.dX); 602 psMetadataAddF32 (recipe, PS_LIST_TAIL, "PSF.CLUMP.DY", 0, "psf clump size", psfClump.dY); 603 604 options->poissonErrorsPhotLMM = psMetadataLookupBool (&status, header, "ERR_LMM"); 605 options->poissonErrorsPhotLin = psMetadataLookupBool (&status, header, "ERR_LIN"); 606 options->poissonErrorsParams = psMetadataLookupBool (&status, header, "ERR_PAR"); 607 608 options->psfFieldNx = psMetadataLookupS32 (&status, header, "IMAXIS1"); 609 options->psfFieldNy = psMetadataLookupS32 (&status, header, "IMAXIS2"); 610 options->psfFieldXo = psMetadataLookupS32 (&status, header, "IMREF1"); 611 options->psfFieldYo = psMetadataLookupS32 (&status, header, "IMREF2"); 612 613 psImageBinning *binning = psImageBinningAlloc(); 614 binning->nXfine = options->psfFieldNx; 615 binning->nYfine = options->psfFieldNy; 517 616 518 617 // we determine the PSF parameter polynomials from the MD-defined polynomials 519 pmPSF *psf = pmPSFAlloc ( type, poissonErrors, NULL);618 pmPSF *psf = pmPSFAlloc (options); 520 619 521 620 // check the number of expected parameters … … 524 623 psAbort("mismatch model par count"); 525 624 526 // load the dimensions of each parameter625 // load the trend mode and dimensions of each parameter 527 626 for (int i = 0; i < nPar; i++) { 528 627 char name[9]; 529 628 snprintf (name, 9, "PAR%02d_NX", i); 530 int nXorder= psMetadataLookupS32 (&status, header, name);629 binning->nXruff = psMetadataLookupS32 (&status, header, name); 531 630 if (!status) continue; // not all parameters are defined 631 532 632 snprintf (name, 9, "PAR%02d_NY", i); 533 int nYorder= psMetadataLookupS32 (&status, header, name);633 binning->nYruff = psMetadataLookupS32 (&status, header, name); 534 634 if (!status) { 535 psError(PS_ERR_UNKNOWN, true, "inconsistent PSF header: NX defined for PAR %d, but no NY", i);635 psError(PS_ERR_UNKNOWN, true, "inconsistent PSF header: NX defined for PAR %d, but not NY", i); 536 636 return false; 537 637 } 538 psf->params->data[i] = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, nXorder, nYorder); 539 } 638 639 snprintf (name, 9, "PAR%02d_MD", i); 640 char *modeName = psMetadataLookupStr (&status, header, name); 641 if (!status) { 642 psError(PS_ERR_UNKNOWN, true, "inconsistent PSF header: NX & NY defined for PAR %d, but not MD", i); 643 return false; 644 } 645 pmTrend2DMode psfTrendMode = pmTrend2DModeFromString (modeName); 646 if (psfTrendMode == PM_TREND_NONE) { 647 psfTrendMode = PM_TREND_POLY_ORD; 648 } 649 650 psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER); 651 psImageBinningSetSkipByOffset (binning, options->psfFieldXo, options->psfFieldYo); 652 psf->params->data[i] = pmTrend2DNoImageAlloc (psfTrendMode, binning, NULL); 653 } 654 psFree (binning); 540 655 541 656 // other required information describing the PSF … … 557 672 int xPow = psMetadataLookupS32 (&status, row, "X_POWER"); 558 673 int yPow = psMetadataLookupS32 (&status, row, "Y_POWER"); 559 // XXX sanity check here 560 561 psPolynomial2D *poly = psf->params->data[iPar]; 562 if (poly == NULL) { 563 psError(PS_ERR_UNKNOWN, true, "values for parameter %d, but missing NX", iPar); 674 675 pmTrend2D *trend = psf->params->data[iPar]; 676 if (trend == NULL) { 677 psError(PS_ERR_UNKNOWN, true, "parameter %d not available", iPar); 564 678 return false; 565 679 } 566 680 567 poly->coeff[xPow][yPow] = psMetadataLookupF32 (&status, row, "VALUE"); 568 poly->coeffErr[xPow][yPow] = psMetadataLookupF32 (&status, row, "ERROR"); 569 poly->mask[xPow][yPow] = psMetadataLookupU8 (&status, row, "MASK"); 681 if (trend->mode == PM_TREND_MAP) { 682 psImageMap *map = trend->map; 683 assert (map); 684 assert (map->map); 685 assert (map->error); 686 assert (xPow >= 0); 687 assert (yPow >= 0); 688 assert (xPow < map->map->numCols); 689 assert (yPow < map->map->numRows); 690 map->map->data.F32[yPow][xPow] = psMetadataLookupF32 (&status, row, "VALUE"); 691 map->error->data.F32[yPow][xPow] = psMetadataLookupF32 (&status, row, "ERROR"); 692 } else { 693 psPolynomial2D *poly = trend->poly; 694 assert (poly); 695 assert (xPow >= 0); 696 assert (yPow >= 0); 697 assert (xPow <= poly->nX); 698 assert (yPow <= poly->nY); 699 poly->coeff[xPow][yPow] = psMetadataLookupF32 (&status, row, "VALUE"); 700 poly->coeffErr[xPow][yPow] = psMetadataLookupF32 (&status, row, "ERROR"); 701 poly->mask[xPow][yPow] = psMetadataLookupU8 (&status, row, "MASK"); 702 } 570 703 } 571 704 psFree (header); … … 618 751 619 752 // create a psMetadata representation (human-readable) of a psf model 620 psMetadata *pmPSFtoMetadata (psMetadata *metadata, pmPSF *psf) 621 { 622 623 if (metadata == NULL) { 624 metadata = psMetadataAlloc (); 625 } 626 627 char *modelName = pmModelClassGetName (psf->type); 628 psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_MODEL_NAME", PS_DATA_STRING, "PSF model name", modelName); 629 630 int nPar = pmModelClassParameterCount (psf->type) ; 631 psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_MODEL_NPAR", PS_DATA_S32, "PSF model parameter count", nPar); 632 633 for (int i = 0; i < nPar; i++) { 634 psPolynomial2D *poly = psf->params->data[i]; 635 if (poly == NULL) 636 continue; 637 psPolynomial2DtoMetadata (metadata, poly, "PSF_PAR%02d", i); 638 } 639 640 // XXX fix this 641 psWarning ("APTREND is currently missing"); 642 // psPolynomial4DtoMetadata (metadata, psf->ApTrend, "APTREND"); 643 644 psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_AP_RESID", PS_DATA_F32, "aperture residual", psf->ApResid); 645 psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_dAP_RESID", PS_DATA_F32, "aperture residual scatter", psf->dApResid); 646 psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_SKY_BIAS", PS_DATA_F32, "sky bias level", psf->skyBias); 647 648 psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_CHISQ", PS_DATA_F32, "chi-square for fit", psf->chisq); 649 psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_NSTARS", PS_DATA_S32, "number of stars used to measure PSF", psf->nPSFstars); 650 psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_POISSON_ERRORS", PS_DATA_BOOL, "Poisson errors for fits", psf->poissonErrors); 651 652 return metadata; 653 } 654 655 // parse a psMetadata representation (human-readable) of a psf model 656 pmPSF *pmPSFfromMetadata (psMetadata *metadata) 657 { 658 659 bool status; 660 char keyword[80]; 661 662 char *modelName = psMetadataLookupPtr (&status, metadata, "PSF_MODEL_NAME"); 663 pmModelType type = pmModelClassGetType (modelName); 664 665 bool poissonErrors = psMetadataLookupPtr (&status, metadata, "PSF_POISSON_ERRORS"); 666 if (!status) 667 poissonErrors = true; 668 669 // we determine the PSF parameter polynomials from the MD-defined polynomials 670 pmPSF *psf = pmPSFAlloc (type, poissonErrors, NULL); 671 672 int nPar = psMetadataLookupS32 (&status, metadata, "PSF_MODEL_NPAR"); 673 if (nPar != pmModelClassParameterCount (psf->type)) 674 psAbort("mismatch model par count"); 675 676 // un-fitted terms, not in the Metadata, are left NULL 677 // XXX add a double-check of the expected number? 678 for (int i = 0; i < nPar; i++) { 679 sprintf (keyword, "PSF_PAR%02d", i); 680 psMetadata *folder = psMetadataLookupPtr (&status, metadata, keyword); 681 if (!status) 682 continue; 683 psPolynomial2D *poly = psPolynomial2DfromMetadata (folder); 684 psFree (psf->params->data[i]); 685 psf->params->data[i] = poly; 686 } 687 688 // load the APTREND data 689 // XXX fix this to work with pmTrend2D 690 psWarning ("APTREND is not being read"); 691 # if (0) 692 sprintf (keyword, "APTREND"); 693 psMetadata *folder = psMetadataLookupPtr (&status, metadata, keyword); 694 psPolynomial4D *poly = psPolynomial4DfromMetadata (folder); 695 psFree (psf->ApTrend); 696 psf->ApTrend = poly; 697 # endif 698 699 psf->ApResid = psMetadataLookupF32 (&status, metadata, "PSF_AP_RESID"); 700 psf->dApResid = psMetadataLookupF32 (&status, metadata, "PSF_dAP_RESID"); 701 psf->skyBias = psMetadataLookupF32 (&status, metadata, "PSF_SKY_BIAS"); 702 703 psf->chisq = psMetadataLookupF32 (&status, metadata, "PSF_CHISQ"); 704 psf->nPSFstars = psMetadataLookupS32 (&status, metadata, "PSF_NSTARS"); 705 706 psFree (metadata); 707 return (psf); 708 } 753 // XXX pmPSF to/from Metadata functions were defined for 1.22 and earlier, but were dropped -
trunk/psModules/src/objects/pmPSF_IO.h
r14652 r15000 6 6 * @author EAM, IfA 7 7 * 8 * @version $Revision: 1. 9$ $Name: not supported by cvs2svn $9 * @date $Date: 2007-0 8-24 00:11:02$8 * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2007-09-24 21:27:58 $ 10 10 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii 11 11 */ … … 33 33 bool pmPSFmodelCheckDataStatusForChip (const pmChip *chip); 34 34 35 psMetadata *pmPSFtoMetadata (psMetadata *metadata, pmPSF *psf);36 pmPSF *pmPSFfromMetadata (psMetadata *metadata);37 38 35 /// @} 39 36 # endif -
trunk/psModules/src/objects/pmPSFtry.c
r14937 r15000 5 5 * @author EAM, IfA 6 6 * 7 * @version $Revision: 1.4 5$ $Name: not supported by cvs2svn $8 * @date $Date: 2007-09-2 1 00:06:57$7 * @version $Revision: 1.46 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2007-09-24 21:27:58 $ 9 9 * 10 10 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 44 44 static void pmPSFtryFree (pmPSFtry *test) 45 45 { 46 47 if (test == NULL) 48 return; 46 if (test == NULL) return; 49 47 50 48 psFree (test->psf); … … 58 56 59 57 // allocate a pmPSFtry based on the desired sources and the model (identified by name) 60 pmPSFtry *pmPSFtryAlloc (psArray *sources, char *modelName, bool poissonErrors, psPolynomial2D *psfTrendMask)58 pmPSFtry *pmPSFtryAlloc (psArray *sources, pmPSFOptions *options) 61 59 { 62 63 // validate the requested model name64 pmModelType type = pmModelClassGetType (modelName);65 if (type == -1) {66 psError (PS_ERR_UNKNOWN, true, "invalid model name %s", modelName);67 return NULL;68 }69 70 60 pmPSFtry *test = (pmPSFtry *) psAlloc(sizeof(pmPSFtry)); 71 72 test->psf = pmPSFAlloc (type, poissonErrors, psfTrendMask); 61 psMemSetDeallocator(test, (psFreeFunc) pmPSFtryFree); 62 63 test->psf = pmPSFAlloc (options); 73 64 test->metric = psVectorAlloc (sources->n, PS_TYPE_F32); 74 65 test->metricErr = psVectorAlloc (sources->n, PS_TYPE_F32); … … 76 67 test->mask = psVectorAlloc (sources->n, PS_TYPE_U8); 77 68 69 psVectorInit (test->mask, 0); 70 psVectorInit (test->metric, 0.0); 71 psVectorInit (test->metricErr, 0.0); 72 psVectorInit (test->fitMag, 0.0); 73 78 74 test->sources = psArrayAlloc (sources->n); 75 79 76 for (int i = 0; i < sources->n; i++) { 80 77 test->sources->data[i] = pmSourceCopy (sources->data[i]); 81 test->mask->data.U8[i] = 0; 82 test->metric->data.F32[i] = 0; 83 test->metricErr->data.F32[i] = 0; 84 test->fitMag->data.F32[i] = 0; 85 } 86 87 psMemSetDeallocator(test, (psFreeFunc) pmPSFtryFree); 78 } 79 88 80 return (test); 89 81 } … … 99 91 100 92 // generate a pmPSFtry with a copy of the test PSF sources 101 pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, float RADIUS, bool poissonErrors, psPolynomial2D *psfTrendMask, bool applyWeights, psMaskType maskVal, psMaskType mark)93 pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, pmPSFOptions *options, psMaskType maskVal, psMaskType mark) 102 94 { 103 95 bool status; … … 105 97 int Npsf = 0; 106 98 107 pmPSFtry *psfTry = pmPSFtryAlloc (sources, modelName, poissonErrors, psfTrendMask); 99 // validate the requested model name 100 options->type = pmModelClassGetType (modelName); 101 if (options->type == -1) { 102 psError (PS_ERR_UNKNOWN, true, "invalid model name %s", modelName); 103 return NULL; 104 } 105 106 pmPSFtry *psfTry = pmPSFtryAlloc (sources, options); 108 107 if (psfTry == NULL) { 109 108 psError (PS_ERR_UNKNOWN, false, "failed to allocate psf model"); … … 111 110 } 112 111 113 // stage 1: fit an independent model (freeModel) to allsources112 // stage 1: fit an EXT model to all candidates PSF sources 114 113 psTimerStart ("fit"); 115 114 for (int i = 0; i < psfTry->sources->n; i++) { … … 123 122 } 124 123 125 // set object mask to define valid pixels 126 psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, RADIUS, "OR", mark);124 // set object mask to define valid pixels -- XXX not unset? 125 psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", mark); 127 126 128 127 // fit model as EXT, not PSF … … 140 139 141 140 // stage 2: construct a psf (pmPSF) from this collection of model fits 142 if (!pmPSFFromPSFtry (psfTry , applyWeights)) {141 if (!pmPSFFromPSFtry (psfTry)) { 143 142 psError(PS_ERR_UNKNOWN, false, "failed to construct a psf model from collection of sources"); 144 143 psFree(psfTry); … … 162 161 continue; 163 162 } 164 source->modelPSF->radiusFit = RADIUS;165 166 // set object mask to define valid pixels 167 psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, RADIUS, "OR", mark);163 source->modelPSF->radiusFit = options->radius; 164 165 // set object mask to define valid pixels -- XXX not unset? 166 psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", mark); 168 167 169 168 // fit the PSF model to the source … … 241 240 242 241 // XXX this function wants aperture radius for pmSourcePhotometry 243 pmPSFtryMetric (psfTry, RADIUS);242 pmPSFtryMetric (psfTry, options->radius); 244 243 psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f : sky bias: %f\n", 245 244 modelName, psfTry->psf->ApResid, psfTry->psf->dApResid, psfTry->psf->skyBias); … … 339 338 /***************************************************************************** 340 339 pmPSFFromPSFtry (psfTry): build a PSF model from a collection of 341 source->modelEXT ent ires. The PSF ignores the first 4 (independent) model340 source->modelEXT entries. The PSF ignores the first 4 (independent) model 342 341 parameters and constructs a polynomial fit to the remaining as a function of 343 342 image coordinate. … … 345 344 Note: some of the array entries may be NULL (failed fits); ignore them. 346 345 *****************************************************************************/ 347 bool pmPSFFromPSFtry (pmPSFtry *psfTry , bool applyWeight)346 bool pmPSFFromPSFtry (pmPSFtry *psfTry) 348 347 { 349 348 pmPSF *psf = psfTry->psf; … … 355 354 356 355 psVector *dz = NULL; 357 if ( applyWeight) {356 if (psf->poissonErrorsParams) { 358 357 dz = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); 359 358 } … … 377 376 // more deviant than three sigma. 378 377 // mask is currently updated for each pass. this is inconsistent? 379 380 psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);381 378 382 379 // The shape parameters (SXX, SXY, SYY) are strongly coupled. We have to handle them very … … 389 386 // threshold. 390 387 388 // XXX this function needs to check the fit residuals and modify Nx,Ny as needed 389 391 390 // convert the measured source shape paramters to polarization terms 392 391 psVector *e0 = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); … … 406 405 // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each. 407 406 // This way, the parameters masked by one of the fits will be applied to the others 408 for (int i = 0; i < stats->clipIter; i++) { 409 psVectorClipFitPolynomial2D (psf->params->data[PM_PAR_E0], stats, psfTry->mask, 0xff, e0, dz, x, y); 410 psTrace ("psModules.pmPSFtry", 4, "clipped E0 : keeping %ld of %ld\n", stats->clippedNvalues, e0->n); 411 psVectorClipFitPolynomial2D (psf->params->data[PM_PAR_E1], stats, psfTry->mask, 0xff, e1, dz, x, y); 412 psTrace ("psModules.pmPSFtry", 4, "clipped E1 : keeping %ld of %ld\n", stats->clippedNvalues, e1->n); 413 psVectorClipFitPolynomial2D (psf->params->data[PM_PAR_E2], stats, psfTry->mask, 0xff, e2, dz, x, y); 414 psTrace ("psModules.pmPSFtry", 4, "clipped E2 : keeping %ld of %ld\n", stats->clippedNvalues, e2->n); 415 } 416 417 // XXX temporary dump of the psf parameters 407 for (int i = 0; i < psf->psfTrendStats->clipIter; i++) { 408 409 // XXX we are using the same stats structure on each pass: do we need to re-init it? 410 pmTrend2DFit (psf->params->data[PM_PAR_E0], psfTry->mask, 0xff, x, y, e0, dz); 411 psTrace ("psModules.pmPSFtry", 4, "clipped E0 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e0->n); 412 pmTrend2DFit (psf->params->data[PM_PAR_E1], psfTry->mask, 0xff, x, y, e1, dz); 413 psTrace ("psModules.pmPSFtry", 4, "clipped E1 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e1->n); 414 pmTrend2DFit (psf->params->data[PM_PAR_E2], psfTry->mask, 0xff, x, y, e2, dz); 415 psTrace ("psModules.pmPSFtry", 4, "clipped E2 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e2->n); 416 } 417 418 // test dump of the psf parameters 418 419 if (psTraceGetLevel("psModules.objects") >= 4) { 419 420 FILE *f = fopen ("pol.dat", "w"); 421 fprintf (f, "# x y : e0obs e1obs e2obs : e0fit e1fit e2fit : mask\n"); 420 422 for (int i = 0; i < e0->n; i++) { 421 fprintf (f, "%f %f : %f %f %f : % d\n",423 fprintf (f, "%f %f : %f %f %f : %f %f %f : %d\n", 422 424 x->data.F32[i], y->data.F32[i], 423 e0->data.F32[i], e1->data.F32[i], e2->data.F32[i], psfTry->mask->data.U8[i]); 425 e0->data.F32[i], e1->data.F32[i], e2->data.F32[i], 426 pmTrend2DEval (psf->params->data[PM_PAR_E0], x->data.F32[i], y->data.F32[i]), 427 pmTrend2DEval (psf->params->data[PM_PAR_E1], x->data.F32[i], y->data.F32[i]), 428 pmTrend2DEval (psf->params->data[PM_PAR_E2], x->data.F32[i], y->data.F32[i]), 429 psfTry->mask->data.U8[i]); 424 430 } 425 431 fclose (f); … … 455 461 // fit the collection of measured parameters to the PSF 2D model 456 462 // the mask is carried from previous steps and updated with this operation 457 // the weight is either the flux error or NULL, depending on ' applyWeights'458 if (!p sVectorClipFitPolynomial2D(psf->params->data[i], stats, psfTry->mask, 0xff, z, NULL, x, y)) {463 // the weight is either the flux error or NULL, depending on 'psf->poissonErrorParams' 464 if (!pmTrend2DFit (psf->params->data[i], psfTry->mask, 0xff, x, y, e0, NULL)) { 459 465 psError(PS_ERR_UNKNOWN, false, "failed to build psf model for parameter %d", i); 460 psFree(stats);461 466 psFree(x); 462 467 psFree(y); … … 485 490 486 491 for (int i = 0; i < psf->params->n; i++) { 487 if (psf->params->data[i] == NULL) 488 continue; 492 if (psf->params->data[i] == NULL) continue; 489 493 fprintf (f, "%f %f : ", model->params->data.F32[i], modelPSF->params->data.F32[i]); 490 494 } … … 496 500 } 497 501 498 psFree (stats);499 502 psFree (x); 500 503 psFree (y); 501 504 psFree (z); 502 505 psFree (dz); 503 return (true);506 return true; 504 507 } 505 508 -
trunk/psModules/src/objects/pmPSFtry.h
r14505 r15000 6 6 * @author EAM, IfA 7 7 * 8 * @version $Revision: 1.1 4$ $Name: not supported by cvs2svn $9 * @date $Date: 2007-0 8-15 20:21:18 $8 * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2007-09-24 21:27:58 $ 10 10 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii 11 11 */ … … 78 78 * 79 79 */ 80 pmPSFtry *pmPSFtryAlloc (psArray *sources, char *modelName, bool poissonErrors, psPolynomial2D *psfTrendMask);80 pmPSFtry *pmPSFtryAlloc (psArray *sources, pmPSFOptions *options); 81 81 82 82 /** pmPSFtryModel() … … 87 87 * 88 88 */ 89 pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, float RADIUS, bool poissonErrors, psPolynomial2D *psfTrendMask, bool applyWeights, psMaskType maskVal, psMaskType mark);89 pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, pmPSFOptions *options, psMaskType maskVal, psMaskType mark); 90 90 91 91 /** pmPSFtryMetric() … … 123 123 * 124 124 */ 125 bool pmPSFFromPSFtry (pmPSFtry *psfTry , bool applyWeights);125 bool pmPSFFromPSFtry (pmPSFtry *psfTry); 126 126 127 127 /// @} -
trunk/psModules/src/objects/pmSourcePhotometry.c
r14962 r15000 3 3 * @author EAM, IfA; GLG, MHPCC 4 4 * 5 * @version $Revision: 1.3 1$ $Name: not supported by cvs2svn $6 * @date $Date: 2007-09-2 1 02:46:46$5 * @version $Revision: 1.32 $ $Name: not supported by cvs2svn $ 6 * @date $Date: 2007-09-24 21:27:58 $ 7 7 * 8 8 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 81 81 // we must have a valid model 82 82 model = pmSourceGetModel (&isPSF, source); 83 if (model == NULL) 84 return false; 83 if (model == NULL) { 84 psTrace ("psModules.objects", 3, "fail mag : no valid model"); 85 return false; 86 } 85 87 86 88 if (model->dparams->data.F32[PM_PAR_I0] > 0) { … … 102 104 source->psfMag = NAN; 103 105 } 104 fprintf (stderr, ".");105 106 } else { 106 107 status = pmSourcePhotometryModel (&source->psfMag, source->modelPSF); … … 109 110 // measure EXT model photometry 110 111 status = pmSourcePhotometryModel (&source->extMag, source->modelEXT); 112 111 113 // for PSFs, correct both apMag and psfMag to same system, consistent with infinite flux star in aperture RADIUS 112 if ((mode & PM_SOURCE_PHOT_APCORR) && isPSF && psf) { 113 // convert to the equivalent 2D model? 114 if ((mode & PM_SOURCE_PHOT_APCORR) && isPSF && psf && psf->ApTrend) { 114 115 source->psfMag += pmTrend2DEval (psf->ApTrend, x, y); 115 116 } 116 117 117 if (!isfinite(SN) || (SN < AP_MIN_SN)) 118 return false; 118 if (!isfinite(SN) || (SN < AP_MIN_SN)) { 119 psTrace ("psModules.objects", 3, "fail mag : bad SN: %f (limit: %f)", SN, AP_MIN_SN); 120 return false; 121 } 119 122 120 123 // replace source flux … … 209 212 status = pmSourcePhotometryAper (&source->apMag, model, flux, source->maskObj, maskVal); 210 213 if (!status) { 214 psTrace ("psModules.objects", 3, "fail mag : bad Ap Mag"); 211 215 psErrorCode last = psErrorCodeLast(); 212 216 if (last == PM_ERR_PHOTOM) { -
trunk/psModules/src/objects/pmTrend2D.c
r14938 r15000 3 3 * @author EAM, IfA 4 4 * 5 * @version $Revision: 1. 2$ $Name: not supported by cvs2svn $6 * @date $Date: 2007-09-2 1 00:09:18 $5 * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ 6 * @date $Date: 2007-09-24 21:27:58 $ 7 7 * 8 8 * Copyright 2004 Institute for Astronomy, University of Hawaii … … 14 14 #endif 15 15 16 # include <strings.h> 16 17 # include <pslib.h> 17 18 # include "pmTrend2D.h" … … 31 32 { 32 33 assert (image); 33 assert (stats);34 34 35 35 pmTrend2D *trend = (pmTrend2D *) psAlloc(sizeof(pmTrend2D)); … … 44 44 case PM_TREND_POLY_ORD: 45 45 trend->poly = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, nXtrend, nYtrend); 46 // set masking somehow 47 for (int nx = 0; nx < trend->poly->nX + 1; nx++) { 48 for (int ny = 0; ny < trend->poly->nY + 1; ny++) { 49 if (nx + ny >= PS_MAX (trend->poly->nX, trend->poly->nY) + 1) { 50 trend->poly->mask[nx][ny] = 1; 51 } else { 52 trend->poly->mask[nx][ny] = 0; 53 } 54 } 55 } 46 56 break; 47 57 … … 69 79 } 70 80 71 pmTrend2D *pmTrend2D FieldAlloc (pmTrend2DMode mode, int nXfield, int nYfield, int nXtrend, int nYtrend, psStats *stats)81 pmTrend2D *pmTrend2DNoImageAlloc (pmTrend2DMode mode, psImageBinning *binning, psStats *stats) 72 82 { 73 83 pmTrend2D *trend = (pmTrend2D *) psAlloc(sizeof(pmTrend2D)); … … 81 91 switch (mode) { 82 92 case PM_TREND_POLY_ORD: 93 trend->poly = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, binning->nXruff, binning->nYruff); 94 // set masking somehow 95 for (int nx = 0; nx < trend->poly->nX + 1; nx++) { 96 for (int ny = 0; ny < trend->poly->nY + 1; ny++) { 97 if (nx + ny >= PS_MAX (trend->poly->nX, trend->poly->nY) + 1) { 98 trend->poly->mask[nx][ny] = 1; 99 } else { 100 trend->poly->mask[nx][ny] = 0; 101 } 102 } 103 } 104 break; 105 106 case PM_TREND_POLY_CHEB: 107 trend->poly = psPolynomial2DAlloc (PS_POLYNOMIAL_CHEB, binning->nXruff, binning->nYruff); 108 break; 109 110 case PM_TREND_MAP: { 111 // binning defines the map scale relationship 112 trend->map = psImageMapNoImageAlloc (binning, stats); 113 break; 114 } 115 116 default: 117 psAbort ("error"); 118 } 119 return (trend); 120 } 121 122 pmTrend2D *pmTrend2DFieldAlloc (pmTrend2DMode mode, int nXfield, int nYfield, int nXtrend, int nYtrend, psStats *stats) 123 { 124 pmTrend2D *trend = (pmTrend2D *) psAlloc(sizeof(pmTrend2D)); 125 psMemSetDeallocator(trend, (psFreeFunc) pmTrend2DFree); 126 127 trend->map = NULL; 128 trend->poly = NULL; 129 trend->stats = psMemIncrRefCounter (stats); 130 trend->mode = mode; 131 132 switch (mode) { 133 case PM_TREND_POLY_ORD: 83 134 trend->poly = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, nXtrend, nYtrend); 135 // set masking somehow 136 for (int nx = 0; nx < trend->poly->nX + 1; nx++) { 137 for (int ny = 0; ny < trend->poly->nY + 1; ny++) { 138 if (nx + ny >= PS_MAX (trend->poly->nX, trend->poly->nY) + 1) { 139 trend->poly->mask[nx][ny] = 1; 140 } else { 141 trend->poly->mask[nx][ny] = 0; 142 } 143 } 144 } 84 145 break; 85 146 … … 111 172 bool status; 112 173 174 assert (trend); 175 assert (x); 176 assert (y); 177 assert (f); 178 113 179 switch (trend->mode) { 114 180 case PM_TREND_POLY_ORD: … … 136 202 double result; 137 203 138 if (!trend) return 0.0;204 assert (trend); 139 205 140 206 switch (trend->mode) { … … 158 224 psVector *result; 159 225 226 assert (trend); 227 assert (x); 228 assert (y); 229 160 230 switch (trend->mode) { 161 231 case PM_TREND_POLY_ORD: … … 173 243 return result; 174 244 } 245 246 psString pmTrend2DModeToString (pmTrend2DMode mode) { 247 248 psString name; 249 250 switch (mode) { 251 case PM_TREND_NONE: 252 name = psStringCopy ("NONE"); 253 break; 254 case PM_TREND_POLY_ORD: 255 name = psStringCopy ("POLY_ORD"); 256 break; 257 case PM_TREND_POLY_CHEB: 258 name = psStringCopy ("POLY_CHEB"); 259 break; 260 case PM_TREND_MAP: 261 name = psStringCopy ("MAP"); 262 break; 263 default: 264 psAbort ("invalid mode %d", mode); 265 } 266 return name; 267 } 268 269 pmTrend2DMode pmTrend2DModeFromString (psString name) { 270 271 if (!name) return PM_TREND_NONE; 272 273 if (!strcasecmp (name, "NONE")) { 274 return PM_TREND_NONE; 275 } 276 if (!strcasecmp (name, "POLY_ORD")) { 277 return PM_TREND_POLY_ORD; 278 } 279 if (!strcasecmp (name, "POLY_CHEB")) { 280 return PM_TREND_POLY_CHEB; 281 } 282 if (!strcasecmp (name, "MAP")) { 283 return PM_TREND_MAP; 284 } 285 psError (PS_ERR_UNKNOWN, true, "Unknown pmTrend2D mode %s\n", name); 286 return PM_TREND_NONE; 287 } -
trunk/psModules/src/objects/pmTrend2D.h
r14938 r15000 5 5 * @author EAM, IfA 6 6 * 7 * @version $Revision: 1. 2$ $Name: not supported by cvs2svn $8 * @date $Date: 2007-09-2 1 00:09:18 $7 * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2007-09-24 21:27:58 $ 9 9 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii 10 10 */ … … 17 17 18 18 typedef enum { 19 PM_TREND_NONE, 19 20 PM_TREND_POLY_ORD, 20 21 PM_TREND_POLY_CHEB, … … 33 34 pmTrend2D *pmTrend2DAlloc (pmTrend2DMode mode, psImage *image, int nXtrend, int nYtrend, psStats *stats); 34 35 36 pmTrend2D *pmTrend2DNoImageAlloc (pmTrend2DMode mode, psImageBinning *binning, psStats *stats); 37 35 38 // allocate a pmTrend2D tied to an abstract field with size nXfield,nYfield 36 39 pmTrend2D *pmTrend2DFieldAlloc (pmTrend2DMode mode, int nXfield, int nYfield, int nXtrend, int nYtrend, psStats *stats); … … 41 44 psVector *pmTrend2DEvalVector (pmTrend2D *trend, psVector *x, psVector *y); 42 45 46 psString pmTrend2DModeToString (pmTrend2DMode mode); 47 pmTrend2DMode pmTrend2DModeFromString (psString name); 48 43 49 /// @} 44 50 # endif
Note:
See TracChangeset
for help on using the changeset viewer.
