Index: /trunk/psphot/src/Makefile.am =================================================================== --- /trunk/psphot/src/Makefile.am (revision 13034) +++ /trunk/psphot/src/Makefile.am (revision 13035) @@ -31,5 +31,5 @@ psphotChoosePSF.c \ psphotGuessModels.c \ - psphotEnsemblePSF.c \ + psphotFitSourcesLinear.c \ psphotBlendFit.c \ psphotReplaceUnfit.c \ @@ -47,8 +47,6 @@ psphotModelTest.c \ psphotFitSet.c \ - psphotWeightBias.c \ psphotSourceFreePixels.c \ psphotSummaryPlots.c \ - psphotTestPSF.c \ psphotMergeSources.c \ psphotReadoutCleanup.c \ @@ -58,5 +56,4 @@ psphotMosaicSubimage.c \ psphotMakeResiduals.c \ - psphotTestSourceOutput.c \ psphotAddNoise.c Index: /trunk/psphot/src/psphot.h =================================================================== --- /trunk/psphot/src/psphot.h (revision 13034) +++ /trunk/psphot/src/psphot.h (revision 13035) @@ -41,4 +41,5 @@ bool psphotMomentsStats (pmReadout *readout, psMetadata *recipe, psArray *sources); bool psphotEnsemblePSF (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final); +bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final); bool psphotGuessModels (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf); bool psphotBlendFit (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf); Index: /trunk/psphot/src/psphotAddNoise.c =================================================================== --- /trunk/psphot/src/psphotAddNoise.c (revision 13034) +++ /trunk/psphot/src/psphotAddNoise.c (revision 13035) @@ -67,9 +67,6 @@ PAR[PM_PAR_SXY] = newshape.sxy; - if (add) { - pmModelAdd (source->weight, source->mask, model, false, false); - } else { - pmModelSub (source->weight, source->mask, model, false, false); - } + // XXX if we use pmSourceOp, the size (and possibly Io) will not be respected + pmSourceOp (source, PM_MODEL_OP_FULL, add); // restore original values Index: /trunk/psphot/src/psphotBlendFit.c =================================================================== --- /trunk/psphot/src/psphotBlendFit.c (revision 13034) +++ /trunk/psphot/src/psphotBlendFit.c (revision 13035) @@ -43,7 +43,9 @@ // limit selection to some SN limit + // XXX this should use peak? if (source->moments == NULL) continue; if (source->moments->SN < FIT_SN_LIM) continue; + // XXX this should use peak? if (source->moments->x < AnalysisRegion.x0) continue; if (source->moments->y < AnalysisRegion.y0) continue; @@ -65,9 +67,10 @@ // replace object in image if (source->mode & PM_SOURCE_MODE_SUBTRACTED) { - pmModelAdd (source->pixels, source->mask, source->modelPSF, false, false); + pmSourceAdd (source, PM_MODEL_OP_FULL); } Nfit ++; // try fitting PSFs, then try extended sources + // these functions subtract the resulting fitted source (XXX and update the modelFlux?) if (psphotFitBlend (readout, source, psf)) { psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->moments->x, source->moments->y); @@ -84,6 +87,7 @@ Nfail ++; - // re-subtract PSF for object, leave local sky - pmModelSub (source->pixels, source->mask, source->modelPSF, false, false); + // re-subtract the object, leave local sky + pmSourceCacheModel (source); + pmSourceSub (source, PM_MODEL_OP_FULL); source->mode |= PM_SOURCE_MODE_SUBTRACTED; source->mode |= PM_SOURCE_MODE_TEMPSUB; Index: /trunk/psphot/src/psphotChoosePSF.c =================================================================== --- /trunk/psphot/src/psphotChoosePSF.c (revision 13034) +++ /trunk/psphot/src/psphotChoosePSF.c (revision 13035) @@ -239,7 +239,10 @@ // set the mask and subtract the PSF model - psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PM_MASK_MARK); - pmModelSub (source->pixels, source->mask, source->modelPSF, false, false); - psImageKeepCircle (source->mask, x, y, RADIUS, "AND", PS_NOT_U8(PM_MASK_MARK)); + // XXX should we be using maskObj? should we be unsetting the mask? + // use pmModelSub because modelFlux has not been generated + assert (source->maskObj); + psImageKeepCircle (source->maskObj, x, y, RADIUS, "OR", PM_MASK_MARK); + pmModelSub (source->pixels, source->maskObj, source->modelPSF, PM_MODEL_OP_FULL); + psImageKeepCircle (source->maskObj, x, y, RADIUS, "AND", PS_NOT_U8(PM_MASK_MARK)); } Index: /trunk/psphot/src/psphotEnsemblePSF.c =================================================================== --- /trunk/psphot/src/psphotEnsemblePSF.c (revision 13034) +++ /trunk/psphot/src/psphotEnsemblePSF.c (revision 13035) @@ -91,7 +91,10 @@ // fill in the model pixel values + // XXX this function should just work with the (normalized) modelFlux pixels + // XXX we should not have to create a complete copy + // XXX review the use of the object mask: can we use a copied mask? psImageInit (fitSource->pixels, 0.0); psImageKeepCircle (fitSource->mask, x, y, model->radiusFit, "OR", PM_MASK_MARK); - pmModelAdd (fitSource->pixels, fitSource->mask, model, false, false); + pmModelAdd (fitSource->pixels, fitSource->mask, model, PM_MODEL_OP_FULL); // save source in list @@ -218,5 +221,5 @@ // subtract object - pmModelSub (Ri->pixels, Ri->mask, model, false, false); + pmModelSub (Ri->pixels, Ri->mask, model, PM_MODEL_OP_FULL); Ri->mode |= PM_SOURCE_MODE_SUBTRACTED; if (!final) Ri->mode |= PM_SOURCE_MODE_TEMPSUB; Index: /trunk/psphot/src/psphotFitSet.c =================================================================== --- /trunk/psphot/src/psphotFitSet.c (revision 13034) +++ /trunk/psphot/src/psphotFitSet.c (revision 13035) @@ -1,5 +1,5 @@ # include "psphotInternal.h" -// XXX this is not used in main psphot code +// This is not used in main psphot code (only in psphotModelTest.c) bool psphotFitSet (pmSource *source, pmModel *oneModel, char *fitset, pmSourceFitMode mode) { @@ -24,4 +24,5 @@ } + // XXX pmSourceFitSet must cache the modelFlux? pmSourceFitSet (source, modelSet, mode); @@ -32,5 +33,5 @@ for (int i = 0; i < modelSet->n; i++) { pmModel *model = modelSet->data[i]; - pmModelSub (source->pixels, source->mask, model, false, false); + pmModelSub (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL); fprintf (stderr, "output parameters (obj %d):\n", i); @@ -42,5 +43,5 @@ // write out psphotSaveImage (NULL, source->pixels, "resid.fits"); - psphotSaveImage (NULL, source->mask, "mask.fits"); + psphotSaveImage (NULL, source->maskObj, "mask.fits"); return true; } Index: /trunk/psphot/src/psphotFitSourcesLinear.c =================================================================== --- /trunk/psphot/src/psphotFitSourcesLinear.c (revision 13035) +++ /trunk/psphot/src/psphotFitSourcesLinear.c (revision 13035) @@ -0,0 +1,281 @@ +# 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 + +static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER); + +bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final) { + + bool status; + float x; + float y; + float f; + // float r; + + psTimerStart ("psphot"); + + // source analysis is done in spatial order + sources = psArraySort (sources, pmSourceSortByY); + + // 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 = + 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; + } + + // select the sources which will be used for the fitting analysis + 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; + if (final) { + if (source->mode & PM_SOURCE_MODE_SUBTRACTED) continue; + } else { + if (source->mode & PM_SOURCE_MODE_BLEND) continue; + } + + // save the original coords + x = source->peak->xf; + y = source->peak->yf; + + // 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; + + psArrayAdd (fitSources, 100, source); + } + psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "built fitSources: %f (%ld objects)\n", psTimerMark ("psphot"), sources->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); + int nBorder = (SKY_FIT_ORDER == 0) ? 1 : 3; + psSparseBorder *border = psSparseBorderAlloc (sparse, nBorder); + + // fill out the sparse matrix elements and border elements (B) + // SRCi is the current source of interest + // SRCj is a possibly overlapping source + for (int i = 0; i < fitSources->n; i++) { + pmSource *SRCi = fitSources->data[i]; + + // diagonal elements of the sparse matrix (auto-cross-product) + f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS); + psSparseMatrixElement (sparse, i, i, f); + + // the formal error depends on the weighting scheme + if (CONSTANT_PHOTOMETRIC_WEIGHTS) { + float var = pmSourceModelDotModel (SRCi, SRCi, false); + errors->data.F32[i] = 1.0 / sqrt(var); + } else { + errors->data.F32[i] = 1.0 / sqrt(f); + } + + + // find the image x model value + f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS); + psSparseVectorElement (sparse, i, f); + + // add the per-source weights (border region) + switch (SKY_FIT_ORDER) { + case 1: + f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS); + psSparseBorderElementB (border, i, 1, f); + f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS); + psSparseBorderElementB (border, i, 2, f); + + case 0: + f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS); + psSparseBorderElementB (border, i, 0, f); + break; + + default: + psAbort("Invalid SKY_FIT_ORDER %d\n", SKY_FIT_ORDER); + break; + } + + // loop over all other stars following this one + for (int j = i + 1; j < fitSources->n; j++) { + pmSource *SRCj = fitSources->data[j]; + + // 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; + + // got an overlap; calculate cross-product and add to output array + f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS); + psSparseMatrixElement (sparse, j, i, f); + } + } + + psSparseResort (sparse); + psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "built matrix: %f (%d elements)\n", psTimerMark ("psphot"), sparse->Nelem); + + // set the sky, sky_x, sky_y components of border matrix + SetBorderMatrixElements (border, readout, fitSources, CONSTANT_PHOTOMETRIC_WEIGHTS, SKY_FIT_ORDER); + + 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; + if (SKY_FIT_LINEAR) { + psSparseBorderSolve (&norm, &skyfit, constraint, border, 5); + fprintf (stderr, "skyfit: %f\n", skyfit->data.F32[0]); + } else { + norm = psSparseSolve (NULL, constraint, sparse, 5); + skyfit = NULL; + } + + // adjust I0 for fitSources and subtract + for (int i = 0; i < fitSources->n; i++) { + pmSource *source = fitSources->data[i]; + pmModel *model = pmSourceGetModel (NULL, source); + + // assign linearly-fitted normalization + if (isnan(norm->data.F32[i])) { + psAbort("linear fitted source is nan"); + } + 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 + pmSourceSub (source, PM_MODEL_OP_FULL); + source->mode |= PM_SOURCE_MODE_SUBTRACTED; + if (!final) source->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++) { + pmSource *source = fitSources->data[i]; + pmModel *model = pmSourceGetModel (NULL, source); + pmSourceChisq (model, source->pixels, source->maskObj, source->weight); + } + + // psFree (index); + psFree (sparse); + 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")); + 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, int SKY_FIT_ORDER) { + + // 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", PM_MASK_MARK); + + // 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_U8(PM_MASK_MARK)); + } + + // accumulate the image statistics from the masked regions + 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, 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 = weight[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_U8(PM_MASK_MARK)); + + // 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; +} Index: /trunk/psphot/src/psphotGrowthCurve.c =================================================================== --- /trunk/psphot/src/psphotGrowthCurve.c (revision 13034) +++ /trunk/psphot/src/psphotGrowthCurve.c (revision 13035) @@ -53,6 +53,7 @@ radius = psf->growth->radius->data.F32[i]; + // NOTE: we use pmModelAdd not pmSourceAdd because we are not working with a normal source psImageKeepCircle (mask, xc, yc, radius, "OR", PM_MASK_MARK); - pmModelAdd (image, mask, model, false, false); + pmModelAdd (image, mask, model, PM_MODEL_OP_FULL); pmSourcePhotometryAper (&apMag, model, image, mask); psImageKeepCircle (mask, xc, yc, radius, "AND", PS_NOT_U8(PM_MASK_MARK)); Index: /trunk/psphot/src/psphotGuessModels.c =================================================================== --- /trunk/psphot/src/psphotGuessModels.c (revision 13034) +++ /trunk/psphot/src/psphotGuessModels.c (revision 13035) @@ -56,4 +56,6 @@ source->modelPSF = modelPSF; source->modelPSF->residuals = psf->residuals; + + pmSourceCacheModel (source); } psLogMsg ("psphot.models", 4, "built models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot")); Index: /trunk/psphot/src/psphotMakeResiduals.c =================================================================== --- /trunk/psphot/src/psphotMakeResiduals.c (revision 13034) +++ /trunk/psphot/src/psphotMakeResiduals.c (revision 13035) @@ -76,8 +76,8 @@ if (model == NULL) continue; // model must be defined - psImage *image = psImageCopy (NULL, source->pixels, PS_TYPE_F32); - psImage *weight = psImageCopy (NULL, source->weight, PS_TYPE_F32); - psImage *mask = psImageCopy (NULL, source->mask, PS_TYPE_U8); - pmModelSub (image, mask, model, false, false); + psImage *image = psImageCopy (NULL, source->pixels, PS_TYPE_F32); + psImage *weight = psImageCopy (NULL, source->weight, PS_TYPE_F32); + psImage *mask = psImageCopy (NULL, source->maskView, PS_TYPE_U8); + pmModelSub (image, mask, model, PM_MODEL_OP_FUNC); // re-normalize image and weight @@ -253,4 +253,6 @@ # endif + psLogMsg ("psphot.pspsf", PS_LOG_INFO, "generate residuals for %ld objects: %f sec\n", input->n, psTimerMark ("residuals")); + psFree (xC); psFree (yC); @@ -272,6 +274,4 @@ psphotSaveImage (NULL, resid->mask, "resid.mk.fits"); - psLogMsg ("psphot.pspsf", PS_LOG_INFO, "generate residuals for %ld objects: %f sec\n", input->n, psTimerMark ("residuals")); - psf->residuals = resid; return true; Index: /trunk/psphot/src/psphotModelTest.c =================================================================== --- /trunk/psphot/src/psphotModelTest.c (revision 13034) +++ /trunk/psphot/src/psphotModelTest.c (revision 13035) @@ -2,4 +2,5 @@ static char DEFAULT_MODE[] = "EXT"; +// XXX consider this function : add more test information? bool psphotModelTest (pmReadout *readout, psMetadata *recipe) { @@ -158,5 +159,5 @@ // define the pixels used for the fit - psImageKeepCircle (source->mask, xObj, yObj, RADIUS, "OR", PM_MASK_MARK); + psImageKeepCircle (source->maskObj, xObj, yObj, RADIUS, "OR", PM_MASK_MARK); char *fitset = psMetadataLookupStr (&status, recipe, "TEST_FIT_SET"); @@ -170,5 +171,5 @@ // measure the source mags pmSourcePhotometryModel (&fitMag, model); - pmSourcePhotometryAper (&obsMag, model, source->pixels, source->mask); + pmSourcePhotometryAper (&obsMag, model, source->pixels, source->maskObj); fprintf (stderr, "ap: %f, fit: %f, apmifit: %f\n", obsMag, fitMag, obsMag - fitMag); @@ -177,5 +178,5 @@ // subtract object, leave local sky - pmModelSub (source->pixels, source->mask, model, false, false); + pmModelSub (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL); fprintf (stderr, "output parameters: \n"); @@ -186,5 +187,5 @@ // write out psphotSaveImage (NULL, source->pixels, "resid.fits"); - psphotSaveImage (NULL, source->mask, "mask.fits"); + psphotSaveImage (NULL, source->maskObj, "mask.fits"); exit (0); Index: /trunk/psphot/src/psphotOutput.c =================================================================== --- /trunk/psphot/src/psphotOutput.c (revision 13034) +++ /trunk/psphot/src/psphotOutput.c (revision 13035) @@ -59,5 +59,5 @@ for (int j = 0; j < source->pixels->numCols; j++) { // skip masked points - if (source->mask->data.U8[i][j]) { + if (source->maskObj->data.U8[i][j]) { continue; } @@ -72,5 +72,5 @@ source->pixels->data.F32[i][j], 1.0 / source->weight->data.F32[i][j], - source->mask->data.U8[i][j]); + source->maskObj->data.U8[i][j]); } } @@ -104,5 +104,5 @@ for (int i = 0; (sources != NULL) && (i < sources->n); i++) { pmSource *source = (pmSource *) sources->data[i]; - pmModel *model = pmSourceSelectModel (source); + pmModel *model = pmSourceGetModel (NULL, source); if (model == NULL) continue; @@ -154,68 +154,2 @@ return header; } - -# if (0) -// output functions: we have several fixed modes (sx, obj, cmp) -void psphotOutput (pmReadout *readout, psMetadata *arguments) { - - bool status; - char *outputFile = NULL; - - psMetadata *header = pmReadoutGetHeader (readout); - - pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); - // sample PSF images?? - if (psfSample != NULL) psphotSamplePSFs (psf, readout->image, psfSample); - if (psfSample != NULL) psphotSamplePSFs (psf, readout->image, psfSample); - - if (psfFile != NULL) { - psMetadata *psfData = pmPSFtoMD (NULL, psf); - psMetadataConfigWrite (psfData, psfFile); - psFree (psfData); - } -} - -psImage *pmModelPSFatXY (psImage *image, pmModel *modelEXT, pmPSF *psf, int x, int y, int dx, int dy) { - - psRegion region = {x - dx, x + dx, y - dy, y + dy}; - psImage *view = psImageSubset (image, region); - psImage *sample = psImageCopy (NULL, view, PS_TYPE_F32); - psImageInit (sample, 0); - modelEXT->params->data.F32[2] = x; - modelEXT->params->data.F32[3] = y; - pmModel *modelPSF = pmModelFromPSF (modelEXT, psf); - pmModelAdd (sample, NULL, modelPSF, false, false); - psFree (modelPSF); - return (sample); -} - -bool psphotSamplePSFs (pmPSF *psf, psImage *image, char *output) { - - // make sample PSFs for 4 corners and the center - psImage *sample; - - // optional dump of all rough source data - if (output[0] == 0) return false; - - pmModel *modelEXT = pmModelAlloc (psf->type); - modelEXT->params->data.F32[0] = 0; - modelEXT->params->data.F32[1] = 1; - - psFits *fits = psFitsOpen (output, "w"); - - // the centers are in parent coordinates; they do not need to correspond to valid pixels... - sample = pmModelPSFatXY (image, modelEXT, psf, 25, 25, 25, 25); - psFitsWriteImage (fits, NULL, sample, 0); - sample = pmModelPSFatXY (image, modelEXT, psf, image->numCols - 25, image->numRows - 25, 25, 25); - psFitsWriteImage (fits, NULL, sample, 0); - sample = pmModelPSFatXY (image, modelEXT, psf, image->numCols - 25, 25, 25, 25); - psFitsWriteImage (fits, NULL, sample, 0); - sample = pmModelPSFatXY (image, modelEXT, psf, 25, image->numRows - 25, 25, 25); - psFitsWriteImage (fits, NULL, sample, 0); - sample = pmModelPSFatXY (image, modelEXT, psf, image->numCols / 2, image->numRows / 2, 25, 25); - psFitsWriteImage (fits, NULL, sample, 0); - - psFitsClose (fits); - return (TRUE); -} -# endif Index: /trunk/psphot/src/psphotRadialPlot.c =================================================================== --- /trunk/psphot/src/psphotRadialPlot.c (revision 13034) +++ /trunk/psphot/src/psphotRadialPlot.c (revision 13035) @@ -60,5 +60,5 @@ for (int iy = 0; iy < source->pixels->numRows; iy++) { for (int ix = 0; ix < source->pixels->numCols; ix++) { - if (source->mask->data.U8[iy][ix]) { + if (source->maskObj->data.U8[iy][ix]) { rb->data.F32[nb] = hypot (ix - Xo, iy - Yo) ; fb->data.F32[nb] = log10(source->pixels->data.F32[iy][ix]); Index: /trunk/psphot/src/psphotRadiusChecks.c =================================================================== --- /trunk/psphot/src/psphotRadiusChecks.c (revision 13034) +++ /trunk/psphot/src/psphotRadiusChecks.c (revision 13035) @@ -25,7 +25,8 @@ bool psphotCheckRadiusPSF (pmReadout *readout, pmSource *source, pmModel *model) { + psF32 *PAR = model->params->data.F32; + + // XXX do we have a better value for the sky noise level? not really... pmMoments *moments = source->moments; - // do we have a better value for the sky noise level? - // not really... // set the fit radius based on the object flux limit and the model @@ -46,9 +47,14 @@ } - bool status = pmSourceRedefinePixels (source, readout, model->params->data.F32[2], model->params->data.F32[3], model->radiusFit); + bool status = pmSourceRedefinePixels (source, readout, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->radiusFit); + + // set the mask to flag the excluded pixels + psImageKeepCircle (source->maskObj, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->radiusFit, "OR", PM_MASK_MARK); return status; } bool psphotCheckRadiusPSFBlend (pmReadout *readout, pmSource *source, pmModel *model, float dR) { + + psF32 *PAR = model->params->data.F32; pmMoments *moments = source->moments; @@ -63,5 +69,8 @@ } - bool status = pmSourceRedefinePixels (source, readout, model->params->data.F32[2], model->params->data.F32[3], model->radiusFit); + bool status = pmSourceRedefinePixels (source, readout, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->radiusFit); + + // set the mask to flag the excluded pixels + psImageKeepCircle (source->maskObj, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->radiusFit, "OR", PM_MASK_MARK); return status; } @@ -86,4 +95,6 @@ bool psphotCheckRadiusEXT (pmReadout *readout, pmSource *source, pmModel *model) { + psF32 *PAR = model->params->data.F32; + pmMoments *moments = source->moments; if (moments == NULL) return false; @@ -94,5 +105,8 @@ // redefine the pixels if needed - bool status = pmSourceRedefinePixels (source, readout, model->params->data.F32[2], model->params->data.F32[3], model->radiusFit); + bool status = pmSourceRedefinePixels (source, readout, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->radiusFit); + + // set the mask to flag the excluded pixels + psImageKeepCircle (source->maskObj, PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], model->radiusFit, "OR", PM_MASK_MARK); return status; } Index: /trunk/psphot/src/psphotReadout.c =================================================================== --- /trunk/psphot/src/psphotReadout.c (revision 13034) +++ /trunk/psphot/src/psphotReadout.c (revision 13035) @@ -126,5 +126,5 @@ // linear PSF fit to peaks - psphotEnsemblePSF (readout, sources, recipe, psf, FALSE); + psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE); if (dump) psphotSaveImage (NULL, readout->image, "image.v1.fits"); if (!strcasecmp (breakPt, "ENSEMBLE")) { @@ -141,5 +141,5 @@ // linear PSF fit to remaining peaks - psphotEnsemblePSF (readout, sources, recipe, psf, TRUE); + psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE); if (dump) psphotSaveImage (NULL, readout->image, "image.v4.fits"); if (!strcasecmp (breakPt, "PASS1")) { @@ -192,5 +192,5 @@ // linear PSF fit to remaining peaks - psphotEnsemblePSF (readout, sources, recipe, psf, TRUE); + psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE); if (dump) psphotSaveImage (NULL, readout->image, "image.v6.fits"); Index: /trunk/psphot/src/psphotReplaceUnfit.c =================================================================== --- /trunk/psphot/src/psphotReplaceUnfit.c (revision 13034) +++ /trunk/psphot/src/psphotReplaceUnfit.c (revision 13035) @@ -1,4 +1,5 @@ # include "psphotInternal.h" +// replace the flux for sources which failed bool psphotReplaceUnfit (psArray *sources) { @@ -15,10 +16,5 @@ replace: - if (source->modelPSF == NULL) continue; - - psTrace ("psphot", 3, "replacing object at %f,%f\n", - source->modelPSF->params->data.F32[PM_PAR_XPOS], source->modelPSF->params->data.F32[PM_PAR_YPOS]); - - pmModelAdd (source->pixels, source->mask, source->modelPSF, false, false); + pmSourceAdd (source, PM_MODEL_OP_FULL); source->mode &= ~PM_SOURCE_MODE_SUBTRACTED; source->mode &= ~PM_SOURCE_MODE_TEMPSUB; @@ -40,12 +36,5 @@ if (!(source->mode & PM_SOURCE_MODE_SUBTRACTED)) continue; - // select appropriate model - pmModel *model = pmSourceGetModel (NULL, source); - if (model == NULL) continue; // model must be defined - - psTrace ("psphot", 3, "replacing object at %f,%f\n", - model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]); - - pmModelAdd (source->pixels, source->mask, model, false, false); + pmSourceAdd (source, PM_MODEL_OP_FULL); source->mode &= ~PM_SOURCE_MODE_SUBTRACTED; } @@ -54,5 +43,5 @@ } -// add or sub source replace or if the source has +// add source, if the source has been subtracted (or if we ignore the state) bool psphotAddWithTest (pmSource *source, bool useState) { @@ -60,20 +49,13 @@ bool state = !(source->mode & PM_SOURCE_MODE_SUBTRACTED); if (state && useState) return true; - - // select appropriate model - pmModel *model = pmSourceGetModel (NULL, source); - if (model == NULL) return false; // model must be defined - - psTrace ("psphot", 3, "replacing object at %f,%f\n", - model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]); // replace the model if 1) state says it is missing or 2) useState is false (just do it) if (!state || !useState) { - pmModelAdd (source->pixels, source->mask, model, false, false); + pmSourceAdd (source, PM_MODEL_OP_FULL); } return true; } -// add or sub source replace or if the source has +// sub source, if the source has been added (or if we ignore the state) bool psphotSubWithTest (pmSource *source, bool useState) { @@ -82,14 +64,7 @@ if (state && useState) return true; - // select appropriate model - pmModel *model = pmSourceGetModel (NULL, source); - if (model == NULL) return false; // model must be defined - - psTrace ("psphot", 3, "replacing object at %f,%f\n", - model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]); - // replace the model if 1) state says it is missing or 2) useState is false (just do it) if (!state || !useState) { - pmModelSub (source->pixels, source->mask, model, false, false); + pmSourceSub (source, PM_MODEL_OP_FULL); } return true; @@ -102,17 +77,10 @@ bool newState = !(source->mode & PM_SOURCE_MODE_SUBTRACTED); if (curState == newState) return true; - - // select appropriate model - pmModel *model = pmSourceGetModel (NULL, source); - if (model == NULL) return false; // model must be defined - - psTrace ("psphot", 3, "replacing object at %f,%f\n", - model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]); if (curState && !newState) { - pmModelSub (source->pixels, source->mask, model, false, false); + pmSourceSub (source, PM_MODEL_OP_FULL); } if (newState && !curState) { - pmModelAdd (source->pixels, source->mask, model, false, false); + pmSourceAdd (source, PM_MODEL_OP_FULL); } return true; Index: /trunk/psphot/src/psphotSourceFits.c =================================================================== --- /trunk/psphot/src/psphotSourceFits.c (revision 13034) +++ /trunk/psphot/src/psphotSourceFits.c (revision 13035) @@ -67,5 +67,5 @@ // these should never be invalid values // XXX drop these tests eventually - if (isnan(model->params->data.F32[PM_PAR_I0])) psAbort("nan in blend fit"); + if (isnan(model->params->data.F32[PM_PAR_I0])) psAbort("nan in blend fit"); if (isnan(model->params->data.F32[PM_PAR_XPOS])) psAbort("nan in blend fit"); if (isnan(model->params->data.F32[PM_PAR_YPOS])) psAbort("nan in blend fit"); @@ -83,7 +83,5 @@ // fit PSF model (set/unset the pixel mask) - psImageKeepCircle (source->mask, x, y, PSF->radiusFit, "OR", PM_MASK_MARK); pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF); - psImageKeepCircle (source->mask, x, y, PSF->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK)); // correct model chisq for flux trend @@ -110,5 +108,8 @@ blend->modelPSF = psMemIncrRefCounter (model); psTrace ("psphot", 5, "fitted blend as PSF\n"); - pmModelSub (source->pixels, source->mask, model, false, false); + + // build cached model and subtract + pmSourceCacheModel (blend); + pmSourceSub (blend, PM_MODEL_OP_FULL); blend->mode |= PM_SOURCE_MODE_SUBTRACTED; blend->mode &= ~PM_SOURCE_MODE_TEMPSUB; @@ -127,8 +128,12 @@ psFree (sourceSet); - psTrace ("psphot", 5, "fitted primary as PSF\n"); - pmModelSub (source->pixels, source->mask, PSF, false, false); + // save the new, successful model psFree (source->modelPSF); source->modelPSF = PSF; + psTrace ("psphot", 5, "fitted primary as PSF\n"); + + // build cached model and subtract + pmSourceCacheModel (source); + pmSourceSub (source, PM_MODEL_OP_FULL); source->mode |= PM_SOURCE_MODE_SUBTRACTED; source->mode &= ~PM_SOURCE_MODE_TEMPSUB; @@ -138,5 +143,4 @@ bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf) { - float x, y; double chiTrend; @@ -150,11 +154,6 @@ psphotCheckRadiusPSF (readout, source, PSF); - x = PSF->params->data.F32[PM_PAR_XPOS]; - y = PSF->params->data.F32[PM_PAR_YPOS]; - // fit PSF model (set/unset the pixel mask) - psImageKeepCircle (source->mask, x, y, PSF->radiusFit, "OR", PM_MASK_MARK); pmSourceFitModel (source, PSF, PM_SOURCE_FIT_PSF); - psImageKeepCircle (source->mask, x, y, PSF->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK)); // correct model chisq for flux trend @@ -168,11 +167,12 @@ } - psTrace ("psphot", 5, "fitted as PSF\n"); - pmModelSub (source->pixels, source->mask, PSF, false, false); - // free old model, save new model psFree (source->modelPSF); source->modelPSF = PSF; - + psTrace ("psphot", 5, "fitted as PSF\n"); + + // build cached model and subtract + pmSourceCacheModel (source); + pmSourceSub (source, PM_MODEL_OP_FULL); source->mode |= PM_SOURCE_MODE_SUBTRACTED; source->mode &= ~PM_SOURCE_MODE_TEMPSUB; @@ -264,9 +264,13 @@ // sub EXT psFree (DBL); - pmModelSub (source->pixels, source->mask, EXT, false, false); - psTrace ("psphot", 5, "blob as EXT: %f %f\n", EXT->params->data.F32[PM_PAR_XPOS], EXT->params->data.F32[PM_PAR_YPOS]); // save new model source->modelEXT = EXT; + + // build cached model and subtract + pmSourceCacheModel (source); + pmSourceSub (source, PM_MODEL_OP_FULL); + psTrace ("psphot", 5, "blob as EXT: %f %f\n", EXT->params->data.F32[PM_PAR_XPOS], EXT->params->data.F32[PM_PAR_YPOS]); + source->mode |= PM_SOURCE_MODE_SUBTRACTED; source->mode &= ~PM_SOURCE_MODE_TEMPSUB; @@ -276,7 +280,4 @@ // sub DLB psFree (EXT); - pmModelSub (source->pixels, source->mask, (pmModel *) DBL->data[0], false, false); - pmModelSub (source->pixels, source->mask, (pmModel *) DBL->data[1], false, false); - psTrace ("psphot", 5, "blob as DBL: %f %f\n", ONE->params->data.F32[PM_PAR_XPOS], ONE->params->data.F32[PM_PAR_YPOS]); // drop old model, save new second model... @@ -288,13 +289,15 @@ // copy most data from the primary source (modelEXT, blends stay NULL) - pmSource *newSrc = pmSourceAlloc (); - newSrc->peak = psMemIncrRefCounter (source->peak); - newSrc->pixels = psMemIncrRefCounter (source->pixels); - newSrc->weight = psMemIncrRefCounter (source->weight); - newSrc->mask = psMemIncrRefCounter (source->mask); - newSrc->moments = psMemIncrRefCounter (source->moments); + // XXX use pmSourceCopy? + pmSource *newSrc = pmSourceCopy (source); newSrc->modelPSF = psMemIncrRefCounter (DBL->data[1]); - newSrc->type = source->type; - newSrc->mode = source->mode; + + // build cached models and subtract + pmSourceCacheModel (source); + pmSourceSub (source, PM_MODEL_OP_FULL); + pmSourceCacheModel (newSrc); + pmSourceSub (newSrc, PM_MODEL_OP_FULL); + psTrace ("psphot", 5, "blob as DBL: %f %f\n", ONE->params->data.F32[PM_PAR_XPOS], ONE->params->data.F32[PM_PAR_YPOS]); + psArrayAdd (sources, 100, newSrc); psFree (newSrc); @@ -306,5 +309,5 @@ psArray *psphotFitDBL (pmReadout *readout, pmSource *source) { - float x, y, dx, dy; + float dx, dy; pmModel *DBL; pmModel *PSF; @@ -344,18 +347,10 @@ modelSet->data[1] = DBL; - x = source->moments->x; - y = source->moments->y; - // fit PSF model (set/unset the pixel mask) - psImageKeepCircle (source->mask, x, y, PSF->radiusFit, "OR", PM_MASK_MARK); pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF); - psImageKeepCircle (source->mask, x, y, PSF->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK)); - return (modelSet); } pmModel *psphotFitEXT (pmReadout *readout, pmSource *source) { - - float x, y; NfitEXT ++; @@ -369,7 +364,4 @@ psphotCheckRadiusEXT (readout, source, EXT); - x = EXT->params->data.F32[PM_PAR_XPOS]; - y = EXT->params->data.F32[PM_PAR_YPOS]; - if ((source->moments->Sx < 1e-3) || (source->moments->Sx < 1e-3)) { psTrace ("psphot", 5, "problem source: moments: %f %f\n", source->moments->Sx, source->moments->Sy); @@ -377,8 +369,5 @@ // fit EXT (not PSF) model (set/unset the pixel mask) - psImageKeepCircle (source->mask, x, y, EXT->radiusFit, "OR", PM_MASK_MARK); pmSourceFitModel (source, EXT, PM_SOURCE_FIT_EXT); - psImageKeepCircle (source->mask, x, y, EXT->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK)); - return (EXT); } Index: /trunk/psphot/src/psphotTestPSF.c =================================================================== --- /trunk/psphot/src/psphotTestPSF.c (revision 13034) +++ /trunk/psphot/src/psphotTestPSF.c (revision 13035) @@ -124,5 +124,5 @@ // subtract model flux - pmModelSub (source->pixels, source->mask, model, false, false); + pmModelSub (source->pixels, source->mask, model, PM_MODEL_OP_FULL); } fclose (f); Index: /trunk/psphot/src/psphotWeightBias.c =================================================================== --- /trunk/psphot/src/psphotWeightBias.c (revision 13034) +++ /trunk/psphot/src/psphotWeightBias.c (revision 13035) @@ -5,4 +5,5 @@ // only allow the normalization to vary // XXX eventually, we should be able to do this with linear fitting... +// XXX UNUSED bool psphotWeightBias (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) { @@ -59,5 +60,5 @@ // replace object in image - pmModelAdd (source->pixels, source->mask, source->modelPSF, false, false); + pmModelAdd (source->pixels, source->mask, source->modelPSF, PM_MODEL_OP_FULL); // make a temporary model (we don't keep the result of this analysis) @@ -76,5 +77,5 @@ // re-subtract PSF for object, leave local sky - pmModelSub (source->pixels, source->mask, source->modelPSF, false, false); + pmModelSub (source->pixels, source->mask, source->modelPSF, PM_MODEL_OP_FULL); PARp = source->modelPSF->params->data.F32;