Index: trunk/psphot/src/psphotFitSourcesLinearStack.c =================================================================== --- trunk/psphot/src/psphotFitSourcesLinearStack.c (revision 27547) +++ trunk/psphot/src/psphotFitSourcesLinearStack.c (revision 27657) @@ -1,54 +1,16 @@ # include "psphotInternal.h" +// XXX need array of covar factors for each image +// XXX define the 'good' / 'bad' flags? -// 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; +# define COVAR_FACTOR 1.0 - 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 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); - - // 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"); @@ -65,36 +27,17 @@ 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); - // 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 = psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS"); + psAssert (status, "You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS"); - 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"); - } - 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; + // 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 @@ -104,6 +47,6 @@ 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; + // XXX 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++) { @@ -123,5 +66,5 @@ } } - 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) { @@ -143,10 +86,10 @@ // diagonal elements of the sparse matrix (auto-cross-product) - f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor); + 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, covarFactor); + float var = pmSourceModelDotModel (SRCi, SRCi, false, COVAR_FACTOR); errors->data.F32[i] = 1.0 / sqrt(var); } else { @@ -155,5 +98,5 @@ // find the image x model value - f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor); + f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR); psSparseVectorElement (sparse, i, f); @@ -162,15 +105,15 @@ pmSource *SRCj = fitSources->data[j]; - // XXX I need to know if this source is on the same image as SRCi -- - if (!sameImge) { continue; } + // 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, covarFactor); + f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR); psSparseMatrixElement (sparse, j, i, f); } @@ -190,10 +133,4 @@ 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 for (int i = 0; i < fitSources->n; i++) { @@ -208,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 @@ -221,5 +157,5 @@ if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) continue; pmModel *model = pmSourceGetModel (NULL, source); - pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal, covarFactor); + 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); @@ -229,100 +165,23 @@ psFree (fitSources); psFree (norm); - psFree (skyfit); psFree (errors); - psFree (border); psLogMsg ("psphot.ensemble", PS_LOG_INFO, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot.linear")); - - psphotVisualShowResidualImage (readout); - 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... +// sort by X (ascending) +int pmPhotObjSortByX (const void **a, const void **b) +{ + pmPhotObj *objA = *(pmPhotObj **)a; + pmPhotObj *objB = *(pmPhotObj **)b; -// 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) { + psF32 fA = objA->x; + psF32 fB = objB->x; - // 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; + psF32 diff = fA - fB; + if (diff > FLT_EPSILON) return (+1); + if (diff < FLT_EPSILON) return (-1); + return (0); }