Index: /trunk/psphot/src/Makefile.am =================================================================== --- /trunk/psphot/src/Makefile.am (revision 12688) +++ /trunk/psphot/src/Makefile.am (revision 12689) @@ -51,4 +51,5 @@ psphotSourcePlots.c \ psphotRadialPlot.c \ + psphotDeblendSatstars.c \ psphotMosaicSubimage.c \ psphotAddNoise.c \ Index: /trunk/psphot/src/psphot.h =================================================================== --- /trunk/psphot/src/psphot.h (revision 12688) +++ /trunk/psphot/src/psphot.h (revision 12689) @@ -111,2 +111,3 @@ bool psphotSubWithTest (pmSource *source, bool useState); bool psphotSetState (pmSource *source, bool curState); +bool psphotDeblendSatstars (psArray *sources, psMetadata *recipe); Index: /trunk/psphot/src/psphotAddNoise.c =================================================================== --- /trunk/psphot/src/psphotAddNoise.c (revision 12688) +++ /trunk/psphot/src/psphotAddNoise.c (revision 12689) @@ -4,9 +4,13 @@ bool status = false; - psTimerStart ("psphot"); + psEllipseShape oldshape; + psEllipseShape newshape; + psEllipseAxes axes; PS_ASSERT (readout, false); PS_ASSERT (readout->parent, false); PS_ASSERT (readout->parent->concepts, false); + + psTimerStart ("psphot"); // increase variance by factor*(object noise): @@ -15,4 +19,6 @@ float FACTOR = psMetadataLookupF32 (&status, recipe, "NOISE.FACTOR"); + PS_ASSERT (status, false); + float SIZE = psMetadataLookupF32 (&status, recipe, "NOISE.SIZE"); PS_ASSERT (status, false); @@ -39,5 +45,22 @@ } - model->params->data.F32[PM_PAR_I0] *= FACTOR; + psF32 *PAR = model->params->data.F32; + + // save original values + float oldI0 = PAR[PM_PAR_I0]; + oldshape.sx = PAR[PM_PAR_SXX]; + oldshape.sy = PAR[PM_PAR_SYY]; + oldshape.sxy = PAR[PM_PAR_SXY]; + + // increase size and height of source + axes = psEllipseShapeToAxes (oldshape, 20.0); + axes.major *= SIZE; + axes.minor *= SIZE; + newshape = psEllipseAxesToShape (axes); + PAR[PM_PAR_I0] = FACTOR*PS_SQR(oldI0); + PAR[PM_PAR_SXX] = newshape.sx; + PAR[PM_PAR_SYY] = newshape.sy; + PAR[PM_PAR_SXY] = newshape.sxy; + if (add) { pmModelAdd (source->weight, source->mask, model, false, false); @@ -45,4 +68,10 @@ pmModelSub (source->weight, source->mask, model, false, false); } + + // restore original values + PAR[PM_PAR_I0] = oldI0; + PAR[PM_PAR_SXX] = oldshape.sx; + PAR[PM_PAR_SYY] = oldshape.sy; + PAR[PM_PAR_SXY] = oldshape.sxy; } if (add) { Index: /trunk/psphot/src/psphotBasicDeblend.c =================================================================== --- /trunk/psphot/src/psphotBasicDeblend.c (revision 12688) +++ /trunk/psphot/src/psphotBasicDeblend.c (revision 12689) @@ -1,5 +1,4 @@ # include "psphot.h" -// 2006.02.07 : no leaks bool psphotBasicDeblend (psArray *sources, psMetadata *recipe) { Index: /trunk/psphot/src/psphotDeblendSatstars.c =================================================================== --- /trunk/psphot/src/psphotDeblendSatstars.c (revision 12689) +++ /trunk/psphot/src/psphotDeblendSatstars.c (revision 12689) @@ -0,0 +1,159 @@ +# include "psphot.h" + +bool psphotDeblendSatstars (psArray *sources, psMetadata *recipe) { + + int N; + pmSource *source; + + psTimerStart ("psphot"); + + int Nblend = 0; + float SAT_TEST_LEVEL = 50000; + float SAT_MIN_RADIUS = 5.0; + + // we need sources spatially-sorted to find overlaps + sources = psArraySort (sources, psphotSortByY); + + // source analysis is done in peak order (brightest first) + // we use an index for this so the spatial sorting is kept + psVector *SN = psVectorAlloc (sources->n, PS_DATA_F32); + for (int i = 0; i < SN->n; i++) { + source = sources->data[i]; + SN->data.F32[i] = source->peak->SN; + } + psVector *index = psVectorSortIndex (NULL, SN); + // this results in an index of increasing SN + + // examine sources in decreasing SN order + for (int i = sources->n - 1; i >= 0; i--) { + N = index->data.U32[i]; + source = sources->data[N]; + + // 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; + + // save these for reference below + int xPeak = source->peak->x; + int yPeak = source->peak->y; + + // generate a basic contour (set of x,y coordinates at-or-below flux level) + psArray *contour = pmSourceContour (source->pixels, xPeak, yPeak, SAT_TEST_LEVEL); + if (contour == NULL) continue; + + // contour consists of a set of X,Y coords giving the boundary + psVector *xVec = contour->data[0]; + psVector *yVec = contour->data[1]; + + // XXX should we filter based on the number of pixels in the contour? + + // find the center of the contour (let's just use mid[x,y]) + int xMin = xVec->data.F32[0]; + int xMax = xVec->data.F32[0]; + int yMin = yVec->data.F32[0]; + int yMax = yVec->data.F32[0]; + for (int j = 0; j < xVec->n; j++) { + xMin = PS_MIN (xMin, xVec->data.F32[j]); + xMax = PS_MAX (xMax, xVec->data.F32[j]); + yMin = PS_MIN (yMin, yVec->data.F32[j]); + yMax = PS_MAX (yMax, yVec->data.F32[j]); + } + int xCenter = 0.5*(xMin + xMax); + int yCenter = 0.5*(yMin + yMax); + psFree (contour); + + // reset the peak for this source to the value of the center pixel + source->peak->x = xCenter; + source->peak->y = yCenter; + source->peak->xf = xCenter; + source->peak->yf = yCenter; + + // temporary array for overlapping objects we find + psArray *overlap = psArrayAllocEmpty (100); + + // search backwards for overlapping sources + for (int j = N - 1; j >= 0; j--) { + pmSource *testSource = sources->data[j]; + if (testSource->peak->x < source->pixels->col0) continue; + if (testSource->peak->x >= source->pixels->col0 + source->pixels->numCols) continue; + if (testSource->peak->y < source->pixels->row0) break; + if (testSource->peak->y >= source->pixels->row0 + source->pixels->numRows) { + fprintf (stderr, "warning: invalid condition\n"); + continue; + } + psArrayAdd (overlap, 100, testSource); + } + + // search forwards for overlapping sources + for (int j = N + 1; j < sources->n; j++) { + pmSource *testSource = sources->data[j]; + if (testSource->peak->x < source->pixels->col0) continue; + if (testSource->peak->x >= source->pixels->col0 + source->pixels->numCols) continue; + if (testSource->peak->y < source->pixels->row0) { + fprintf (stderr, "warning: invalid condition\n"); + continue; + } + if (testSource->peak->y >= source->pixels->row0 + source->pixels->numRows) break; + psArrayAdd (overlap, 100, testSource); + } + + if (overlap->n == 0) { + psFree (overlap); + continue; + } + + // now find the contour which is at 0.5*SAT_TEST_LEVEL (xPeak, yPeak is valid high pixel) + contour = pmSourceContour (source->pixels, xPeak, yPeak, 0.5*SAT_TEST_LEVEL); + if (contour == NULL) { + psFree (overlap); + continue; + } + + // any peaks within this contour should be dropped (mark as blend, drop after loop is done) + // also drop any peaks which are too close to this peal + psVector *xv = contour->data[0]; + psVector *yv = contour->data[1]; + for (int k = 0; k < overlap->n; k++) { + pmSource *testSource = overlap->data[k]; + float radius = hypot((testSource->peak->x - xCenter), (testSource->peak->y - yCenter)); + if (radius < SAT_MIN_RADIUS) { + testSource->mode |= PM_SOURCE_MODE_BLEND; + Nblend ++; + continue; + } + for (int j = 0; j < xv->n; j+=2) { + if (fabs(yv->data.F32[j] - testSource->peak->y) > 0.5) continue; + if (xv->data.F32[j+0] > testSource->peak->x) break; + if (xv->data.F32[j+1] < testSource->peak->x) break; + testSource->mode |= PM_SOURCE_MODE_BLEND; + Nblend ++; + j = xv->n; // skip rest of contour + } + } + psFree (overlap); + psFree (contour); + } + + // drop the sources marked as BLEND + for (int i = 0; i < sources->n;) { + pmSource *source = sources->data[i]; + + if (!(source->mode & PM_SOURCE_MODE_BLEND)) { + i++; + continue; + } + + // shuffle the remaining sources forward + for (int j = i; j < sources->n - 1; j++) { + sources->data[j] = sources->data[j+1]; + } + psFree (source); + sources->n --; + } + + psFree (SN); + psFree (index); + + psLogMsg ("psphot", PS_LOG_INFO, "found %d satstar blend peaks, leaving %ld sources: %f sec\n", Nblend, sources->n, psTimerMark ("psphot")); + return true; +} Index: /trunk/psphot/src/psphotEvalPSF.c =================================================================== --- /trunk/psphot/src/psphotEvalPSF.c (revision 12688) +++ /trunk/psphot/src/psphotEvalPSF.c (revision 12689) @@ -73,4 +73,7 @@ case PM_MODEL_UNTRIED: source->mode &= ~PM_SOURCE_MODE_FITTED; + if (source->mode & PM_SOURCE_MODE_SATSTAR) { + psTrace ("psphot", 5, "satstar was not fitted"); + } return false; case PM_MODEL_BADARGS: @@ -79,4 +82,7 @@ default: source->mode |= PM_SOURCE_MODE_FAIL; + if (source->mode & PM_SOURCE_MODE_SATSTAR) { + psTrace ("psphot", 5, "satstar failed fit"); + } return false; } @@ -103,4 +109,7 @@ if (model->params->data.F32[PM_PAR_I0] <= 0) { source->mode |= PM_SOURCE_MODE_FAIL; + if (source->mode & PM_SOURCE_MODE_SATSTAR) { + psTrace ("psphot", 5, "satstar failed fit (peak below 0)"); + } return false; } @@ -130,4 +139,11 @@ keep &= (SN > PSF_MIN_SN); keep &= (Chi < PSF_MAX_CHI); + + if (source->mode & PM_SOURCE_MODE_SATSTAR) { + psTrace ("psphot", 5, "satstar fit results: %f, %f %d : %f %f : %f %f : %f %f\n", + model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS], keep, + dSX, nSx, dSY, nSy, SN, Chi); + } + if (keep) { if (source->mode & PM_SOURCE_MODE_PSFSTAR) { Index: /trunk/psphot/src/psphotFindPeaks.c =================================================================== --- /trunk/psphot/src/psphotFindPeaks.c (revision 12688) +++ /trunk/psphot/src/psphotFindPeaks.c (revision 12689) @@ -39,6 +39,9 @@ // optionally save example images under trace if (psTraceGetLevel("psphot") > 5) { - psphotSaveImage (NULL, smooth_im, "imsmooth.fits"); - psphotSaveImage (NULL, smooth_wt, "wtsmooth.fits"); + char name[64]; + sprintf (name, "imsmooth.v%d.fits", pass); + psphotSaveImage (NULL, smooth_im, name); + sprintf (name, "wtsmooth.v%d.fits", pass); + psphotSaveImage (NULL, smooth_wt, name); } @@ -58,5 +61,7 @@ // optionally save example images under trace if (psTraceGetLevel("psphot") > 5) { - psphotSaveImage (NULL, smooth_im, "snsmooth.fits"); + char name[64]; + sprintf (name, "snsmooth.v%d.fits", pass); + psphotSaveImage (NULL, smooth_im, name); } Index: /trunk/psphot/src/psphotReadout.c =================================================================== --- /trunk/psphot/src/psphotReadout.c (revision 12688) +++ /trunk/psphot/src/psphotReadout.c (revision 12689) @@ -61,4 +61,5 @@ return psphotReadoutCleanup (config, readout, recipe, NULL, NULL); } + pmPeaksWriteText (peaks, "oldpeaks.dat"); // construct sources and measure basic stats @@ -70,4 +71,6 @@ return psphotReadoutCleanup(config, readout, recipe, NULL, sources); } + + psphotDeblendSatstars (sources, recipe); // mark blended peaks PS_SOURCE_BLEND @@ -102,24 +105,24 @@ psphotGuessModels (readout, sources, recipe, psf); + psphotSaveImage (NULL, readout->image, "image.v0.fits"); + // linear PSF fit to peaks psphotEnsemblePSF (readout, sources, recipe, psf, FALSE); + psphotSaveImage (NULL, readout->image, "image.v1.fits"); if (!strcasecmp (breakPt, "ENSEMBLE")) { goto finish; } - // plot positive + negative sources (replace, then remove) - psphotSourcePlots (readout, sources, recipe); - // non-linear PSF and EXT fit to brighter sources psphotBlendFit (readout, sources, recipe, psf); + psphotSaveImage (NULL, readout->image, "image.v2.fits"); // replace all sources psphotReplaceAll (sources); - - // plot positive sources - // psphotSourcePlots (readout, sources, recipe); + psphotSaveImage (NULL, readout->image, "image.v3.fits"); // linear PSF fit to remaining peaks psphotEnsemblePSF (readout, sources, recipe, psf, TRUE); + psphotSaveImage (NULL, readout->image, "image.v4.fits"); if (!strcasecmp (breakPt, "PASS1")) { goto finish; @@ -137,4 +140,5 @@ // find the peaks in the image psArray *newPeaks = psphotFindPeaks (readout, recipe, 2); + pmPeaksWriteText (newPeaks, "newpeaks.dat"); // remove noise for subtracted objects @@ -153,4 +157,5 @@ // replace all sources psphotReplaceAll (sources); + psphotSaveImage (NULL, readout->image, "image.v5.fits"); // merge the newly selected peaks into the existing list @@ -160,6 +165,10 @@ // linear PSF fit to remaining peaks psphotEnsemblePSF (readout, sources, recipe, psf, TRUE); + psphotSaveImage (NULL, readout->image, "image.v6.fits"); finish: + + // plot positive sources + psphotSourcePlots (readout, sources, recipe); // measure aperture photometry corrections