Index: /trunk/psphot/src/psphotEnsemblePSF.c =================================================================== --- /trunk/psphot/src/psphotEnsemblePSF.c (revision 10051) +++ /trunk/psphot/src/psphotEnsemblePSF.c (revision 10052) @@ -1,8 +1,19 @@ # include "psphot.h" -// fit all reasonable sources with the linear PSF model -bool psphotEnsemblePSF (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final) { - - bool status; +// 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 + +static bool SetBorderT (psSparseBorder *border, pmReadout *readout, psArray *refSources, bool fitSlope); + +bool psphotEnsemblePSF (pmReadout *readout, psArray *refSources, psMetadata *recipe, pmPSF *psf, bool final) { + + bool isPSF; + bool status; float x; float y; @@ -10,21 +21,12 @@ float r; - // psRegion allArray = psRegionSet (0, 0, 0, 0); - psTimerStart ("psphot"); // source analysis is done in spatial order - sources = psArraySort (sources, psphotSortByY); - - // radius for source moments if star is saturated - // XXX we already do this in psphotSourceStats, is it still needed? - float INNER_RADIUS = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS"); - - // setup the PSF fit details - psphotInitRadiusPSF (recipe, psf->type); - - // pre-calculate all model pixels - psArray *models = psArrayAllocEmpty (sources->n); - psVector *index = psVectorAllocEmpty (sources->n, PS_TYPE_U32); + refSources = psArraySort (refSources, psphotSortByY); + + // storage arrays for fitSources and sequence index + psArray *fitSources = psArrayAllocEmpty (refSources->n); + psVector *index = psVectorAllocEmpty (refSources->n, PS_TYPE_U32); // option to limit analysis to a specific region @@ -32,5 +34,4 @@ psRegion AnalysisRegion = psRegionFromString (region); AnalysisRegion = psRegionForImage (readout->image, AnalysisRegion); - // psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region)); if (psRegionIsNaN (AnalysisRegion)) psAbort ("psphot", "analysis region mis-defined"); @@ -41,61 +42,20 @@ } - for (int i = 0; i < sources->n; i++) { - pmSource *inSource = sources->data[i]; + // construct source table with stand-alone fitSources + for (int i = 0; i < refSources->n; i++) { + pmSource *refSource = refSources->data[i]; // skip non-astronomical objects (very likely defects) - // XXX EAM : should we try these anyway? - if (inSource->type == PM_SOURCE_TYPE_DEFECT) continue; - if (inSource->type == PM_SOURCE_TYPE_SATURATED) continue; + if (refSource->type == PM_SOURCE_TYPE_DEFECT) continue; + if (refSource->type == PM_SOURCE_TYPE_SATURATED) continue; if (final) { - if (inSource->mode & PM_SOURCE_MODE_SUBTRACTED) continue; + if (refSource->mode & PM_SOURCE_MODE_SUBTRACTED) continue; } else { - if (inSource->mode & PM_SOURCE_MODE_BLEND) continue; + if (refSource->mode & PM_SOURCE_MODE_BLEND) continue; } - if (inSource->moments->x < AnalysisRegion.x0) continue; - if (inSource->moments->y < AnalysisRegion.y0) continue; - if (inSource->moments->x > AnalysisRegion.x1) continue; - if (inSource->moments->y > AnalysisRegion.y1) continue; - - pmSource *otSource = pmSourceAlloc (); - - // really saturated stars should be re-measured for a better centroid - // XXX EAM : move this to a 'clear satstar function' - // XXX : place this in pmSourceMoments or wherever it is called? - if (inSource->mode & PM_SOURCE_MODE_SATSTAR) { - status = pmSourceMoments (inSource, INNER_RADIUS); - } - - // XXX EAM : add option to use EXT or PSF form - // use the source moments, etc to guess basic model parameters - pmModel *modelEXT = pmSourceModelGuess (inSource, psf->type); - if (inSource->mode & PM_SOURCE_MODE_SATSTAR) { - modelEXT->params->data.F32[PM_PAR_XPOS] = inSource->moments->x; - modelEXT->params->data.F32[PM_PAR_YPOS] = inSource->moments->y; - } else { - // peak-up on peak (for non-sat objects) - - // ix,iy must land on inSource->pixels - int ix = PS_MAX (inSource->pixels->col0 + 1, PS_MIN (inSource->pixels->col0 + inSource->pixels->numCols - 2, inSource->peak->x)); - int iy = PS_MAX (inSource->pixels->row0 + 1, PS_MIN (inSource->pixels->row0 + inSource->pixels->numRows - 2, inSource->peak->y)); - - psPolynomial2D *bicube = psImageBicubeFit (inSource->pixels, ix, iy); - psPlane min = psImageBicubeMin (bicube); - - psTrace ("psphot", 5, "peak coord: %f %f -> %f %f\n", - modelEXT->params->data.F32[PM_PAR_XPOS], modelEXT->params->data.F32[PM_PAR_YPOS], min.x + ix, min.y + iy); - - // if min point is too deviant, keep the old value - if ((fabs(min.x) < 1.5) && (fabs(min.y) < 1.5)) { - modelEXT->params->data.F32[PM_PAR_XPOS] = min.x + ix; - modelEXT->params->data.F32[PM_PAR_YPOS] = min.y + iy; - } - psFree (bicube); - } - - // set PSF parameters for this model - pmModel *model = pmModelFromPSF (modelEXT, psf); - psFree (modelEXT); + // which model to use? + pmModel *model = pmSourceGetModel (&isPSF, refSource); + if (model == NULL) continue; // model must be defined // save the original coords @@ -103,59 +63,64 @@ y = model->params->data.F32[PM_PAR_YPOS]; - // set the fit radius based on the object flux limit and the model - psphotCheckRadiusPSF (readout, inSource, model); + // is the source in the region of interest? + if (x < AnalysisRegion.x0) continue; + if (y < AnalysisRegion.y0) continue; + if (x > AnalysisRegion.x1) continue; + if (y > AnalysisRegion.y1) continue; + + pmSource *fitSource = pmSourceAlloc (); // make temporary copies of the image pixels and mask // we need to have a copy which will not be modified by changes to its neighbor - otSource->mask = psImageCopy (NULL, inSource->mask, PS_TYPE_U8); - otSource->pixels = psImageCopy (NULL, inSource->pixels, PS_TYPE_F32); - otSource->weight = psImageCopy (NULL, inSource->weight, PS_TYPE_F32); - - // build the model image - psImage *flux = otSource->pixels; - psImage *mask = otSource->mask; + fitSource->mask = psImageCopy (NULL, refSource->mask, PS_TYPE_U8); + fitSource->pixels = psImageCopy (NULL, refSource->pixels, PS_TYPE_F32); + fitSource->weight = psImageCopy (NULL, refSource->weight, PS_TYPE_F32); // set model to unit peak, zero sky (we assume sky is subtracted) + // XXX do this with the actual model or use keep the original guess? model->params->data.F32[PM_PAR_SKY] = 0.0; model->params->data.F32[PM_PAR_I0] = 1.0; // fill in the model pixel values - psImageInit (flux, 0.0); - psImageKeepCircle (mask, x, y, model->radiusFit, "OR", PM_MASK_MARK); - pmModelAdd (flux, mask, model, false, false); - - // calculate nDOF (nPix - 1) - // int Nmaskpix = psImageCountPixelMask (mask, allArray, PM_MASK_SAT); - // model->nDOF = mask->numCols*mask->numRows - Nmaskpix - 1; + psImageInit (fitSource->pixels, 0.0); + psImageKeepCircle (fitSource->mask, x, y, model->radiusFit, "OR", PM_MASK_MARK); + pmModelAdd (fitSource->pixels, fitSource->mask, model, false, false); // save source in list - otSource->modelPSF = model; - index->data.U32[models->n] = i; - psArrayAdd (models, 100, otSource); - psFree (otSource); - } - psLogMsg ("psphot.emsemble", 4, "built models: %f (%ld objects)\n", psTimerMark ("psphot"), sources->n); - - // fill out the sparse matrix - psSparse *sparse = psSparseAlloc (models->n, 100); - psVector *weight = psVectorAlloc (models->n, PS_TYPE_F32); - psVector *errors = psVectorAlloc (models->n, PS_TYPE_F32); + if (isPSF) { + fitSource->modelPSF = psMemIncrRefCounter (model); + } else { + fitSource->modelEXT = psMemIncrRefCounter (model); + } + index->data.U32[fitSources->n] = i; + psArrayAdd (fitSources, 100, fitSource); + psFree (fitSource); + } + psLogMsg (__func__, 4, "built fitSources: %f (%ld objects)\n", psTimerMark ("psphot"), refSources->n); + + // vectors to store stats for each object + psVector *weight = psVectorAlloc (fitSources->n, PS_TYPE_F32); + 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 + psSparse *sparse = psSparseAlloc (fitSources->n, 100); psSparseBorder *border = psSparseBorderAlloc (sparse, 1); - for (int i = 0; i < models->n; i++) { + // fill out the sparse matrix elements and border elements (B) + // Ri is the current refSource of interest (sci pixels) + // Fi is the current fitSource of interest (fit pixels) + for (int i = 0; i < fitSources->n; i++) { int N = index->data.U32[i]; - pmSource *Fi = sources->data[N]; - pmSource *Mi = models->data[i]; + pmSource *Ri = refSources->data[N]; + pmSource *Fi = fitSources->data[i]; // scale by diagonal element (auto-cross-product) - r = pmSourceCrossProduct (Mi, Mi, CONSTANT_PHOTOMETRIC_WEIGHTS); + r = pmSourceCrossProduct (Fi, Fi, CONSTANT_PHOTOMETRIC_WEIGHTS); weight->data.F32[i] = r; // the formal error depends on the weighting scheme if (CONSTANT_PHOTOMETRIC_WEIGHTS) { - float var = pmSourceCrossProduct (Mi, Mi, false); + float var = pmSourceCrossProduct (Fi, Fi, false); errors->data.F32[i] = 1.0 / sqrt(var); } else { @@ -166,11 +131,11 @@ // find the image x model value - f = pmSourceCrossProduct (Fi, Mi, CONSTANT_PHOTOMETRIC_WEIGHTS); + f = pmSourceCrossProduct (Ri, Fi, CONSTANT_PHOTOMETRIC_WEIGHTS); psSparseVectorElement (sparse, i, f / r); // add the per-source weights (border region) - p = pmSourceWeight (Mi, CONST, CONSTANT_PHOTOMETRIC_WEIGHTS); - // px = pmSourceWeight (Mi, XTERM, CONSTANT_PHOTOMETRIC_WEIGHTS); - // py = pmSourceWeight (Mi, YTERM, CONSTANT_PHOTOMETRIC_WEIGHTS); + float p = pmSourceWeight (Fi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS); + // float px = pmSourceWeight (Fi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS); + // float py = pmSourceWeight (Fi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS); psSparseBorderElementB (border, i, 0, p); // psSparseBorderElementB (border, i, 1, px); @@ -178,47 +143,132 @@ // loop over all other stars following this one - for (int j = i + 1; j < models->n; j++) { - pmSource *Mj = models->data[j]; + for (int j = i + 1; j < fitSources->n; j++) { + pmSource *Fj = fitSources->data[j]; // skip over disjoint source images, break after last possible overlap - if (Mi->pixels->row0 + Mi->pixels->numRows < Mj->pixels->row0) break; - if (Mj->pixels->row0 + Mj->pixels->numRows < Mi->pixels->row0) continue; - if (Mi->pixels->col0 + Mi->pixels->numCols < Mj->pixels->col0) continue; - if (Mj->pixels->col0 + Mj->pixels->numCols < Mi->pixels->col0) continue; + if (Fi->pixels->row0 + Fi->pixels->numRows < Fj->pixels->row0) break; + if (Fj->pixels->row0 + Fj->pixels->numRows < Fi->pixels->row0) continue; + if (Fi->pixels->col0 + Fi->pixels->numCols < Fj->pixels->col0) continue; + if (Fj->pixels->col0 + Fj->pixels->numCols < Fi->pixels->col0) continue; // got an overlap; calculate cross-product and add to output array - f = pmSourceCrossProduct (Mi, Mj, CONSTANT_PHOTOMETRIC_WEIGHTS); + f = pmSourceCrossProduct (Fi, Fj, CONSTANT_PHOTOMETRIC_WEIGHTS); psSparseMatrixElement (sparse, j, i, f / r); } } - psLogMsg ("psphot.emsemble", 4, "built matrix: %f (%d elements)\n", psTimerMark ("psphot"), sparse->Nelem); + psSparseResort (sparse); + psLogMsg (__func__, 4, "built matrix: %f (%d elements)\n", psTimerMark ("psphot"), sparse->Nelem); + + // set the sky, sky_x, sky_y components of border matrix + // XXX ignore sky slope terms for now + SetBorderT (border, readout, refSources, false); + + psSparseConstraint constraint; + constraint.paramMin = 0.0; + constraint.paramMax = 1e8; + constraint.paramDelta = 1e8; + + // solve for normalization terms (need include local sky?) + psVector *norm = NULL; + psVector *skyfit = NULL; + psSparseBorderSolve (&norm, &skyfit, constraint, border, 3); + + // adjust fitSources, set refSources and subtract + for (int i = 0; i < fitSources->n; i++) { + int N = index->data.U32[i]; + pmSource *Ri = refSources->data[N]; + pmModel *model = pmSourceGetModel (NULL, Ri); + + // assign linearly-fitted normalization + if (isnan(norm->data.F32[i])) { + psAbort ("psphot", "ensemble source is nan"); + } + model->params->data.F32[PM_PAR_I0] = norm->data.F32[i]; + // XXX note: is the value of 'errors' modified by the sky fit? + model->dparams->data.F32[PM_PAR_I0] = errors->data.F32[i]; + + // subtract object + pmModelSub (Ri->pixels, Ri->mask, model, false, false); + Ri->mode |= PM_SOURCE_MODE_SUBTRACTED; + if (!final) Ri->mode |= PM_SOURCE_MODE_TEMPSUB; + // XXX not sure about the use of TEMPSUB + } + + // measure chisq for each source + for (int i = 0; final && (i < fitSources->n); i++) { + int N = index->data.U32[i]; + pmSource *Ri = refSources->data[N]; + pmModel *model = pmSourceGetModel (NULL, Ri); + + x = model->params->data.F32[PM_PAR_XPOS]; + y = model->params->data.F32[PM_PAR_YPOS]; + + psImageKeepCircle (Ri->mask, x, y, model->radiusFit, "OR", PM_MASK_MARK); + pmSourceChisq (model, Ri->pixels, Ri->mask, Ri->weight); + psImageKeepCircle (Ri->mask, x, y, model->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK)); + } + + psFree (index); + psFree (sparse); + psFree (fitSources); + psFree (norm); + psFree (weight); + psFree (errors); + + psLogMsg (__func__, 3, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot")); + return true; +} + +// calculate the weight terms for the sky fit component of the matrix +static bool SetBorderT (psSparseBorder *border, pmReadout *readout, psArray *refSources, bool fitSlope) { // generate the image-wide weight terms - // first, mask all object pixels - for (int i = 0; i < models->n; i++) { - int N = index->data.U32[i]; - pmSource *Fi = sources->data[N]; - pmSource *Mi = models->data[i]; - - // XXX need to select the correct mode (check this against the selections below - // - which model is appropriate for the given source? - // - which radius is appropriate for the given source? - // - where do x,y come from? - psImageKeepCircle (Mi->mask, x, y, model->radiusFit, "OR", PM_MASK_MARK); + // turn on MARK for all image pixels + psRegion fullArray = psRegionSet (0, 0, 0, 0); + fullArray = psRegionForImage (readout->mask, fullArray); + psImageMaskRegion (readout->mask, fullArray, "OR", PM_MASK_MARK); + + // turn off MARK for all object pixels + // we must use the refSources here since their pixels point at the science image + for (int i = 0; i < refSources->n; i++) { + pmSource *refSource = refSources->data[i]; + pmModel *model = pmSourceGetModel (NULL, refSource); + float x = model->params->data.F32[PM_PAR_XPOS]; + float y = model->params->data.F32[PM_PAR_YPOS]; + psImageMaskCircle (refSource->mask, x, y, model->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK)); } // accumulate the image statistics from the masked regions - for (int j = 0; j < numRows; j++) { - for (int i = 0; i < numCols; i++) { - w += image->data.F32[j][i]/weight->data.F32[j][i]; - x += i*image->data.F32[j][i]/weight->data.F32[j][i]; - y += j*image->data.F32[j][i]/weight->data.F32[j][i]; - x2 += i*i*image->data.F32[j][i]/weight->data.F32[j][i]; - xy += i*j*image->data.F32[j][i]/weight->data.F32[j][i]; - y2 += j*j*image->data.F32[j][i]/weight->data.F32[j][i]; + psF32 **image = readout->image->data.F32; + psF32 **weight = readout->weight->data.F32; + psU8 **mask = readout->mask->data.U8; + + double w, x, y, x2, xy, y2, xc, yc; + w = x = y = x2 = xy = y2 = 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; + w += image[j][i]/weight[j][i]; + if (!fitSlope) continue; + xc = i + col0; + yc = j + row0; + x += xc*image[j][i]/weight[j][i]; + y += yc*image[j][i]/weight[j][i]; + x2 += xc*xc*image[j][i]/weight[j][i]; + xy += xc*yc*image[j][i]/weight[j][i]; + y2 += yc*yc*image[j][i]/weight[j][i]; } } + + // turn off MARK for all image pixels + psImageMaskRegion (readout->mask, fullArray, "AND", PS_NOT_U8(PM_MASK_MARK)); + + // set the Border T elements psSparseBorderElementT (border, 0, 0, w); - if (0) { + if (fitSlope) { psSparseBorderElementT (border, 1, 0, x); psSparseBorderElementT (border, 2, 0, y); @@ -230,66 +280,4 @@ psSparseBorderElementT (border, 2, 2, y2); } - - psSparseConstraint constraint; - constraint.paramMin = 0.0; - constraint.paramMax = 1e8; - constraint.paramDelta = 1e8; - - // solve for normalization terms (need include local sky?) - psSparseResort (sparse); - psSparseBorderSolve (&xFit, &yFit, constraint, border, 3); - - // psVector *norm = psSparseSolve (NULL, constraint, sparse, 3); - - // adjust models, set sources and subtract - for (int i = 0; i < models->n; i++) { - int N = index->data.U32[i]; - pmSource *Fi = sources->data[N]; - pmSource *Mi = models->data[i]; - - // if we already have a PSF model, free it. - psFree (Fi->modelPSF); - - // need to increment counter so we can free models here and sources above - Fi->modelPSF = psMemIncrRefCounter (Mi->modelPSF); - - // assign linearly-fitted normalization - if (isnan(norm->data.F32[i])) { - psAbort ("psphot", "ensemble source is nan"); - } - Fi->modelPSF->params->data.F32[PM_PAR_I0] = norm->data.F32[i]; - // XXX I was using this before, but it is not strictly accurate - // XXX EAM : this factor of sqrt(2) makes the errors consistent, but I don't understand it - // Fi->modelPSF->dparams->data.F32[PM_PAR_I0] = sqrt(sqrt(2) * norm->data.F32[i] / (sparse->Bfj->data.F32[i] * weight->data.F32[i])); - Fi->modelPSF->dparams->data.F32[PM_PAR_I0] = errors->data.F32[i]; - - // subtract object - pmModelSub (Fi->pixels, Fi->mask, Fi->modelPSF, false, false); - Fi->mode |= PM_SOURCE_MODE_SUBTRACTED; - if (!final) Fi->mode |= PM_SOURCE_MODE_TEMPSUB; - } - - // measure chisq for each source - for (int i = 0; final && (i < models->n); i++) { - int N = index->data.U32[i]; - pmSource *Fi = sources->data[N]; - pmModel *model = Fi->modelPSF; - - x = model->params->data.F32[PM_PAR_XPOS]; - y = model->params->data.F32[PM_PAR_YPOS]; - - psImageKeepCircle (Fi->mask, x, y, model->radiusFit, "OR", PM_MASK_MARK); - pmSourceChisq (model, Fi->pixels, Fi->mask, Fi->weight); - psImageKeepCircle (Fi->mask, x, y, model->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK)); - } - - psFree (index); - psFree (sparse); - psFree (models); - psFree (norm); - psFree (weight); - psFree (errors); - - psLogMsg ("psphot.emsemble", 3, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot")); return true; } Index: /trunk/psphot/src/psphotGuessModels.c =================================================================== --- /trunk/psphot/src/psphotGuessModels.c (revision 10052) +++ /trunk/psphot/src/psphotGuessModels.c (revision 10052) @@ -0,0 +1,41 @@ +# include "psphot.h" + +// construct an initial PSF model for each object +bool psphotGuessModels (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) { + + // setup the PSF fit radius details + psphotInitRadiusPSF (recipe, psf->type); + + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + + // skip non-astronomical objects (very likely defects) + if (source->type == PM_SOURCE_TYPE_DEFECT) continue; + if (source->type == PM_SOURCE_TYPE_SATURATED) continue; + + // use the source moments, etc to guess basic model parameters + pmModel *modelEXT = pmSourceModelGuess (source, psf->type); + if (source->mode & PM_SOURCE_MODE_SATSTAR) { + modelEXT->params->data.F32[PM_PAR_XPOS] = source->moments->x; + modelEXT->params->data.F32[PM_PAR_YPOS] = source->moments->y; + } 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); + psFree (modelEXT); + + // XXX need to define the guess flux? + // set the fit radius based on the object flux limit and the model + psphotCheckRadiusPSF (readout, source, modelPSF); + + // set the source PSF model + source->modelPSF = modelPSF; + } + psLogMsg ("psphot.emsemble", 4, "built models: %f (%ld objects)\n", psTimerMark ("psphot"), sources->n); + return true; +} + +// XXX do we always know which model is supposed to be used? Index: /trunk/psphot/src/psphotReadout.c =================================================================== --- /trunk/psphot/src/psphotReadout.c (revision 10051) +++ /trunk/psphot/src/psphotReadout.c (revision 10052) @@ -78,4 +78,6 @@ } + psphotGuessModels (readout, sources, recipe, psf); + // linear PSF fit to peaks psphotEnsemblePSF (readout, sources, recipe, psf, FALSE);