Changeset 35769 for trunk/psphot/src/psphotExtendedSourceFits.c
- Timestamp:
- Jul 3, 2013, 2:43:13 PM (13 years ago)
- Location:
- trunk/psphot
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/psphot
- Property svn:mergeinfo changed
/branches/eam_branches/ipp-20130509/psphot (added) merged: 35594,35614,35630-35634,35640-35641,35650-35652,35654-35655,35659,35661,35701,35751
- Property svn:mergeinfo changed
-
trunk/psphot/src
- Property svn:mergeinfo changed
/branches/eam_branches/ipp-20130509/psphot/src (added) merged: 35594,35614,35630-35634,35640,35651,35654,35659,35661,35751
- Property svn:mergeinfo changed
-
trunk/psphot/src/psphotExtendedSourceFits.c
r35559 r35769 1 1 # include "psphotInternal.h" 2 bool psphotSetWindowTrail (float *fitRadius, float *windowRadius, pmReadout *readout, pmSource *source, psImageMaskType markVal, float newRadius); 2 3 3 4 // for now, let's store the detections on the readout->analysis for each readout … … 52 53 psTimerStart ("psphot.extended"); 53 54 55 psphotFitInitExtended(); 56 54 57 // find the currently selected readout 55 58 pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest … … 70 73 if (!sources->n) { 71 74 psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size"); 75 psphotSersicModelClassCleanup(); 72 76 return true; 73 77 } … … 80 84 // do not thread if we are trying to study the fitting process 81 85 if (psTraceGetLevel ("psphot.psphotFitEXT") >= 6) { 82 86 nThreads = 0; 83 87 } 84 88 … … 107 111 } 108 112 113 bool chisqConvergence = psMetadataLookupBool (&status, recipe, "LMM_FIT_CHISQ_CONVERGENCE"); // Fit tolerance 114 if (!status) { 115 // default to the old method (chisqConvergence) 116 chisqConvergence = true; 117 } 118 119 int gainFactorMode = psMetadataLookupS32 (&status, recipe, "LMM_FIT_GAIN_FACTOR_MODE"); // Fit tolerance 120 if (!status) { 121 // default to the old method (chisqConvergence) 122 gainFactorMode = 0; 123 } 124 125 // Define source fitting parameters for extended source fits 126 pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc(); 127 fitOptions->mode = PM_SOURCE_FIT_EXT; 128 fitOptions->saveCovariance = true; // XXX make this a user option? 129 fitOptions->covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix 130 fitOptions->nIter = fitIter; 131 fitOptions->minTol = fitMinTol; 132 fitOptions->maxTol = fitMaxTol; 133 134 fitOptions->gainFactorMode = gainFactorMode; 135 fitOptions->chisqConvergence = chisqConvergence; 136 109 137 // maskVal is used to test for rejected pixels, and must include markVal 110 138 maskVal |= markVal; 111 139 112 140 // select the collection of desired models 113 psMetadata * models = psMetadataLookupMetadata (&status, recipe, "EXTENDED_SOURCE_MODELS");141 psMetadata *allModels = psMetadataLookupMetadata (&status, recipe, "EXTENDED_SOURCE_MODELS"); 114 142 if (!status) { 115 143 psWarning ("extended source model fits requested but model model is missing (EXTENDED_SOURCE_MODELS)\n"); 144 psFree (fitOptions); 145 psphotSersicModelClassCleanup(); 116 146 return true; 117 147 } 148 if (allModels->list->n == 0) { 149 psWarning ("extended source model fits requested but no models are specified\n"); 150 psFree (fitOptions); 151 psphotSersicModelClassCleanup(); 152 return true; 153 } 154 155 psMetadata *models = NULL; 156 char *modelSelection = psMetadataLookupStr (&status, recipe, "EXTENDED_SOURCE_MODELS_SELECTION"); 157 if (!modelSelection || !status) { 158 models = psMemIncrRefCounter(allModels); 159 goto got_models; 160 } 161 162 if (!strcasecmp(modelSelection, "ALL")) { 163 models = psMemIncrRefCounter(allModels); 164 goto got_models; 165 } 166 167 if (!strcasecmp(modelSelection, "NONE")) { 168 psWarning ("extended source model fits requested but no models are selected (EXTENDED_SOURCE_MODELS_SELECTION = NONE)\n"); 169 psFree (fitOptions); 170 psphotSersicModelClassCleanup(); 171 return true; 172 } 173 174 psArray *selection = psStringSplitArray (modelSelection, ",", false); 175 if (selection->n == 0) { 176 psWarning ("extended source model fits requested but model selection string is empty (EXTENDED_SOURCE_MODELS_SELECTION)\n"); 177 psFree (fitOptions); 178 psFree (selection); 179 psphotSersicModelClassCleanup(); 180 return true; 181 } 182 models = psMetadataAlloc(); 183 for (int i = 0; i < selection->n; i++) { 184 psMetadata *model = psMetadataLookupMetadata (&status, allModels, selection->data[i]); 185 if (!model) { 186 psWarning ("extended source model selection string includes an invalid name %s\n", (char *) selection->data[i]); 187 continue; 188 } 189 psMetadataAddMetadata (models, PS_LIST_TAIL, selection->data[i], PS_META_REPLACE, "extended model", model); 190 } 191 psFree (selection); 192 118 193 if (models->list->n == 0) { 119 psWarning ("extended source model fits requested but no models are specified\n"); 120 return true; 121 } 194 psWarning ("extended source model fits requested but no valid models in selection string (EXTENDED_SOURCE_MODELS_SELECTION = %s)\n", modelSelection); 195 psFree (fitOptions); 196 psFree (models); 197 psphotSersicModelClassCleanup(); 198 return true; 199 } 200 201 got_models: 122 202 123 203 psphotInitRadiusEXT (recipe, readout); … … 128 208 while ((item = psMetadataGetAndIncrement (iter)) != NULL) { 129 209 130 if (item->type != PS_DATA_METADATA) { 131 // XXX we could cull the bad entries or build a validated model folder 132 psAbort ("Invalid type for EXTENDED_SOURCE_MODEL entry %s, not a metadata folder", item->name); 133 } 134 135 psMetadata *model = (psMetadata *) item->data.md; 136 137 // check on the model type 138 char *modelName = psMetadataLookupStr (&status, model, "MODEL"); 139 int modelType = pmModelClassGetType (modelName); 140 if (modelType < 0) { 141 psAbort ("Unknown model class for EXTENDED_SOURCE_MODEL entry %s: %s", item->name, modelName); 142 } 143 psMetadataAddS32 (model, PS_LIST_TAIL, "MODEL_TYPE", PS_META_REPLACE, "", modelType); 144 145 // check on the SNLIM, set a float value 146 char *SNword = psMetadataLookupStr (&status, model, "SNLIM"); 147 if (!status) { 148 psAbort("SNLIM not defined for extended source model %s\n", item->name); 149 } 150 float SNlim = atof (SNword); 151 psMetadataAddF32 (model, PS_LIST_TAIL, "SNLIM_VALUE", PS_META_REPLACE, "", SNlim); 152 153 // check on the PSF-Convolution status 154 char *convolvedWord = psMetadataLookupStr (&status, model, "PSF_CONVOLVED"); 155 if (!status || (strcasecmp (convolvedWord, "true") && strcasecmp (convolvedWord, "false"))) { 156 psAbort ("PSF_CONVOLVED entry invalid or missing for EXTENDED_SOURCE_MODEL entry %s", item->name); 157 } 158 bool convolved = !strcasecmp (convolvedWord, "true"); 159 psMetadataAddBool (model, PS_LIST_TAIL, "PSF_CONVOLVED_VALUE", PS_META_REPLACE, "", convolved); 210 if (item->type != PS_DATA_METADATA) { 211 // XXX we could cull the bad entries or build a validated model folder 212 psAbort ("Invalid type for EXTENDED_SOURCE_MODEL entry %s, not a metadata folder", item->name); 213 } 214 215 psMetadata *model = (psMetadata *) item->data.md; 216 217 // check on the model type 218 char *modelName = psMetadataLookupStr (&status, model, "MODEL"); 219 int modelType = pmModelClassGetType (modelName); 220 if (modelType < 0) { 221 psAbort ("Unknown model class for EXTENDED_SOURCE_MODEL entry %s: %s", item->name, modelName); 222 } 223 psMetadataAddS32 (model, PS_LIST_TAIL, "MODEL_TYPE", PS_META_REPLACE, "", modelType); 224 225 // check on the SNLIM, set a float value 226 char *SNword = psMetadataLookupStr (&status, model, "SNLIM"); 227 if (!status) { 228 psAbort("SNLIM not defined for extended source model %s\n", item->name); 229 } 230 float SNlim = atof (SNword); 231 psMetadataAddF32 (model, PS_LIST_TAIL, "SNLIM_VALUE", PS_META_REPLACE, "", SNlim); 232 233 // check on the PSF-Convolution status 234 char *convolvedWord = psMetadataLookupStr (&status, model, "PSF_CONVOLVED"); 235 if (!status || (strcasecmp (convolvedWord, "true") && strcasecmp (convolvedWord, "false"))) { 236 psAbort ("PSF_CONVOLVED entry invalid or missing for EXTENDED_SOURCE_MODEL entry %s", item->name); 237 } 238 bool convolved = !strcasecmp (convolvedWord, "true"); 239 psMetadataAddBool (model, PS_LIST_TAIL, "PSF_CONVOLVED_VALUE", PS_META_REPLACE, "", convolved); 240 241 if (convolved) { 242 psLogMsg ("psphot", PS_LOG_INFO, "using convolved model class %s (%s)", modelName, item->name); 243 } else { 244 psLogMsg ("psphot", PS_LOG_INFO, "using simple model class %s (%s)", modelName, item->name); 245 } 160 246 } 161 247 psFree (iter); … … 195 281 psMetadataIterator *iter = psMetadataIteratorAlloc (models, PS_LIST_HEAD, NULL); 196 282 psArrayAdd(job->args, 1, iter); 197 psArrayAdd(job->args, 1, AnalysisRegion); // XXX make a pointer 283 psArrayAdd(job->args, 1, AnalysisRegion); 284 psArrayAdd(job->args, 1, fitOptions); 198 285 199 286 PS_ARRAY_ADD_SCALAR(job->args, psfSize, PS_TYPE_S32); 200 PS_ARRAY_ADD_SCALAR(job->args, fitIter, PS_TYPE_S32);201 PS_ARRAY_ADD_SCALAR(job->args, fitMinTol, PS_TYPE_F32);202 PS_ARRAY_ADD_SCALAR(job->args, fitMaxTol, PS_TYPE_F32);203 204 287 PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); 205 288 PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK); … … 218 301 psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); 219 302 psFree(AnalysisRegion); 303 psFree (fitOptions); 304 psFree (models); 305 psphotSersicModelClassCleanup(); 220 306 return false; 221 307 } … … 224 310 psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); 225 311 psFree(AnalysisRegion); 312 psFree (fitOptions); 313 psFree (models); 314 psphotSersicModelClassCleanup(); 226 315 return false; 227 316 } 228 317 psScalar *scalar = NULL; 229 scalar = job->args->data[ 8];318 scalar = job->args->data[9]; 230 319 Next += scalar->data.S32; 231 scalar = job->args->data[ 9];320 scalar = job->args->data[10]; 232 321 Nconvolve += scalar->data.S32; 233 scalar = job->args->data[1 0];322 scalar = job->args->data[11]; 234 323 NconvolvePass += scalar->data.S32; 235 scalar = job->args->data[1 1];324 scalar = job->args->data[12]; 236 325 Nplain += scalar->data.S32; 237 scalar = job->args->data[1 2];326 scalar = job->args->data[13]; 238 327 NplainPass += scalar->data.S32; 239 scalar = job->args->data[1 3];328 scalar = job->args->data[14]; 240 329 Nfaint += scalar->data.S32; 241 scalar = job->args->data[1 4];330 scalar = job->args->data[15]; 242 331 Nfail += scalar->data.S32; 243 332 psFree(job->args->data[3]); // iterator allocated above … … 250 339 psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); 251 340 psFree(AnalysisRegion); 341 psFree (fitOptions); 342 psFree (models); 343 psphotSersicModelClassCleanup(); 252 344 return false; 253 345 } … … 260 352 } else { 261 353 psScalar *scalar = NULL; 262 scalar = job->args->data[ 8];354 scalar = job->args->data[9]; 263 355 Next += scalar->data.S32; 264 scalar = job->args->data[ 9];356 scalar = job->args->data[10]; 265 357 Nconvolve += scalar->data.S32; 266 scalar = job->args->data[1 0];358 scalar = job->args->data[11]; 267 359 NconvolvePass += scalar->data.S32; 268 scalar = job->args->data[1 1];360 scalar = job->args->data[12]; 269 361 Nplain += scalar->data.S32; 270 scalar = job->args->data[1 2];362 scalar = job->args->data[13]; 271 363 NplainPass += scalar->data.S32; 272 scalar = job->args->data[1 3];364 scalar = job->args->data[14]; 273 365 Nfaint += scalar->data.S32; 274 scalar = job->args->data[1 4];366 scalar = job->args->data[15]; 275 367 Nfail += scalar->data.S32; 276 368 psFree(job->args->data[3]); // metadata iterator allocated above … … 281 373 psFree (cellGroups); 282 374 psFree(AnalysisRegion); 375 psFree (fitOptions); 376 psFree (models); 283 377 284 378 psphotSersicModelClassCleanup(); … … 290 384 psLogMsg ("psphot", PS_LOG_INFO, " %d plain models (%d passed)\n", Nplain, NplainPass); 291 385 psLogMsg ("psphot", PS_LOG_INFO, " %d too faint to fit, %d failed\n", Nfaint, Nfail); 386 387 psphotFitSummaryExtended(); 388 292 389 return true; 293 390 } … … 308 405 309 406 // arguments: readout, sources, models, region, psfSize, maskVal, markVal 310 pmReadout *readout = job->args->data[0];311 psArray *sources = job->args->data[1];312 psMetadata *models = job->args->data[2];407 pmReadout *readout = job->args->data[0]; 408 psArray *sources = job->args->data[1]; 409 psMetadata *models = job->args->data[2]; 313 410 psMetadataIterator *iter = job->args->data[3]; 314 psRegion *region = job->args->data[4]; 315 int psfSize = PS_SCALAR_VALUE(job->args->data[5],S32); 316 int fitIter = PS_SCALAR_VALUE(job->args->data[6],S32); 317 float fitMinTol = PS_SCALAR_VALUE(job->args->data[7],F32); 318 float fitMaxTol = PS_SCALAR_VALUE(job->args->data[8],F32); 319 psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[9],PS_TYPE_IMAGE_MASK_DATA); 320 psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[10],PS_TYPE_IMAGE_MASK_DATA); 321 322 // Define source fitting parameters for extended source fits 323 pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc(); 324 fitOptions->mode = PM_SOURCE_FIT_EXT; 325 fitOptions->saveCovariance = true; // XXX make this a user option? 326 fitOptions->covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix 327 fitOptions->nIter = fitIter; 328 fitOptions->minTol = fitMinTol; 329 fitOptions->maxTol = fitMaxTol; 411 psRegion *region = job->args->data[4]; 412 pmSourceFitOptions *fitOptions = job->args->data[5]; 413 414 int psfSize = PS_SCALAR_VALUE(job->args->data[6],S32); 415 psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[7],PS_TYPE_IMAGE_MASK_DATA); 416 psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[8],PS_TYPE_IMAGE_MASK_DATA); 417 418 // psTraceSetLevel ("psLib.math.psMinimizeLMChi2_Alt", 5); 330 419 331 420 // choose the sources of interest … … 348 437 if (source->peak->y > region->y1) continue; 349 438 439 440 // XXX for a test, just do the obvious trail 441 // XXX if (source->peak->xf < 1100) continue; 442 // XXX if (source->peak->xf > 1400) continue; 443 // XXX if (source->peak->yf > 245) continue; 444 350 445 // replace object in image 351 446 if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { … … 384 479 assert (status); 385 480 386 // limit selection to some SN limit387 // assert (source->peak); // how can a source not have a peak?388 481 // limit selection to some SN limit 389 482 bool skipSource = false; … … 424 517 } 425 518 } else { 426 psFree (source->modelFlux); 427 source->modelFlux = NULL; 428 modelFit = psphotFitEXT (readout, source, fitOptions, modelType, maskVal, markVal); 429 if (!modelFit) { 430 psTrace ("psphot", 5, "failed to fit plain model for object at %f, %f", source->moments->Mx, source->moments->My); 431 Nfail ++; 432 continue; 433 } 434 psTrace ("psphot", 4, "fit plain model for %f, %f : %s chisq = %f (npix: %d, niter: %d)\n", 435 source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq, modelFit->nPix, modelFit->nIter); 436 Nplain ++; 437 if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) { 438 NplainPass ++; 439 source->mode |= PM_SOURCE_MODE_EXTENDED_FIT; 440 } 519 bool doneFits = false; 520 while (!doneFits) { 521 psFree (source->modelFlux); 522 source->modelFlux = NULL; 523 modelFit = psphotFitEXT (modelFit, readout, source, fitOptions, modelType, maskVal, markVal); 524 if (!modelFit) { 525 psTrace ("psphot", 5, "failed to fit plain model for object at %f, %f", source->moments->Mx, source->moments->My); 526 Nfail ++; 527 doneFits = true; 528 continue; 529 } 530 psTrace ("psphot", 4, "fit plain model for %f, %f : %s chisq = %f (npix: %d, niter: %d)\n", source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq, modelFit->nPix, modelFit->nIter); 531 Nplain ++; 532 if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) { 533 NplainPass ++; 534 source->mode |= PM_SOURCE_MODE_EXTENDED_FIT; 535 } 536 doneFits = true; 537 538 if (modelType == pmModelClassGetType("PS_MODEL_TRAIL")) { 539 // calculate the object end points: 540 float *PAR = modelFit->params->data.F32; 541 float Xmin = PAR[PM_PAR_XPOS] - 0.5*PAR[PM_PAR_LENGTH]*sin(PAR[PM_PAR_THETA]); 542 float Ymin = PAR[PM_PAR_YPOS] - 0.5*PAR[PM_PAR_LENGTH]*cos(PAR[PM_PAR_THETA]); 543 float Xmax = PAR[PM_PAR_XPOS] + 0.5*PAR[PM_PAR_LENGTH]*sin(PAR[PM_PAR_THETA]); 544 float Ymax = PAR[PM_PAR_YPOS] + 0.5*PAR[PM_PAR_LENGTH]*cos(PAR[PM_PAR_THETA]); 545 546 if (false && (source->peak->xf > 1100) && 547 (source->peak->xf < 1400) && 548 (source->peak->yf < 245)) { 549 fprintf (stderr, "src vs fit : %d %d - %d %d | %f %f - %f %f\n", 550 source->pixels->col0, source->pixels->row0, 551 source->pixels->col0 + source->pixels->numCols, source->pixels->row0 + source->pixels->numRows, 552 Xmin, Ymin, Xmax, Ymax); 553 } 554 if (PAR[PM_PAR_LENGTH] > 0.9*fitRadius) { 555 doneFits = false; 556 fprintf (stderr, "update window : %f %f : %f -> %f\n", source->peak->xf, source->peak->yf, fitRadius, 2*fitRadius); 557 psphotSetWindowTrail (&fitRadius, &windowRadius, readout, source, markVal, fitRadius*2.0); 558 } 559 } 560 } 441 561 } 442 443 // XXX deprecate? 444 // XXX pmSourceExtFitPars *extFitPars = pmSourceExtFitParsAlloc(); 445 // XXX extFitPars->Mxx = source->moments->Mxx; 446 // XXX extFitPars->Mxy = source->moments->Mxy; 447 // XXX extFitPars->Myy = source->moments->Myy; 448 // XXX extFitPars->Mrf = source->moments->Mrf; 449 // XXX extFitPars->Mrh = source->moments->Mrh; 450 // XXX extFitPars->peakMag = (source->peak->rawFlux > 0) ? -2.5*log10(source->peak->rawFlux) : NAN; 451 452 // save kron mag, but assign apMag & psfMag on output (not yet calculated) 453 // XXX extFitPars->krMag = -2.5*log10(source->moments->KronFlux); 454 455 // psArrayAdd (source->extFitPars, 4, extFitPars); 456 // psFree (extFitPars); 562 // XXX really need to do this in a cleaner way: 563 if (!modelFit) continue; 457 564 458 565 // test for fit quality / result … … 519 626 psTrace ("psphot", 5, "extended source model for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); 520 627 } 521 psFree (fitOptions);522 628 523 629 // change the value of a scalar on the array (wrap this and put it in psArray.h) 630 scalar = job->args->data[9]; 631 scalar->data.S32 = Next; 632 633 scalar = job->args->data[10]; 634 scalar->data.S32 = Nconvolve; 635 524 636 scalar = job->args->data[11]; 525 scalar->data.S32 = N ext;637 scalar->data.S32 = NconvolvePass; 526 638 527 639 scalar = job->args->data[12]; 528 scalar->data.S32 = N convolve;640 scalar->data.S32 = Nplain; 529 641 530 642 scalar = job->args->data[13]; 531 scalar->data.S32 = N convolvePass;643 scalar->data.S32 = NplainPass; 532 644 533 645 scalar = job->args->data[14]; 534 scalar->data.S32 = N plain;646 scalar->data.S32 = Nfaint; 535 647 536 648 scalar = job->args->data[15]; 537 scalar->data.S32 = NplainPass;538 539 scalar = job->args->data[16];540 scalar->data.S32 = Nfaint;541 542 scalar = job->args->data[17];543 649 scalar->data.S32 = Nfail; 544 650 545 651 return true; 546 652 } 653 654 # define PAD_WINDOW 3.0 655 bool psphotSetWindowTrail (float *fitRadius, float *windowRadius, pmReadout *readout, pmSource *source, psImageMaskType markVal, float newRadius) { 656 657 psRegion newRegion; 658 659 psAssert (source, "source not defined??"); 660 psAssert (source->moments, "moments not defined??"); 661 662 *fitRadius = newRadius; 663 *windowRadius = *fitRadius + PAD_WINDOW; 664 665 // check to see if new region is completely contained within old region 666 newRegion = psRegionForSquare (source->peak->xf, source->peak->yf, *windowRadius); 667 newRegion = psRegionForImage (readout->image, newRegion); 668 669 // redefine the pixels to match 670 pmSourceRedefinePixelsByRegion (source, readout, newRegion); 671 672 // clear the circular mask 673 psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); 674 675 // set the mask to flag the excluded pixels 676 psImageKeepCircle (source->maskObj, source->peak->xf, source->peak->yf, *fitRadius, "OR", markVal); 677 678 return true; 679 }
Note:
See TracChangeset
for help on using the changeset viewer.
