Index: trunk/psphot/src/psphotAddNoise.c =================================================================== --- trunk/psphot/src/psphotAddNoise.c (revision 21458) +++ trunk/psphot/src/psphotAddNoise.c (revision 21519) @@ -50,5 +50,5 @@ // skip sources which were not subtracted - if (!(source->mode & PM_SOURCE_MODE_SUBTRACTED)) continue; + if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) continue; // select appropriate model Index: trunk/psphot/src/psphotApResid.c =================================================================== --- trunk/psphot/src/psphotApResid.c (revision 21458) +++ trunk/psphot/src/psphotApResid.c (revision 21519) @@ -1,5 +1,3 @@ # include "psphotInternal.h" - -bool psphotApResidMags_Unthreaded (int *nskip, int *nfail, psArray *sources, pmPSF *psf, pmSourcePhotometryMode photMode, psImageMaskType maskVal); # define SKIPSTAR(MSG) { psTrace ("psphot", 3, "invalid : %s", MSG); continue; } @@ -95,23 +93,23 @@ for (int j = 0; j < cells->n; j++) { - if (nThreads) { - // allocate a job -- if threads are not defined, this just runs the job - psThreadJob *job = psThreadJobAlloc ("PSPHOT_APRESID_MAGS"); - - psArrayAdd(job->args, 1, cells->data[j]); // sources - psArrayAdd(job->args, 1, psf); - PS_ARRAY_ADD_SCALAR(job->args, photMode, PS_TYPE_S32); - PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); - - PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nskip - PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail - - if (!psThreadJobAddPending(job)) { - psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - psFree (job); - return false; - } - psFree(job); - } else { + // allocate a job -- if threads are not defined, this just runs the job + psThreadJob *job = psThreadJobAlloc ("PSPHOT_APRESID_MAGS"); + + psArrayAdd(job->args, 1, cells->data[j]); // sources + psArrayAdd(job->args, 1, psf); + PS_ARRAY_ADD_SCALAR(job->args, photMode, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); + + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nskip + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail + + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + psFree (job); + return false; + } + psFree(job); + +# if (0) int nskip = 0; int nfail = 0; @@ -123,29 +121,27 @@ Nskip += nskip; Nfail += nfail; - } +# endif } - if (nThreads) { - // wait for the threads to finish and manage results - if (!psThreadPoolWait (false)) { - psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - return false; + // wait for the threads to finish and manage results + if (!psThreadPoolWait (false)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + return false; + } + + // we have only supplied one type of job, so we can assume the types here + psThreadJob *job = NULL; + while ((job = psThreadJobGetDone()) != NULL) { + if (job->args->n < 1) { + fprintf (stderr, "error with job\n"); + } else { + psScalar *scalar = NULL; + scalar = job->args->data[4]; + Nskip += scalar->data.S32; + scalar = job->args->data[5]; + Nfail += scalar->data.S32; } - - // we have only supplied one type of job, so we can assume the types here - psThreadJob *job = NULL; - while ((job = psThreadJobGetDone()) != NULL) { - if (job->args->n < 1) { - fprintf (stderr, "error with job\n"); - } else { - psScalar *scalar = NULL; - scalar = job->args->data[4]; - Nskip += scalar->data.S32; - scalar = job->args->data[5]; - Nfail += scalar->data.S32; - } - psFree(job); - } + psFree(job); } } @@ -497,4 +493,5 @@ continue; } + source->mode |= PM_SOURCE_MODE_AP_MAGS; } @@ -509,4 +506,5 @@ } +# if (0) bool psphotApResidMags_Unthreaded (int *nskip, int *nfail, psArray *sources, pmPSF *psf, pmSourcePhotometryMode photMode, psImageMaskType maskVal) { @@ -542,2 +540,3 @@ return true; } +# endif Index: trunk/psphot/src/psphotBlendFit.c =================================================================== --- trunk/psphot/src/psphotBlendFit.c (revision 21458) +++ trunk/psphot/src/psphotBlendFit.c (revision 21519) @@ -1,5 +1,3 @@ # include "psphotInternal.h" - -bool psphotBlendFit_Unthreaded (int *nfit, int *npsf, int *next, int *nfail, pmReadout *readout, psMetadata *recipe, psArray *sources, pmPSF *psf, psArray *newSources); // XXX I don't like this name @@ -62,28 +60,29 @@ for (int j = 0; j < cells->n; j++) { - if (nThreads) { - // allocate a job -- if threads are not defined, this just runs the job - psThreadJob *job = psThreadJobAlloc ("PSPHOT_BLEND_FIT"); - psArray *newSources = psArrayAllocEmpty(16); - - psArrayAdd(job->args, 1, readout); - psArrayAdd(job->args, 1, recipe); - psArrayAdd(job->args, 1, cells->data[j]); // sources - psArrayAdd(job->args, 1, psf); - psArrayAdd(job->args, 1, newSources); // return for new sources - psFree (newSources); - - PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfit - PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Npsf - PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Next - PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail - - if (!psThreadJobAddPending(job)) { - psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - psFree (job); - return NULL; - } - psFree(job); - } else { + // allocate a job -- if threads are not defined, this just runs the job + psThreadJob *job = psThreadJobAlloc ("PSPHOT_BLEND_FIT"); + psArray *newSources = psArrayAllocEmpty(16); + + psArrayAdd(job->args, 1, readout); + psArrayAdd(job->args, 1, recipe); + psArrayAdd(job->args, 1, cells->data[j]); // sources + psArrayAdd(job->args, 1, psf); + psArrayAdd(job->args, 1, newSources); // return for new sources + psFree (newSources); + + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfit + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Npsf + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Next + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail + + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + psFree (job); + return NULL; + } + psFree(job); + +# if (0) + { int nfit = 0; int npsf = 0; @@ -107,38 +106,37 @@ psFree (newSources); } - } - - if (nThreads) { - // wait for the threads to finish and manage results - if (!psThreadPoolWait (false)) { - psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - return NULL; - } - - // we have only supplied one type of job, so we can assume the types here - psThreadJob *job = NULL; - while ((job = psThreadJobGetDone()) != NULL) { - if (job->args->n < 1) { - fprintf (stderr, "error with job\n"); - } else { - psScalar *scalar = NULL; - scalar = job->args->data[5]; - Nfit += scalar->data.S32; - scalar = job->args->data[6]; - Npsf += scalar->data.S32; - scalar = job->args->data[7]; - Next += scalar->data.S32; - scalar = job->args->data[8]; - Nfail += scalar->data.S32; - - // add these back onto sources - psArray *newSources = job->args->data[4]; - for (int j = 0; j < newSources->n; j++) { - psArrayAdd (sources, 16, newSources->data[j]); - } - } - psFree(job); - } - } +# endif + } + + // wait for the threads to finish and manage results + if (!psThreadPoolWait (false)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + return NULL; + } + + // we have only supplied one type of job, so we can assume the types here + psThreadJob *job = NULL; + while ((job = psThreadJobGetDone()) != NULL) { + if (job->args->n < 1) { + fprintf (stderr, "error with job\n"); + } else { + psScalar *scalar = NULL; + scalar = job->args->data[5]; + Nfit += scalar->data.S32; + scalar = job->args->data[6]; + Npsf += scalar->data.S32; + scalar = job->args->data[7]; + Next += scalar->data.S32; + scalar = job->args->data[8]; + Nfail += scalar->data.S32; + + // add these back onto sources + psArray *newSources = job->args->data[4]; + for (int j = 0; j < newSources->n; j++) { + psArrayAdd (sources, 16, newSources->data[j]); + } + } + psFree(job); + } } psFree (cellGroups); @@ -225,5 +223,118 @@ // replace object in image - if (source->mode & PM_SOURCE_MODE_SUBTRACTED) { + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + } + Nfit ++; + + // try fitting PSFs or extended sources depending on source->mode + // these functions subtract the resulting fitted source + if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { + if (psphotFitBlob (readout, source, newSources, psf, maskVal, markVal)) { + source->type = PM_SOURCE_TYPE_EXTENDED; + psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf); + Next ++; + source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT; + continue; + } + } else { + if (psphotFitBlend (readout, source, psf, maskVal, markVal)) { + source->type = PM_SOURCE_TYPE_STAR; + psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf); + Npsf ++; + source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT; + continue; + } + } + + psTrace ("psphot", 5, "source at %7.1f, %7.1f failed", source->peak->xf, source->peak->yf); + Nfail ++; + + // re-subtract the object, leave local sky + pmSourceCacheModel (source, maskVal); + pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED; + } + + // change the value of a scalar on the array (wrap this and put it in psArray.h) + scalar = job->args->data[5]; + scalar->data.S32 = Nfit; + + scalar = job->args->data[6]; + scalar->data.S32 = Npsf; + + scalar = job->args->data[7]; + scalar->data.S32 = Next; + + scalar = job->args->data[8]; + scalar->data.S32 = Nfail; + + return true; +} + +# if (0) +bool psphotBlendFit_Unthreaded (int *nfit, int *npsf, int *next, int *nfail, pmReadout *readout, psMetadata *recipe, psArray *sources, pmPSF *psf, psArray *newSources) { + + bool status = false; + int Nfit = 0; + int Npsf = 0; + int Next = 0; + int Nfail = 0; + + // bit-masks to test for good/bad pixels + psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); + assert (maskVal); + + // bit-mask to mark pixels not used in analysis + psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); + assert (markVal); + + // maskVal is used to test for rejected pixels, and must include markVal + maskVal |= markVal; + + // S/N limit to perform full non-linear fits + float FIT_SN_LIM = psMetadataLookupF32 (&status, recipe, "FULL_FIT_SN_LIM"); + + // option to limit analysis to a specific region + char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION"); + psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region)); + if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined"); + + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + + // skip non-astronomical objects (very likely defects) + if (source->mode & PM_SOURCE_MODE_BLEND) continue; + if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue; + if (source->type == PM_SOURCE_TYPE_DEFECT) continue; + if (source->type == PM_SOURCE_TYPE_SATURATED) continue; + + // skip DBL second sources (ie, added by psphotFitBlob) + if (source->mode & PM_SOURCE_MODE_PAIR) continue; + + // limit selection to some SN limit + if (source->peak->SN < FIT_SN_LIM) continue; + + // exclude sources outside optional analysis region + if (source->peak->xf < AnalysisRegion.x0) continue; + if (source->peak->yf < AnalysisRegion.y0) continue; + if (source->peak->xf > AnalysisRegion.x1) continue; + if (source->peak->yf > AnalysisRegion.y1) continue; + + // if model is NULL, we don't have a starting guess + if (source->modelPSF == NULL) continue; + + // skip sources which are insignificant flux? + // XXX this is somewhat ad-hoc + if (source->modelPSF->params->data.F32[1] < 0.1) { + psTrace ("psphot", 5, "skipping near-zero source: %f, %f : %f\n", + source->modelPSF->params->data.F32[1], + source->modelPSF->params->data.F32[2], + source->modelPSF->params->data.F32[3]); + continue; + } + + // replace object in image + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); } @@ -254,115 +365,5 @@ pmSourceCacheModel (source, maskVal); pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - source->mode |= PM_SOURCE_MODE_SUBTRACTED; - } - - // change the value of a scalar on the array (wrap this and put it in psArray.h) - scalar = job->args->data[5]; - scalar->data.S32 = Nfit; - - scalar = job->args->data[6]; - scalar->data.S32 = Npsf; - - scalar = job->args->data[7]; - scalar->data.S32 = Next; - - scalar = job->args->data[8]; - scalar->data.S32 = Nfail; - - return true; -} - -bool psphotBlendFit_Unthreaded (int *nfit, int *npsf, int *next, int *nfail, pmReadout *readout, psMetadata *recipe, psArray *sources, pmPSF *psf, psArray *newSources) { - - bool status = false; - int Nfit = 0; - int Npsf = 0; - int Next = 0; - int Nfail = 0; - - // bit-masks to test for good/bad pixels - psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); - assert (maskVal); - - // bit-mask to mark pixels not used in analysis - psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); - assert (markVal); - - // maskVal is used to test for rejected pixels, and must include markVal - maskVal |= markVal; - - // S/N limit to perform full non-linear fits - float FIT_SN_LIM = psMetadataLookupF32 (&status, recipe, "FULL_FIT_SN_LIM"); - - // option to limit analysis to a specific region - char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION"); - psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region)); - if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined"); - - for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - - // skip non-astronomical objects (very likely defects) - if (source->mode & PM_SOURCE_MODE_BLEND) continue; - if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue; - if (source->type == PM_SOURCE_TYPE_DEFECT) continue; - if (source->type == PM_SOURCE_TYPE_SATURATED) continue; - - // skip DBL second sources (ie, added by psphotFitBlob) - if (source->mode & PM_SOURCE_MODE_PAIR) continue; - - // limit selection to some SN limit - if (source->peak->SN < FIT_SN_LIM) continue; - - // exclude sources outside optional analysis region - if (source->peak->xf < AnalysisRegion.x0) continue; - if (source->peak->yf < AnalysisRegion.y0) continue; - if (source->peak->xf > AnalysisRegion.x1) continue; - if (source->peak->yf > AnalysisRegion.y1) continue; - - // if model is NULL, we don't have a starting guess - if (source->modelPSF == NULL) continue; - - // skip sources which are insignificant flux? - // XXX this is somewhat ad-hoc - if (source->modelPSF->params->data.F32[1] < 0.1) { - psTrace ("psphot", 5, "skipping near-zero source: %f, %f : %f\n", - source->modelPSF->params->data.F32[1], - source->modelPSF->params->data.F32[2], - source->modelPSF->params->data.F32[3]); - continue; - } - - // replace object in image - if (source->mode & PM_SOURCE_MODE_SUBTRACTED) { - pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); - } - Nfit ++; - - // try fitting PSFs or extended sources depending on source->mode - // these functions subtract the resulting fitted source - if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { - if (psphotFitBlob (readout, source, newSources, psf, maskVal, markVal)) { - source->type = PM_SOURCE_TYPE_EXTENDED; - psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf); - Next ++; - continue; - } - } else { - if (psphotFitBlend (readout, source, psf, maskVal, markVal)) { - source->type = PM_SOURCE_TYPE_STAR; - psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf); - Npsf ++; - continue; - } - } - - psTrace ("psphot", 5, "source at %7.1f, %7.1f failed", source->peak->xf, source->peak->yf); - Nfail ++; - - // re-subtract the object, leave local sky - pmSourceCacheModel (source, maskVal); - pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - source->mode |= PM_SOURCE_MODE_SUBTRACTED; + source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED; } @@ -375,2 +376,3 @@ return true; } +# endif Index: trunk/psphot/src/psphotDeblendSatstars.c =================================================================== --- trunk/psphot/src/psphotDeblendSatstars.c (revision 21458) +++ trunk/psphot/src/psphotDeblendSatstars.c (revision 21519) @@ -111,5 +111,5 @@ // any peaks within this contour should be dropped (mark as blend, drop after loop is done) - // also drop any peaks which are too close to this peal + // also drop any peaks which are too close to this peak psVector *xv = contour->data[0]; psVector *yv = contour->data[1]; Index: trunk/psphot/src/psphotDetectReadout.c =================================================================== --- trunk/psphot/src/psphotDetectReadout.c (revision 21458) +++ trunk/psphot/src/psphotDetectReadout.c (revision 21519) @@ -22,7 +22,4 @@ // Generate the mask and weight images, including the user-defined analysis region of interest psphotSetMaskAndWeight (config, readout, recipe); - if (!strcasecmp (breakPt, "NOTHING")) { - return psphotReadoutCleanup(config, readout, recipe, NULL, NULL, NULL); - } // load the psf model, if suppled. FWHM_X,FWHM_Y,etc are saved in the recipe @@ -38,22 +35,5 @@ // include externally-supplied sources (supplied as PSPHOT.INPUT.CMF) - pmDetections *detections = psphotDetectionsFromSources (config, inSources); - if (!detections || !detections->peaks) { - psError(PSPHOT_ERR_ARGUMENTS, true, "Can't find PSF stars"); - return psphotReadoutCleanup(config, readout, recipe, detections, NULL, NULL); - } - - // construct sources and measure basic stats - psArray *sources = psphotSourceStats (config, readout, detections); - if (!sources) return false; - if (!strcasecmp (breakPt, "PEAKS")) { - return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources); - } - - // classify sources based on moments, brightness - if (!psphotRoughClass (readout, sources, recipe, havePSF)) { - psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); - } + psphotSetSourceParams (config, sources, psf); // calculate source magnitudes @@ -63,2 +43,12 @@ return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources); } + +// deep in pmSourcePixelWeight, we need the following items for each location of interest: +// pmModel *model (should be a realized version of the PSF, can have unity normalization) +// psImage *mask (local subimage of mask for source) +// psImageMaskType maskVal (from recipe & mask header) + +{ + pmModel *model = pmModelFromPSFforXY (psf, x, y, 1.0); +} + Index: trunk/psphot/src/psphotExtendedSourceAnalysis.c =================================================================== --- trunk/psphot/src/psphotExtendedSourceAnalysis.c (revision 21458) +++ trunk/psphot/src/psphotExtendedSourceAnalysis.c (revision 21519) @@ -61,5 +61,5 @@ // replace object in image - if (source->mode & PM_SOURCE_MODE_SUBTRACTED) { + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); } @@ -73,7 +73,8 @@ psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - source->mode |= PM_SOURCE_MODE_SUBTRACTED; + source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED; continue; } + source->mode |= PM_SOURCE_MODE_RADIAL_FLUX; } @@ -85,4 +86,5 @@ psTrace ("psphot", 5, "measured isophotal mags for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); Nisophot ++; + source->mode |= PM_SOURCE_MODE_EXTENDED_STATS; } } @@ -95,4 +97,5 @@ psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); Npetro ++; + source->mode |= PM_SOURCE_MODE_EXTENDED_STATS; } } @@ -105,4 +108,5 @@ psTrace ("psphot", 5, "measure kron mags for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); Nkron ++; + source->mode |= PM_SOURCE_MODE_EXTENDED_STATS; } } @@ -116,9 +120,10 @@ psTrace ("psphot", 5, "measured annuli for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); Nannuli ++; + source->mode |= PM_SOURCE_MODE_EXTENDED_STATS; } // re-subtract the object, leave local sky pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - source->mode |= PM_SOURCE_MODE_SUBTRACTED; + source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED; } Index: trunk/psphot/src/psphotExtendedSourceFits.c =================================================================== --- trunk/psphot/src/psphotExtendedSourceFits.c (revision 21458) +++ trunk/psphot/src/psphotExtendedSourceFits.c (revision 21519) @@ -110,5 +110,5 @@ // replace object in image - if (source->mode & PM_SOURCE_MODE_SUBTRACTED) { + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); } @@ -170,4 +170,5 @@ if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) { NconvolvePass ++; + source->mode |= PM_SOURCE_MODE_EXTENDED_FIT; } } else { @@ -184,4 +185,5 @@ if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) { NplainPass ++; + source->mode |= PM_SOURCE_MODE_EXTENDED_FIT; } } @@ -224,5 +226,5 @@ pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - source->mode |= PM_SOURCE_MODE_SUBTRACTED; + source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED; psFree (modelFluxes); @@ -251,5 +253,5 @@ // subtract the best fit from the object, leave local sky pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - source->mode |= PM_SOURCE_MODE_SUBTRACTED; + source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED; // the initial model flux is no longer needed Index: trunk/psphot/src/psphotExtendedSources.c =================================================================== --- trunk/psphot/src/psphotExtendedSources.c (revision 21458) +++ (revision ) @@ -1,137 +1,0 @@ -# include "psphot.h" - -bool psphotExtendedSources (pmReadout *readout, psArray *sources, psMetadata *recipe) { - - bool status; - int Next = 0; - int Nconvolve = 0; - int Npetro = 0; - int Nisophot = 0; - int Nannuli = 0; - int Nkron = 0; - - // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) - psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels - assert (maskVal); - - // perform full non-linear fits / extended source analysis? - if (!psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_FITS")) { - psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source measurements\n"); - return true; - } - - // option to limit analysis to a specific region - char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION"); - psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region)); - if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined"); - - // S/N limit to perform full non-linear fits - float SN_LIM = psMetadataLookupF32 (&status, recipe, "EXTENDED_SOURCE_SN_LIM"); - - // which extended source analyses should we perform? - bool doPSFConvModel = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PSF_CONVOLVED_MODEL"); - bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN"); - bool doIsophotal = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL"); - bool doAnnuli = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI"); - bool doKron = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON"); - - // source analysis is done in S/N order (brightest first) - sources = psArraySort (sources, pmSourceSortBySN); - - // XXX some init functions for the extended source recipe options? - - // choose the sources of interest - for (int i = 0; i < sources->n; i++) { - - pmSource *source = sources->data[i]; - - // skip PSF-like and non-astronomical objects - if (source->type == PM_SOURCE_TYPE_STAR) continue; - if (source->type == PM_SOURCE_TYPE_DEFECT) continue; - if (source->type == PM_SOURCE_TYPE_SATURATED) continue; - - // limit selection to some SN limit - assert (source->peak); // how can a source not have a peak? - if (source->peak->SN < SN_LIM) continue; - - // XXX this should use peak? - if (source->peak->x < AnalysisRegion.x0) continue; - if (source->peak->y < AnalysisRegion.y0) continue; - if (source->peak->x > AnalysisRegion.x1) continue; - if (source->peak->y > AnalysisRegion.y1) continue; - - // if model is NULL, we don't have a starting guess - // XXX this is probably not needed for objects where we only measure aperture-like terms - if (source->modelEXT == NULL) continue; - - // replace object in image - if (source->mode & PM_SOURCE_MODE_SUBTRACTED) { - pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); - } - Next ++; - - if (doPSFConvModel && !psphotPSFConvModel (source, recipe, maskVal)) { - psError(PSPHOT_ERR_UNKNOWN, false, "failure in PSF Convolved Model fit"); - return false; - } else { - // XXX why am I caching the model? - pmSourceCacheModel (source, maskVal); // XXX put this in the source model function? - psTrace ("psphot", 5, "psf-convolved model for source at %7.1f, %7.1f", source->moments->x, source->moments->y); - Nconvolve ++; - } - - // all of these below require the radial profile - // XXX push this as a test and call in each of the functions below? - // XXX this constructs a pmSourceExtendedParameters element - if (doPetrosian || doIsophotal || doAnnuli || doKron) { - if (!psphotRadialProfile (source, recipe, maskVal)) { - psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate radial profile"); - return false; - } - } - - if (doPetrosian && !psphotPetrosian (source, recipe, maskVal)) { - psError(PSPHOT_ERR_UNKNOWN, false, "failure in Petrosian analysis"); - return false; - } else { - psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->x, source->moments->y); - Npetro ++; - } - - if (doIsophotal && !psphotIsophotal (source, recipe, maskVal)) { - psError(PSPHOT_ERR_UNKNOWN, false, "failure in Isophotal analysis"); - return false; - } else { - psTrace ("psphot", 5, "psf-convolved model for source at %7.1f, %7.1f", source->moments->x, source->moments->y); - Nisophot ++; - } - - if (doAnnuli && !psphotAnnuli (source, recipe, maskVal)) { - psError(PSPHOT_ERR_UNKNOWN, false, "failure in Annuli analysis"); - return false; - } else { - psTrace ("psphot", 5, "psf-convolved model for source at %7.1f, %7.1f", source->moments->x, source->moments->y); - Nannuli ++; - } - - if (doKron && !psphotKron (source, recipe, maskVal)) { - psError(PSPHOT_ERR_UNKNOWN, false, "failure in Kron analysis"); - return false; - } else { - psTrace ("psphot", 5, "psf-convolved model for source at %7.1f, %7.1f", source->moments->x, source->moments->y); - Nkron ++; - } - - // re-subtract the object, leave local sky - pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - source->mode |= PM_SOURCE_MODE_SUBTRACTED; - } - - psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot"), Next); - psLogMsg ("psphot", PS_LOG_INFO, " %d convolved models\n", Nconvolve); - psLogMsg ("psphot", PS_LOG_INFO, " %d petrosian\n", Npetro); - psLogMsg ("psphot", PS_LOG_INFO, " %d isophotal\n", Nisophot); - psLogMsg ("psphot", PS_LOG_INFO, " %d annuli\n", Nannuli); - psLogMsg ("psphot", PS_LOG_INFO, " %d kron\n", Nkron); - return true; -} Index: trunk/psphot/src/psphotFindPeaks.c =================================================================== --- trunk/psphot/src/psphotFindPeaks.c (revision 21458) +++ trunk/psphot/src/psphotFindPeaks.c (revision 21519) @@ -21,7 +21,8 @@ } - // correct the peak values to S/N = sqrt(significance) - // get the peak flux from the unsmoothed image - // the peak pixel coords are guaranteed to be on the image + // Convert the peak values to S/N = sqrt(significance). + // Get the peak flux from the unsmoothed image. + // Rescale the peak position errors using the peak variance + // The peak pixel coords are guaranteed to be on the image int row0 = readout->image->row0; int col0 = readout->image->col0; @@ -30,4 +31,10 @@ peak->SN = sqrt(peak->value); peak->flux = readout->image->data.F32[peak->y-row0][peak->x-col0]; + if (readout->variance && isfinite (peak->dx)) { + peak->dx *= sqrt(readout->variance->data.F32[peak->y-row0][peak->x-col0]); + } + if (readout->variance && isfinite (peak->dy)) { + peak->dy *= sqrt(readout->variance->data.F32[peak->y-row0][peak->x-col0]); + } } Index: trunk/psphot/src/psphotFitSourcesLinear.c =================================================================== --- trunk/psphot/src/psphotFitSourcesLinear.c (revision 21458) +++ trunk/psphot/src/psphotFitSourcesLinear.c (revision 21519) @@ -63,15 +63,18 @@ pmSource *source = sources->data[i]; + // turn this bit off and turn it on again if we pass this test + source->mode &= ~PM_SOURCE_MODE_LINEAR_FIT; + // skip non-astronomical objects (very likely defects) if (source->type == PM_SOURCE_TYPE_DEFECT) continue; if (source->type == PM_SOURCE_TYPE_SATURATED) continue; - // if (source->type == PM_SOURCE_TYPE_STAR && + // do not include CRs in the full ensemble fit if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue; if (final) { - if (source->mode & PM_SOURCE_MODE_SUBTRACTED) continue; + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) continue; } else { - if (source->mode & PM_SOURCE_MODE_BLEND) continue; + if (source->mode & PM_SOURCE_MODE_BLEND) continue; } @@ -91,4 +94,5 @@ if (y > AnalysisRegion.y1) continue; + source->mode |= PM_SOURCE_MODE_LINEAR_FIT; psArrayAdd (fitSources, 100, source); } @@ -185,4 +189,10 @@ psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "solve matrix: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem); + // XXXX **** philosophical question: + // we measure bright objects in three passes: 1) linear fit; 2) non-linear fit; 3) linear fit: + // should retain the chisq and errors from the intermediate non-linear fit? + // the non-linear fit provides better values for the position errors, and for + // extended sources, the shape errors + // adjust I0 for fitSources and subtract for (int i = 0; i < fitSources->n; i++) { @@ -194,4 +204,5 @@ psAbort("linear fitted source is nan"); } + model->params->data.F32[PM_PAR_I0] = norm->data.F32[i]; model->dparams->data.F32[PM_PAR_I0] = errors->data.F32[i]; @@ -200,5 +211,5 @@ // subtract object pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - source->mode |= PM_SOURCE_MODE_SUBTRACTED; + source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED; } psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "sub models: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem); @@ -207,4 +218,5 @@ for (int i = 0; final && (i < fitSources->n); i++) { pmSource *source = fitSources->data[i]; + if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) continue; pmModel *model = pmSourceGetModel (NULL, source); pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal); Index: trunk/psphot/src/psphotGuessModels.c =================================================================== --- trunk/psphot/src/psphotGuessModels.c (revision 21458) +++ trunk/psphot/src/psphotGuessModels.c (revision 21519) @@ -22,6 +22,4 @@ } -bool psphotGuessModel_Unthreaded (pmReadout *readout, psArray *sources, pmPSF *psf, psImageMaskType maskVal, psImageMaskType markVal); - // construct an initial PSF model for each object bool psphotGuessModels (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf) { @@ -67,47 +65,45 @@ for (int j = 0; j < cells->n; j++) { - if (nThreads) { - // allocate a job -- if threads are not defined, this just runs the job - psThreadJob *job = psThreadJobAlloc ("PSPHOT_GUESS_MODEL"); - psArrayAdd(job->args, 1, readout); - psArrayAdd(job->args, 1, cells->data[j]); // sources - psArrayAdd(job->args, 1, psf); - - // XXX change these to use abstract mask type info - PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); - PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK); - - if (!psThreadJobAddPending(job)) { - psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - psFree (job); - return false; - } - psFree(job); - } else { + // allocate a job -- if threads are not defined, this just runs the job + psThreadJob *job = psThreadJobAlloc ("PSPHOT_GUESS_MODEL"); + psArrayAdd(job->args, 1, readout); + psArrayAdd(job->args, 1, cells->data[j]); // sources + psArrayAdd(job->args, 1, psf); + + // XXX change these to use abstract mask type info + PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); + PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK); + + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + psFree (job); + return false; + } + psFree(job); + +# if (0) if (!psphotGuessModel_Unthreaded (readout, cells->data[j], psf, maskVal, markVal)) { psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); return false; } - } - } - - if (nThreads) { - // wait for the threads to finish and manage results - // wait here for the threaded jobs to finish - // fprintf (stderr, "wait for threads (%d, %d)\n", jx, jy); - if (!psThreadPoolWait (false)) { - psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - return false; - } - - // we have only supplied one type of job, so we can assume the types here - psThreadJob *job = NULL; - while ((job = psThreadJobGetDone()) != NULL) { - // we have no returned data from this operation - if (job->args->n < 1) { - fprintf (stderr, "error with job\n"); - } - psFree(job); - } +# endif + } + + // wait for the threads to finish and manage results + // wait here for the threaded jobs to finish + // fprintf (stderr, "wait for threads (%d, %d)\n", jx, jy); + if (!psThreadPoolWait (false)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + return false; + } + + // we have only supplied one type of job, so we can assume the types here + psThreadJob *job = NULL; + while ((job = psThreadJobGetDone()) != NULL) { + // we have no returned data from this operation + if (job->args->n < 1) { + fprintf (stderr, "error with job\n"); + } + psFree(job); } } @@ -116,14 +112,11 @@ int nMiss = 0; for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { - source->mode &= ~PM_SOURCE_MODE_EXT_LIMIT; + if (source->tmpFlags & PM_SOURCE_TMPF_MODEL_GUESS) { continue; } - nMiss ++; } - psLogMsg ("psphot.models", 4, "failed to build models for %d objects\n", nMiss); + psAssert (nMiss == 0, "failed to attempt to build models for %d objects\n", nMiss); psFree (cellGroups); @@ -148,7 +141,7 @@ pmSource *source = sources->data[i]; - // XXXX this is just for a test: use this to mark sources for which the model is measured - // check later that all are used. - source->mode |= PM_SOURCE_MODE_EXT_LIMIT; + // this is used to mark sources for which the model is measured. We check later that + // all are used. + source->tmpFlags |= PM_SOURCE_TMPF_MODEL_GUESS; // skip non-astronomical objects (very likely defects) @@ -185,20 +178,16 @@ pmModel *modelPSF = pmModelFromPSF (modelEXT, psf); // ALLOC X5 if (modelPSF == NULL) { - psError(PSPHOT_ERR_PSF, false, - "Failed to determine PSF model at r,c = (%d,%d); trying centre of image", + psWarning ("Failed to determine PSF model at r,c = (%d,%d); trying centre of image", source->peak->y, source->peak->x); - // - // Try the centre of the image - // + + // Try the center of the image modelEXT->params->data.F32[PM_PAR_XPOS] = 0.5*readout->image->numCols; modelEXT->params->data.F32[PM_PAR_YPOS] = 0.5*readout->image->numRows; modelPSF = pmModelFromPSF (modelEXT, psf); if (modelPSF == NULL) { - psError(PSPHOT_ERR_PSF, false, - "Failed to determine PSF model at centre of image"); + psError(PSPHOT_ERR_PSF, false, "Failed to determine PSF model at center of image"); psFree(modelEXT); return false; } - source->mode |= PM_SOURCE_MODE_BADPSF; } @@ -221,4 +210,5 @@ } +# if (0) // construct models only for sources in the specified region bool psphotGuessModel_Unthreaded (pmReadout *readout, psArray *sources, pmPSF *psf, psImageMaskType maskVal, psImageMaskType markVal) { @@ -301,2 +291,3 @@ return true; } +# endif Index: trunk/psphot/src/psphotMagnitudes.c =================================================================== --- trunk/psphot/src/psphotMagnitudes.c (revision 21458) +++ trunk/psphot/src/psphotMagnitudes.c (revision 21519) @@ -1,5 +1,3 @@ # include "psphotInternal.h" - -bool psphotMagnitudes_Unthreaded (int *nap, psArray *sources, pmPSF *psf, psImageBinning *binning, pmReadout *backModel, pmReadout *backStdev, pmSourcePhotometryMode photMode, psImageMaskType maskVal); bool psphotMagnitudes(pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources, pmPSF *psf) { @@ -59,25 +57,25 @@ for (int j = 0; j < cells->n; j++) { - if (nThreads) { - // allocate a job -- if threads are not defined, this just runs the job - psThreadJob *job = psThreadJobAlloc ("PSPHOT_MAGNITUDES"); - - psArrayAdd(job->args, 1, cells->data[j]); // sources - psArrayAdd(job->args, 1, psf); - psArrayAdd(job->args, 1, binning); - psArrayAdd(job->args, 1, backModel); - psArrayAdd(job->args, 1, backStdev); - - PS_ARRAY_ADD_SCALAR(job->args, photMode, PS_TYPE_S32); - PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); - PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for nAp - - if (!psThreadJobAddPending(job)) { - psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - psFree (job); - return false; - } - psFree(job); - } else { + // allocate a job -- if threads are not defined, this just runs the job + psThreadJob *job = psThreadJobAlloc ("PSPHOT_MAGNITUDES"); + + psArrayAdd(job->args, 1, cells->data[j]); // sources + psArrayAdd(job->args, 1, psf); + psArrayAdd(job->args, 1, binning); + psArrayAdd(job->args, 1, backModel); + psArrayAdd(job->args, 1, backStdev); + + PS_ARRAY_ADD_SCALAR(job->args, photMode, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for nAp + + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + psFree (job); + return false; + } + psFree(job); + +# if (0) int nap = 0; if (!psphotMagnitudes_Unthreaded (&nap, cells->data[j], psf, binning, backModel, backStdev, photMode, maskVal)) { @@ -86,25 +84,23 @@ } Nap += nap; - } - } - - if (nThreads) { - // wait for the threads to finish and manage results - if (!psThreadPoolWait (false)) { - psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - return false; - } - - // we have only supplied one type of job, so we can assume the types here - psThreadJob *job = NULL; - while ((job = psThreadJobGetDone()) != NULL) { - if (job->args->n < 1) { - fprintf (stderr, "error with job\n"); - } else { - psScalar *scalar = job->args->data[7]; - Nap += scalar->data.S32; - } - psFree(job); - } +# endif + } + + // wait for the threads to finish and manage results + if (!psThreadPoolWait (false)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + return false; + } + + // we have only supplied one type of job, so we can assume the types here + psThreadJob *job = NULL; + while ((job = psThreadJobGetDone()) != NULL) { + if (job->args->n < 1) { + fprintf (stderr, "error with job\n"); + } else { + psScalar *scalar = job->args->data[7]; + Nap += scalar->data.S32; + } + psFree(job); } } @@ -162,4 +158,5 @@ } +# if (0) bool psphotMagnitudes_Unthreaded (int *nap, psArray *sources, pmPSF *psf, psImageBinning *binning, pmReadout *backModel, pmReadout *backStdev, pmSourcePhotometryMode photMode, psImageMaskType maskVal) { @@ -198,4 +195,5 @@ return true; } +# endif bool psphotPSFWeights(pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources) { @@ -295,5 +293,5 @@ } - status = pmSourcePixelWeight (&source->pixWeight, model, source->pixels, source->maskObj, maskVal); + status = pmSourcePixelWeight (&source->pixWeight, model, source->maskObj, maskVal); if (!status) { psTrace ("psphot", 3, "fail to measure pixel weight"); Index: trunk/psphot/src/psphotMergeSources.c =================================================================== --- trunk/psphot/src/psphotMergeSources.c (revision 21458) +++ trunk/psphot/src/psphotMergeSources.c (revision 21519) @@ -122,2 +122,42 @@ return detections; } + +// generate the detection structure for the supplied array of sources +bool psphotSetSourceParams (pmConfig *config, psArray *sources, pmPSF *psf) { + + psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE); + if (!recipe) { + psError(PSPHOT_ERR_CONFIG, false, "missing recipe %s", PSPHOT_RECIPE); + return false; + } + + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + pmModel *model = source->modelPSF; + + if (source->mode & PSF_SOURCE_MASK || !isfinite(source->psfMag)) { + continue; + } + + float flux = powf(10.0, -0.4 * source->psfMag); + float xpos = model->params->data.F32[PM_PAR_XPOS]; + float ypos = model->params->data.F32[PM_PAR_YPOS]; + + pmPeak *peak = pmPeakAlloc(xpos, ypos, 1.0, PM_PEAK_LONE); + peak->xf = xpos; + peak->yf = ypos; + peak->flux = flux; // this are being set wrong, but does it matter? + + if (isfinite (source->errMag) && (source->errMag > 0.0)) { + peak->SN = 1.0 / source->errMag; + } else { + peak->SN = 0.0; + } + + source->peak = peak; + } + + psLogMsg ("psphot", 3, "%ld PSF sources loaded", sources->n); + + return true; +} Index: trunk/psphot/src/psphotReplaceUnfit.c =================================================================== --- trunk/psphot/src/psphotReplaceUnfit.c (revision 21458) +++ trunk/psphot/src/psphotReplaceUnfit.c (revision 21519) @@ -17,5 +17,5 @@ replace: pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); - source->mode &= ~PM_SOURCE_MODE_SUBTRACTED; + source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED; } psLogMsg ("psphot.replace", 3, "replace unfitted models: %f sec (%ld objects)\n", psTimerMark ("psphot.replace"), sources->n); @@ -38,8 +38,8 @@ // replace other sources? - if (!(source->mode & PM_SOURCE_MODE_SUBTRACTED)) continue; + if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) continue; pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); - source->mode &= ~PM_SOURCE_MODE_SUBTRACTED; + source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED; } psLogMsg ("psphot.replace", PS_LOG_INFO, "replaced models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.replace")); @@ -62,8 +62,8 @@ // replace other sources? - if (source->mode & PM_SOURCE_MODE_SUBTRACTED) continue; + if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) continue; pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); - source->mode |= PM_SOURCE_MODE_SUBTRACTED; + source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED; } psLogMsg ("psphot.replace", PS_LOG_INFO, "replaced models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.replace")); @@ -75,5 +75,5 @@ // what is current state? (true : add; false : sub) - bool state = !(source->mode & PM_SOURCE_MODE_SUBTRACTED); + bool state = !(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); if (state && useState) return true; @@ -86,5 +86,5 @@ // what is current state? (true : sub; false : add) - bool state = (source->mode & PM_SOURCE_MODE_SUBTRACTED); + bool state = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); if (state && useState) return true; @@ -97,5 +97,5 @@ // what is desired state? (true : add; false : sub) - bool newState = !(source->mode & PM_SOURCE_MODE_SUBTRACTED); + bool newState = !(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); if (curState == newState) return true; Index: trunk/psphot/src/psphotSourceFits.c =================================================================== --- trunk/psphot/src/psphotSourceFits.c (revision 21458) +++ trunk/psphot/src/psphotSourceFits.c (revision 21519) @@ -120,5 +120,6 @@ pmSourceCacheModel (blend, maskVal); pmSourceSub (blend, PM_MODEL_OP_FULL, maskVal); - blend->mode |= PM_SOURCE_MODE_SUBTRACTED; + blend->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED; + blend->mode |= PM_SOURCE_MODE_BLEND_FIT; } NfitBlend += modelSet->n; @@ -143,5 +144,6 @@ pmSourceCacheModel (source, maskVal); pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - source->mode |= PM_SOURCE_MODE_SUBTRACTED; + source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED; + source->mode |= PM_SOURCE_MODE_BLEND_FIT; return true; } @@ -184,5 +186,5 @@ pmSourceCacheModel (source, maskVal); pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - source->mode |= PM_SOURCE_MODE_SUBTRACTED; + source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED; return true; } @@ -285,5 +287,5 @@ pmSourceCacheModel (source, maskVal); pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - source->mode |= PM_SOURCE_MODE_SUBTRACTED; + source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED; psTrace ("psphot", 5, "blob as EXT: %f %f\n", EXT->params->data.F32[PM_PAR_XPOS], EXT->params->data.F32[PM_PAR_YPOS]); return true; @@ -296,6 +298,6 @@ psFree (source->modelPSF); source->modelPSF = psMemIncrRefCounter (DBL->data[0]); - source->mode |= PM_SOURCE_MODE_SUBTRACTED; - source->mode |= PM_SOURCE_MODE_PAIR; + source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED; + source->mode |= PM_SOURCE_MODE_PAIR; // copy most data from the primary source (modelEXT, blends stay NULL) Index: trunk/psphot/src/psphotSourceSize.c =================================================================== --- trunk/psphot/src/psphotSourceSize.c (revision 21458) +++ trunk/psphot/src/psphotSourceSize.c (revision 21519) @@ -56,5 +56,6 @@ // source must have been subtracted - if (!(source->mode & PM_SOURCE_MODE_SUBTRACTED)) { + if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) { + source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED; psTrace("psphot", 7, "Not calculating extNsigma,crNsigma since source is not subtracted\n"); continue; @@ -63,5 +64,5 @@ psF32 **resid = source->pixels->data.F32; psF32 **variance = source->variance->data.F32; - psImageMaskType **mask = source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA; + psImageMaskType **mask = source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA; // check for extendedness: measure the delta flux significance at the 1 sigma contour @@ -82,4 +83,5 @@ if (xPeak < 1 || xPeak > source->pixels->numCols - 2 || yPeak < 1 || yPeak > source->pixels->numRows - 2) { + source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED; psTrace("psphot", 7, "Not calculating crNsigma due to edge\n"); continue; @@ -98,4 +100,5 @@ if (!keep) { psTrace("psphot", 7, "Not calculating crNsigma due to masked pixels\n"); + source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED; continue; } @@ -156,5 +159,7 @@ } source->crNsigma = (nCR > 0) ? fCR / nCR : 0.0; - if (!isfinite(source->crNsigma)) continue; + if (!isfinite(source->crNsigma)) { + continue; + } // this source is thought to be a cosmic ray. flag the detection and mask the pixels @@ -264,5 +269,5 @@ // replace the source flux pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); - source->mode &= ~PM_SOURCE_MODE_SUBTRACTED; + source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED; // flag this as a CR Index: trunk/psphot/src/psphotSourceStats.c =================================================================== --- trunk/psphot/src/psphotSourceStats.c (revision 21458) +++ trunk/psphot/src/psphotSourceStats.c (revision 21519) @@ -1,5 +1,3 @@ # include "psphotInternal.h" - -bool psphotSourceStats_Unthreaded (int *nfail, int *nmoments, psArray *sources, psMetadata *recipe); psArray *psphotSourceStats (pmConfig *config, pmReadout *readout, pmDetections *detections) { @@ -80,20 +78,20 @@ for (int j = 0; j < cells->n; j++) { - if (nThreads) { - // allocate a job -- if threads are not defined, this just runs the job - psThreadJob *job = psThreadJobAlloc ("PSPHOT_SOURCE_STATS"); - - psArrayAdd(job->args, 1, cells->data[j]); // sources - psArrayAdd(job->args, 1, recipe); - PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nmoments - PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail - - if (!psThreadJobAddPending(job)) { - psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - psFree (job); - return NULL; - } - psFree(job); - } else { + // allocate a job -- if threads are not defined, this just runs the job + psThreadJob *job = psThreadJobAlloc ("PSPHOT_SOURCE_STATS"); + + psArrayAdd(job->args, 1, cells->data[j]); // sources + psArrayAdd(job->args, 1, recipe); + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nmoments + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail + + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + psFree (job); + return NULL; + } + psFree(job); + +# if (0) int nfail = 0; int nmoments = 0; @@ -104,28 +102,26 @@ Nfail += nfail; Nmoments += nmoments; +# endif + } + + // wait for the threads to finish and manage results + if (!psThreadPoolWait (false)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + return NULL; + } + + // we have only supplied one type of job, so we can assume the types here + psThreadJob *job = NULL; + while ((job = psThreadJobGetDone()) != NULL) { + if (job->args->n < 1) { + fprintf (stderr, "error with job\n"); + } else { + psScalar *scalar = NULL; + scalar = job->args->data[2]; + Nmoments += scalar->data.S32; + scalar = job->args->data[3]; + Nfail += scalar->data.S32; } - } - - if (nThreads) { - // wait for the threads to finish and manage results - if (!psThreadPoolWait (false)) { - psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - return NULL; - } - - // we have only supplied one type of job, so we can assume the types here - psThreadJob *job = NULL; - while ((job = psThreadJobGetDone()) != NULL) { - if (job->args->n < 1) { - fprintf (stderr, "error with job\n"); - } else { - psScalar *scalar = NULL; - scalar = job->args->data[2]; - Nmoments += scalar->data.S32; - scalar = job->args->data[3]; - Nfail += scalar->data.S32; - } - psFree(job); - } + psFree(job); } } @@ -148,4 +144,94 @@ psArray *sources = job->args->data[0]; psMetadata *recipe = job->args->data[1]; + + float INNER = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS"); + if (!status) return false; + float MIN_SN = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN"); + if (!status) return false; + float RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS"); + if (!status) return false; + + // bit-masks to test for good/bad pixels + psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); + assert (maskVal); + + // bit-mask to mark pixels not used in analysis + psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); + assert (markVal); + + // maskVal is used to test for rejected pixels, and must include markVal + maskVal |= markVal; + + // threaded measurement of the sources moments + int Nfail = 0; + int Nmoments = 0; + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + + // skip faint sources for moments measurement + if (source->peak->SN < MIN_SN) { + source->mode |= PM_SOURCE_MODE_BELOW_MOMENTS_SN; + continue; + } + + // measure a local sky value + // the local sky is now ignored; kept here for reference only + status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal); + if (!status) { + source->mode |= PM_SOURCE_MODE_SKY_FAILURE; + psErrorClear(); // XXX re-consider the errors raised here + Nfail ++; + continue; + } + + // measure the local sky variance (needed if noise is not sqrt(signal)) + // XXX EAM : this should use ROBUST not SAMPLE median, but it is broken + status = pmSourceLocalSkyVariance (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal); + if (!status) { + source->mode |= PM_SOURCE_MODE_SKYVAR_FAILURE; + Nfail ++; + psErrorClear(); + continue; + } + + // measure basic source moments + status = pmSourceMoments (source, RADIUS); + if (status) { + Nmoments ++; + continue; + } + + // if no valid pixels, or massive swing, likely saturated source, + // try a much larger box + BIG_RADIUS = PS_MIN (INNER, 3*RADIUS); + psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y); + status = pmSourceMoments (source, BIG_RADIUS); + if (status) { + source->mode |= PM_SOURCE_MODE_BIG_RADIUS; + Nmoments ++; + continue; + } + + source->mode |= PM_SOURCE_MODE_MOMENTS_FAILURE; + Nfail ++; + psErrorClear(); + continue; + } + + // change the value of a scalar on the array (wrap this and put it in psArray.h) + scalar = job->args->data[2]; + scalar->data.S32 = Nmoments; + + scalar = job->args->data[3]; + scalar->data.S32 = Nfail; + + return true; +} + +# if (0) +bool psphotSourceStats_Unthreaded (int *nfail, int *nmoments, psArray *sources, psMetadata *recipe) { + + bool status = false; + float BIG_RADIUS; float INNER = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS"); @@ -219,88 +305,4 @@ // change the value of a scalar on the array (wrap this and put it in psArray.h) - scalar = job->args->data[2]; - scalar->data.S32 = Nmoments; - - scalar = job->args->data[3]; - scalar->data.S32 = Nfail; - - return true; -} - -bool psphotSourceStats_Unthreaded (int *nfail, int *nmoments, psArray *sources, psMetadata *recipe) { - - bool status = false; - float BIG_RADIUS; - - float INNER = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS"); - if (!status) return false; - float MIN_SN = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN"); - if (!status) return false; - float RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS"); - if (!status) return false; - - // bit-masks to test for good/bad pixels - psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); - assert (maskVal); - - // bit-mask to mark pixels not used in analysis - psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); - assert (markVal); - - // maskVal is used to test for rejected pixels, and must include markVal - maskVal |= markVal; - - // threaded measurement of the sources moments - int Nfail = 0; - int Nmoments = 0; - for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - - // skip faint sources for moments measurement - if (source->peak->SN < MIN_SN) { - continue; - } - - // measure a local sky value - // the local sky is now ignored; kept here for reference only - status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal); - if (!status) { - psErrorClear(); // XXX re-consider the errors raised here - Nfail ++; - continue; - } - - // measure the local sky variance (needed if noise is not sqrt(signal)) - // XXX EAM : this should use ROBUST not SAMPLE median, but it is broken - status = pmSourceLocalSkyVariance (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal); - if (!status) { - Nfail ++; - psErrorClear(); - continue; - } - - // measure basic source moments - status = pmSourceMoments (source, RADIUS); - if (status) { - Nmoments ++; - continue; - } - - // if no valid pixels, or massive swing, likely saturated source, - // try a much larger box - BIG_RADIUS = PS_MIN (INNER, 3*RADIUS); - psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y); - status = pmSourceMoments (source, BIG_RADIUS); - if (status) { - Nmoments ++; - continue; - } - - Nfail ++; - psErrorClear(); - continue; - } - - // change the value of a scalar on the array (wrap this and put it in psArray.h) *nmoments = Nmoments; *nfail = Nfail; @@ -308,2 +310,3 @@ return true; } +# endif Index: trunk/psphot/src/psphotVisual.c =================================================================== --- trunk/psphot/src/psphotVisual.c (revision 21458) +++ trunk/psphot/src/psphotVisual.c (revision 21519) @@ -1101,5 +1101,5 @@ Graphdata graphdata; - bool state = !(source->mode & PM_SOURCE_MODE_SUBTRACTED); + bool state = !(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); psphotAddWithTest (source, true, maskVal); // replace source if subtracted