Index: branches/eam_branches/20090715/psphot/src/psphot.h =================================================================== --- branches/eam_branches/20090715/psphot/src/psphot.h (revision 24799) +++ branches/eam_branches/20090715/psphot/src/psphot.h (revision 25022) @@ -157,5 +157,5 @@ bool psphotSubWithTest (pmSource *source, bool useState, psImageMaskType maskVal); bool psphotSetState (pmSource *source, bool curState, psImageMaskType maskVal); -bool psphotDeblendSatstars (psArray *sources, psMetadata *recipe); +bool psphotDeblendSatstars (pmReadout *readout, psArray *sources, psMetadata *recipe); bool psphotSourceSize (pmConfig *config, pmReadout *readout, psArray *sources, psMetadata *recipe, long first); Index: branches/eam_branches/20090715/psphot/src/psphotApResid.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotApResid.c (revision 24799) +++ branches/eam_branches/20090715/psphot/src/psphotApResid.c (revision 25022) @@ -168,4 +168,7 @@ if (source->mode & PM_SOURCE_MODE_POOR) SKIPSTAR ("POOR STAR"); + if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) SKIPSTAR ("EXTENDED"); + if (source->mode & PM_SOURCE_MODE_CR_LIMIT) SKIPSTAR ("COSMIC RAY"); + if (!isfinite(source->apMag) || !isfinite(source->psfMag)) { continue; @@ -317,4 +320,5 @@ */ +// XXX this still sucks... need a better way to estimate the error floor... bool psphotMagErrorScale (float *errorScale, float *errorFloor, psVector *dMag, psVector *dap, psVector *mask, int nGroup) { @@ -385,4 +389,5 @@ for (int i = 0; i < dSo->n; i++) { *errorFloor = dSo->data.F32[i]; + if (fabs(*errorFloor) <= FLT_EPSILON) continue; if (isfinite(*errorFloor)) break; } Index: branches/eam_branches/20090715/psphot/src/psphotCullPeaks.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotCullPeaks.c (revision 24799) +++ branches/eam_branches/20090715/psphot/src/psphotCullPeaks.c (revision 25022) @@ -18,25 +18,14 @@ nsigma_min = 0; } + float fPadding = psMetadataLookupF32(&status, recipe, "FOOTPRINT_CULL_NSIGMA_PAD"); + if (!status) { + fPadding = 0; + } const float skyStdev = psMetadataLookupF32(NULL, recipe, "SKY_STDEV"); - - return pmFootprintArrayCullPeaks(image, weight, footprints, - nsigma_delta, nsigma_min*skyStdev); -} - - -/* - * Cull an entire psArray of pmFootprints - * XXX drop this intermediate level function? - */ -psErrorCode -pmFootprintArrayCullPeaks(const psImage *img, // the image wherein lives the footprint - const psImage *weight, // corresponding variance image - psArray *footprints, // array of pmFootprints - const float nsigma_delta, // how many sigma above local background a peak - // needs to be to survive - const float min_threshold) { // minimum permitted coll height + const float min_threshold = nsigma_min*skyStdev; + for (int i = 0; i < footprints->n; i++) { pmFootprint *fp = footprints->data[i]; - if (pmFootprintCullPeaks(img, weight, fp, nsigma_delta, min_threshold) != PS_ERR_NONE) { + if (pmFootprintCullPeaks(image, weight, fp, nsigma_delta, fPadding, min_threshold) != PS_ERR_NONE) { return psError(PS_ERR_UNKNOWN, false, "Culling pmFootprint %d", fp->id); } @@ -45,220 +34,2 @@ return PS_ERR_NONE; } - - /* - * Examine the peaks in a pmFootprint, and throw away the ones that are not sufficiently - * isolated. More precisely, for each peak find the highest coll that you'd have to traverse - * to reach a still higher peak --- and if that coll's more than nsigma DN below your - * starting point, discard the peak. - */ -psErrorCode pmFootprintCullPeaks_OLD(const psImage *img, // the image wherein lives the footprint - const psImage *weight, // corresponding variance image - pmFootprint *fp, // Footprint containing mortal peaks - const float nsigma_delta, // how many sigma above local background a peak - // needs to be to survive - const float min_threshold) { // minimum permitted coll height - assert (img != NULL); assert (img->type.type == PS_TYPE_F32); - assert (weight != NULL); assert (weight->type.type == PS_TYPE_F32); - assert (img->row0 == weight->row0 && img->col0 == weight->col0); - assert (fp != NULL); - - if (fp->peaks == NULL || fp->peaks->n == 0) { // nothing to do - return PS_ERR_NONE; - } - - psRegion subRegion; // desired subregion; 1 larger than bounding box (grr) - subRegion.x0 = fp->bbox.x0; subRegion.x1 = fp->bbox.x1 + 1; - subRegion.y0 = fp->bbox.y0; subRegion.y1 = fp->bbox.y1 + 1; - const psImage *subImg = psImageSubset((psImage *)img, subRegion); - const psImage *subWt = psImageSubset((psImage *)weight, subRegion); - assert (subImg != NULL && subWt != NULL); - // - // We need a psArray of peaks brighter than the current peak. We'll fake this - // by reusing the fp->peaks but lying about n. - // - // We do this for efficiency (otherwise I'd need two peaks lists), and we are - // rather too chummy with psArray in consequence. But it works. - // - psArray *brightPeaks = psArrayAlloc(0); - psFree(brightPeaks->data); - brightPeaks->data = psMemIncrRefCounter(fp->peaks->data);// use the data from fp->peaks - // - // The brightest peak is always safe; go through other peaks trying to cull them - // - for (int i = 1; i < fp->peaks->n; i++) { // n.b. fp->peaks->n can change within the loop - const pmPeak *peak = fp->peaks->data[i]; - int x = peak->x - subImg->col0; - int y = peak->y - subImg->row0; - // - // Find the level nsigma below the peak that must separate the peak - // from any of its friends - // - assert (x >= 0 && x < subImg->numCols && y >= 0 && y < subImg->numRows); - const float stdev = sqrt(subWt->data.F32[y][x]); - float threshold = subImg->data.F32[y][x] - nsigma_delta*stdev; - if (isnan(threshold) || threshold < min_threshold) { -#if 1 // min_threshold is assumed to be below the detection threshold, - // so all the peaks are pmFootprint, and this isn't the brightest - // XXX mark peak to be dropped - (void)psArrayRemoveIndex(fp->peaks, i); - i--; // we moved everything down one - continue; -#else -#error n.b. We will be running LOTS of checks at this threshold, so only find the footprint once - threshold = min_threshold; -#endif - } - - // XXX EAM : if stdev >= 0, i'm not sure how this can ever be true? - if (threshold > subImg->data.F32[y][x]) { - threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON; - } - - // XXX this is a bit expensive: psImageAlloc for every peak contained in this footprint - // perhaps this should alloc a single ID image above and pass it in to be set. - - const int peak_id = 1; // the ID for the peak of interest - brightPeaks->n = i; // only stop at a peak brighter than we are - - // XXX optionally use the faster pmFootprintsFind if the subimage size is large (eg, M31) - - pmFootprint *peakFootprint = pmFootprintsFindAtPoint(subImg, threshold, brightPeaks, peak->y, peak->x); - brightPeaks->n = 0; // don't double free - psImage *idImg = pmSetFootprintID(NULL, peakFootprint, peak_id); - psFree(peakFootprint); - - // Check if any of the previous (brighter) peaks are within the footprint of this peak - // If so, the current peak is bogus; drop it. - int j; - for (j = 0; j < i; j++) { - const pmPeak *peak2 = fp->peaks->data[j]; - int x2 = peak2->x - subImg->col0; - int y2 = peak2->y - subImg->row0; - const int peak2_id = idImg->data.S32[y2][x2]; // the ID for some other peak - - if (peak2_id == peak_id) { // There's a brighter peak within the footprint above - ; // threshold; so cull our initial peak - (void)psArrayRemoveIndex(fp->peaks, i); - i--; // we moved everything down one - break; - } - } - if (j == i) { - j++; - } - - psFree(idImg); - } - - brightPeaks->n = 0; psFree(brightPeaks); - psFree((psImage *)subImg); - psFree((psImage *)subWt); - - return PS_ERR_NONE; -} - - /* - * Examine the peaks in a pmFootprint, and throw away the ones that are not sufficiently - * isolated. More precisely, for each peak find the highest coll that you'd have to traverse - * to reach a still higher peak --- and if that coll's more than nsigma DN below your - * starting point, discard the peak. - */ - -# define IN_PEAK 1 -psErrorCode pmFootprintCullPeaks(const psImage *img, // the image wherein lives the footprint - const psImage *weight, // corresponding variance image - pmFootprint *fp, // Footprint containing mortal peaks - const float nsigma_delta, // how many sigma above local background a peak - // needs to be to survive - const float min_threshold) { // minimum permitted coll height - assert (img != NULL); assert (img->type.type == PS_TYPE_F32); - assert (weight != NULL); assert (weight->type.type == PS_TYPE_F32); - assert (img->row0 == weight->row0 && img->col0 == weight->col0); - assert (fp != NULL); - - if (fp->peaks == NULL || fp->peaks->n == 0) { // nothing to do - return PS_ERR_NONE; - } - - psRegion subRegion; // desired subregion; 1 larger than bounding box (grr) - subRegion.x0 = fp->bbox.x0; subRegion.x1 = fp->bbox.x1 + 1; - subRegion.y0 = fp->bbox.y0; subRegion.y1 = fp->bbox.y1 + 1; - - psImage *subImg = psImageSubset((psImage *)img, subRegion); - psImage *subWt = psImageSubset((psImage *)weight, subRegion); - assert (subImg != NULL && subWt != NULL); - - psImage *idImg = psImageAlloc(subImg->numCols, subImg->numRows, PS_TYPE_S32); - - // We need a psArray of peaks brighter than the current peak. - // We reject peaks which either: - // 1) are below the local threshold - // 2) have a brighter peak within their threshold - - // allocate the full-sized array. if the final array is much smaller, we can realloc - // at that point. - psArray *brightPeaks = psArrayAllocEmpty(fp->peaks->n); - psArrayAdd (brightPeaks, 128, fp->peaks->data[0]); - - // The brightest peak is always safe; go through other peaks trying to cull them - for (int i = 1; i < fp->peaks->n; i++) { // n.b. fp->peaks->n can change within the loop - const pmPeak *peak = fp->peaks->data[i]; - int x = peak->x - subImg->col0; - int y = peak->y - subImg->row0; - // - // Find the level nsigma below the peak that must separate the peak - // from any of its friends - // - assert (x >= 0 && x < subImg->numCols && y >= 0 && y < subImg->numRows); - const float stdev = sqrt(subWt->data.F32[y][x]); - float threshold = subImg->data.F32[y][x] - nsigma_delta*stdev; - if (isnan(threshold) || threshold < min_threshold) { - // min_threshold is assumed to be below the detection threshold, - // so all the peaks are pmFootprint, and this isn't the brightest - continue; - } - - // XXX EAM : if stdev >= 0, i'm not sure how this can ever be true? - if (threshold > subImg->data.F32[y][x]) { - threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON; - } - - // XXX this is a bit expensive: psImageAlloc for every peak contained in this footprint - // perhaps this should alloc a single ID image above and pass it in to be set. - - // XXX optionally use the faster pmFootprintsFind if the subimage size is large (eg, M31) - - // at this point brightPeaks only has the peaks brighter than the current - pmFootprint *peakFootprint = pmFootprintsFindAtPoint(subImg, threshold, brightPeaks, peak->y, peak->x); - - // XXX need to supply the image here - // we set the IDs to either 1 (in peak) or 0 (not in peak) - pmSetFootprintID (idImg, peakFootprint, IN_PEAK); - psFree(peakFootprint); - - // Check if any of the previous (brighter) peaks are within the footprint of this peak - // If so, the current peak is bogus; drop it. - bool keep = true; - for (int j = 0; keep && (j < brightPeaks->n); j++) { - const pmPeak *peak2 = fp->peaks->data[j]; - int x2 = peak2->x - subImg->col0; - int y2 = peak2->y - subImg->row0; - if (idImg->data.S32[y2][x2] == IN_PEAK) - // There's a brighter peak within the footprint above threshold; so cull our initial peak - keep = false; - } - if (!keep) continue; - - psArrayAdd (brightPeaks, 128, fp->peaks->data[i]); - } - - psFree (fp->peaks); - fp->peaks = brightPeaks; - - psFree(idImg); - psFree(subImg); - psFree(subWt); - - return PS_ERR_NONE; -} - Index: branches/eam_branches/20090715/psphot/src/psphotDeblendSatstars.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotDeblendSatstars.c (revision 24799) +++ branches/eam_branches/20090715/psphot/src/psphotDeblendSatstars.c (revision 25022) @@ -1,5 +1,5 @@ # include "psphotInternal.h" -bool psphotDeblendSatstars (psArray *sources, psMetadata *recipe) { +bool psphotDeblendSatstars (pmReadout *readout, psArray *sources, psMetadata *recipe) { int N; @@ -9,6 +9,10 @@ int Nblend = 0; - float SAT_TEST_LEVEL = 50000; float SAT_MIN_RADIUS = 5.0; + + bool status; + pmCell *cell = readout->parent; + float SATURATION = 0.75*psMetadataLookupF32 (&status, cell->concepts, "CELL.SATURATION"); + float SAT_TEST_LEVEL = 0.5*SATURATION; // we need sources spatially-sorted to find overlaps @@ -32,5 +36,5 @@ // XXX filter? if (source->mode & PM_SOURCE_MODE_SATSTAR) continue; if (source->mode & PM_SOURCE_MODE_BLEND) continue; - if (source->peak->flux < SAT_TEST_LEVEL) continue; + if (source->peak->flux < SATURATION) continue; // save these for reference below @@ -45,6 +49,5 @@ psVector *xVec = contour->data[0]; psVector *yVec = contour->data[1]; - - // XXX should we filter based on the number of pixels in the contour? + if (xVec->n < 5) continue; // find the center of the contour (let's just use mid[x,y]) @@ -62,4 +65,9 @@ int yCenter = 0.5*(yMin + yMax); psFree (contour); + + psAssert (xCenter >= source->pixels->col0, "invalid shift in object center"); + psAssert (xCenter < source->pixels->col0 + source->pixels->numCols, "invalid shift in object center"); + psAssert (yCenter >= source->pixels->row0, "invalid shift in object center"); + psAssert (yCenter < source->pixels->row0 + source->pixels->numRows, "invalid shift in object center"); // reset the peak for this source to the value of the center pixel Index: branches/eam_branches/20090715/psphot/src/psphotGuessModels.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotGuessModels.c (revision 24799) +++ branches/eam_branches/20090715/psphot/src/psphotGuessModels.c (revision 25022) @@ -6,19 +6,4 @@ // 2) loop over the sources once and associate them with their cell // 3) define the threaded function to work with sources for a given cell - -// A guess for when the moments aren't available -static pmModel *wildGuess(pmSource *source, // Source for which to guess - pmPSF *psf // The point-spread function - ) -{ - pmModel *model = pmModelAlloc(psf->type); - psF32 *PAR = model->params->data.F32; - PAR[PM_PAR_SKY] = 0; - // XXX get this from the image pixels - PAR[PM_PAR_I0] = source->peak->flux; - PAR[PM_PAR_XPOS] = source->peak->xf; - PAR[PM_PAR_YPOS] = source->peak->yf; - return model; -} // construct an initial PSF model for each object @@ -81,11 +66,4 @@ } psFree(job); - -# if (0) - if (!psphotGuessModel_Unthreaded (readout, cells->data[j], psf, maskVal, markVal)) { - psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); - return false; - } -# endif } @@ -152,47 +130,46 @@ nSrc ++; - // XXX if a source is faint, it will not have moments measured. - // it must be modelled as a PSF. In this case, we need to use - // the peak centroid to get the coordinates and get the peak flux - // from the image? - pmModel *modelEXT; - if (!source->moments) { - modelEXT = wildGuess(source, psf); + // the guess central intensity comes from the peak: + float Io = source->peak->flux; + + // We have two options to get a guess for the object position: the position from the + // peak and the position from the moments. Use the peak position if (a) there are no + // moments and (b) the sources is not saturated + + bool useMoments = false; + useMoments = (source->mode & PM_SOURCE_MODE_SATSTAR); // we only want to try if SATSTAR is set, but.. + useMoments = (useMoments && source->moments); // can't if there are no moments + useMoments = (useMoments && source->moments->nPixels); // can't if the moments were not measured + useMoments = (useMoments && !(source->mode && PM_SOURCE_MODE_MOMENTS_FAILURE)); // can't if the moments failed... + + float Xo, Yo; + if (useMoments) { + Xo = source->moments->Mx; + Yo = source->moments->My; } else { - // use the source moments, etc to guess basic model parameters - modelEXT = pmSourceModelGuess (source, psf->type); // ALLOC X5 - if (!modelEXT) { - modelEXT = wildGuess(source, psf); - } - // these valuse are set in pmSourceModelGuess, should this rule be in there as well? - if (source->mode & PM_SOURCE_MODE_SATSTAR) { - modelEXT->params->data.F32[PM_PAR_XPOS] = source->moments->Mx; - modelEXT->params->data.F32[PM_PAR_YPOS] = source->moments->My; - } else { - modelEXT->params->data.F32[PM_PAR_XPOS] = source->peak->xf; - modelEXT->params->data.F32[PM_PAR_YPOS] = source->peak->yf; - } + Xo = source->peak->xf; + Yo = source->peak->yf; } - // set PSF parameters for this model (apply 2D shape model) - pmModel *modelPSF = pmModelFromPSF (modelEXT, psf); // ALLOC X5 + // set PSF parameters for this model (apply 2D shape model to coordinates Xo, Yo) + pmModel *modelPSF = pmModelFromPSFforXY(psf, Xo, Yo, Io); + if (modelPSF == NULL) { - psWarning ("Failed to determine PSF model at r,c = (%d,%d); trying centre of image", - source->peak->y, source->peak->x); + psWarning ("Failed to determine PSF model at (%f,%f); trying image center", Xo, Yo); - // Try the center of the image - modelEXT->params->data.F32[PM_PAR_XPOS] = 0.5*readout->image->numCols; - modelEXT->params->data.F32[PM_PAR_YPOS] = 0.5*readout->image->numRows; - modelPSF = pmModelFromPSF (modelEXT, psf); + float Xc = 0.5*readout->image->numCols; + float Yc = 0.5*readout->image->numRows; + pmModel *modelPSF = pmModelFromPSFforXY(psf, Xc, Yc, Io); if (modelPSF == NULL) { psError(PSPHOT_ERR_PSF, false, "Failed to determine PSF model at center of image"); - psFree(modelEXT); return false; } + + // Now set the object position at the expected location: + modelPSF->params->data.F32[PM_PAR_XPOS] = Xo; + modelPSF->params->data.F32[PM_PAR_YPOS] = Yo; source->mode |= PM_SOURCE_MODE_BADPSF; } - psFree (modelEXT); // FREE (x3) - // XXX need to define the guess flux? // set the fit radius based on the object flux limit and the model // this function affects the mask pixels @@ -209,85 +186,2 @@ return true; } - -# if (0) -// construct models only for sources in the specified region -bool psphotGuessModel_Unthreaded (pmReadout *readout, psArray *sources, pmPSF *psf, psImageMaskType maskVal, psImageMaskType markVal) { - - int nSrc = 0; - - for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - - // XXXX this is just for a test: use this to mark sources for which the model is measured - // check later that all are used. - source->mode |= PM_SOURCE_MODE_EXT_LIMIT; - - // skip non-astronomical objects (very likely defects) - if (source->type == PM_SOURCE_TYPE_DEFECT) continue; - if (source->type == PM_SOURCE_TYPE_SATURATED) continue; - if (!source->peak) continue; - - nSrc ++; - - // XXX if a source is faint, it will not have moments measured. - // it must be modelled as a PSF. In this case, we need to use - // the peak centroid to get the coordinates and get the peak flux - // from the image? - pmModel *modelEXT; - if (!source->moments) { - modelEXT = wildGuess(source, psf); - } else { - // use the source moments, etc to guess basic model parameters - modelEXT = pmSourceModelGuess (source, psf->type); // ALLOC - if (!modelEXT) { - modelEXT = wildGuess(source, psf); - } - // these valuse are set in pmSourceModelGuess, should this rule be in there as well? - if (source->mode & PM_SOURCE_MODE_SATSTAR) { - modelEXT->params->data.F32[PM_PAR_XPOS] = source->moments->Mx; - modelEXT->params->data.F32[PM_PAR_YPOS] = source->moments->My; - } 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); // ALLOC - if (modelPSF == NULL) { - psError(PSPHOT_ERR_PSF, false, - "Failed to determine PSF model at r,c = (%d,%d); trying centre of image", - source->peak->y, source->peak->x); - // - // Try the centre of the image - // - modelEXT->params->data.F32[PM_PAR_XPOS] = 0.5*readout->image->numCols; - modelEXT->params->data.F32[PM_PAR_YPOS] = 0.5*readout->image->numRows; - modelPSF = pmModelFromPSF (modelEXT, psf); - if (modelPSF == NULL) { - psError(PSPHOT_ERR_PSF, false, - "Failed to determine PSF model at centre of image"); - psFree(modelEXT); - return false; - } - - source->mode |= PM_SOURCE_MODE_BADPSF; - } - psFree (modelEXT); - - // XXX need to define the guess flux? - // set the fit radius based on the object flux limit and the model - // this function affects the mask pixels - psphotCheckRadiusPSF (readout, source, modelPSF, markVal); - - // set the source PSF model - source->modelPSF = modelPSF; - source->modelPSF->residuals = psf->residuals; - - pmSourceCacheModel (source, maskVal); - - } - - return true; -} -# endif Index: branches/eam_branches/20090715/psphot/src/psphotReadout.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotReadout.c (revision 24799) +++ branches/eam_branches/20090715/psphot/src/psphotReadout.c (revision 25022) @@ -87,5 +87,5 @@ // find blended neighbors of very saturated stars // XXX merge this with Basic Deblend? - psphotDeblendSatstars (sources, recipe); + psphotDeblendSatstars (readout, sources, recipe); // mark blended peaks PS_SOURCE_BLEND @@ -236,2 +236,3 @@ return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources); } + Index: branches/eam_branches/20090715/psphot/src/psphotReadoutMinimal.c =================================================================== --- branches/eam_branches/20090715/psphot/src/psphotReadoutMinimal.c (revision 24799) +++ branches/eam_branches/20090715/psphot/src/psphotReadoutMinimal.c (revision 25022) @@ -59,5 +59,5 @@ // find blended neighbors of very saturated stars // XXX merge this with Basic Deblend? - psphotDeblendSatstars (sources, recipe); + psphotDeblendSatstars (readout, sources, recipe); // mark blended peaks PS_SOURCE_BLEND