Index: branches/simmosaic_branches/psphot/src/psphotGuessModels.c =================================================================== --- branches/simmosaic_branches/psphot/src/psphotGuessModels.c (revision 24860) +++ branches/simmosaic_branches/psphot/src/psphotGuessModels.c (revision 27839) @@ -3,29 +3,58 @@ // XXX : the threading here is not great. this may be due to blocks between elements, but // the selection of the objects in a cell is not optimal. To fix: -// 1) define the boundaries of the cells up front +// 1) define the boundaries of the cells up front // 2) loop over the sources once and associate them with their cell // 3) define the threaded function to work with sources for a given cell -// A guess for when the moments aren't available -static pmModel *wildGuess(pmSource *source, // Source for which to guess - pmPSF *psf // The point-spread function - ) +// for now, let's store the detections on the readout->analysis for each readout +bool psphotGuessModels (pmConfig *config, const pmFPAview *view) { - pmModel *model = pmModelAlloc(psf->type); - psF32 *PAR = model->params->data.F32; - PAR[PM_PAR_SKY] = 0; - // XXX get this from the image pixels - PAR[PM_PAR_I0] = source->peak->flux; - PAR[PM_PAR_XPOS] = source->peak->xf; - PAR[PM_PAR_YPOS] = source->peak->yf; - return model; + bool status = true; + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // skip the chisq image (optionally?) + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; // skip chisq image + if (!psphotGuessModelsReadout (config, view, "PSPHOT.INPUT", i)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on to guess models for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; } -// construct an initial PSF model for each object -bool psphotGuessModels (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf) { +// construct an initial PSF model for each object (new sources only) +bool psphotGuessModelsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index) { bool status; psTimerStart ("psphot.models"); + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->newSources; + psAssert (sources, "missing sources?"); + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping model guess"); + return true; + } + + pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); + psAssert (psf, "missing psf?"); // select the appropriate recipe information @@ -36,5 +65,5 @@ int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads if (!status) { - nThreads = 0; + nThreads = 0; } @@ -51,5 +80,5 @@ // setup the PSF fit radius details - psphotInitRadiusPSF (recipe, psf->type); + psphotInitRadiusPSF (recipe, readout->analysis, psf->type); // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts) @@ -61,50 +90,43 @@ 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_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; - } -# 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); - } + 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_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); + } + + // 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); + } } @@ -112,9 +134,9 @@ int nMiss = 0; for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - if (source->tmpFlags & PM_SOURCE_TMPF_MODEL_GUESS) { - continue; - } - nMiss ++; + pmSource *source = sources->data[i]; + if (source->tmpFlags & PM_SOURCE_TMPF_MODEL_GUESS) { + continue; + } + nMiss ++; } psAssert (nMiss == 0, "failed to attempt to build models for %d objects\n", nMiss); @@ -132,5 +154,5 @@ psArray *sources = job->args->data[1]; pmPSF *psf = job->args->data[2]; - + psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[3],PS_TYPE_IMAGE_MASK_DATA); psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA); @@ -139,155 +161,71 @@ for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - - // 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) - 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 X5 - 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 X5 - if (modelPSF == NULL) { - psWarning ("Failed to determine PSF model at r,c = (%d,%d); trying centre of image", - source->peak->y, source->peak->x); - - // 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 center of image"); - psFree(modelEXT); - return false; - } - source->mode |= PM_SOURCE_MODE_BADPSF; - } - psFree (modelEXT); // FREE (x3) - - // 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); // ALLOC x14 (!) - + pmSource *source = sources->data[i]; + + // 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) + if (source->type == PM_SOURCE_TYPE_DEFECT) continue; + if (source->type == PM_SOURCE_TYPE_SATURATED) continue; + if (!source->peak) continue; + + nSrc ++; + + // the guess central intensity comes from the peak: + float Io = source->peak->flux; + + // We have two options to get a guess for the object position: the position from the + // peak and the position from the moments. Use the peak position if (a) there are no + // moments and (b) the sources is not saturated + + bool useMoments = false; + useMoments = (source->mode & PM_SOURCE_MODE_SATSTAR); // we only want to try if SATSTAR is set, but.. + useMoments = (useMoments && source->moments); // can't if there are no moments + useMoments = (useMoments && source->moments->nPixels); // can't if the moments were not measured + useMoments = (useMoments && !(source->mode && PM_SOURCE_MODE_MOMENTS_FAILURE)); // can't if the moments failed... + + float Xo, Yo; + if (useMoments) { + Xo = source->moments->Mx; + Yo = source->moments->My; + } else { + Xo = source->peak->xf; + Yo = source->peak->yf; + } + + // set PSF parameters for this model (apply 2D shape model to coordinates Xo, Yo) + pmModel *modelPSF = pmModelFromPSFforXY(psf, Xo, Yo, Io); + + if (modelPSF == NULL) { + psWarning ("Failed to determine PSF model at (%f,%f); trying image center", Xo, Yo); + + float Xc = 0.5*readout->image->numCols; + float Yc = 0.5*readout->image->numRows; + pmModel *modelPSF = pmModelFromPSFforXY(psf, Xc, Yc, Io); + if (modelPSF == NULL) { + psError(PSPHOT_ERR_PSF, false, "Failed to determine PSF model at center of image"); + return false; + } + + // Now set the object position at the expected location: + modelPSF->params->data.F32[PM_PAR_XPOS] = Xo; + modelPSF->params->data.F32[PM_PAR_YPOS] = Yo; + source->mode |= PM_SOURCE_MODE_BADPSF; + } + + // 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 + psAssert (source->modelPSF == NULL, "failed to free one of the models?"); + source->modelPSF = modelPSF; + source->modelPSF->residuals = psf->residuals; + + pmSourceCacheModel (source, maskVal); // ALLOC x14 (!) } return true; } - -# if (0) -// 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; -} -# endif