Index: trunk/psphot/src/psphotStackUpdateReadout.c =================================================================== --- trunk/psphot/src/psphotStackUpdateReadout.c (revision 38353) +++ trunk/psphot/src/psphotStackUpdateReadout.c (revision 38383) @@ -58,5 +58,7 @@ } - // the intitial run. + // Load the psf from the previous run. + // Note this must be done before we load the sources because it is used. + // I believe only to set the residuals if (!psphotLoadPSF (config, view, STACK_RAW)) { // this only happens if we had a programming error in psphotLoadPSF @@ -78,8 +80,19 @@ } +#ifdef MEASURE_PSF + if (!psphotChoosePSF (config, view, STACK_RAW, true)) { + psLogMsg ("psphot", 3, "failure to construct a psf model"); + return psphotReadoutCleanup (config, view, STACK_RAW); + } + if (!strcasecmp (breakPt, "PSFMODEL")) { + return psphotReadoutCleanup (config, view, STACK_RAW); + } +#endif + if (!psphotDeblendSatstars (config, view, STACK_RAW)) { psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis"); return psphotReadoutCleanup (config, view, STACK_RAW); } + if (!psphotMergeSources (config, view, STACK_RAW)) { psError (PSPHOT_ERR_UNKNOWN, false, "failed to merge sources"); @@ -218,4 +231,5 @@ COPY_PARAM( "MSKY_NX" ); COPY_PARAM( "MSKY_NY" ); + COPY_PARAM( "CHIP_SEEING" ); } @@ -294,16 +308,38 @@ psphotInitRadiusEXT (recipe, readout); + // not really fitting but need to get some things set up to instantiate the pcm models + + pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc(); + fitOptions->mode = PM_SOURCE_FIT_EXT_AND_SKY; + fitOptions->covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix + // XXX: change all of this status checking to asserts. All should be defined with the current recipe. + float fitNsigmaConv = psMetadataLookupF32 (&status, recipe, "EXT_FIT_NSIGMA_CONV"); // number of sigma for the convolution + if (!status || !isfinite(fitNsigmaConv) || fitNsigmaConv <= 0) { + fitNsigmaConv = 5.0; + } + fitOptions->nsigma = fitNsigmaConv; + // Poisson or Constant weight for chisq tests? + fitOptions->poissonErrors = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_FITS_POISSON"); + if (!status) fitOptions->poissonErrors = true; + int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE"); + assert (status); + for (int i = 0; i < sources->n; i++) { pmSource *source = sources->data[i]; + source->imageID = index; // XXX is this necessary? pmModel *modelPSF = source->modelPSF; - if (!modelPSF) { - continue; - } // free any previous radial aperture measurements psFree(source->radialAper); // Guess Models ususally does this - psphotCheckRadiusPSF (readout, source, modelPSF, markVal); - source->modelPSF->residuals = psf->residuals; + if (modelPSF) { + psphotCheckRadiusPSF (readout, source, modelPSF, markVal); + source->modelPSF->residuals = psf->residuals; + } + if (!(source->mode & PM_SOURCE_MODE_PSFMODEL)) { + // The cmf loaders unconditionallly create psf models even if the source did not have one. + fprintf(stderr, "source %d %5d should not have a psf model 0x%08X%08X\n", index, source->seq, source->mode2, source->mode); + psFree(source->modelPSF) ; + } // make sure that the window radius is large enough. Should we do this regardless of whether @@ -318,39 +354,85 @@ } + // If we find a good model we will cache it below. bool goodModel = false; - // the cmf readers do not set PM_PAR_I0 for the models. - // Calculate it here. + + // the cmf readers do not set PM_PAR_I0 for the models so we calculate it here. // We may want to put this code in the cmf readers but for now I don't want to change psphot operation - // for anything except psphotStack -updatemode - if (isfinite(source->psfFlux)) { + // for anything except psphotStack -updatemode. + if (modelPSF && isfinite(source->psfFlux)) { + // The cmf reader is incorrectly setting the sky parameter for psf models + modelPSF->params->data.F32[PM_PAR_SKY] = 0.0; + modelPSF->dparams->data.F32[PM_PAR_SKY] = 0.0; + // Calculate flux for normalized model + float saveI0 = modelPSF->params->data.F32[PM_PAR_I0]; modelPSF->params->data.F32[PM_PAR_I0] = 1.0; float normFlux = modelPSF->class->modelFlux(modelPSF->params); - modelPSF->params->data.F32[PM_PAR_I0] = source->psfFlux / normFlux; - goodModel = true; - } - if (source->modelFits && source->modelFits->n) { - for (int m = 0; m < source->modelFits->n; m++) { - pmModel *model = source->modelFits->data[m]; - if (isfinite(model->mag)) { - // calculate flux from model magnitude - float modelFlux = pow(10., -0.4 * model->mag); - model->params->data.F32[PM_PAR_I0] = 1.0; - float normFlux = modelPSF->class->modelFlux(modelPSF->params); - model->params->data.F32[PM_PAR_I0] = modelFlux / normFlux; - if (model == source->modelEXT) { - goodModel = true; - } - } else { - // nan model magnitude no way to calculate the flux. - // The flux was probably negative. - // If this is the extended model for this source we can't use the model. - // Set source type to star to prevent source subtraction from crashing. - // This is somewhat rare .2% of extened sources. Is there another way that we can estimate - // the flux, or more to the point the value for PM_PAR_I0? - if (model == source->modelEXT && source->type == PM_SOURCE_TYPE_EXTENDED) { - source->type = PM_SOURCE_TYPE_STAR; - goodModel = false; - } + float psfFlux = source->psfFlux; +#define APPLY_APCORR +#ifdef APPLY_APCORR + double apTrend = pmTrend2DEval (psf->ApTrend, (float)source->peak->x, (float)source->peak->y); + double adjustment = pow(10.0, -0.4*apTrend); + if (isfinite(adjustment) && adjustment != 0.0) { + psfFlux = psfFlux / adjustment; + } else { + fprintf(stderr, "apTrend adjustment for source %5d invalid: %f\n", source->seq, adjustment); + } +#endif + float newI0 = psfFlux / normFlux; + if (isfinite (newI0) ) { + // psf model looks to be good. + modelPSF->params->data.F32[PM_PAR_I0] = newI0; +#ifdef USE_SAVED_I0 + // I am testing with a modified cmf that saves PM_PAR_I0 + modelPSF->params->data.F32[PM_PAR_I0] = saveI0; +#endif + // fprintf (stderr, "%5d %12f %12f %12f\n", source->seq, saveI0, newI0, saveI0 - newI0); + goodModel = true; + } else { + psAssert(isfinite (modelPSF->params->data.F32[PM_PAR_I0]), "got infinite I0 for psf model"); + } + } + + // It is rather expensive initializing all of the extended models. Lets only initialize + // the extended model that will be used if any. + + bool isPSF = false; + pmModel *model = pmSourceGetModel (&isPSF, source); + if (model && !isPSF) { + // selected model is extended + goodModel = false; + // calculate flux from the previously measured model magnitude. It will be nan if the + // flux was negative unfortunately. + if (isfinite(model->mag)) { + float modelFlux = pow(10., -0.4 * model->mag); + + // XXX: We are assuming here that we only use PCM models. Get from recipe. + pmPCMdata *pcm = pmPCMinit (source, fitOptions, model, maskVal, psfSize); + if (!pcm) { + psError (PSPHOT_ERR_UNKNOWN, false, "failed to to initialize PCM for model"); + return false; } + model->params->data.F32[PM_PAR_I0] = 1.0; + pmPCMMakeModel (source, model, pcm->nsigma, maskVal, psfSize); + float normFlux = model->class->modelFlux(model->params); + float I0 = modelFlux / normFlux; + if (isfinite(I0)) { + model->params->data.F32[PM_PAR_I0] = I0; + goodModel = true; + // Usually this it is set when doing the fit. Since we are instantiating the model we + // set it explicitly. + model->isPCM = true; + } + psFree(pcm); + } + if (!goodModel) { + // We do not have a usable extended model for this source. + // The original flux was probably negative. + // This is somewhat rare. I saw it in about .2% of extended sources in my first test. + // Set source type to star to prevent source subtraction from crashing. + // this may cause model psf to be used in places. + // XXX However we aren't going to cache it below. + source->type = PM_SOURCE_TYPE_STAR; + model = NULL; } } @@ -373,13 +455,20 @@ } - // Should we use pmPCMCacheModel for extended sources? - // XXX: are the extended models ready? - if (goodModel) { - pmSourceCacheModel (source, maskVal); // ALLOC x14 (!) - } - } - - // as a test drop psf so we can try and remeasure it -// psMetadataRemoveKey(readout->analysis, "PSPHOT.PSF"); + if (goodModel && model) { + if (model->isPCM) { + pmPCMCacheModel (source, maskVal, psfSize, fitNsigmaConv); + } else { + pmSourceCacheModel (source, maskVal); + } + } + } + + psFree(fitOptions); + +#ifdef MEASURE_PSF + // need to drop the psf model in order for it to be re-measured + // Hmm, will this invalidate the results from above where we used the loaded one? + psMetadataRemoveKey(readout->analysis, "PSPHOT.PSF"); +#endif return true;