Changeset 18847
- Timestamp:
- Aug 1, 2008, 8:39:03 AM (18 years ago)
- File:
-
- 1 edited
-
trunk/psphot/src/psphotExtendedSourceFits.c (modified) (8 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psphot/src/psphotExtendedSourceFits.c
r18555 r18847 24 24 // perform full extended source non-linear fits? 25 25 if (!psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_FITS")) { 26 psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source measurements\n");27 return true;26 psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source measurements\n"); 27 return true; 28 28 } 29 29 … … 31 31 psMetadata *models = psMetadataLookupMetadata (&status, recipe, "EXTENDED_SOURCE_MODELS"); 32 32 if (!status) { 33 psWarning ("extended source model fits requested but model model is missing (EXTENDED_SOURCE_MODELS)\n");34 return true;33 psWarning ("extended source model fits requested but model model is missing (EXTENDED_SOURCE_MODELS)\n"); 34 return true; 35 35 } 36 36 if (models->list->n == 0) { 37 psWarning ("extended source model fits requested but no models are specified\n");38 return true;37 psWarning ("extended source model fits requested but no models are specified\n"); 38 return true; 39 39 } 40 40 … … 45 45 46 46 if (item->type != PS_DATA_METADATA) { 47 psAbort ("Invalid type for EXTENDED_SOURCE_MODEL entry %s, not a metadata folder", item->name);48 // XXX we could cull the bad entries or build a validated model folder47 psAbort ("Invalid type for EXTENDED_SOURCE_MODEL entry %s, not a metadata folder", item->name); 48 // XXX we could cull the bad entries or build a validated model folder 49 49 } 50 50 … … 55 55 int modelType = pmModelClassGetType (modelName); 56 56 if (modelType < 0) { 57 psAbort ("Unknown model class for EXTENDED_SOURCE_MODEL entry %s: %s", item->name, modelName);58 } 57 psAbort ("Unknown model class for EXTENDED_SOURCE_MODEL entry %s: %s", item->name, modelName); 58 } 59 59 psMetadataAddS32 (model, PS_LIST_TAIL, "MODEL_TYPE", PS_META_REPLACE, "", modelType); 60 60 … … 62 62 char *SNword = psMetadataLookupStr (&status, model, "SNLIM"); 63 63 if (!status) { 64 psAbort("SNLIM not defined for extended source model %s\n", item->name);64 psAbort("SNLIM not defined for extended source model %s\n", item->name); 65 65 } 66 66 float SNlim = atof (SNword); 67 67 psMetadataAddF32 (model, PS_LIST_TAIL, "SNLIM_VALUE", PS_META_REPLACE, "", SNlim); 68 68 69 69 // check on the PSF-Convolution status 70 70 char *convolvedWord = psMetadataLookupStr (&status, model, "PSF_CONVOLVED"); 71 71 if (!status || (strcasecmp (convolvedWord, "true") && strcasecmp (convolvedWord, "false"))) { 72 psAbort ("PSF_CONVOLVED entry invalid or missing for EXTENDED_SOURCE_MODEL entry %s", item->name);73 } 72 psAbort ("PSF_CONVOLVED entry invalid or missing for EXTENDED_SOURCE_MODEL entry %s", item->name); 73 } 74 74 bool convolved = !strcasecmp (convolvedWord, "true"); 75 75 psMetadataAddBool (model, PS_LIST_TAIL, "PSF_CONVOLVED_VALUE", PS_META_REPLACE, "", convolved); … … 106 106 107 107 // if model is NULL, we don't have a starting guess 108 // XXX use the parameters guessed from moments108 // XXX use the parameters guessed from moments 109 109 // if (source->modelEXT == NULL) continue; 110 110 … … 115 115 Next ++; 116 116 117 // save the modelFlux here in case we need to subtract it (for failure) 118 psImage *modelFluxStart = psMemIncrRefCounter (source->modelFlux); 119 120 if (savePics) { 121 psphotSaveImage (NULL, readout->image, "image.xp.fits"); 122 } 123 124 // array to store the pointers to the model flux images while the models are being fitted 125 psArray *modelFluxes = psArrayAllocEmpty (4); 126 127 // allocate the array to store the model fits 128 if (source->modelFits == NULL) { 129 source->modelFits = psArrayAllocEmpty (4); 130 } 131 132 // loop here over the models chosen for each source (exclude by S/N) 133 psMetadataIterator *iter = psMetadataIteratorAlloc (models, PS_LIST_HEAD, NULL); 134 psMetadataItem *item = NULL; 135 while ((item = psMetadataGetAndIncrement (iter)) != NULL) { 136 137 // XXX this should have been forced above 138 assert (item->type == PS_DATA_METADATA); 139 psMetadata *model = (psMetadata *) item->data.md; 140 141 // check the SNlim and skip model if source is too faint 142 float SNlim = psMetadataLookupF32 (&status, model, "SNLIM_VALUE"); 143 assert (status); 144 145 // limit selection to some SN limit 146 assert (source->peak); // how can a source not have a peak? 147 if (source->peak->SN < SNlim) continue; 148 149 // check on the model type 150 pmModelType modelType = psMetadataLookupS32 (&status, model, "MODEL_TYPE"); 151 assert (status); 152 153 // check on the PSF-Convolution status 154 bool convolved = psMetadataLookupBool (&status, model, "PSF_CONVOLVED_VALUE"); 155 assert (status); 156 157 // XXX psTraceSetLevel ("psLib.math.psMinimizeLMChi2", 6); 158 // XXX psTraceSetLevel ("psphot.psphotModelWithPSF_LMM", 6); 159 160 // fit the model as convolved or not 161 pmModel *modelFit = NULL; 162 if (convolved) { 163 modelFit = psphotPSFConvModel (readout, source, modelType, maskVal, markVal, psfSize); 164 if (!modelFit) { 165 psTrace ("psphot", 5, "failed to fit psf-conv model for object at %f, %f", source->moments->x, source->moments->y); 166 continue; 167 } 168 psTrace ("psphot", 4, "fit psf-conv model for %f, %f : %s chisq = %f\n", source->moments->x, source->moments->y, pmModelClassGetName (modelFit->type), modelFit->chisq); 169 Nconvolve ++; 170 if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) { 171 NconvolvePass ++; 172 } 173 } else { 174 psFree (source->modelFlux); 175 source->modelFlux = NULL; 176 modelFit = psphotFitEXT (readout, source, modelType, maskVal, markVal); 177 if (!modelFit) { 178 psTrace ("psphot", 5, "failed to fit plain model for object at %f, %f", source->moments->x, source->moments->y); 179 continue; 180 } 181 pmSourceCacheModel (source, maskVal); 182 psTrace ("psphot", 4, "fit plain model for %f, %f : %s chisq = %f\n", source->moments->x, source->moments->y, pmModelClassGetName (modelFit->type), modelFit->chisq); 183 Nplain ++; 184 if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) { 185 NplainPass ++; 186 } 187 } 188 189 // save each of the model flux images and store the best 190 psArrayAdd (modelFluxes, 4, source->modelFlux); 191 192 // test for fit quality / result 193 psArrayAdd (source->modelFits, 4, modelFit); 194 195 psFree (modelFit); 196 } 197 psFree (iter); 198 199 // evaluate the relative quality of the models, choose one 200 float minChisq = NAN; 201 int minModel = -1; 202 for (int i = 0; i < source->modelFits->n; i++) { 203 pmModel *model = source->modelFits->data[i]; 204 205 if (!(model->flags & PM_MODEL_STATUS_FITTED)) continue; 206 207 if (model->flags & (PM_MODEL_STATUS_BADARGS)) continue; 208 if (model->flags & (PM_MODEL_STATUS_NONCONVERGE)) continue; 209 if (model->flags & (PM_MODEL_STATUS_OFFIMAGE)) continue; 210 211 if ((minModel < 0) || (model->chisq < minChisq)) { 212 minChisq = model->chisq; 213 minModel = i; 214 } 215 } 216 217 if (minModel == -1) { 218 // re-subtract the object, leave local sky 219 psTrace ("psphot", 5, "failed to fit extended source model to object at %f, %f", source->moments->x, source->moments->y); 220 221 // replace original model, subtract it 222 psFree (source->modelFlux); 223 source->modelFlux = modelFluxStart; 224 225 pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); 226 source->mode |= PM_SOURCE_MODE_SUBTRACTED; 227 228 psFree (modelFluxes); 229 230 if (savePics) { 231 psphotSaveImage (NULL, readout->image, "image.xp.fits"); 232 char key[10]; 233 fprintf (stdout, "continue? "); 234 fgets (key, 8, stdin); 235 if (key[0] == 'n') { 236 savePics = false; 237 } 238 } 239 continue; 240 } 241 242 // save the best extended model in modelEXT 243 psFree (source->modelEXT); 244 source->modelEXT = psMemIncrRefCounter (source->modelFits->data[minModel]); 245 246 // save the modelFlux for the best model 247 psFree (source->modelFlux); 248 source->modelFlux = psMemIncrRefCounter (modelFluxes->data[minModel]); 117 // save the modelFlux here in case we need to subtract it (for failure) 118 psImage *modelFluxStart = psMemIncrRefCounter (source->modelFlux); 119 120 if (savePics) { 121 psphotSaveImage (NULL, readout->image, "image.xp.fits"); 122 } 123 124 // array to store the pointers to the model flux images while the models are being fitted 125 psArray *modelFluxes = psArrayAllocEmpty (4); 126 127 // allocate the array to store the model fits 128 if (source->modelFits == NULL) { 129 source->modelFits = psArrayAllocEmpty (4); 130 } 131 132 // loop here over the models chosen for each source (exclude by S/N) 133 psMetadataIterator *iter = psMetadataIteratorAlloc (models, PS_LIST_HEAD, NULL); 134 psMetadataItem *item = NULL; 135 while ((item = psMetadataGetAndIncrement (iter)) != NULL) { 136 137 // XXX this should have been forced above 138 assert (item->type == PS_DATA_METADATA); 139 psMetadata *model = (psMetadata *) item->data.md; 140 141 // check the SNlim and skip model if source is too faint 142 float SNlim = psMetadataLookupF32 (&status, model, "SNLIM_VALUE"); 143 assert (status); 144 145 // limit selection to some SN limit 146 assert (source->peak); // how can a source not have a peak? 147 if (source->peak->SN < SNlim) continue; 148 149 // check on the model type 150 pmModelType modelType = psMetadataLookupS32 (&status, model, "MODEL_TYPE"); 151 assert (status); 152 153 // check on the PSF-Convolution status 154 bool convolved = psMetadataLookupBool (&status, model, "PSF_CONVOLVED_VALUE"); 155 assert (status); 156 157 // XXX psTraceSetLevel ("psLib.math.psMinimizeLMChi2", 6); 158 // XXX psTraceSetLevel ("psphot.psphotModelWithPSF_LMM", 6); 159 160 // fit the model as convolved or not 161 pmModel *modelFit = NULL; 162 if (convolved) { 163 modelFit = psphotPSFConvModel (readout, source, modelType, maskVal, markVal, psfSize); 164 if (!modelFit) { 165 psTrace ("psphot", 5, "failed to fit psf-conv model for object at %f, %f", source->moments->x, source->moments->y); 166 continue; 167 } 168 psTrace ("psphot", 4, "fit psf-conv model for %f, %f : %s chisq = %f\n", source->moments->x, source->moments->y, pmModelClassGetName (modelFit->type), modelFit->chisq); 169 Nconvolve ++; 170 if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) { 171 NconvolvePass ++; 172 } 173 } else { 174 psFree (source->modelFlux); 175 source->modelFlux = NULL; 176 modelFit = psphotFitEXT (readout, source, modelType, maskVal, markVal); 177 if (!modelFit) { 178 psTrace ("psphot", 5, "failed to fit plain model for object at %f, %f", source->moments->x, source->moments->y); 179 continue; 180 } 181 pmSourceCacheModel (source, maskVal); 182 psTrace ("psphot", 4, "fit plain model for %f, %f : %s chisq = %f\n", source->moments->x, source->moments->y, pmModelClassGetName (modelFit->type), modelFit->chisq); 183 Nplain ++; 184 if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) { 185 NplainPass ++; 186 } 187 } 188 189 // save each of the model flux images and store the best 190 psArrayAdd (modelFluxes, 4, source->modelFlux); 191 192 // test for fit quality / result 193 psArrayAdd (source->modelFits, 4, modelFit); 194 195 psFree (modelFit); 196 } 197 psFree (iter); 198 199 // evaluate the relative quality of the models, choose one 200 float minChisq = NAN; 201 int minModel = -1; 202 for (int i = 0; i < source->modelFits->n; i++) { 203 pmModel *model = source->modelFits->data[i]; 204 205 if (!(model->flags & PM_MODEL_STATUS_FITTED)) continue; 206 207 if (model->flags & (PM_MODEL_STATUS_BADARGS)) continue; 208 if (model->flags & (PM_MODEL_STATUS_NONCONVERGE)) continue; 209 if (model->flags & (PM_MODEL_STATUS_OFFIMAGE)) continue; 210 211 if ((minModel < 0) || (model->chisq < minChisq)) { 212 minChisq = model->chisq; 213 minModel = i; 214 } 215 } 216 217 if (minModel == -1) { 218 // re-subtract the object, leave local sky 219 psTrace ("psphot", 5, "failed to fit extended source model to object at %f, %f", source->moments->x, source->moments->y); 220 221 // replace original model, subtract it 222 psFree (source->modelFlux); 223 source->modelFlux = modelFluxStart; 224 225 pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); 226 source->mode |= PM_SOURCE_MODE_SUBTRACTED; 227 228 psFree (modelFluxes); 229 230 if (savePics) { 231 psphotSaveImage (NULL, readout->image, "image.xp.fits"); 232 char key[10]; 233 fprintf (stdout, "continue? "); 234 if (!fgets (key, 8, stdin)) { 235 psWarning("Couldn't read anything."); 236 } else if (key[0] == 'n') { 237 savePics = false; 238 } 239 } 240 continue; 241 } 242 243 // save the best extended model in modelEXT 244 psFree (source->modelEXT); 245 source->modelEXT = psMemIncrRefCounter (source->modelFits->data[minModel]); 246 247 // save the modelFlux for the best model 248 psFree (source->modelFlux); 249 source->modelFlux = psMemIncrRefCounter (modelFluxes->data[minModel]); 249 250 250 251 // subtract the best fit from the object, leave local sky … … 252 253 source->mode |= PM_SOURCE_MODE_SUBTRACTED; 253 254 254 // the initial model flux is no longer needed255 psFree (modelFluxStart);256 psFree (modelFluxes);257 258 psTrace ("psphot", 4, "best ext model for %f, %f : %s chisq = %f\n", source->moments->x, source->moments->y, pmModelClassGetName (source->modelEXT->type), source->modelEXT->chisq);259 psTrace ("psphot", 5, "extended source model for source at %7.1f, %7.1f", source->moments->x, source->moments->y);260 261 if (savePics) {262 psphotSaveImage (NULL, readout->image, "image.xm.fits");263 char key[10];264 fprintf (stdout, "continue? ");265 fgets (key, 8, stdin);266 if (key[0] == 'n') {267 savePics = false;268 }269 }255 // the initial model flux is no longer needed 256 psFree (modelFluxStart); 257 psFree (modelFluxes); 258 259 psTrace ("psphot", 4, "best ext model for %f, %f : %s chisq = %f\n", source->moments->x, source->moments->y, pmModelClassGetName (source->modelEXT->type), source->modelEXT->chisq); 260 psTrace ("psphot", 5, "extended source model for source at %7.1f, %7.1f", source->moments->x, source->moments->y); 261 262 if (savePics) { 263 psphotSaveImage (NULL, readout->image, "image.xm.fits"); 264 char key[10]; 265 fprintf (stdout, "continue? "); 266 fgets (key, 8, stdin); 267 if (key[0] == 'n') { 268 savePics = false; 269 } 270 } 270 271 } 271 272
Note:
See TracChangeset
for help on using the changeset viewer.
