Index: trunk/psphot/src/psphot.h =================================================================== --- trunk/psphot/src/psphot.h (revision 38514) +++ trunk/psphot/src/psphot.h (revision 38515) @@ -513,5 +513,5 @@ bool psphotSourceChildren (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc); bool psphotSourceChildrenReadout (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc, int index); -psArray *psphotSourceChildrenByObject (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objectsSrc, bool sourcesSubtracted); +psArray *psphotSourceChildrenByObject (pmConfig *config, const pmFPAview *view, const char *fileruleOut, const char *fileruleSrc, psArray *objectsSrc, bool sourcesSubtracted); bool psphotSourceParents (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc); Index: trunk/psphot/src/psphotExtendedSourceFits.c =================================================================== --- trunk/psphot/src/psphotExtendedSourceFits.c (revision 38514) +++ trunk/psphot/src/psphotExtendedSourceFits.c (revision 38515) @@ -149,4 +149,8 @@ fitOptions->poissonErrors = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_FITS_POISSON"); if (!status) fitOptions->poissonErrors = true; + + // Save pointer to fitOptions for possible use later reinstantiating the pcm models on the output readout + psMetadataAddPtr (readout->analysis, PS_LIST_TAIL, "PCM_FIT_OPTIONS", PS_DATA_UNKNOWN | PS_META_REPLACE, "pcm fit options", + fitOptions); // maskVal is used to test for rejected pixels, and must include markVal Index: trunk/psphot/src/psphotMergeSources.c =================================================================== --- trunk/psphot/src/psphotMergeSources.c (revision 38514) +++ trunk/psphot/src/psphotMergeSources.c (revision 38515) @@ -853,22 +853,37 @@ // array containing the child sources. XXX currently, this is only used by psphotStackReadout // (sources go on allSources so that psphotChoosePSF can be called repeatedly) -psArray *psphotSourceChildrenByObject (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objectsSrc, bool sourcesSubtracted) { +psArray *psphotSourceChildrenByObject (pmConfig *config, const pmFPAview *view, const char *fileruleOut, const char *fileruleSrc, psArray *objectsSrc, bool sourcesSubtracted) { bool status; - int nImages = psphotFileruleCount(config, filerule); + int nImages = psphotFileruleCount(config, fileruleOut); // generate look-up arrays for detections and readouts psArray *detArrays = psArrayAlloc(nImages); psArray *readouts = psArrayAlloc(nImages); + psArray *fitOptionsArray = psArrayAlloc(nImages); + + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + assert (recipe); + psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); + assert (maskVal); + int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE"); + assert (status); for (int i = 0; i < nImages; i++) { // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest + pmFPAfile *file = pmFPAfileSelectSingle(config->files, fileruleOut, i); // File of interest psAssert (file, "missing file?"); pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); psAssert (readout, "missing readout?"); + + pmFPAfile *fileSrc = pmFPAfileSelectSingle(config->files, fileruleSrc, i); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readoutSrc = pmFPAviewThisReadout(view, fileSrc->fpa); + psAssert (readoutSrc, "missing readout?"); + // create DETECTIONS containers for each image, in case one lacks it @@ -886,8 +901,12 @@ psAssert (detections, "missing detections?"); } + pmSourceFitOptions *fitOptions = psMetadataLookupPtr (&status, readoutSrc->analysis, "PCM_FIT_OPTIONS"); + psAssert (fitOptions, "missing pcm fit options"); + psMetadataAddPtr (readout->analysis, PS_LIST_TAIL, "PCM_FIT_OPTIONS", PS_DATA_UNKNOWN | PS_META_REPLACE, "pcm fit options", fitOptions); // we need to save the new sources on the detection arrays of the appropriate image detArrays->data[i] = psMemIncrRefCounter(detections); readouts->data[i] = psMemIncrRefCounter(readout); + fitOptionsArray->data[i] = psMemIncrRefCounter(fitOptions); } @@ -932,12 +951,22 @@ // does this copy all model data? (NO) sourceOut->modelPSF = pmModelCopy(sourceSrc->modelPSF); - sourceOut->modelEXT = pmModelCopy(sourceSrc->modelEXT); - + + bool foundModelEXT = false; if (sourceSrc->modelFits) { sourceOut->modelFits = psArrayAlloc(sourceSrc->modelFits->n); for (int j = 0; j < sourceSrc->modelFits->n; j++) { - sourceOut->modelFits->data[j] = pmModelCopy(sourceSrc->modelFits->data[j]); - } + pmModel *modelSrc = sourceSrc->modelFits->data[j]; + pmModel *modelOut = sourceOut->modelFits->data[j] = pmModelCopy(modelSrc); + if (modelSrc == sourceSrc->modelEXT) { + foundModelEXT = true; + sourceOut->modelEXT = psMemIncrRefCounter (modelOut); + } + modelOut->isPCM = modelSrc->isPCM; + } } + if (!foundModelEXT && sourceSrc->modelEXT) { + // Will this ever happen? + sourceOut->modelEXT = pmModelCopy(sourceSrc->modelEXT); + } // drop the references to the original image pixels: @@ -949,8 +978,10 @@ pmReadout *readout = readouts->data[index]; + pmSourceFitOptions *fitOptions = fitOptionsArray->data[index]; + // allocate image, weight, mask for the new image for each peak if (sourceOut->modelPSF) { pmSourceRedefinePixels (sourceOut, readout, sourceOut->peak->x, sourceOut->peak->y, - sourceOut->modelPSF->fitRadius); + sourceSrc->windowRadius); } else { // if we have no pixels we can't use it to determine the psf so make sure this bit is off @@ -965,5 +996,23 @@ if (!sourcesSubtracted) { sourceOut->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED; - } + } else { + if (sourceSrc->modelFlux) { + bool isPSF = false; + pmModel *model = pmSourceGetModel (&isPSF, sourceOut); + if (model->isPCM) { + pmPCMdata *pcm = pmPCMinit (sourceOut, fitOptions, model, maskVal, psfSize); + if (pcm) { + // pmPCMMakeModel (sourceOut, model, pcm->nsigma, maskVal, psfSize); + pmPCMCacheModel (sourceOut, maskVal, psfSize, pcm->nsigma); + psFree(pcm); + } else { + // What to do here? + psAssert (pcm, "pmPCMinit failed!"); + } + } else { + pmSourceCacheModel (sourceOut, maskVal); + } + } + } // set the output detections: @@ -981,4 +1030,5 @@ psFree (detArrays); psFree (readouts); + psFree (fitOptionsArray); return objectsOut; Index: trunk/psphot/src/psphotReplaceUnfit.c =================================================================== --- trunk/psphot/src/psphotReplaceUnfit.c (revision 38514) +++ trunk/psphot/src/psphotReplaceUnfit.c (revision 38515) @@ -316,4 +316,6 @@ } + pmSourceFitOptions *fitOptions = psMetadataLookupPtr (&status, readout->analysis, "PCM_FIT_OPTIONS"); + psAssert (fitOptions, "missing pcm fit options"); // XXX the sources have already been copied (merge into here?) @@ -368,8 +370,9 @@ source->modelPSF = modelPSF; source->modelPSF->fitRadius = radius; - model = source->modelPSF; + // model = source->modelPSF; } if (model->isPCM) { +#if 0 psAssert(false, "this section is not complete"); @@ -389,10 +392,19 @@ psImageConvolveFFT (source->modelFlux, rawModelFlux, NULL, 0, psfKernel); - psFree (psfKernel); psFree (rawModelFlux); - } - - pmSourceCacheModel (source, maskVal); // ALLOC x14 (!) +#endif + pmPCMdata *pcm = pmPCMinit (source, fitOptions, model, maskVal, psfSize); + if (pcm) { + // pmPCMMakeModel (source, model, pcm->nsigma, maskVal, psfSize); + pmPCMCacheModel (source, maskVal, psfSize, pcm->nsigma); + psFree(pcm); + } else { + psAssert (pcm, "pmPCMinit failed!"); + } + + } else { + pmSourceCacheModel (source, maskVal); // ALLOC x14 (!) + } } Index: trunk/psphot/src/psphotStackUpdateReadout.c =================================================================== --- trunk/psphot/src/psphotStackUpdateReadout.c (revision 38514) +++ trunk/psphot/src/psphotStackUpdateReadout.c (revision 38515) @@ -116,5 +116,7 @@ // XXX NOTE : if we use the pre-20130914 psphotStackReadout code, we need to use 'false' for the // sourcesSubtracted argument - psArray *objectsOut = psphotSourceChildrenByObject (config, view, STACK_OUT, objects, true); + // XXX: do this here so that the variance is availble to rebuild the models + psphotStackMatchPSFsetup (config, view, STACK_OUT, STACK_RAW); + psArray *objectsOut = psphotSourceChildrenByObject (config, view, STACK_OUT, STACK_RAW, objects, true); if (!objectsOut) { psFree(objects); @@ -135,6 +137,6 @@ // this forces photometry on the undetected sources from other images - // set up the FWHM vector - psphotStackMatchPSFsetup (config, view, STACK_OUT, STACK_RAW); +// Moved this up above +// psphotStackMatchPSFsetup (config, view, STACK_OUT, STACK_RAW); psphotDumpImages (config, view, STACK_RAW, "raw.t0"); psphotDumpImages (config, view, STACK_OUT, "out.t0"); @@ -177,4 +179,9 @@ // and to subtract the sources psphotFitSourcesLinear (config, view, STACK_OUT, false, false); + +#ifdef notyet + psphotRemoveAllSources (config, view, STACK_OUT, false); +#endif + snprintf (line, 256, "%s.%d", "out.t4", entry); psphotDumpImages (config, view, STACK_OUT, line); @@ -221,5 +228,5 @@ COPY_PARAM( "MSKY_NX" ); COPY_PARAM( "MSKY_NY" ); - COPY_PARAM( "CHIP_SEEING" ); + // COPY_PARAM( "CHIP_SEEING" ); } @@ -312,6 +319,10 @@ fitOptions->poissonErrors = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_FITS_POISSON"); if (!status) fitOptions->poissonErrors = true; + int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE"); assert (status); + + // Remember these data for later + psMetadataAddPtr (readout->analysis, PS_LIST_TAIL, "PCM_FIT_OPTIONS", PS_DATA_UNKNOWN | PS_META_REPLACE, "pcm fit options", fitOptions); for (int i = 0; i < sources->n; i++) { @@ -331,11 +342,12 @@ // 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) ; + modelPSF = NULL; } // make sure that the window radius is large enough. Should we do this regardless of whether // or not the source is extended? - if (source->modelEXT) { - float fitRadius = NAN; - float windowRadius = NAN; + float fitRadius = NAN; + float windowRadius = NAN; + if (1 || source->modelEXT) { if (!psphotSetRadiusMoments (&fitRadius, &windowRadius, readout, source, markVal)) { psError (PSPHOT_ERR_UNKNOWN, false, "failed to set radius for ext source"); @@ -387,5 +399,8 @@ float modelFlux = pow(10., -0.4 * model->mag); - // XXX: We are assuming here that we only use PCM models. Should get that information from the recipe. + model->fitRadius = fitRadius; + + // XXX: We are assuming here that we only use PCM models. + // I Should be getting that information from the recipe. pmPCMdata *pcm = pmPCMinit (source, fitOptions, model, maskVal, psfSize); if (!pcm) { @@ -394,5 +409,6 @@ } model->params->data.F32[PM_PAR_I0] = 1.0; - pmPCMMakeModel (source, model, pcm->nsigma, maskVal, psfSize); + // XXX: this duplicates work that pmPCMCacheModel does + // pmPCMMakeModel (source, model, pcm->nsigma, maskVal, psfSize); float normFlux = model->class->modelFlux(model->params); float I0 = modelFlux / normFlux;