Index: trunk/psphot/src/psphotFitSourcesLinearStack.c =================================================================== --- trunk/psphot/src/psphotFitSourcesLinearStack.c (revision 27532) +++ trunk/psphot/src/psphotFitSourcesLinearStack.c (revision 27547) @@ -1,14 +1,23 @@ # include "psphotInternal.h" -// fit flux (and optionally sky model) to all reasonable sources -// with the linear fitting process. sources must have an associated -// model with selected pixels, and the fit radius must be defined - -// given the set of sources, each of which points to the pixels in the -// science image, we construct a set of simulated sources with their own pixels. -// these are used to determine the simultaneous linear fit of fluxes. -// the analysis is performed wrt the simulated pixel values - -bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final) { +// for now, let's store the detections on the readout->analysis for each readout +bool psphotFitSourcesLinearStack (pmConfig *config, const pmFPAview *view, bool final) +{ + bool status = true; + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotFitSourcesLinearReadoutStack (config, view, "PSPHOT.INPUT", i, final)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +} + +bool psphotFitSourcesLinearReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, bool final) { bool status; @@ -18,4 +27,29 @@ // float r; + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + assert (recipe); + + // 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->allSources; + psAssert (sources, "missing sources?"); + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping linear fit"); + return true; + } + + pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); + psAssert (sources, "missing psf?"); + psTimerStart ("psphot.linear"); @@ -36,4 +70,10 @@ // storage array for fitSources psArray *fitSources = psArrayAllocEmpty (sources->n); + + // option to limit analysis to a specific region + char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION"); + psRegion AnalysisRegion = psRegionFromString (region); + AnalysisRegion = psRegionForImage (readout->image, AnalysisRegion); + if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined"); bool CONSTANT_PHOTOMETRIC_WEIGHTS = @@ -42,32 +82,44 @@ psAbort("You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS"); } + int SKY_FIT_ORDER = psMetadataLookupS32(&status, recipe, "SKY_FIT_ORDER"); + if (!status) { + SKY_FIT_ORDER = 0; + } + bool SKY_FIT_LINEAR = psMetadataLookupBool(&status, recipe, "SKY_FIT_LINEAR"); + if (!status) { + SKY_FIT_LINEAR = false; + } + + // XXX test: choose a larger-than expected radius: + float covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix + psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "covariance factor: %f\n", covarFactor); + + // XXX do not apply covarFactor for the moment... + // covarFactor = 1.0; // select the sources which will be used for the fitting analysis - for (int i = 0; i < sources->n; i++) { - 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; - - // do not include CRs in the full ensemble fit - if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue; - - if (final) { - if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) continue; - } else { - if (source->mode & PM_SOURCE_MODE_BLEND) continue; - } - - // generate model for sources without, or skip if we can't - if (!source->modelFlux) { + for (int i = 0; i < objects->n; i++) { + pmPhotObj *object = objects->data[i]; + if (!object) continue; + if (!object->sources) continue; + + // check an element of the group to see if we should use it + if (!object->flags & PM_PHOT_OBJ_BAD) continue; + + for (int j = 0; j < object->sources->n; j++) { + pmSource *source = object->sources->data[j]; + if (!source) continue; + + // turn this bit off and turn it on again if we keep this source + source->mode &= ~PM_SOURCE_MODE_LINEAR_FIT; + + // generate model for sources without, or skip if we can't + if (!source->modelFlux) { if (!pmSourceCacheModel (source, maskVal)) continue; - } - - source->mode |= PM_SOURCE_MODE_LINEAR_FIT; - psArrayAdd (fitSources, 100, source); + } + + source->mode |= PM_SOURCE_MODE_LINEAR_FIT; + psArrayAdd (fitSources, 100, source); + } } psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "built fitSources: %f sec (%ld objects)\n", psTimerMark ("psphot.linear"), sources->n); @@ -81,8 +133,6 @@ psVector *errors = psVectorAlloc (fitSources->n, PS_TYPE_F32); - // create the border matrix (includes the sparse matrix) - // for just sky: 1 row; for x,y terms: 3 rows + // create the sparse matrix psSparse *sparse = psSparseAlloc (fitSources->n, 100); - psSparseBorder *border = psSparseBorderAlloc (sparse, 1); // fill out the sparse matrix elements and border elements (B) @@ -93,10 +143,10 @@ // diagonal elements of the sparse matrix (auto-cross-product) - f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS); + f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor); psSparseMatrixElement (sparse, i, i, f); // the formal error depends on the weighting scheme if (CONSTANT_PHOTOMETRIC_WEIGHTS) { - float var = pmSourceModelDotModel (SRCi, SRCi, false); + float var = pmSourceModelDotModel (SRCi, SRCi, false, covarFactor); errors->data.F32[i] = 1.0 / sqrt(var); } else { @@ -104,15 +154,14 @@ } - // find the image x model value - f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS); + f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor); psSparseVectorElement (sparse, i, f); - - f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS); - psSparseBorderElementB (border, i, 0, f); // loop over all other stars following this one for (int j = i + 1; j < fitSources->n; j++) { pmSource *SRCj = fitSources->data[j]; + + // XXX I need to know if this source is on the same image as SRCi -- + if (!sameImge) { continue; } // skip over disjoint source images, break after last possible overlap @@ -123,5 +172,5 @@ // got an overlap; calculate cross-product and add to output array - f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS); + f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor); psSparseMatrixElement (sparse, j, i, f); } @@ -130,7 +179,4 @@ psSparseResort (sparse); psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "built matrix: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem); - - // set the sky, sky_x, sky_y components of border matrix - SetBorderMatrixElements (border, readout, fitSources, CONSTANT_PHOTOMETRIC_WEIGHTS, SKY_FIT_ORDER, markVal); psSparseConstraint constraint; @@ -140,6 +186,6 @@ // solve for normalization terms (need include local sky?) - psVector *norm = psSparseSolve (NULL, constraint, sparse, 5); - + psVector *norm = NULL; + norm = psSparseSolve (NULL, constraint, sparse, 5); psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "solve matrix: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem); @@ -170,9 +216,10 @@ // measure chisq for each source - for (int i = 0; final && (i < fitSources->n); i++) { + // for (int i = 0; final && (i < fitSources->n); i++) { + for (int i = 0; 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); + pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal, covarFactor); } psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "get chisqs: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem); @@ -182,4 +229,5 @@ psFree (fitSources); psFree (norm); + psFree (skyfit); psFree (errors); psFree (border); @@ -188,6 +236,93 @@ psphotVisualShowResidualImage (readout); - psphotVisualShowFlags (sources); + psphotVisualPlotChisq (sources); + // psphotVisualShowFlags (sources); + + // We have to place this visualization here because the models are not realized until + // psphotGuessModels or fitted until psphotFitSourcesLinear. + psphotVisualShowPSFStars (recipe, psf, sources); return true; } + +// XXX do we need this? +// XXX disallow the simultaneous sky fit and remove this code... + +// Calculate the weight terms for the sky fit component of the matrix. This function operates +// on the pixels which correspond to all of the sources of interest. These elements fill in +// the border matrix components in the sparse matrix equation. +static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER, psImageMaskType markVal) { + + // generate the image-wide weight terms + // turn on MARK for all image pixels + psRegion fullArray = psRegionSet (0, 0, 0, 0); + fullArray = psRegionForImage (readout->mask, fullArray); + psImageMaskRegion (readout->mask, fullArray, "OR", markVal); + + // turn off MARK for all object pixels + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + pmModel *model = pmSourceGetModel (NULL, source); + if (model == NULL) continue; + float x = model->params->data.F32[PM_PAR_XPOS]; + float y = model->params->data.F32[PM_PAR_YPOS]; + psImageMaskCircle (source->maskView, x, y, model->fitRadius, "AND", PS_NOT_IMAGE_MASK(markVal)); + } + + // accumulate the image statistics from the masked regions + psF32 **image = readout->image->data.F32; + psF32 **variance = readout->variance->data.F32; + psImageMaskType **mask = readout->mask->data.PS_TYPE_IMAGE_MASK_DATA; + + double w, x, y, x2, xy, y2, xc, yc, wt, f, fo, fx, fy; + w = x = y = x2 = xy = y2 = fo = fx = fy = 0; + + int col0 = readout->image->col0; + int row0 = readout->image->row0; + + for (int j = 0; j < readout->image->numRows; j++) { + for (int i = 0; i < readout->image->numCols; i++) { + if (mask[j][i]) continue; + if (constant_weights) { + wt = 1.0; + } else { + wt = variance[j][i]; + } + f = image[j][i]; + w += 1/wt; + fo += f/wt; + if (SKY_FIT_ORDER == 0) continue; + + xc = i + col0; + yc = j + row0; + x += xc/wt; + y += yc/wt; + x2 += xc*xc/wt; + xy += xc*yc/wt; + y2 += yc*yc/wt; + fx += f*xc/wt; + fy += f*yc/wt; + } + } + + // turn off MARK for all image pixels + psImageMaskRegion (readout->mask, fullArray, "AND", PS_NOT_IMAGE_MASK(markVal)); + + // set the Border T elements + psSparseBorderElementG (border, 0, fo); + psSparseBorderElementT (border, 0, 0, w); + if (SKY_FIT_ORDER > 0) { + psSparseBorderElementG (border, 0, fx); + psSparseBorderElementG (border, 0, fy); + psSparseBorderElementT (border, 1, 0, x); + psSparseBorderElementT (border, 2, 0, y); + psSparseBorderElementT (border, 0, 1, x); + psSparseBorderElementT (border, 1, 1, x2); + psSparseBorderElementT (border, 2, 1, xy); + psSparseBorderElementT (border, 0, 2, y); + psSparseBorderElementT (border, 1, 2, xy); + psSparseBorderElementT (border, 2, 2, y2); + } + + return true; +}