Index: branches/eam_branches/ipp-20130509/psphot/src/psphot.h =================================================================== --- branches/eam_branches/ipp-20130509/psphot/src/psphot.h (revision 35634) +++ branches/eam_branches/ipp-20130509/psphot/src/psphot.h (revision 35640) @@ -67,4 +67,5 @@ bool psphotModelBackground (pmConfig *config, const pmFPAview *view, const char *filerule); bool psphotModelBackgroundReadoutFileIndex (pmConfig *config, const pmFPAview *view, const char *filerule, int index); +bool psphotModelBackground_Threaded (psThreadJob *job); bool psphotMaskBackground (pmConfig *config, const pmFPAview *view, const char *filerule); @@ -119,4 +120,5 @@ bool psphotAddOrSubNoise (pmConfig *config, const pmFPAview *view, const char *filerule, bool add); bool psphotAddOrSubNoiseReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool add); +bool psphotAddOrSubNoise_Threaded (psThreadJob *job); bool psphotExtendedSourceAnalysis (pmConfig *config, const pmFPAview *view, const char *filerule); @@ -243,4 +245,6 @@ bool psphotFitInit (int nThreads); bool psphotFitSummary (void); +bool psphotFitSummaryExtended (void); +bool psphotFitInitExtended (void); bool psphotLoadPSF (pmConfig *config, const pmFPAview *view, const char *filerule); Index: branches/eam_branches/ipp-20130509/psphot/src/psphotAddNoise.c =================================================================== --- branches/eam_branches/ipp-20130509/psphot/src/psphotAddNoise.c (revision 35634) +++ branches/eam_branches/ipp-20130509/psphot/src/psphotAddNoise.c (revision 35640) @@ -1,4 +1,5 @@ # include "psphotInternal.h" +bool psphotAddOrSubNoise_Threaded (psThreadJob *job); bool psphotMaskSource(pmSource *source, bool add, psImageMaskType maskVal); @@ -32,10 +33,10 @@ } -static int Nmasked = 0; - // the return state indicates if any sources were actually replaced bool psphotAddOrSubNoiseReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool add) { bool status = false; + + psTimerStart ("psphot.noise"); // find the currently selected readout @@ -55,5 +56,9 @@ if (!sources) return false; - psTimerStart ("psphot.noise"); + // determine the number of allowed threads + int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads + if (!status) { + nThreads = 0; + } // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) @@ -84,4 +89,79 @@ psphotVisualShowImage (readout); + + // 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_ADD_NOISE"); + psArrayAdd(job->args, 1, cells->data[j]); // sources + PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); + PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK); + PS_ARRAY_ADD_SCALAR(job->args, FACTOR, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR(job->args, SIZE, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR(job->args, add, PS_TYPE_U8); + +# if (1) + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to add/sub noise."); + return false; + } +# else + if (!psphotAddOrSubNoise_Threaded(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to add/sub noise."); + return false; + } + psFree(job); +# endif + } + + // wait for the threads to finish and manage results + if (!psThreadPoolWait (false, true)) { + psFree(cellGroups); + psError(PS_ERR_UNKNOWN, false, "Unable to add/sub noise."); + 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); + } + } + + if (add) { + psLogMsg ("psphot.noise", PS_LOG_WARN, "add noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.noise")); + } else { + psLogMsg ("psphot.noise", PS_LOG_WARN, "sub noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.noise")); + } + + psFree (cellGroups); + + psphotVisualShowImage (readout); + + return true; +} + +bool psphotAddOrSubNoise_Threaded (psThreadJob *job) { + + psArray *sources = job->args->data[0]; + psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[1], PS_TYPE_IMAGE_MASK_DATA); + psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[2], PS_TYPE_IMAGE_MASK_DATA); + psF32 FACTOR = PS_SCALAR_VALUE(job->args->data[3], F32); + psF32 SIZE = PS_SCALAR_VALUE(job->args->data[4], F32); + psBool add = PS_SCALAR_VALUE(job->args->data[5], U8); // loop over all source @@ -104,13 +184,4 @@ psphotMaskSource (source, add, markVal); } - if (add) { - psLogMsg ("psphot.noise", PS_LOG_WARN, "add noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.noise")); - } else { - psLogMsg ("psphot.noise", PS_LOG_WARN, "sub noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.noise")); - } - fprintf (stderr, "masked %d objects\n", Nmasked); - - psphotVisualShowImage (readout); - return true; } @@ -142,7 +213,5 @@ } } - Nmasked ++; - - return true; -} - + return true; +} + Index: branches/eam_branches/ipp-20130509/psphot/src/psphotBlendFit.c =================================================================== --- branches/eam_branches/ipp-20130509/psphot/src/psphotBlendFit.c (revision 35634) +++ branches/eam_branches/ipp-20130509/psphot/src/psphotBlendFit.c (revision 35640) @@ -86,4 +86,16 @@ if (!status || !isfinite(fitMaxTol) || fitMaxTol <= 0) { fitMaxTol = 1.0; + } + + bool chisqConvergence = psMetadataLookupBool (&status, recipe, "LMM_FIT_CHISQ_CONVERGENCE"); // Fit tolerance + if (!status) { + // default to the old method (chisqConvergence) + chisqConvergence = true; + } + + int gainFactorMode = psMetadataLookupS32 (&status, recipe, "LMM_FIT_GAIN_FACTOR_MODE"); // Fit tolerance + if (!status) { + // default to the old method (chisqConvergence) + gainFactorMode = 0; } @@ -119,4 +131,7 @@ fitOptions->mode = PM_SOURCE_FIT_PSF; fitOptions->covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix + + fitOptions->gainFactorMode = gainFactorMode; + fitOptions->chisqConvergence = chisqConvergence; psphotInitLimitsPSF (recipe, readout); @@ -211,4 +226,5 @@ psLogMsg ("psphot.psphotBlendFit", PS_LOG_WARN, "fit models: %f sec for %d objects (%d psf, %d ext, %d failed, %ld skipped)\n", psTimerMark ("psphot.fit.nonlinear"), Nfit, Npsf, Next, Nfail, sources->n - Nfit); + psphotFitSummary (); psphotVisualShowResidualImage (readout, false); @@ -271,32 +287,32 @@ // 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; + if (source->mode & PM_SOURCE_MODE_BLEND) goto skip_blend; + if (source->mode & PM_SOURCE_MODE_CR_LIMIT) goto skip_cr; + if (source->type == PM_SOURCE_TYPE_DEFECT) goto skip_defect; + if (source->type == PM_SOURCE_TYPE_SATURATED) goto skip_sat; // skip saturated stars modeled with a radial profile - if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue; + if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) goto skip_sat; // skip DBL second sources (ie, added by psphotFitBlob - if (source->mode & PM_SOURCE_MODE_PAIR) continue; + if (source->mode & PM_SOURCE_MODE_PAIR) goto skip_blend; // do not include MOMENTS_FAILURES in the fit - if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue; + if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) goto skip_generic; // limit selection to some SN limit if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { - if (source->moments->KronFlux < FIT_SN_LIM * source->moments->KronFluxErr) continue; + if (source->moments->KronFlux < FIT_SN_LIM * source->moments->KronFluxErr) goto skip_generic; } else { - if (sqrt(source->peak->detValue) < FIT_SN_LIM) continue; + if (sqrt(source->peak->detValue) < FIT_SN_LIM) goto skip_generic; } // 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 (source->peak->xf < AnalysisRegion.x0) goto skip_generic; + if (source->peak->yf < AnalysisRegion.y0) goto skip_generic; + if (source->peak->xf > AnalysisRegion.x1) goto skip_generic; + if (source->peak->yf > AnalysisRegion.y1) goto skip_generic; // if model is NULL, we don't have a starting guess - if (source->modelPSF == NULL) continue; + if (source->modelPSF == NULL) goto skip_generic; // skip sources which are insignificant flux? @@ -307,5 +323,5 @@ source->modelPSF->params->data.F32[2], source->modelPSF->params->data.F32[3]); - continue; + goto skip_generic; } @@ -357,4 +373,21 @@ pmSourceCacheModel (source, maskVal); pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + continue; + + skip_blend: + psTrace ("psphot", 5, "source at %7.1f, %7.1f skipped", source->peak->xf, source->peak->yf); + continue; + skip_cr: + psTrace ("psphot", 5, "source at %7.1f, %7.1f skipped", source->peak->xf, source->peak->yf); + continue; + skip_defect: + psTrace ("psphot", 5, "source at %7.1f, %7.1f skipped", source->peak->xf, source->peak->yf); + continue; + skip_sat: + psTrace ("psphot", 5, "source at %7.1f, %7.1f skipped", source->peak->xf, source->peak->yf); + continue; + skip_generic: + psTrace ("psphot", 5, "source at %7.1f, %7.1f skipped", source->peak->xf, source->peak->yf); + continue; } Index: branches/eam_branches/ipp-20130509/psphot/src/psphotChoosePSF.c =================================================================== --- branches/eam_branches/ipp-20130509/psphot/src/psphotChoosePSF.c (revision 35634) +++ branches/eam_branches/ipp-20130509/psphot/src/psphotChoosePSF.c (revision 35640) @@ -158,4 +158,16 @@ } float maxChisqDOF = psMetadataLookupF32 (&status, recipe, "PSF_FIT_MAX_CHISQ"); // Fit tolerance + + bool chisqConvergence = psMetadataLookupBool (&status, recipe, "LMM_FIT_CHISQ_CONVERGENCE"); // Fit tolerance + if (!status) { + // default to the old method (chisqConvergence) + chisqConvergence = true; + } + + int gainFactorMode = psMetadataLookupS32 (&status, recipe, "LMM_FIT_GAIN_FACTOR_MODE"); // Fit tolerance + if (!status) { + // default to the old method (chisqConvergence) + gainFactorMode = 0; + } // options which modify the behavior of the model fitting @@ -169,6 +181,8 @@ options->fitOptions->mode = PM_SOURCE_FIT_PSF; options->fitOptions->covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix - + options->fitOptions->gainFactorMode = gainFactorMode; + options->fitOptions->chisqConvergence = chisqConvergence; + psArray *stars = psArrayAllocEmpty (sources->n); Index: branches/eam_branches/ipp-20130509/psphot/src/psphotExtendedSourceFits.c =================================================================== --- branches/eam_branches/ipp-20130509/psphot/src/psphotExtendedSourceFits.c (revision 35634) +++ branches/eam_branches/ipp-20130509/psphot/src/psphotExtendedSourceFits.c (revision 35640) @@ -53,4 +53,6 @@ psTimerStart ("psphot.extended"); + psphotFitInitExtended(); + // find the currently selected readout pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest @@ -107,4 +109,28 @@ fitMaxTol = 1.0; } + + bool chisqConvergence = psMetadataLookupBool (&status, recipe, "LMM_FIT_CHISQ_CONVERGENCE"); // Fit tolerance + if (!status) { + // default to the old method (chisqConvergence) + chisqConvergence = true; + } + + int gainFactorMode = psMetadataLookupS32 (&status, recipe, "LMM_FIT_GAIN_FACTOR_MODE"); // Fit tolerance + if (!status) { + // default to the old method (chisqConvergence) + gainFactorMode = 0; + } + + // Define source fitting parameters for extended source fits + pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc(); + fitOptions->mode = PM_SOURCE_FIT_EXT; + fitOptions->saveCovariance = true; // XXX make this a user option? + fitOptions->covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix + fitOptions->nIter = fitIter; + fitOptions->minTol = fitMinTol; + fitOptions->maxTol = fitMaxTol; + + fitOptions->gainFactorMode = gainFactorMode; + fitOptions->chisqConvergence = chisqConvergence; // maskVal is used to test for rejected pixels, and must include markVal @@ -196,11 +222,8 @@ psMetadataIterator *iter = psMetadataIteratorAlloc (models, PS_LIST_HEAD, NULL); psArrayAdd(job->args, 1, iter); - psArrayAdd(job->args, 1, AnalysisRegion); // XXX make a pointer + psArrayAdd(job->args, 1, AnalysisRegion); + psArrayAdd(job->args, 1, fitOptions); PS_ARRAY_ADD_SCALAR(job->args, psfSize, PS_TYPE_S32); - PS_ARRAY_ADD_SCALAR(job->args, fitIter, PS_TYPE_S32); - PS_ARRAY_ADD_SCALAR(job->args, fitMinTol, PS_TYPE_F32); - PS_ARRAY_ADD_SCALAR(job->args, fitMaxTol, PS_TYPE_F32); - PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK); @@ -291,4 +314,7 @@ psLogMsg ("psphot", PS_LOG_INFO, " %d plain models (%d passed)\n", Nplain, NplainPass); psLogMsg ("psphot", PS_LOG_INFO, " %d too faint to fit, %d failed\n", Nfaint, Nfail); + + psphotFitSummaryExtended(); + return true; } @@ -314,19 +340,9 @@ psMetadataIterator *iter = job->args->data[3]; psRegion *region = job->args->data[4]; - int psfSize = PS_SCALAR_VALUE(job->args->data[5],S32); - int fitIter = PS_SCALAR_VALUE(job->args->data[6],S32); - float fitMinTol = PS_SCALAR_VALUE(job->args->data[7],F32); - float fitMaxTol = PS_SCALAR_VALUE(job->args->data[8],F32); - psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[9],PS_TYPE_IMAGE_MASK_DATA); - psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[10],PS_TYPE_IMAGE_MASK_DATA); - - // Define source fitting parameters for extended source fits - pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc(); - fitOptions->mode = PM_SOURCE_FIT_EXT; - fitOptions->saveCovariance = true; // XXX make this a user option? - fitOptions->covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix - fitOptions->nIter = fitIter; - fitOptions->minTol = fitMinTol; - fitOptions->maxTol = fitMaxTol; + pmSourceFitOptions *fitOptions = job->args->data[5]; + + int psfSize = PS_SCALAR_VALUE(job->args->data[6],S32); + psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[7],PS_TYPE_IMAGE_MASK_DATA); + psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[8],PS_TYPE_IMAGE_MASK_DATA); // psTraceSetLevel ("psLib.math.psMinimizeLMChi2_Alt", 5); @@ -540,26 +556,25 @@ psTrace ("psphot", 5, "extended source model for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My); } - psFree (fitOptions); // change the value of a scalar on the array (wrap this and put it in psArray.h) + scalar = job->args->data[9]; + scalar->data.S32 = Next; + + scalar = job->args->data[10]; + scalar->data.S32 = Nconvolve; + scalar = job->args->data[11]; - scalar->data.S32 = Next; + scalar->data.S32 = NconvolvePass; scalar = job->args->data[12]; - scalar->data.S32 = Nconvolve; + scalar->data.S32 = Nplain; scalar = job->args->data[13]; - scalar->data.S32 = NconvolvePass; + scalar->data.S32 = NplainPass; scalar = job->args->data[14]; - scalar->data.S32 = Nplain; + scalar->data.S32 = Nfaint; scalar = job->args->data[15]; - scalar->data.S32 = NplainPass; - - scalar = job->args->data[16]; - scalar->data.S32 = Nfaint; - - scalar = job->args->data[17]; scalar->data.S32 = Nfail; Index: branches/eam_branches/ipp-20130509/psphot/src/psphotModelBackground.c =================================================================== --- branches/eam_branches/ipp-20130509/psphot/src/psphotModelBackground.c (revision 35634) +++ branches/eam_branches/ipp-20130509/psphot/src/psphotModelBackground.c (revision 35640) @@ -28,4 +28,6 @@ } +// track background cell failures for log reporting purposes +static int nFailures = 0; // Generate the background model @@ -33,5 +35,5 @@ // linear interpolation to generate full-scale model // -// NOTE that the 'analysis' metedata pass in here is used to store the binning information. +// NOTE that the 'analysis' metedata passed in here is used to store the binning information. // This may be the analysis for this readout, but it may be the analysis for the pmFPAfile // corresponding to the model. Other information about the background model is saved on the @@ -109,38 +111,33 @@ const psStatsOptions statsOption = statsOptionLocation | statsOptionWidth; - psStats *stats = psStatsAlloc (statsOption); + psStats *statsDefaults = psStatsAlloc (statsOption); // set range for old-version of sky statistic if (statsOptionLocation & PS_STAT_ROBUST_QUARTILE) { - stats->min = 0.25; - stats->max = 0.75; + statsDefaults->min = 0.25; + statsDefaults->max = 0.75; } // set user-option for number of pixels per region - stats->nSubsample = psMetadataLookupF32 (&status, recipe, "IMSTATS_NPIX"); + statsDefaults->nSubsample = psMetadataLookupF32 (&status, recipe, "IMSTATS_NPIX"); if (!status) { - stats->nSubsample = 1000; + statsDefaults->nSubsample = 1000; } // optionally set the binsize - stats->binsize = psMetadataLookupF32 (&status, recipe, "SKY_HISTOGRAM_BINS"); + statsDefaults->binsize = psMetadataLookupF32 (&status, recipe, "SKY_HISTOGRAM_BINS"); if (status) { - stats->options |= PS_STAT_USE_BINSIZE; + statsDefaults->options |= PS_STAT_USE_BINSIZE; } // optionally set the sigma clipping - stats->clipSigma = psMetadataLookupF32 (&status, recipe, "SKY_CLIP_SIGMA"); + statsDefaults->clipSigma = psMetadataLookupF32 (&status, recipe, "SKY_CLIP_SIGMA"); if (!status) { - if (stats->options & PS_STAT_FITTED_MEAN) { - stats->clipSigma = 1.0; + if (statsDefaults->options & PS_STAT_FITTED_MEAN) { + statsDefaults->clipSigma = 1.0; } else { - stats->clipSigma = 3.0; + statsDefaults->clipSigma = 3.0; } } - - // stats is not initialized by psStats??? use this to save the input options - // XXX re-work this to use the new psStatsInit function - psStats *statsDefaults = psStatsAlloc (statsOption); - *statsDefaults = *stats; // we save the binning structure for use in psphotMagnitudes @@ -154,7 +151,6 @@ // XXXX we can thread this here by running blocks in parallel - - int nFailures = 0; psImageBackgroundInit(); + nFailures = 0; // reset for this pass // we have Nx * Ny model points, but we can use a window which is larger (or smaller) than @@ -163,95 +159,71 @@ // measure clipped median for subimages - psRegion ruffRegion = psRegionSet (0,0,0,0); - psRegion fineRegion = psRegionSet (0,0,0,0); for (int iy = 0; iy < model->numRows; iy++) { for (int ix = 0; ix < model->numCols; ix++) { - - // convert the ruff grid cell to the equivalent fine grid cell - ruffRegion = psRegionSet (ix + 0.5 - 0.5*dXsample, ix + 0.5 + 0.5*dXsample, - iy + 0.5 - 0.5*dYsample, iy + 0.5 + 0.5*dYsample); - fineRegion = psImageBinningSetFineRegion (binning, ruffRegion); - fineRegion = psRegionForImage (image, fineRegion); - - psImage *subset = psImageSubset (image, fineRegion); - if (!subset->numCols || !subset->numRows) { - psFree (subset); - continue; + + // allocate a job -- if threads are not defined, this just runs the job + psThreadJob *job = psThreadJobAlloc ("PSPHOT_MODEL_BACKGROUND"); + + psArrayAdd(job->args, 1, image); + psArrayAdd(job->args, 1, mask); + psArrayAdd(job->args, 1, binning); + psArrayAdd(job->args, 1, rng); + psArrayAdd(job->args, 1, statsDefaults); + + psArrayAdd(job->args, 1, modelData); + psArrayAdd(job->args, 1, modelStdevData); + + PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); + + PS_ARRAY_ADD_SCALAR(job->args, ix, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, iy, PS_TYPE_S32); + + PS_ARRAY_ADD_SCALAR(job->args, dXsample, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR(job->args, dYsample, PS_TYPE_F32); + + PS_ARRAY_ADD_SCALAR(job->args, statsOptionLocation, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, statsOptionWidth, PS_TYPE_S32); + + PS_ARRAY_ADD_SCALAR(job->args, NAN, PS_TYPE_F32); // this is used as a return value for dQvalue + +# if (1) + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + return NULL; } - psImage *submask = psImageSubset (mask, fineRegion); - - // reset the default values - *stats = *statsDefaults; - - // Use the selected background statistic for the first pass - // If it fails, fall back on the "ROBUST_MEDIAN" version - // If both fail, set the pixel to NAN and (below) interpolate - // XXX psImageBackground will probably be renamed psImageStats - // XXX don't bother trying if there are no valid pixels... - - psVector *sample = NULL; - - // turn on stats tracing in desired cells - # if (0) - psMetadata *plots = psMetadataLookupPtr (&status, recipe, "DIAGNOSTIC.PLOTS"); - assert (plots); - - int xPlot = psMetadataLookupS32 (&status, plots, "IMAGE.BACKGROUND.CELL.HISTOGRAM.X"); - assert (status); - int yPlot = psMetadataLookupS32 (&status, plots, "IMAGE.BACKGROUND.CELL.HISTOGRAM.Y"); - assert (status); - - bool gotX = (xPlot < 0) || (xPlot == ix); - bool gotY = (yPlot < 0) || (yPlot == iy); - - if (gotX && gotY) { - (void) psTraceSetLevel ("psLib.math.vectorFittedStats", 6); - (void) psTraceSetLevel ("psLib.math.vectorRobustStats", 6); - } else { - (void) psTraceSetLevel ("psLib.math.vectorFittedStats", 0); - (void) psTraceSetLevel ("psLib.math.vectorRobustStats", 0); +# else + if (!psphotModelBackground_Threaded(job)) { + psError(PS_ERR_UNKNOWN, false, "Failure to model background."); + return NULL; } - # endif - - if (psImageBackground(stats, &sample, subset, submask, maskVal, rng)) { - if (stats->options & PS_STAT_ROBUST_QUARTILE) { - modelData[iy][ix] = stats->robustMedian; - } else { - modelData[iy][ix] = psStatsGetValue(stats, statsOptionLocation); - } - modelStdevData[iy][ix] = psStatsGetValue(stats, statsOptionWidth); - // fprintf (stderr, "dQ : %f - %f - %f = %f\n", stats->robustLQ, stats->robustMedian, stats->robustUQ, stats->robustUQ + stats->robustLQ - 2*stats->robustMedian); - psVectorAppend (dQ, stats->robustUQ + stats->robustLQ - 2*stats->robustMedian); - - // supply sample to plotting routing - psphotDiagnosticPlots (config, "IMAGE.BACKGROUND.CELL.HISTOGRAM", ix, iy, modelData[iy][ix], modelStdevData[iy][ix], sample); - } else { - psStatsOptions currentOptions = stats->options; - stats->options = PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV; - if (!psImageBackground(stats, &sample, subset, submask, maskVal, rng)) { - if ((nFailures < 3) || (nFailures % 100 == 0)) { - psLogMsg ("psphot", PS_LOG_WARN, "Failed to estimate background using ROBUST_MEDIAN for " - "(%dx%d, (row0,col0) = (%d,%d)", - subset->numRows, subset->numCols, subset->row0, subset->col0); - } - nFailures ++; - modelData[iy][ix] = modelStdevData[iy][ix] = NAN; - } else { - modelData[iy][ix] = psStatsGetValue (stats, PS_STAT_ROBUST_MEDIAN); - modelStdevData[iy][ix] = psStatsGetValue(stats, PS_STAT_ROBUST_STDEV); - - // supply sample to plotting routing - psphotDiagnosticPlots (config, "IMAGE.BACKGROUND.CELL.HISTOGRAM", ix, iy, modelData[iy][ix], modelStdevData[iy][ix], sample); - } - // drop errors caused by psImageBackground failures - // XXX we probably should trap and exit on serious failures - psErrorClear(); - stats->options = currentOptions; - } - psFree(sample); - modelData[iy][ix] += SKY_BIAS; - psFree (subset); - psFree (submask); + if (job->args->n < 1) { + fprintf (stderr, "error with job\n"); + } else { + psScalar *scalar = job->args->data[14]; + float dQvalue = scalar->data.F32; + psVectorAppend (dQ, dQvalue); + } + psFree(job); +# endif } + } + + // wait for the threads to finish and manage results + if (!psThreadPoolWait (false, true)) { + 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 = job->args->data[14]; + float dQvalue = scalar->data.F32; + psVectorAppend (dQ, dQvalue); + } + psFree(job); } @@ -302,5 +274,4 @@ if (Count == 0) { psError (PSPHOT_ERR_DATA, true, "failed to build background image"); - psFree(stats); psFree(statsDefaults); psFree(binning); @@ -318,4 +289,11 @@ modelData[iy][ix] = Value; modelStdevData[iy][ix] = ValueStdev; + } + } + + // apply artificial offset if desired + for (int iy = 0; iy < model->numRows; iy++) { + for (int ix = 0; ix < model->numCols; ix++) { + modelData[iy][ix] += SKY_BIAS; } } @@ -364,5 +342,4 @@ psFree(dQ); - psFree(stats); psFree(statsDefaults); psFree(binning); @@ -433,2 +410,127 @@ return true; } + +bool psphotModelBackground_Threaded (psThreadJob *job) { + + psImage *image = job->args->data[0]; + psImage *mask = job->args->data[1]; + + psImageBinning *binning = job->args->data[2]; + + psRandom *rng = job->args->data[3]; + psStats *statsDefaults = job->args->data[4]; + + psF32 **modelData = job->args->data[5]; + psF32 **modelStdevData = job->args->data[6]; + + psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[7],PS_TYPE_IMAGE_MASK_DATA); + + int ix = PS_SCALAR_VALUE(job->args->data[8],S32); + int iy = PS_SCALAR_VALUE(job->args->data[9],S32); + + float dXsample = PS_SCALAR_VALUE(job->args->data[10],F32); + float dYsample = PS_SCALAR_VALUE(job->args->data[11],F32); + + psStatsOptions statsOptionLocation = PS_SCALAR_VALUE(job->args->data[12],S32); + psStatsOptions statsOptionWidth = PS_SCALAR_VALUE(job->args->data[13],S32); + + // convert the ruff grid cell to the equivalent fine grid cell + psRegion ruffRegion = psRegionSet (ix + 0.5 - 0.5*dXsample, ix + 0.5 + 0.5*dXsample, iy + 0.5 - 0.5*dYsample, iy + 0.5 + 0.5*dYsample); + psRegion fineRegion = psImageBinningSetFineRegion (binning, ruffRegion); + fineRegion = psRegionForImage (image, fineRegion); + + psImage *subset = psImageSubset (image, fineRegion); + if (!subset->numCols || !subset->numRows) { + psFree (subset); + return false; // XXX do we / should we fail on this? + } + psImage *submask = psImageSubset (mask, fineRegion); + + psStats *stats = psStatsAlloc (PS_STAT_NONE); + *stats = *statsDefaults; + + // Use the selected background statistic for the first pass + // If it fails, fall back on the "ROBUST_MEDIAN" version + // If both fail, set the pixel to NAN and (later) interpolate + + psVector *sample = NULL; + float dQvalue = NAN; + + if (psImageBackground(stats, &sample, subset, submask, maskVal, rng)) { + if (stats->options & PS_STAT_ROBUST_QUARTILE) { + modelData[iy][ix] = stats->robustMedian; + } else { + modelData[iy][ix] = psStatsGetValue(stats, statsOptionLocation); + } + modelStdevData[iy][ix] = psStatsGetValue(stats, statsOptionWidth); + + // fprintf (stderr, "dQ : %f - %f - %f = %f\n", stats->robustLQ, stats->robustMedian, stats->robustUQ, stats->robustUQ + stats->robustLQ - 2*stats->robustMedian); + // XXX this operation is not thread safe -- move out somehow... + // psVectorAppend (dQ, stats->robustUQ + stats->robustLQ - 2*stats->robustMedian); + dQvalue = stats->robustUQ + stats->robustLQ - 2*stats->robustMedian; + // return dQvalue to main thread + + // supply sample to plotting routing + // only allow in a non-threaded context -- can plot more than one cell + // psphotDiagnosticPlots (config, "IMAGE.BACKGROUND.CELL.HISTOGRAM", ix, iy, modelData[iy][ix], modelStdevData[iy][ix], sample); + } else { + // psStatsOptions currentOptions = stats->options; + stats->options = PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV; + if (!psImageBackground(stats, &sample, subset, submask, maskVal, rng)) { + if ((nFailures < 3) || (nFailures % 100 == 0)) { + psLogMsg ("psphot", PS_LOG_WARN, "Failed to estimate background using ROBUST_MEDIAN for " + "(%dx%d, (row0,col0) = (%d,%d)", + subset->numRows, subset->numCols, subset->row0, subset->col0); + } + nFailures ++; // static + modelData[iy][ix] = modelStdevData[iy][ix] = NAN; + } else { + modelData[iy][ix] = psStatsGetValue (stats, PS_STAT_ROBUST_MEDIAN); + modelStdevData[iy][ix] = psStatsGetValue(stats, PS_STAT_ROBUST_STDEV); + + // supply sample to plotting routing + // only allow in a non-threaded context -- can plot more than one cell + // psphotDiagnosticPlots (config, "IMAGE.BACKGROUND.CELL.HISTOGRAM", ix, iy, modelData[iy][ix], modelStdevData[iy][ix], sample); + } + // drop errors caused by psImageBackground failures + // NOTE : psStats raises errors when it cannot find a solution; this + // probably should not be errors but exit stats info only, but c'est la code + // XXX we probably should trap and exit on serious failures + psErrorClear(); + // stats->options = currentOptions; // this is not needed (set by init above) + } + psFree (stats); + psFree (sample); + psFree (subset); + psFree (submask); + + // return the dQvalue to the calling thread + psScalar *scalar = job->args->data[14]; + scalar->data.F32 = dQvalue; + + return true; +} + +// code for in-line plotting from above +// turn on stats tracing in desired cells +// XXX this code may need to be re-worked for a threaded context +# if (0) + psMetadata *plots = psMetadataLookupPtr (&status, recipe, "DIAGNOSTIC.PLOTS"); + assert (plots); + + int xPlot = psMetadataLookupS32 (&status, plots, "IMAGE.BACKGROUND.CELL.HISTOGRAM.X"); + assert (status); + int yPlot = psMetadataLookupS32 (&status, plots, "IMAGE.BACKGROUND.CELL.HISTOGRAM.Y"); + assert (status); + + bool gotX = (xPlot < 0) || (xPlot == ix); + bool gotY = (yPlot < 0) || (yPlot == iy); + + if (gotX && gotY) { + (void) psTraceSetLevel ("psLib.math.vectorFittedStats", 6); + (void) psTraceSetLevel ("psLib.math.vectorRobustStats", 6); + } else { + (void) psTraceSetLevel ("psLib.math.vectorFittedStats", 0); + (void) psTraceSetLevel ("psLib.math.vectorRobustStats", 0); + } +# endif Index: branches/eam_branches/ipp-20130509/psphot/src/psphotSetThreads.c =================================================================== --- branches/eam_branches/ipp-20130509/psphot/src/psphotSetThreads.c (revision 35634) +++ branches/eam_branches/ipp-20130509/psphot/src/psphotSetThreads.c (revision 35640) @@ -5,6 +5,18 @@ psThreadTask *task = NULL; + pmPSFThreads (); + + task = psThreadTaskAlloc("PSPHOT_MODEL_BACKGROUND", 15); + task->function = &psphotModelBackground_Threaded; + psThreadTaskAdd(task); + psFree(task); + task = psThreadTaskAlloc("PSPHOT_GUESS_MODEL", 5); task->function = &psphotGuessModel_Threaded; + psThreadTaskAdd(task); + psFree(task); + + task = psThreadTaskAlloc("PSPHOT_ADD_NOISE", 6); + task->function = &psphotAddOrSubNoise_Threaded; psThreadTaskAdd(task); psFree(task); @@ -40,5 +52,5 @@ psFree(task); - task = psThreadTaskAlloc("PSPHOT_EXTENDED_FIT", 18); + task = psThreadTaskAlloc("PSPHOT_EXTENDED_FIT", 16); task->function = &psphotExtendedSourceFits_Threaded; psThreadTaskAdd(task); Index: branches/eam_branches/ipp-20130509/psphot/src/psphotSourceFits.c =================================================================== --- branches/eam_branches/ipp-20130509/psphot/src/psphotSourceFits.c (revision 35634) +++ branches/eam_branches/ipp-20130509/psphot/src/psphotSourceFits.c (revision 35640) @@ -3,9 +3,20 @@ // given a source with an existing modelPSF, attempt a full PSF fit, subtract if successful +static int NfitBlend = 0; + static int NfitPSF = 0; -static int NfitBlend = 0; +static int NfitIterPSF = 0; + static int NfitDBL = 0; +static int NfitIterDBL = 0; +static int NfitPixDBL = 0; + static int NfitEXT = 0; +static int NfitIterEXT = 0; +static int NfitPixEXT = 0; + static int NfitPCM = 0; +static int NfitIterPCM = 0; +static int NfitPixPCM = 0; bool psphotFitInit (int nThreads) { @@ -17,6 +28,23 @@ bool psphotFitSummary () { - psLogMsg ("psphot.pspsf", PS_LOG_INFO, "fitted %5d psf, %5d blend, %5d ext, %5d dbl : %6.2f sec\n", - NfitPSF, NfitBlend, NfitEXT, NfitDBL, psTimerMark ("psphot.fits")); + psLogMsg ("psphot.pspsf", PS_LOG_INFO, "fitted %5d psf (%d iter), %5d blend: %5d ext (%d iter, %d pix), %5d dbl (%d iter, %d pix) : %6.2f sec\n", + NfitPSF, NfitIterPSF, NfitBlend, NfitEXT, NfitIterEXT, NfitPixEXT, NfitDBL, NfitIterDBL, NfitPixDBL, psTimerMark ("psphot.fits")); + return true; +} + +bool psphotFitSummaryExtended () { + + psLogMsg ("psphot.pspsf", PS_LOG_INFO, "fitted %5d ext (%d iter, %d pix), %5d pcm (%d iter, %d pix)\n", + NfitEXT, NfitIterEXT, NfitPixEXT, NfitPCM, NfitIterPCM, NfitPixPCM); + return true; +} + +bool psphotFitInitExtended () { + NfitEXT = 0; + NfitIterEXT = 0; + NfitPixEXT = 0; + NfitPCM = 0; + NfitIterPCM = 0; + NfitPixPCM = 0; return true; } @@ -169,4 +197,5 @@ options.mode = PM_SOURCE_FIT_PSF; pmSourceFitModel (source, PSF, &options, maskVal); + NfitIterPSF += PSF->nIter; if (!isfinite(PSF->params->data.F32[PM_PAR_I0])) psAbort("nan in fit"); @@ -440,4 +469,7 @@ options.mode = PM_SOURCE_FIT_PSF; pmSourceFitSet (source, modelSet, &options, maskVal); + NfitIterDBL += DBL->nIter; + NfitPixDBL += DBL->nDOF; + return (modelSet); } @@ -507,4 +539,6 @@ pmSourceFitModel (source, model, &options, maskVal); + NfitIterEXT += model->nIter; + NfitPixEXT += model->nDOF; # if (PS_TRACE_ON) @@ -587,4 +621,6 @@ // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5); pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize); // NOTE : 1687 allocs in here + NfitIterPCM += pcm->modelConv->nIter; + NfitPixPCM += pcm->modelConv->nDOF; if (TIMING) { t5 = psTimerMark ("psphotFitPCM"); }