Index: /trunk/psphot/Makefile =================================================================== --- /trunk/psphot/Makefile (revision 5671) +++ /trunk/psphot/Makefile (revision 5672) @@ -29,4 +29,6 @@ $(SRC)/psphotApplyPSF.$(ARCH).o \ $(SRC)/psphotFixedPSF.$(ARCH).o \ +$(SRC)/psphotEnsemblePSF.$(ARCH).o \ +$(SRC)/psphotReapplyPSF.$(ARCH).o \ $(SRC)/psphotFitGalaxies.$(ARCH).o \ $(SRC)/psphotOutput.$(ARCH).o \ @@ -36,8 +38,10 @@ $(SRC)/psphotDefinePixels.$(ARCH).o\ $(SRC)/psphotMagnitudes.$(ARCH).o \ +$(SRC)/psphotImageBackground.$(ARCH).o \ $(SRC)/psLine.$(ARCH).o \ $(SRC)/psModulesUtils.$(ARCH).o \ $(SRC)/pmPeaksSigmaLimit.$(ARCH).o \ $(SRC)/pmSourceFitFixed.$(ARCH).o \ +$(SRC)/psSparse.$(ARCH).o \ $(SRC)/psImageData.$(ARCH).o @@ -69,4 +73,9 @@ TEST = \ $(SRC)/psphotTest.$(ARCH).o \ +$(SRC)/psphotTestArguments.$(ARCH).o \ +$(SRC)/pmSourceContour.$(ARCH).o \ +$(SRC)/psImageData.$(ARCH).o \ +$(SRC)/psphotSetup.$(ARCH).o \ +$(SRC)/psModulesUtils.$(ARCH).o \ $(SRC)/psSparse.$(ARCH).o @@ -85,5 +94,5 @@ $(TEST): $(SRC)/psphot.h -INSTALL = psphot psphotModelTest +INSTALL = psphot psphotTest psphotModelTest # dependancy rules for binary code ######################### Index: /trunk/psphot/src/modelTestFitSource.c =================================================================== --- /trunk/psphot/src/modelTestFitSource.c (revision 5671) +++ /trunk/psphot/src/modelTestFitSource.c (revision 5672) @@ -7,5 +7,5 @@ bool status; int modelType; - float sky, obsMag, fitMag, value; + float obsMag, fitMag, value; char name[64]; @@ -96,8 +96,5 @@ // subtract object, leave local sky - sky = model->params->data.F32[0]; - model->params->data.F32[0] = 0; - pmSourceSubModel (source->pixels, source->mask, model, false); - model->params->data.F32[0] = sky; + pmSourceSubModel (source->pixels, source->mask, model, false, false); // write out Index: /trunk/psphot/src/pmPeaksSigmaLimit.c =================================================================== --- /trunk/psphot/src/pmPeaksSigmaLimit.c (revision 5671) +++ /trunk/psphot/src/pmPeaksSigmaLimit.c (revision 5672) @@ -23,5 +23,7 @@ // set peak threshold NSIGMA = psMetadataLookupF32 (&status, config, "PEAKS_NSIGMA_LIMIT"); - threshold = NSIGMA*sky->sampleStdev + sky->sampleMean; + + // threshold = NSIGMA*sky->sampleStdev + sky->sampleMean; + threshold = NSIGMA*sky->sampleStdev; psLogMsg ("psphot", 3, "threshold: %f DN\n", threshold); Index: /trunk/psphot/src/pmSourceContour.c =================================================================== --- /trunk/psphot/src/pmSourceContour.c (revision 5672) +++ /trunk/psphot/src/pmSourceContour.c (revision 5672) @@ -0,0 +1,210 @@ +# include "psphot.h" + +/****************************************************************************** +findValue(source, level, row, col, dir): a private function which determines +the column coordinate of the model function which has the value "level". If +dir equals 0, then you loop leftwards from the peak pixel, otherwise, +rightwards. + +XXX: reverse order of row,col args? + +XXX: Input row/col are in image coords. + +XXX: The result is returned in image coords. +*****************************************************************************/ +# define LEFT false +# define RIGHT true + +// return the first coordinate at or below the threshold in the requested direction +static int findContourNeg( + psImage *image, + float threshold, + int x, + int y, + bool right) +{ + + psTrace(__func__, 4, "---- %s() begin ----\n", __func__); + + // We define variables incr and lastColumn so that we can use the same loop + // whether we are stepping leftwards, or rightwards. + + int incr; + int subCol; + int lastColumn; + if (right) { + incr = 1; + lastColumn = image->numCols - 1; + } else { + incr = -1; + lastColumn = 0; + } + + subCol = x; + + while (subCol != lastColumn) { + float value = image->data.F32[y][subCol]; + if (value <= threshold) { + psTrace(__func__, 4, "---- %s() end ----\n", __func__); + return (subCol); + } + subCol += incr; + } + psTrace(__func__, 4, "---- %s() end ----\n", __func__); + return (lastColumn); +} + +// return the last coordinate at or below the threshold in the requested direction +static int findContourPos( + psImage *image, + float threshold, + int x, + int y, + bool right, + int xEnd) +{ + + psTrace(__func__, 4, "---- %s() begin ----\n", __func__); + + // We define variables incr and lastColumn so that we can use the same loop + // whether we are stepping leftwards, or rightwards. + + int incr; + int subCol; + int lastColumn; + if (right) { + incr = 1; + lastColumn = PS_MIN (image->numCols - 1, xEnd); + } else { + incr = -1; + lastColumn = PS_MAX (0, xEnd); + } + + subCol = x; + while (subCol != lastColumn) { + float value = image->data.F32[y][subCol]; + if (value >= threshold) { + psTrace(__func__, 4, "---- %s() end ----\n", __func__); + if (subCol == x) { + return (subCol); + } + return (subCol); + } + subCol += incr; + } + psTrace(__func__, 4, "---- %s() end ----\n", __func__); + return (lastColumn); +} + +/****************************************************************************** +new implementation of source contour function +*****************************************************************************/ +psArray *pmSourceContour_EAM (psImage *image, int x, int y, float threshold) { + + int xR, yR, x0, x1, x0s, x1s; + + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(image, false); + + // the requested point must be within the contour + if (image->data.F32[y][x] < threshold) return NULL; + + // Ensure that the starting column is allowable. + if (x < 0) return NULL; + if (y < 0) return NULL; + if (x >= image->numCols) return NULL; + if (y >= image->numRows) return NULL; + + // Allocate data for x/y pairs. + psVector *xVec = psVectorAlloc(100, PS_TYPE_F32); + psVector *yVec = psVectorAlloc(100, PS_TYPE_F32); + + // First row: find the left and right end-points + int Npt = 0; + + x0 = findContourNeg (image, threshold, x, y, LEFT); + x1 = findContourNeg (image, threshold, x, y, RIGHT); + xVec->data.F32[Npt + 0] = image->col0 + x0; + xVec->data.F32[Npt + 1] = image->col0 + x1; + yVec->data.F32[Npt + 0] = image->row0 + y; + yVec->data.F32[Npt + 1] = image->row0 + y; + Npt += 2; + + x0s = x0; + x1s = x1; + + // look for contour outline above row + xR = x0s; + yR = y + 1; + while (yR < image->numRows) { + if (image->data.F32[yR][xR] < threshold) { + x0 = findContourPos (image, threshold, xR, yR, RIGHT, x1); + if (x0 == x1) { + fprintf (stderr, "top: %d (%d - %d)\n", yR, xR, x1); + goto pt1; + } + x1 = findContourNeg (image, threshold, x0, yR, RIGHT); + x0--; + } else { + x0 = findContourNeg (image, threshold, xR, yR, LEFT); + x1 = findContourNeg (image, threshold, xR, yR, RIGHT); + } + fprintf (stderr, "pos: %d (%d - %d)\n", yR, x0, x1); + + xVec->data.F32[Npt + 0] = image->col0 + x0; + xVec->data.F32[Npt + 1] = image->col0 + x1; + yVec->data.F32[Npt + 0] = image->row0 + yR; + yVec->data.F32[Npt + 1] = image->row0 + yR; + Npt += 2; + + if (Npt >= xVec->nalloc - 1) { + psVectorRealloc (xVec, xVec->nalloc + 100); + psVectorRealloc (yVec, yVec->nalloc + 100); + } + yR ++; + xR = x0; + } + +pt1: + // look for contour outline below row + xR = x0s; + x1 = x1s; + yR = y - 1; + while (yR >= 0) { + if (image->data.F32[yR][xR] < threshold) { + x0 = findContourPos (image, threshold, xR, yR, RIGHT, x1); + if (x0 == x1) { + fprintf (stderr, "top: %d (%d - %d)\n", yR, xR, x1); + goto pt2; + } + x1 = findContourNeg (image, threshold, x0, yR, RIGHT); + x0--; + } else { + x0 = findContourNeg (image, threshold, xR, yR, LEFT); + x1 = findContourNeg (image, threshold, xR, yR, RIGHT); + } + fprintf (stderr, "neg: %d (%d - %d)\n", yR, x0, x1); + + xVec->data.F32[Npt + 0] = image->col0 + x0; + xVec->data.F32[Npt + 1] = image->col0 + x1; + yVec->data.F32[Npt + 0] = image->row0 + yR; + yVec->data.F32[Npt + 1] = image->row0 + yR; + Npt += 2; + + if (Npt >= xVec->nalloc - 1) { + psVectorRealloc (xVec, xVec->nalloc + 100); + psVectorRealloc (yVec, yVec->nalloc + 100); + } + yR --; + } +pt2: + xVec->n = Npt; + yVec->n = Npt; + + fprintf (stderr, "done\n"); + psArray *tmpArray = psArrayAlloc(2); + tmpArray->data[0] = (psPtr *) xVec; + tmpArray->data[1] = (psPtr *) yVec; + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(tmpArray); +} Index: /trunk/psphot/src/psphot.c =================================================================== --- /trunk/psphot/src/psphot.c (revision 5671) +++ /trunk/psphot/src/psphot.c (revision 5672) @@ -11,4 +11,6 @@ pmPSFClump psfClump; bool status; + + psphotModelGroupInit (); config = psphotArguments (&argc, argv); @@ -25,4 +27,7 @@ // measure image stats for initial guess sky = psphotImageStats (imdata, config); + + // psPolynomial2D *skyModel = psphotImageBackground (imdata, config, sky); + psphotImageBackground (imdata, config, sky); // find the peaks in the image @@ -53,4 +58,13 @@ switch (FITMODE) { + case -2: + psphotEnsemblePSF (imdata, config, sources, psf, sky); + psphotReapplyPSF (imdata, config, sources, psf, sky); + break; + + case -1: + psphotEnsemblePSF (imdata, config, sources, psf, sky); + break; + case 0: psphotFixedPSF (imdata, config, sources, psf, sky); @@ -63,5 +77,5 @@ case 2: - // fit extended objects with galaxy models + // fit extended objects with galaxy models psphotApplyPSF (imdata, config, sources, psf, sky); psphotFitGalaxies (imdata, config, sources, sky); @@ -69,5 +83,5 @@ case 3: - // fit extended objects with galaxy models + // fit extended objects with galaxy models psphotFixedPSF (imdata, config, sources, psf, sky); psphotFitGalaxies (imdata, config, sources, sky); Index: /trunk/psphot/src/psphot.h =================================================================== --- /trunk/psphot/src/psphot.h (revision 5671) +++ /trunk/psphot/src/psphot.h (revision 5672) @@ -3,11 +3,13 @@ # include // for unlink # include +# include +# include +# include +# include # include "psLibUtils.h" -# include "pmObjects_EAM.h" # include "psModulesUtils.h" -# include "pmPSF.h" -# include "pmPSFtry.h" -# include "pmModelGroup.h" # include "psSparse.h" + +# define psMemCopy(A)(psMemIncrRefCounter((A))) typedef struct { @@ -32,4 +34,5 @@ bool psphotSubtractPSF (pmSource *source); int psphotSortBySN (const void **a, const void **b); +int psphotSortByY (const void **a, const void **b); int psphotSaveImage (psMetadata *header, psImage *image, char *filename); bool psphotDefinePixels (pmSource *mySource, const eamReadout *imdata, psF32 x, psF32 y, psF32 Radius); @@ -61,2 +64,10 @@ psMetadata *modelTestArguments (int *argc, char **argv); bool modelTestFitSource (eamReadout *imdata, psMetadata *config); +bool psphotEnsemblePSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); +float psphotCrossProduct (pmSource *Mi, pmSource *Mj); +psPolynomial2D *psphotImageBackground (eamReadout *imdata, psMetadata *config, psStats *sky); +bool psphotReapplyPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky); + +pmModel *pmModelCopy (pmModel *model); +psArray *pmSourceContour_EAM (psImage *image, int x, int y, float threshold); +psMetadata *psphotTestArguments (int *argc, char **argv); Index: /trunk/psphot/src/psphotBasicDeblend.c =================================================================== --- /trunk/psphot/src/psphotBasicDeblend.c (revision 5672) +++ /trunk/psphot/src/psphotBasicDeblend.c (revision 5672) @@ -0,0 +1,77 @@ +# include "psphot.h" + +// try PSF models and select best option + +pmPSF *psphotBasicDeblend (psMetadata *config, psArray *sources, psStats *skystats) +{ + + // sources must be sorted by Y + // index must be sorted by SN + + // search for sources which fall within current box + + // generate source contour (1/4 peak counts) + + // find sources from overlap which fall in contour + // mark source + + // select the candidate PSF stars (pointers to original sources) + for (int i = 0; (i < sources->n) && (stars->n < NSTARS); i++) { + pmSource *source = sources->data[i]; + if (source->type != PM_SOURCE_PSFSTAR) continue; + psArrayAdd (stars, 200, source); + } + psTrace (".psphot.pspsf", 3, "selected candidate %d PSF objects\n", stars->n); + + // get the fixed PSF fit radius + // XXX EAM : check that PSF_FIT_RADIUS < SKY_OUTER_RADIUS + float RADIUS = psMetadataLookupF32 (&status, config, "PSF_FIT_RADIUS"); + + // get the list pointers for the PSF_MODEL entries + psMetadataItem *mdi = psMetadataLookup (config, "PSF_MODEL"); + if (mdi == NULL) psAbort ("psphotChoosePSF", "missing PSF_MODEL selection"); + + psList *list = (psList *) mdi->data.list; + psListIterator *iter = psListIteratorAlloc (list, PS_LIST_HEAD, FALSE); + + // set up an array to store the results + psArray *models = psArrayAlloc (list->n); + + // try each model option listed in config + for (int i = 0; i < models->n; i++) { + + item = psListGetAndIncrement (iter); + modelName = item->data.V; + + models->data[i] = pmPSFtryModel (stars, modelName, RADIUS); + psFree (modelName); + psFree (item); + } + psFree (iter); + // psFree (list); XXX EAM - is list freed with iter? + psFree (stars); + + // select the best of the models + // here we are using the clippedStdev on the metric as the indicator + try = models->data[0]; + int bestN = 0; + float bestM = try->psf->dApResid; + for (int i = 1; i < models->n; i++) { + try = models->data[i]; + float M = try->psf->dApResid; + if (M < bestM) { + bestM = M; + bestN = i; + } + } + + // keep only the selected model: + try = models->data[bestN]; + pmPSF *psf = psMemCopy(try->psf); + psFree (models); // keep only the pmPSF resulting from this analysis + + modelName = pmModelGetType (psf->type); + psLogMsg ("psphot.pspsf", 3, "selected psf model %s, ApResid: %f +/- %f\n", modelName, psf->ApResid, psf->dApResid); + + return (psf); +} Index: /trunk/psphot/src/psphotEnsemblePSF.c =================================================================== --- /trunk/psphot/src/psphotEnsemblePSF.c (revision 5672) +++ /trunk/psphot/src/psphotEnsemblePSF.c (revision 5672) @@ -0,0 +1,173 @@ +# include "psphot.h" + +bool psphotEnsemblePSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) { + + bool status; + float x; + float y; + float Sky; + + psTimerStart ("psphot"); + + // source analysis is done in spatial order + sources = psArraySort (sources, psphotSortByY); + + float OUTER_RADIUS = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS"); + float PSF_FIT_NSIGMA = psMetadataLookupF32 (&status, config, "PSF_FIT_NSIGMA"); + float PSF_FIT_PADDING = psMetadataLookupF32 (&status, config, "PSF_FIT_PADDING"); + + // set the object surface-brightness limit for fitted pixels + float FLUX_LIMIT = PSF_FIT_NSIGMA * sky->sampleStdev; + psLogMsg ("psphot.ensemble", 4, "fitting pixels with at least %f object counts\n", FLUX_LIMIT); + + // pre-calculate all model pixels + psArray *models = psArrayAlloc (sources->n); + for (int i = 0; i < sources->n; i++) { + pmSource *inSource = sources->data[i]; + pmSource *otSource = pmSourceAlloc (); + + // XXX EAM : add option to use FLT or PSF form + // use the source moments, etc to guess basic model parameters + pmModel *modelFLT = pmSourceModelGuess (inSource, psf->type); + + // set PSF parameters for this model + pmModel *model = pmModelFromPSF (modelFLT, psf); + psFree (modelFLT); + + // save the original coords + x = model->params->data.F32[2]; + y = model->params->data.F32[3]; + + // get function which specifies the radius at this the model hits the given flux + pmModelRadius modelRadius = pmModelRadius_GetFunction (psf->type); + + // set the fit radius based on the object flux limit and the model + model->radius = modelRadius (model->params, FLUX_LIMIT) + PSF_FIT_PADDING; + if (isnan(model->radius)) psAbort ("psphotEnsemblePSF", "error in radius"); + + // if needed, ask for more object pixels + if (model->radius > OUTER_RADIUS) { + psphotDefinePixels (inSource, imdata, x, y, model->radius); + } + + // make temporary copies of the image pixels and mask + otSource->mask = psImageCopy (NULL, inSource->mask, PS_TYPE_U8); + otSource->pixels = psImageCopy (NULL, inSource->pixels, PS_TYPE_F32); + + // build the model image + psImage *flux = otSource->pixels; + psImage *mask = otSource->mask; + + // XXX EAM : set model to unit peak, zero sky, maybe use peak (x,y) + // model->params->data.F32[2] = inSource->peak->x; + // model->params->data.F32[3] = inSource->peak->y; + model->params->data.F32[0] = 0.0; + model->params->data.F32[1] = 1.0; + + // fill in the model pixel values + psImageInit (flux, 0.0); + psImageKeepCircle (mask, x, y, model->radius, "OR", PSPHOT_MASK_MARKED); + pmSourceAddModel (flux, mask, model, false, false); + + // save source in list + otSource->modelPSF = model; + models->data[i] = otSource; + } + psLogMsg ("psphot.emsemble", 4, "built models: %f (%d objects)\n", psTimerMark ("psphot"), sources->n); + + float f; + + // fill out the sparse matrix + psSparse *sparse = psSparseAlloc (sources->n, 100); + for (int i = 0; i < sources->n; i++) { + pmSource *Fi = sources->data[i]; + pmSource *Mi = models->data[i]; + + // diagonal element (auto-cross-product) + f = psphotCrossProduct (Mi, Mi); + psSparseMatrixElement (sparse, i, i, f); + + // find the image x model value + f = psphotCrossProduct (Fi, Mi); + psSparseVectorElement (sparse, i, f); + + // loop over all other stars following this one + for (int j = i + 1; j < sources->n; j++) { + pmSource *Mj = models->data[j]; + + if (Mi->pixels->col0 + Mi->pixels->numCols < Mj->pixels->col0) break; + if (Mj->pixels->col0 + Mj->pixels->numCols < Mi->pixels->col0) continue; + if (Mi->pixels->row0 + Mi->pixels->numRows < Mj->pixels->row0) continue; + if (Mj->pixels->row0 + Mj->pixels->numRows < Mi->pixels->row0) continue; + + // got an overlap; calculate cross-product and add to output array + f = psphotCrossProduct (Mi, Mj); + psSparseMatrixElement (sparse, j, i, f); + } + } + psLogMsg ("psphot.emsemble", 4, "built matrix: %f (%d elements)\n", psTimerMark ("psphot"), sparse->Nelem); + + // solve for normalization terms (need include local sky) + psSparseResort (sparse); + psVector *norm = psSparseSolve (NULL, sparse, 3); + + // adjust models, set sources and subtract + for (int i = 0; i < sources->n; i++) { + pmSource *Fi = sources->data[i]; + pmSource *Mi = models->data[i]; + + Fi->modelPSF = Mi->modelPSF; + Fi->modelPSF->params->data.F32[1] = norm->data.F32[i]; + + // subtract object + pmSourceSubModel (Fi->pixels, Fi->mask, Fi->modelPSF, false, false); + + Fi->modelPSF->params->data.F32[0] = Sky; + // need to set this! + } + + // XXX EAM : need to free up many things here + + psLogMsg ("psphot.emsemble", 4, "apply models: %f\n", psTimerMark ("psphot")); + return true; +} + + +float psphotCrossProduct (pmSource *Mi, pmSource *Mj) { + + int Xs, Xe, Ys, Ye; + int xi, xj, yi, yj; + int xIs, xJs, yIs, yJs; + int xIe, yIe; + float flux; + + psImage *Pi = Mi->pixels; + psImage *Pj = Mj->pixels; + + psImage *Ti = Mi->mask; + psImage *Tj = Mj->mask; + + Xs = PS_MAX (Pi->col0, Pj->col0); + Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols); + + Ys = PS_MAX (Pi->row0, Pj->row0); + Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows); + + xIs = Xs - Pi->col0; + xJs = Xs - Pj->col0; + yIs = Ys - Pi->row0; + yJs = Ys - Pj->row0; + + xIe = Xe - Pi->col0; + yIe = Ye - Pi->row0; + + flux = 0; + for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) { + for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) { + if (Ti->data.U8[yi][xi]) continue; + if (Tj->data.U8[yj][xj]) continue; + flux += Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj]; + } + } + return (flux); +} Index: /trunk/psphot/src/psphotFitGalaxies.c =================================================================== --- /trunk/psphot/src/psphotFitGalaxies.c (revision 5671) +++ /trunk/psphot/src/psphotFitGalaxies.c (revision 5672) @@ -6,5 +6,4 @@ float x; float y; - float sky; int Nfit = 0; int Nfail = 0; @@ -104,9 +103,5 @@ // subtract object, leave local sky - sky = model->params->data.F32[0]; - model->params->data.F32[0] = 0; - pmSourceSubModel (source->pixels, source->mask, model, false); - model->params->data.F32[0] = sky; - + pmSourceSubModel (source->pixels, source->mask, model, false, false); } psLogMsg ("psphot", 3, "fit galaxies: %f sec for %d galaxies (%d failures, %d total iterations)\n", psTimerMark ("psphot"), Nfit, Nfail, Niter); Index: /trunk/psphot/src/psphotFixedPSF.c =================================================================== --- /trunk/psphot/src/psphotFixedPSF.c (revision 5671) +++ /trunk/psphot/src/psphotFixedPSF.c (revision 5672) @@ -8,5 +8,4 @@ float x; float y; - float Sky; int Nfit = 0; @@ -77,8 +76,5 @@ // subtract object, leave local sky - Sky = model->params->data.F32[0]; - model->params->data.F32[0] = 0; - pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false); - model->params->data.F32[0] = Sky; + pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false, false); } Index: /trunk/psphot/src/psphotImageBackground.c =================================================================== --- /trunk/psphot/src/psphotImageBackground.c (revision 5672) +++ /trunk/psphot/src/psphotImageBackground.c (revision 5672) @@ -0,0 +1,73 @@ +# include "psphot.h" + +psPolynomial2D *psphotImageBackground (eamReadout *imdata, psMetadata *config, psStats *sky) +{ + + // fit background to subset of image points within 3 sigma of sky + psTimerStart ("psphot"); + + psImage *image = imdata->image; + psImage *mask = imdata->mask; + + bool status = false; + int Npixels = image->numRows*image->numCols; + int Nsubset = PS_MIN (Npixels, psMetadataLookupF32 (&status, config, "IMSTATS_NPIX")); + psRandom *rnd = psRandomAlloc (PS_RANDOM_TAUS, 0); + + psVector *x = psVectorAlloc (Nsubset, PS_TYPE_F32); + psVector *y = psVectorAlloc (Nsubset, PS_TYPE_F32); + psVector *z = psVectorAlloc (Nsubset, PS_TYPE_F32); + + float Sky = sky->sampleMean; + float dSky = 3*sky->sampleStdev; + + // choose Nsubset points between 0 and Nx*Ny, convert to coords + int Nout = 0; + Npixels = image->numRows*image->numCols; + for (int i = 0; (i < 3*Nsubset) && (Nout < Nsubset); i++) { + double frnd = psRandomUniform (rnd); + int pixel = Npixels * frnd; + int ix = pixel / image->numCols; + int iy = pixel % image->numCols; + float value = image->data.F32[iy][ix]; + if (fabs(value - Sky) > dSky) continue; + z->data.F32[Nout] = value; + x->data.F32[Nout] = ix; + y->data.F32[Nout] = iy; + Nout ++; + } + x->n = y->n = z->n = Nout; + psLogMsg ("psphot", 5, "back: %f sec (select %d points)\n", psTimerMark ("psphot"), Nout); + + psPolynomial2D *skyModel = psPolynomial2DAlloc(1, 1, PS_POLYNOMIAL_ORD); + skyModel->mask[1][1] = 1; + + // psVector *mask = psVectorAlloc (Nout, PS_TYPE_U8); + // psVectorInit (mask, 0); + // psStats *fitstat = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); + // skyModel = psVectorFitPolynomial2D (skyModel, mask, 1, z, NULL, x, y); + + skyModel = psVectorFitPolynomial2D (skyModel, NULL, 0, z, NULL, x, y); + psLogMsg ("psphot", 5, "back: %f sec (fit model)\n", psTimerMark ("psphot")); + fprintf (stderr, "model: %f %f %f %f\n", + skyModel->coeff[0][0], skyModel->coeff[1][0], + skyModel->coeff[0][1], skyModel->coeff[1][1]); + + // this is a very inefficient way to evaluate the function.. + // lame, temporary step to zero masked pixels + for (int j = 0; j < image->numRows; j++) { + for (int i = 0; i < image->numCols; i++) { + if (mask->data.U8[j][i]) { + image->data.F32[j][i] = 0; + } else { + Sky = psPolynomial2DEval (skyModel, i, j); + image->data.F32[j][i] -= Sky; + } + } + } + psLogMsg ("psphot", 5, "back: %f sec (fit model)\n", psTimerMark ("psphot")); + + psphotSaveImage (imdata->header, imdata->image, "backsub.fits"); + + return (skyModel); +} Index: /trunk/psphot/src/psphotMagnitudes.c =================================================================== --- /trunk/psphot/src/psphotMagnitudes.c (revision 5671) +++ /trunk/psphot/src/psphotMagnitudes.c (revision 5672) @@ -24,8 +24,5 @@ // replace source flux - // XXX EAM : need to add 'sky' boolean option here - model->params->data.F32[0] = 0; - pmSourceAddModel (source->pixels, source->mask, model, false); - model->params->data.F32[0] = sky; + pmSourceAddModel (source->pixels, source->mask, model, false, false); // measure object photometry @@ -38,8 +35,5 @@ // subtract object, leave local sky - // XXX EAM : need to add 'sky' boolean option here - model->params->data.F32[0] = 0; - pmSourceSubModel (source->pixels, source->mask, model, false); - model->params->data.F32[0] = sky; + pmSourceSubModel (source->pixels, source->mask, model, false, false); // unmask aperture Index: /trunk/psphot/src/psphotMarkPSF.c =================================================================== --- /trunk/psphot/src/psphotMarkPSF.c (revision 5671) +++ /trunk/psphot/src/psphotMarkPSF.c (revision 5672) @@ -30,9 +30,30 @@ float nSx, nSy, Chi; - if (source->modelPSF == NULL) return (false); + pmModel *model = source->modelPSF; - // if object has fitted peak above saturation, label as SATSTAR + // do we actually have a valid PSF model? + if (model == NULL) return (false); + + // did the model fit fail for one or another reason? + switch (model->status) { + case PM_MODEL_BADARGS: + case PM_MODEL_UNTRIED: + source->type = PM_SOURCE_FAIL_FIT_PSF; // better choice? + return false; + case PM_MODEL_SUCCESS: + break; + case PM_MODEL_NONCONVERGE: + case PM_MODEL_OFFIMAGE: + psLogMsg ("psphot", 5, "PSF fit failed for %f, %f (%d iterations)\n", model->params->data.F32[2], model->params->data.F32[3], model->nIter); + source->type = PM_SOURCE_FAIL_FIT_PSF; // better choice? + return false; + default: + return false; + } + + // if the object has fitted peak above saturation, label as SATSTAR + // this is a valid PSF object, but ignore the other quality tests // remember: fit does not use saturated pixels (masked) - if (source->modelPSF->params->data.F32[1] >= SATURATE) { + if (model->params->data.F32[1] >= SATURATE) { if (source->type == PM_SOURCE_PSFSTAR) { psLogMsg ("psphot", 5, "PSFSTAR marked SATSTAR\n"); @@ -41,4 +62,13 @@ return (true); } + + // if the object has a fitted peak below 0, the fit did not converge cleanly + if (model->params->data.F32[1] < 0) { + source->type = PM_SOURCE_FAIL_FIT_PSF; + return (false); + } + + // if the source was predicted to be a SATSTAR, but it fitted below saturation, + // make a note to the user if (source->type == PM_SOURCE_SATSTAR) { psLogMsg ("psphot", 5, "SATSTAR marked GOODSTAR (fitted peak below saturation)\n"); @@ -46,10 +76,10 @@ } - SN = source->modelPSF->params->data.F32[1]/source->modelPSF->dparams->data.F32[1]; - SX = source->modelPSF->params->data.F32[4]; - SY = source->modelPSF->params->data.F32[5]; - dSX = source->modelPSF->dparams->data.F32[4]; - dSY = source->modelPSF->dparams->data.F32[5]; - Chi = source->modelPSF->chisq / source->modelPSF->nDOF; + SN = model->params->data.F32[1]/model->dparams->data.F32[1]; + SX = model->params->data.F32[4]; + SY = model->params->data.F32[5]; + dSX = model->dparams->data.F32[4]; + dSY = model->dparams->data.F32[5]; + Chi = model->chisq / model->nDOF; // swing of sigma_x,y in sigmas Index: /trunk/psphot/src/psphotModelGroupInit.c =================================================================== --- /trunk/psphot/src/psphotModelGroupInit.c (revision 5672) +++ /trunk/psphot/src/psphotModelGroupInit.c (revision 5672) @@ -0,0 +1,20 @@ +# include "psphot.h" + +# include "models/pmModel_WAUSS.c" + +static pmModelGroup userModels[] = { + {"PS_MODEL_WAUSS", 9, pmModelFunc_WAUSS, pmModelFlux_WAUSS, pmModelRadius_WAUSS, pmModelLimits_WAUSS, pmModelGuess_WAUSS, pmModelFromPSF_WAUSS, pmModelFitStatus_WAUSS}, +}; + + +void psphotModelGroupInit (void) +{ + + pmModelGroupInit (); + + int Nmodels = sizeof (userModels) / sizeof (pmModelGroup); + for (int i = 0; i < Nmodels; i++) { + pmModelGroupAdd (&userModels[i]); + } + return true; +} Index: /trunk/psphot/src/psphotOutput.c =================================================================== --- /trunk/psphot/src/psphotOutput.c (revision 5671) +++ /trunk/psphot/src/psphotOutput.c (revision 5672) @@ -583,5 +583,5 @@ unlink (filename); psFits *fits = psFitsAlloc (filename); - psFitsWriteImage (fits, header, image, 0); + psFitsWriteImage (fits, NULL, image, 0); psFree (fits); return (TRUE); Index: /trunk/psphot/src/psphotReapplyPSF.c =================================================================== --- /trunk/psphot/src/psphotReapplyPSF.c (revision 5672) +++ /trunk/psphot/src/psphotReapplyPSF.c (revision 5672) @@ -0,0 +1,82 @@ +# include "psphot.h" + +bool psphotReapplyPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) +{ + bool status; + float x, y, r; + int Nfit = 0; + int Nsub = 0; + int Niter = 0; + + psTimerStart ("psphot"); + + // source analysis is done in S/N order (brightest first) + sources = psArraySort (sources, psphotSortBySN); + + // we may set this differently here from the value used to mark likely saturated stars + float SATURATION = psMetadataLookupF32 (&status, config, "SATURATION"); + float PSF_MIN_SN = psMetadataLookupF32 (&status, config, "PSF_MIN_SN"); + float PSF_MAX_CHI = psMetadataLookupF32 (&status, config, "PSF_MAX_CHI"); + float PSF_SHAPE_NSIGMA = psMetadataLookupF32 (&status, config, "PSF_SHAPE_NSIGMA"); + + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + + // skip non-astronomical objects (very likely defects) + // XXX EAM : should we try these anyway? + if (source->type == PM_SOURCE_DEFECT) continue; + if (source->type == PM_SOURCE_SATURATED) continue; + + pmModel *model = source->modelPSF; + pmModel *save = pmModelCopy (model); + x = model->params->data.F32[2]; + y = model->params->data.F32[3]; + r = model->radius; + + // mask outside radius to limit the following actions to the source + psImageKeepCircle (source->mask, x, y, r, "OR", PSPHOT_MASK_MARKED); + + // replace object in image + pmSourceAddModel (source->pixels, source->mask, source->modelPSF, false, false); + + // fit PSF model (set/unset the pixel mask) + pmSourceFitModel (source, model, true); + + // model succeeded : tentatively keep it + Niter += model[0].nIter; + Nfit ++; + + // is it a good PSF model? + if (psphotMarkPSF (source, PSF_SHAPE_NSIGMA, PSF_MIN_SN, PSF_MAX_CHI, SATURATION)) { + psFree (save); + Nsub ++; + } else { + psFree (model); + source->modelPSF = save; + } + pmSourceSubModel (source->pixels, source->mask, source->modelPSF, false, false); + + // clear the masked region + psImageKeepCircle (source->mask, x, y, r, "AND", ~PSPHOT_MASK_MARKED); + } + + psLogMsg ("psphot", 3, "fit PSF models: %f sec for %d objects (%d total iterations)\n", psTimerMark ("psphot"), Nfit, Niter); + psLogMsg ("psphot", 4, "subtracted %d PSF objects\n", Nsub); + return (true); +} + +pmModel *pmModelCopy (pmModel *model) { + + pmModel *new = pmModelAlloc (model->type); + + new->chisq = model->chisq; + new->nIter = model->nIter; + + for (int i = 0; i < new->params->n; i++) { + new->params->data.F32[i] = model->params->data.F32[i]; + new->dparams->data.F32[i] = model->dparams->data.F32[i]; + } + + return (new); +} Index: /trunk/psphot/src/psphotSetup.c =================================================================== --- /trunk/psphot/src/psphotSetup.c (revision 5671) +++ /trunk/psphot/src/psphotSetup.c (revision 5672) @@ -84,8 +84,18 @@ // XXX EAM does the mask need to grow? float SATURATION = psMetadataLookupF32 (&status, config, "SATURATION"); - for (int i = 0; i < image->numRows; i++) { + for (int i = 0; status && (i < image->numRows); i++) { for (int j = 0; j < image->numCols; j++) { if (image->data.F32[i][j] >= SATURATION) { mask->data.U8[i][j] |= PSPHOT_MASK_SATURATED; + } + } + } + + // mask the pixels below min threshold + float MIN_VALID = psMetadataLookupF32 (&status, config, "MIN_VALID_PIXEL"); + for (int i = 0; status && (i < image->numRows); i++) { + for (int j = 0; j < image->numCols; j++) { + if (image->data.F32[i][j] < MIN_VALID) { + mask->data.U8[i][j] |= PSPHOT_MASK_INVALID; } } Index: /trunk/psphot/src/psphotSortBySN.c =================================================================== --- /trunk/psphot/src/psphotSortBySN.c (revision 5671) +++ /trunk/psphot/src/psphotSortBySN.c (revision 5672) @@ -18,2 +18,17 @@ } +// sort by Y (ascending) +int psphotSortByY (const void **a, const void **b) +{ + pmSource *A = *(pmSource **)a; + pmSource *B = *(pmSource **)b; + + psF32 fA = (A->moments == NULL) ? 0 : A->moments->y; + psF32 fB = (B->moments == NULL) ? 0 : B->moments->y; + + psF32 diff = fA - fB; + if (diff > FLT_EPSILON) return (+1); + if (diff < FLT_EPSILON) return (-1); + return (0); +} + Index: /trunk/psphot/src/psphotSubtractPSF.c =================================================================== --- /trunk/psphot/src/psphotSubtractPSF.c (revision 5671) +++ /trunk/psphot/src/psphotSubtractPSF.c (revision 5672) @@ -6,7 +6,5 @@ bool psphotSubtractPSF (pmSource *source) { - float sky; pmModel *model; - psImage *pixels; // only subtract successful fits @@ -23,11 +21,6 @@ } - pixels = source->pixels; - // subtract object, leave local sky - sky = model->params->data.F32[0]; - model->params->data.F32[0] = 0; - pmSourceSubModel (pixels, source->mask, model, false); - model->params->data.F32[0] = sky; + pmSourceSubModel (source->pixels, source->mask, model, false, false); // XXX EAM : amplify the weight matrix a la dophot? Index: /trunk/psphot/src/psphotTest.c =================================================================== --- /trunk/psphot/src/psphotTest.c (revision 5671) +++ /trunk/psphot/src/psphotTest.c (revision 5672) @@ -3,5 +3,32 @@ int main (int argc, char **argv) { - psSparseMatrixTest (); + bool status; + psMetadata *config = psphotTestArguments (&argc, argv); + eamReadout *imdata = psphotSetup (config); + + int x = psMetadataLookupF32 (&status, config, "TEST_FIT_X"); + int y = psMetadataLookupF32 (&status, config, "TEST_FIT_Y"); + float fraction = psMetadataLookupF32 (&status, config, "TEST_FRACTION"); + float threshold = psMetadataLookupF32 (&status, config, "TEST_THRESHOLD"); + + if (threshold == -1) { + threshold = fraction * imdata->image->data.F32[y][x]; + } + + psArray *contour = pmSourceContour_EAM (imdata->image, x, y, threshold); + if (contour == NULL) { + fprintf (stderr, "error: invalid contour\n"); + exit (1); + } + + psVector *xV = contour->data[0]; + psVector *yV = contour->data[1]; + + for (int i = 0; i < xV->n; i++) { + fprintf (stdout, "%d %f %f\n", i, xV->data.F32[i], yV->data.F32[i]); + } + + // psSparseMatrixTest (); + exit (0); } Index: /trunk/psphot/src/psphotTestArguments.c =================================================================== --- /trunk/psphot/src/psphotTestArguments.c (revision 5672) +++ /trunk/psphot/src/psphotTestArguments.c (revision 5672) @@ -0,0 +1,70 @@ +# include "psphot.h" +static int usage (); + +psMetadata *psphotTestArguments (int *argc, char **argv) { + + int N; + unsigned int Nfail; + int mode = PS_DATA_STRING | PS_META_REPLACE; + psF32 UseCoords_X, UseCoords_Y; + + // basic pslib options + psLogSetFormat ("M"); + psArgumentVerbosity (argc, argv); + + // identify options in args list - these go on config + bool UseCoords = false; + if ((N = psArgumentGet (*argc, argv, "-coords"))) { + UseCoords = true; + psArgumentRemove (N, argc, argv); + UseCoords_X = atof (argv[N]); + psArgumentRemove (N, argc, argv); + UseCoords_Y = atof (argv[N]); + psArgumentRemove (N, argc, argv); + } + + float threshold = -1; + if ((N = psArgumentGet (*argc, argv, "-threshold"))) { + psArgumentRemove (N, argc, argv); + threshold = atof (argv[N]); + psArgumentRemove (N, argc, argv); + } + + float fraction = 0.5; + if ((N = psArgumentGet (*argc, argv, "-frac"))) { + psArgumentRemove (N, argc, argv); + fraction = atof (argv[N]); + psArgumentRemove (N, argc, argv); + } + + if (*argc != 3) usage (); + + // load config information + psMetadata *config = psMetadataAlloc (); + psMetadataAdd (config, PS_LIST_HEAD, "PSF_MODEL", PS_DATA_METADATA_MULTI, "folder for psf model entries", NULL); + config = psMetadataConfigParse (config, &Nfail, argv[2], FALSE); + + // identify input image + psMetadataAdd (config, PS_LIST_HEAD, "IMAGE", mode, "", argv[1]); + + mode = PS_DATA_F32 | PS_META_REPLACE; + psMetadataAdd (config, PS_LIST_HEAD, "TEST_THRESHOLD", mode, "", threshold); + psMetadataAdd (config, PS_LIST_HEAD, "TEST_FRACTION", mode, "", fraction); + + if (UseCoords) { + psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_X", mode, "", UseCoords_X); + psMetadataAdd (config, PS_LIST_HEAD, "TEST_FIT_Y", mode, "", UseCoords_Y); + } + + return (config); +} + +static int usage () { + + fprintf (stderr, "USAGE: psphotModelTest (image.fits) (config)\n"); + fprintf (stderr, "options: \n"); + fprintf (stderr, " -coords x y : specify object center\n"); + fprintf (stderr, " -threshold (threshold) : (defaults to FRAC of peak\n"); + fprintf (stderr, " -frac (FRAC) : specify how to determine threshold\n"); + exit (2); +}