Index: branches/tap_branches/psphot/src/psphotFitSourcesLinearStack.c =================================================================== --- branches/tap_branches/psphot/src/psphotFitSourcesLinearStack.c (revision 25900) +++ branches/tap_branches/psphot/src/psphotFitSourcesLinearStack.c (revision 27838) @@ -1,22 +1,16 @@ # include "psphotInternal.h" +// XXX need array of covar factors for each image +// XXX define the 'good' / 'bad' flags? -// 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 +# define COVAR_FACTOR 1.0 -// 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 - -static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, psImageMaskType markVal); - -bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final) { +bool psphotFitSourcesLinearStack (pmConfig *config, psArray *objects, bool final) { bool status; - float x; - float y; float f; - // float r; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + assert (recipe); psTimerStart ("psphot.linear"); @@ -33,45 +27,44 @@ maskVal |= markVal; - // source analysis is done in spatial order - sources = psArraySort (sources, pmSourceSortByY); + // analysis is done in spatial order (to speed up overlap search) + // sort by first element in each source list + objects = psArraySort (objects, pmPhotObjSortByX); // storage array for fitSources - psArray *fitSources = psArrayAllocEmpty (sources->n); + psArray *fitSources = psArrayAllocEmpty (objects->n); - bool CONSTANT_PHOTOMETRIC_WEIGHTS = - psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS"); - if (!status) { - psAbort("You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS"); - } + bool CONSTANT_PHOTOMETRIC_WEIGHTS = psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS"); + psAssert (status, "You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS"); + + // XXX store a local static array of covar factors for each of the images (by image ID) + // float covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix + // psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "covariance factor: %f\n", covarFactor); // select the sources which will be used for the fitting analysis - for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; + for (int i = 0; i < objects->n; i++) { + pmPhotObj *object = objects->data[i]; + if (!object) continue; + if (!object->sources) continue; - // turn this bit off and turn it on again if we pass this test - source->mode &= ~PM_SOURCE_MODE_LINEAR_FIT; + // XXX check an element of the group to see if we should use it + // if (!object->flags & PM_PHOT_OBJ_BAD) continue; - // skip non-astronomical objects (very likely defects) - if (source->type == PM_SOURCE_TYPE_DEFECT) continue; - if (source->type == PM_SOURCE_TYPE_SATURATED) continue; + for (int j = 0; j < object->sources->n; j++) { + pmSource *source = object->sources->data[j]; + if (!source) continue; - // do not include CRs in the full ensemble fit - if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue; + // turn this bit off and turn it on again if we keep this source + source->mode &= ~PM_SOURCE_MODE_LINEAR_FIT; - 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) { + if (!pmSourceCacheModel (source, maskVal)) continue; + } - // 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); + psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "built fitSources: %f sec (%ld objects)\n", psTimerMark ("psphot.linear"), objects->n); if (fitSources->n == 0) { @@ -83,8 +76,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) @@ -95,10 +86,10 @@ // diagonal elements of the sparse matrix (auto-cross-product) - f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS); + f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR); 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, COVAR_FACTOR); errors->data.F32[i] = 1.0 / sqrt(var); } else { @@ -106,11 +97,7 @@ } - // find the image x model value - f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS); + f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR); psSparseVectorElement (sparse, i, f); - - f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS); - psSparseBorderElementB (border, i, 0, f); // loop over all other stars following this one @@ -118,12 +105,15 @@ pmSource *SRCj = fitSources->data[j]; + // we only need to generate dot terms for source on the same image + if (SRCj->imageID != SRCi->imageID) { continue; } + // skip over disjoint source images, break after last possible overlap - if (SRCi->pixels->row0 + SRCi->pixels->numRows < SRCj->pixels->row0) break; - if (SRCj->pixels->row0 + SRCj->pixels->numRows < SRCi->pixels->row0) continue; - if (SRCi->pixels->col0 + SRCi->pixels->numCols < SRCj->pixels->col0) continue; - if (SRCj->pixels->col0 + SRCj->pixels->numCols < SRCi->pixels->col0) continue; + if (SRCj->pixels->row0 + SRCj->pixels->numRows < SRCi->pixels->row0) continue; // source(i) is above source(j) + if (SRCi->pixels->row0 + SRCi->pixels->numRows < SRCj->pixels->row0) continue; // source(i) is below source(j) + if (SRCj->pixels->col0 + SRCj->pixels->numCols < SRCi->pixels->col0) continue; // source(i) is right of source(j) + if (SRCi->pixels->col0 + SRCi->pixels->numCols < SRCj->pixels->col0) break; // source(i) is left of source(j) [no other source(j) can overlap source(i)] // 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, COVAR_FACTOR); psSparseMatrixElement (sparse, j, i, f); } @@ -133,7 +123,4 @@ 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; constraint.paramMin = 0.0; @@ -142,13 +129,7 @@ // 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); - - // XXXX **** philosophical question: - // we measure bright objects in three passes: 1) linear fit; 2) non-linear fit; 3) linear fit: - // should retain the chisq and errors from the intermediate non-linear fit? - // the non-linear fit provides better values for the position errors, and for - // extended sources, the shape errors // adjust I0 for fitSources and subtract @@ -164,5 +145,4 @@ model->params->data.F32[PM_PAR_I0] = norm->data.F32[i]; model->dparams->data.F32[PM_PAR_I0] = errors->data.F32[i]; - // XXX is the value of 'errors' modified by the sky fit? // subtract object @@ -172,9 +152,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, COVAR_FACTOR); } psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "get chisqs: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem); @@ -185,67 +166,22 @@ psFree (norm); psFree (errors); - psFree (border); psLogMsg ("psphot.ensemble", PS_LOG_INFO, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot.linear")); - - psphotVisualShowResidualImage (readout); - psphotVisualShowFlags (sources); return true; } -// 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, psImageMaskType markVal) { +// sort by X (ascending) +int pmPhotObjSortByX (const void **a, const void **b) +{ + pmPhotObj *objA = *(pmPhotObj **)a; + pmPhotObj *objB = *(pmPhotObj **)b; - // 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); + psF32 fA = objA->x; + psF32 fB = objB->x; - // 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->radiusFit, "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; - } - } - - // 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); - - return true; + psF32 diff = fA - fB; + if (diff > FLT_EPSILON) return (+1); + if (diff < FLT_EPSILON) return (-1); + return (0); }