Index: trunk/psphot/src/psphotApResid.c =================================================================== --- trunk/psphot/src/psphotApResid.c (revision 21357) +++ trunk/psphot/src/psphotApResid.c (revision 21359) @@ -1,3 +1,5 @@ # 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; } @@ -29,5 +31,4 @@ nThreads = 0; } - // nThreads = 0; // XXX until testing is complete, do not thread this function bool measureAptrend = psMetadataLookupBool (&status, recipe, "MEASURE.APTREND"); @@ -94,43 +95,57 @@ for (int j = 0; j < cells->n; j++) { - // 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)) { + 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 { + int nskip = 0; + int nfail = 0; + + if (!psphotApResidMags_Unthreaded (&nskip, &nfail, sources, psf, photMode, maskVal)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + return false; + } + Nskip += nskip; + Nfail += nfail; + } + + } + + if (nThreads) { + // wait for the threads to finish and manage results + if (!psThreadPoolWait (false)) { psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - psFree (job); return false; } - psFree(job); - - } - - // 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); } } @@ -493,2 +508,36 @@ return true; } + +bool psphotApResidMags_Unthreaded (int *nskip, int *nfail, psArray *sources, pmPSF *psf, pmSourcePhotometryMode photMode, psImageMaskType maskVal) { + + int Nskip = 0; + int Nfail = 0; + + for (int i = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + + if (source->type != PM_SOURCE_TYPE_STAR) SKIPSTAR ("NOT STAR"); + if (source->mode & PM_SOURCE_MODE_SATSTAR) SKIPSTAR ("SATSTAR"); + if (source->mode & PM_SOURCE_MODE_BLEND) SKIPSTAR ("BLEND"); + if (source->mode & PM_SOURCE_MODE_FAIL) SKIPSTAR ("FAIL STAR"); + if (source->mode & PM_SOURCE_MODE_POOR) SKIPSTAR ("POOR STAR"); + + if (!pmSourceMagnitudes (source, psf, photMode, maskVal)) { + Nskip ++; + psTrace ("psphot", 3, "skip : bad source mag"); + continue; + } + + if (!isfinite(source->apMag) || !isfinite(source->psfMag)) { + Nfail ++; + psTrace ("psphot", 3, "fail : nan mags : %f %f", source->apMag, source->psfMag); + continue; + } + } + + // change the value of a scalar on the array (wrap this and put it in psArray.h) + *nskip = Nskip; + *nfail = Nfail; + + return true; +} Index: trunk/psphot/src/psphotBlendFit.c =================================================================== --- trunk/psphot/src/psphotBlendFit.c (revision 21357) +++ trunk/psphot/src/psphotBlendFit.c (revision 21359) @@ -1,3 +1,5 @@ # 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 @@ -21,5 +23,4 @@ nThreads = 0; } - // nThreads = 0; // XXX until testing is complete, do not thread this function // source fitting parameters for extended source fits @@ -63,58 +64,82 @@ for (int j = 0; j < cells->n; j++) { - // 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)) { + 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 { + int nfit = 0; + int npsf = 0; + int next = 0; + int nfail = 0; + psArray *newSources = psArrayAllocEmpty(16); + + if (!psphotBlendFit_Unthreaded (&nfit, &npsf, &next, &nfail, readout, recipe, cells->data[j], psf, newSources)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + return NULL; + } + Nfit += nfit; + Npsf += npsf; + Next += next; + Nfail += nfail; + + // add these back onto sources + for (int k = 0; k < newSources->n; k++) { + psArrayAdd (sources, 16, newSources->data[k]); + } + 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."); - psFree (job); return NULL; } - psFree(job); - - } - - // 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]); + + // 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(job); + } } } @@ -249,2 +274,105 @@ 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; + } + + // change the value of a scalar on the array (wrap this and put it in psArray.h) + *nfit = Nfit; + *npsf = Npsf; + *next = Next; + *nfail = Nfail; + + return true; +} Index: trunk/psphot/src/psphotGuessModels.c =================================================================== --- trunk/psphot/src/psphotGuessModels.c (revision 21357) +++ trunk/psphot/src/psphotGuessModels.c (revision 21359) @@ -22,4 +22,6 @@ } +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) { @@ -38,5 +40,4 @@ nThreads = 0; } - // nThreads = 0; // XXX until testing is complete, do not thread this function // bit-masks to test for good/bad pixels @@ -66,38 +67,47 @@ for (int j = 0; j < cells->n; j++) { - // 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)) { + 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 { + 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."); - psFree (job); return false; } - psFree(job); - } - - // 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); + + // 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); + } } } @@ -211,2 +221,82 @@ } +// construct models only for sources in the specified region +bool psphotGuessModel_Unthreaded (pmReadout *readout, psArray *sources, pmPSF *psf, psImageMaskType maskVal, psImageMaskType markVal) { + + int nSrc = 0; + + for (int i = 0; i < sources->n; i++) { + 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; + + // 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->peak) continue; + + nSrc ++; + + // XXX if a source is faint, it will not have moments measured. + // it must be modelled as a PSF. In this case, we need to use + // the peak centroid to get the coordinates and get the peak flux + // from the image? + pmModel *modelEXT; + if (!source->moments) { + modelEXT = wildGuess(source, psf); + } else { + // use the source moments, etc to guess basic model parameters + modelEXT = pmSourceModelGuess (source, psf->type); // ALLOC + if (!modelEXT) { + modelEXT = wildGuess(source, psf); + } + // these valuse are set in pmSourceModelGuess, should this rule be in there as well? + if (source->mode & PM_SOURCE_MODE_SATSTAR) { + modelEXT->params->data.F32[PM_PAR_XPOS] = source->moments->Mx; + modelEXT->params->data.F32[PM_PAR_YPOS] = source->moments->My; + } else { + modelEXT->params->data.F32[PM_PAR_XPOS] = source->peak->xf; + modelEXT->params->data.F32[PM_PAR_YPOS] = source->peak->yf; + } + } + + // set PSF parameters for this model (apply 2D shape model) + pmModel *modelPSF = pmModelFromPSF (modelEXT, psf); // ALLOC + if (modelPSF == NULL) { + psError(PSPHOT_ERR_PSF, false, + "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 + // + 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"); + psFree(modelEXT); + return false; + } + + source->mode |= PM_SOURCE_MODE_BADPSF; + } + psFree (modelEXT); + + // XXX need to define the guess flux? + // set the fit radius based on the object flux limit and the model + // this function affects the mask pixels + psphotCheckRadiusPSF (readout, source, modelPSF, markVal); + + // set the source PSF model + source->modelPSF = modelPSF; + source->modelPSF->residuals = psf->residuals; + + pmSourceCacheModel (source, maskVal); + + } + + return true; +} Index: trunk/psphot/src/psphotMagnitudes.c =================================================================== --- trunk/psphot/src/psphotMagnitudes.c (revision 21357) +++ trunk/psphot/src/psphotMagnitudes.c (revision 21359) @@ -1,3 +1,5 @@ # 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) { @@ -17,5 +19,4 @@ nThreads = 0; } - // nThreads = 0; // XXX until testing is complete, do not thread this function // if backModel or backStdev are missing, the values of sky and/or skyErr will be set to NAN @@ -58,42 +59,52 @@ for (int j = 0; j < cells->n; j++) { - // 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_U8); // XXX change this to use abstract mask type info - PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for nAp - - if (!psThreadJobAddPending(job)) { + 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 { + int nap = 0; + if (!psphotMagnitudes_Unthreaded (&nap, cells->data[j], psf, binning, backModel, backStdev, photMode, maskVal)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + return false; + } + 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."); - psFree (job); return false; } - psFree(job); - - } - - // 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); + + // 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); + } } } @@ -116,5 +127,5 @@ pmReadout *backStdev = job->args->data[4]; pmSourcePhotometryMode photMode = PS_SCALAR_VALUE(job->args->data[5],S32); - psMaskType maskVal = PS_SCALAR_VALUE(job->args->data[6],U8); + psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA); for (int i = 0; i < sources->n; i++) { @@ -150,2 +161,158 @@ return true; } + +bool psphotMagnitudes_Unthreaded (int *nap, psArray *sources, pmPSF *psf, psImageBinning *binning, pmReadout *backModel, pmReadout *backStdev, pmSourcePhotometryMode photMode, psImageMaskType maskVal) { + + bool status; + int Nap = 0; + + for (int i = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + status = pmSourceMagnitudes (source, psf, photMode, maskVal); + if (status && isfinite(source->apMag)) Nap ++; + + if (backModel) { + psAssert (binning, "if backModel is defined, so should binning be"); + source->sky = psImageUnbinPixel(source->peak->x, source->peak->y, backModel->image, binning); + if (isnan(source->sky) && false) { + psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "error setting pmSource.sky"); + } + } else { + source->sky = NAN; + } + + if (backStdev) { + psAssert (binning, "if backStdev is defined, so should binning be"); + source->skyErr = psImageUnbinPixel(source->peak->x, source->peak->y, backStdev->image, binning); + if (isnan(source->skyErr) && false) { + psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "error setting pmSource.skyErr"); + } + } else { + source->skyErr = NAN; + } + } + + // change the value of a scalar on the array (wrap this and put it in psArray.h) + *nap = Nap; + + return true; +} + +# if (0) +bool psphotPSFWeights(pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources, pmPSF *psf) { + + bool status = false; + int Nap = 0; + + psTimerStart ("psphot.mags"); + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + assert (recipe); + + // determine the number of allowed threads + int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads + if (!status) { + nThreads = 0; + } + nThreads = 0; // XXX until testing is complete, do not thread this function + + // 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; + + pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_APCORR | PM_SOURCE_PHOT_WEIGHT; + if (!IGNORE_GROWTH) photMode |= PM_SOURCE_PHOT_GROWTH; + if (INTERPOLATE_AP) photMode |= PM_SOURCE_PHOT_INTERP; + + // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts) + int Cx = 1, Cy = 1; + psphotChooseCellSizes (&Cx, &Cy, readout, nThreads); + + psArray *cellGroups = psphotAssignSources (Cx, Cy, sources); + + for (int i = 0; i < cellGroups->n; i++) { + + psArray *cells = cellGroups->data[i]; + + for (int j = 0; j < cells->n; j++) { + + // allocate a job -- if threads are not defined, this just runs the job + psThreadJob *job = psThreadJobAlloc ("PSPHOT_PSF_WEIGHTS"); + + psArrayAdd(job->args, 1, cells->data[j]); // sources + psArrayAdd(job->args, 1, psf); + PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); + + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + psFree (job); + return false; + } + psFree(job); + + } + + // 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"); + } + psFree(job); + } + } + + psFree (cellGroups); + + psLogMsg ("psphot", PS_LOG_DETAIL, "measure psf weights : %f sec for %ld objects\n", psTimerMark ("psphot.mags"), sources->n); + return true; +} + +bool psphotPSFWeights_Threaded (psThreadJob *job) { + + bool status; + bool isPSF; + + psArray *sources = job->args->data[0]; + pmPSF *psf = job->args->data[1]; + psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[2],PS_TYPE_IMAGE_MASK_DATA); + + for (int i = 0; i < sources->n; i++) { + pmSource *source = (pmSource *) sources->data[i]; + + // we must have a valid model + pmModel *model = pmSourceGetModel (&isPSF, source); + if (model == NULL) { + psTrace ("psphot", 3, "fail mag : no valid model"); + source->pixWeight = NAN; + continue; + } + + status = pmSourcePixelWeight (&source->pixWeight, model, source->pixels, source->maskObj, maskVal); + if (!status) { + psTrace ("psphot", 3, "fail to measure pixel weight"); + source->pixWeight = NAN; + continue; + } + + } + + return true; +} + +# endif + Index: trunk/psphot/src/psphotSourceStats.c =================================================================== --- trunk/psphot/src/psphotSourceStats.c (revision 21357) +++ trunk/psphot/src/psphotSourceStats.c (revision 21359) @@ -1,3 +1,5 @@ # include "psphotInternal.h" + +bool psphotSourceStats_Unthreaded (int *nfail, int *nmoments, psArray *sources, psMetadata *recipe); psArray *psphotSourceStats (pmConfig *config, pmReadout *readout, pmDetections *detections) { @@ -78,40 +80,52 @@ for (int j = 0; j < cells->n; j++) { - // 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)) { + 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 { + int nfail = 0; + int nmoments = 0; + if (!psphotSourceStats_Unthreaded (&nfail, &nmoments, cells->data[j], recipe)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + return NULL; + } + Nfail += nfail; + Nmoments += nmoments; + } + } + + if (nThreads) { + // wait for the threads to finish and manage results + if (!psThreadPoolWait (false)) { psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - psFree (job); return NULL; } - psFree(job); - - } - - // 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; + + // 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); } } @@ -213,2 +227,83 @@ 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; + + return true; +}