Index: /branches/eam_branches/ipp-20120805/psphot/src/psphot.h =================================================================== --- /branches/eam_branches/ipp-20120805/psphot/src/psphot.h (revision 34367) +++ /branches/eam_branches/ipp-20120805/psphot/src/psphot.h (revision 34368) @@ -281,5 +281,5 @@ bool psphotVisualShowSatStars (psMetadata *recipe, pmPSF *psf, psArray *sources); bool psphotVisualShowPSFModel (pmReadout *readout, pmPSF *psf); -bool psphotVisualPlotRadialProfile (int myKapa, pmSource *source, psImageMaskType maskVal); +bool psphotVisualPlotRadialProfile (int myKapa, pmSource *source, psImageMaskType maskVal, pmSourceMode showmode); bool psphotVisualPlotRadialProfiles (psMetadata *recipe, psArray *sources, pmSourceMode showmode); bool psphotVisualShowFlags (psArray *sources); @@ -295,4 +295,8 @@ bool psphotVisualClose(void); +int psphotKapaChannel (int channel); +void plotline (int myKapa, Graphdata *graphdata, float x0, float y0, float x1, float y1); + + bool psphotPetrosian (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal); bool psphotPetrosianRadialBins (pmSource *source, float radiusMax, float skynoise); Index: /branches/eam_branches/ipp-20120805/psphot/src/psphotDeblendSatstars.c =================================================================== --- /branches/eam_branches/ipp-20120805/psphot/src/psphotDeblendSatstars.c (revision 34367) +++ /branches/eam_branches/ipp-20120805/psphot/src/psphotDeblendSatstars.c (revision 34368) @@ -1,3 +1,19 @@ # include "psphotInternal.h" + +typedef struct { + float min; + float max; + float lower20; + float upper20; + int Npts; +} QuickStats; + +bool psImageQuickStats (QuickStats *stats, psImage *image, psRegion region); + +bool psphotTestRadialModel (pmSource *source, psVector **logRmodelOut, psVector **logFmodelOut, psVector *logFlux, psVector *logRad, float Rmax); +bool psphotTestRadialModelSub (pmSource *source, psVector *logRmodel, psVector *logFmodel, float Xo, float Yo, float Rmax, psImageMaskType maskVal); +bool psphotVisualRadialProfileSatstar (pmSource *source, psImageMaskType maskVal); +float InterpolateValues (float X0, float Y0, float X1, float Y1, float X); +bool psphotVisualScaleImage (int kapaFD, psImage *inImage, psImage *inMask, const char *name, float factor, int channel); // for now, let's store the detections on the readout->analysis for each readout @@ -17,4 +33,140 @@ bool psphotDeblendSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int fileIndex) { + + int N; + pmSource *source; + bool status; + + psTimerStart ("psphot.deblend.sat"); + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, fileIndex); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->newSources; + psAssert (sources, "missing sources?"); + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping satstar blend"); + return true; + } + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (maskVal); + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psImageMaskType maskSat = psMetadataLookupImageMask(&status, recipe, "MASK.SAT"); // Mask value for bad pixels + assert (maskSat); + + float SATURATION = NAN; + { + // XXX do we need to set this differently from the value used to mark saturated pixels? + pmCell *cell = readout->parent; + + // do not completely trust the values in the header... + float CELL_SATURATION = psMetadataLookupF32 (&status, cell->concepts, "CELL.SATURATION"); + float MIN_SATURATION = psMetadataLookupF32 (&status, recipe, "DEBLEND_MIN_SATURATION"); + if (!status || !isfinite(MIN_SATURATION)) { + MIN_SATURATION = 40000.0; + } + if (!isfinite(CELL_SATURATION)) { + SATURATION = MIN_SATURATION; + } else { + SATURATION = PS_MAX(MIN_SATURATION, CELL_SATURATION); + } + } + + // 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->moments->SN; + } + psVector *index = psVectorSortIndex (NULL, SN); + // this results in an index of increasing SN + + // psphotVisualPlotRadialProfiles (recipe, sources, PM_SOURCE_MODE_SATSTAR); + + int BIG_RADIUS = 250; + int BIG_SIGMA = BIG_RADIUS / 4.0; + + int display = psphotKapaChannel (1); + psphotVisualScaleImage (display, (psImage *) source->pixels->parent, NULL, "image", 1.0, 0); + + // examine sources in decreasing SN order + for (int i = sources->n - 1; i >= 0; i--) { + N = index->data.U32[i]; + source = sources->data[N]; + + bool isSat = source->mode & PM_SOURCE_MODE_SATSTAR; + + // for fairly bright stars, check the pixels near the peak again + if (source->moments->Peak > 0.5*SATURATION) { + // pmSourceRoughClass does this analysis, but uses a small (5x5) window. + // here we use a larger window since stacks can have some funny features + psRegion inner; + QuickStats stats; + // grow out the search radius until we have lower20 < 0.1*SATURATION + for (int radius = 2; radius < 30; radius ++) { + inner = psRegionForSquare (source->peak->x, source->peak->y, radius); + inner = psRegionForImage (source->maskView, inner); + psImageQuickStats (&stats, readout->image, inner); + if ((stats.Npts > 1) && (stats.lower20 < 0.1*SATURATION)) break; + } + int Nsatpix = psImageCountPixelMask (source->maskView, inner, maskSat); + // fprintf (stderr, "test object: %d,%d : %d vs %d : %f - %f - %f - %f\n", + // source->peak->x, source->peak->y, Nsatpix, stats.Npts, + // stats.min, stats.lower20, stats.upper20, stats.max); + if (Nsatpix > 1) isSat = true; + } + if (!isSat) continue; + + // For saturated stars, choose a much larger box NOTE this is slightly sleazy, but + // only slightly: pmSourceRedefinePixels uses the readout to pass the pointers to + // the parent image data. I guess the API could be simplified: we could recover + // this from the source in the function + + pmReadout tmpReadout; + tmpReadout.image = (psImage *)source->pixels->parent; + tmpReadout.mask = (psImage *)source->maskView->parent; + tmpReadout.variance = (psImage *)source->variance->parent; + + // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS + pmSourceRedefinePixels (source, &tmpReadout, source->peak->x, source->peak->y, BIG_RADIUS + 2); + + psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y); + status = pmSourceMoments (source, BIG_RADIUS, BIG_SIGMA, 0.0, 5.0, maskVal); + source->mode |= PM_SOURCE_MODE_BIG_RADIUS; + + // XXX visualize, model, and subtract + if (!psphotVisualRadialProfileSatstar (source, maskVal)) { + break; + } + + // generate radial profile, store on the source structure + } + // show the image after object have been subtracted + psphotVisualScaleImage (display, (psImage *) source->pixels->parent, NULL, "image", 1.0, 1); + + psFree (SN); + psFree (index); + + psLogMsg ("psphot", PS_LOG_INFO, "deblend satstar: %f sec\n", psTimerMark ("psphot.deblend.sat")); + return true; +} + +bool psphotDeblendSatstarsReadoutOld (pmConfig *config, const pmFPAview *view, const char *filerule, int fileIndex) { int N; @@ -78,6 +230,4 @@ psVector *index = psVectorSortIndex (NULL, SN); // this results in an index of increasing SN - - psphotVisualPlotRadialProfiles (recipe, sources, PM_SOURCE_MODE_SATSTAR); // examine sources in decreasing SN order @@ -222,205 +372,512 @@ } -bool psphotDeblendSatstarsReadoutOld (pmConfig *config, const pmFPAview *view, const char *filerule, int fileIndex) { - - int N; - pmSource *source; - bool status; - - psTimerStart ("psphot.deblend.sat"); - - int Nblend = 0; - float SAT_MIN_RADIUS = 5.0; - - // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, fileIndex); // File of interest - psAssert (file, "missing file?"); - - pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); - psAssert (readout, "missing readout?"); - - pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); - psAssert (detections, "missing detections?"); - - psArray *sources = detections->newSources; - psAssert (sources, "missing sources?"); - - if (!sources->n) { - psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping satstar blend"); - return true; - } - - // select the appropriate recipe information - psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); - psAssert (recipe, "missing recipe?"); - - pmCell *cell = readout->parent; - - float SATURATION = NAN; - - // do not completely trust the values in the header... - float CELL_SATURATION = psMetadataLookupF32 (&status, cell->concepts, "CELL.SATURATION"); - float MIN_SATURATION = psMetadataLookupF32 (&status, recipe, "DEBLEND_MIN_SATURATION"); - if (!status || !isfinite(MIN_SATURATION)) { - MIN_SATURATION = 40000.0; - } - if (!isfinite(CELL_SATURATION)) { - SATURATION = MIN_SATURATION; +static bool skipDisplay = false; +bool psphotVisualRadialProfileSatstar (pmSource *source, psImageMaskType maskVal) { + + int kapaImage = -1; + int kapaGraph = -1; + Graphdata graphdata; + + if (!pmVisualTestLevel("psphot.satstar", 3)) { + skipDisplay = true; + } + + float Rmax = 320.0; + + bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); + if (subtracted) pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + + int nPts = source->pixels->numRows * source->pixels->numCols; + psVector *rg = psVectorAllocEmpty (nPts, PS_TYPE_F32); + psVector *Rg = psVectorAllocEmpty (nPts, PS_TYPE_F32); + psVector *fg = psVectorAllocEmpty (nPts, PS_TYPE_F32); + psVector *Fg = psVectorAllocEmpty (nPts, PS_TYPE_F32); + psVector *rb = psVectorAllocEmpty (nPts, PS_TYPE_F32); + psVector *Rb = psVectorAllocEmpty (nPts, PS_TYPE_F32); + psVector *fb = psVectorAllocEmpty (nPts, PS_TYPE_F32); + + int ng = 0; + int nb = 0; + + float Xo = NAN; + float Yo = NAN; + + if (source->modelPSF) { + Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS] - source->pixels->col0; + Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS] - source->pixels->row0; } else { - SATURATION = PS_MAX(MIN_SATURATION, CELL_SATURATION); - } - float SAT_TEST_LEVEL = 0.5*SATURATION; - - // we need sources spatially-sorted to find overlaps - sources = psArraySort (sources, pmSourceSortByY); - - // 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->rawFlux; - } - 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->rawFlux < 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]; - if (xVec->n < 5) { - psFree(contour); - continue; - } - - // 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); - - 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 - source->peak->x = xCenter; - source->peak->y = yCenter; - source->peak->xf = xCenter; - source->peak->yf = yCenter; + Xo = source->moments->Mx - source->pixels->col0; + Yo = source->moments->My - source->pixels->row0; + } + + for (int iy = 0; iy < source->pixels->numRows; iy++) { + for (int ix = 0; ix < source->pixels->numCols; ix++) { + if ((source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal)) { + rb->data.F32[nb] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ; + // rb->data.F32[nb] = hypot (ix - Xo, iy - Yo) ; + Rb->data.F32[nb] = log10(rb->data.F32[nb]); + fb->data.F32[nb] = log10(source->pixels->data.F32[iy][ix]); + nb++; + } else { + rg->data.F32[ng] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ; + // rg->data.F32[ng] = hypot (ix - Xo, iy - Yo) ; + Rg->data.F32[ng] = log10(rg->data.F32[ng]); + Fg->data.F32[ng] = source->pixels->data.F32[iy][ix]; + fg->data.F32[ng] = log10(Fg->data.F32[ng]); + ng++; + } + } + } + rg->n = ng; + Rg->n = ng; + fg->n = ng; + Fg->n = ng; + rb->n = nb; + Rb->n = nb; + fb->n = nb; + + if (!skipDisplay) { + kapaImage = psphotKapaChannel (1); + if (kapaImage == -1) return false; + + kapaGraph = psphotKapaChannel (2); + if (kapaGraph == -1) return false; + + KapaSection section; // put the positive profile in one and the residuals in another? + + // first section : mag vs CR nSigma + section.dx = 1.0; + section.dy = 0.5; + section.x = 0.0; + section.y = 0.0; + section.bg = -1; + section.name = NULL; + psStringAppend (§ion.name, "linlog"); + KapaSetSection (kapaGraph, §ion); + psFree (section.name); + + // first section : mag vs CR nSigma + section.dx = 1.0; + section.dy = 0.5; + section.x = 0.0; + section.y = 0.5; + section.bg = -1; + section.name = NULL; + psStringAppend (§ion.name, "loglog"); + KapaSetSection (kapaGraph, §ion); + psFree (section.name); + + KapaInitGraph (&graphdata); + + // ** linlog ** + KapaSelectSection (kapaGraph, "linlog"); + + // examine sources to set data range + graphdata.xmin = -0.05; + graphdata.xmax = Rmax + 0.05; + graphdata.ymin = -0.05; + graphdata.ymax = +8.05; + KapaSetLimits (kapaGraph, &graphdata); + + KapaSetFont (kapaGraph, "helvetica", 14); + KapaBox (kapaGraph, &graphdata); + KapaSendLabel (kapaGraph, "radius (pixels)", KAPA_LABEL_XM); + KapaSendLabel (kapaGraph, "log flux (counts)", KAPA_LABEL_YM); + + graphdata.color = KapaColorByName ("black"); + graphdata.ptype = 2; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (kapaGraph, ng, &graphdata); + KapaPlotVector (kapaGraph, ng, rg->data.F32, "x"); + KapaPlotVector (kapaGraph, ng, fg->data.F32, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 0; + graphdata.size = 0.3; + graphdata.style = 2; + KapaPrepPlot (kapaGraph, nb, &graphdata); + KapaPlotVector (kapaGraph, nb, rb->data.F32, "x"); + KapaPlotVector (kapaGraph, nb, fb->data.F32, "y"); + + // ** loglog ** + KapaSelectSection (kapaGraph, "loglog"); + + // examine sources to set data range + graphdata.xmin = -1.51; + graphdata.xmax = log10(Rmax) + 0.02; + graphdata.ymin = -5.05; + graphdata.ymax = +8.05; + graphdata.color = KapaColorByName ("black"); + KapaSetLimits (kapaGraph, &graphdata); + + KapaSetFont (kapaGraph, "helvetica", 14); + KapaBox (kapaGraph, &graphdata); + KapaSendLabel (kapaGraph, "log radius (pixels)", KAPA_LABEL_XM); + KapaSendLabel (kapaGraph, "log flux (counts)", KAPA_LABEL_YM); + + graphdata.color = KapaColorByName ("black"); + graphdata.ptype = 2; + graphdata.size = 0.5; + graphdata.style = 2; + KapaPrepPlot (kapaGraph, ng, &graphdata); + KapaPlotVector (kapaGraph, ng, Rg->data.F32, "x"); + KapaPlotVector (kapaGraph, ng, fg->data.F32, "y"); + + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 0; + graphdata.size = 0.3; + graphdata.style = 2; + KapaPrepPlot (kapaGraph, nb, &graphdata); + KapaPlotVector (kapaGraph, nb, Rb->data.F32, "x"); + KapaPlotVector (kapaGraph, nb, fb->data.F32, "y"); + } + + // measure the radial profile + psVector *logFmodel = NULL; + psVector *logRmodel = NULL; + psphotTestRadialModel (source, &logRmodel, &logFmodel, Rg, Fg, Rmax); + + if (!skipDisplay) { + graphdata.color = KapaColorByName ("red"); + graphdata.ptype = 2; + graphdata.size = 1.0; + graphdata.style = 2; + KapaPrepPlot (kapaGraph, logRmodel->n, &graphdata); + KapaPlotVector (kapaGraph, logRmodel->n, logRmodel->data.F32, "x"); + KapaPlotVector (kapaGraph, logRmodel->n, logFmodel->data.F32, "y"); + } + + // subtract the model from the images + psphotTestRadialModelSub (source, logRmodel, logFmodel, Xo, Yo, Rmax, maskVal); + + psFree (logRmodel); + psFree (logFmodel); + + if (!skipDisplay && source->modelPSF) { + // generate model profiles (major and minor axis): + // create a model with theta = 0.0 so major and minor axes are equiv to x and y: + psEllipseShape rawShape, rotShape; + + rawShape.sx = source->modelPSF->params->data.F32[PM_PAR_SXX] / M_SQRT2; + rawShape.sy = source->modelPSF->params->data.F32[PM_PAR_SYY] / M_SQRT2; + rawShape.sxy = source->modelPSF->params->data.F32[PM_PAR_SXY]; + + psEllipseAxes axes = psEllipseShapeToAxes (rawShape, 20.0); + + axes.theta = 0.0; + + rotShape = psEllipseAxesToShape (axes); + + psVector *params = psVectorAlloc(source->modelPSF->params->n, PS_TYPE_F32); + for (int i = 0; i < source->modelPSF->params->n; i++) { + params->data.F32[i] = source->modelPSF->params->data.F32[i]; + } + params->data.F32[PM_PAR_SXX] = rotShape.sx * M_SQRT2; + params->data.F32[PM_PAR_SYY] = rotShape.sy * M_SQRT2; + params->data.F32[PM_PAR_SXY] = rotShape.sxy; + params->data.F32[PM_PAR_XPOS] = 0.0; + params->data.F32[PM_PAR_YPOS] = 0.0; + + psVector *rmod = psVectorAlloc(Rmax*10, PS_TYPE_F32); + psVector *fmaj = psVectorAlloc(Rmax*10, PS_TYPE_F32); + psVector *fmin = psVectorAlloc(Rmax*10, PS_TYPE_F32); + + psVector *coord = psVectorAlloc(2, PS_TYPE_F32); + + float r = 0.0; + for (int i = 0; i < rmod->n; i++) { + r = i*0.1; + rmod->data.F32[i] = r; + + coord->data.F32[1] = r; + coord->data.F32[0] = 0.0; + fmaj->data.F32[i] = log10(source->modelPSF->modelFunc (NULL, params, coord)); + + coord->data.F32[0] = r; + coord->data.F32[1] = 0.0; + fmin->data.F32[i] = log10(source->modelPSF->modelFunc (NULL, params, coord)); + } + psFree (coord); + psFree (params); + + float FWHM_MAJOR = 2.0*source->modelPSF->modelRadius (source->modelPSF->params, 0.5*source->modelPSF->params->data.F32[PM_PAR_I0]); + float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major); + if (FWHM_MAJOR < FWHM_MINOR) PS_SWAP (FWHM_MAJOR, FWHM_MINOR); + + psEllipseMoments emoments; + emoments.x2 = source->moments->Mxx; + emoments.xy = source->moments->Mxy; + emoments.y2 = source->moments->Myy; + axes = psEllipseMomentsToAxes (emoments, 20.0); + float MOMENTS_MAJOR = 2.355*axes.major; + float MOMENTS_MINOR = 2.355*axes.minor; + + float logHM = log10(0.5*source->modelPSF->params->data.F32[PM_PAR_I0]); + + // reset source Add/Sub state to recorded + if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); + + // ** linlog ** + KapaSelectSection (kapaGraph, "linlog"); + KapaGetGraphData (kapaGraph, &graphdata); + graphdata.color = KapaColorByName ("blue"); + graphdata.ptype = 0; + graphdata.size = 0.0; + graphdata.style = 0; + KapaPrepPlot (kapaGraph, rmod->n, &graphdata); + KapaPlotVector (kapaGraph, rmod->n, rmod->data.F32, "x"); + KapaPlotVector (kapaGraph, rmod->n, fmin->data.F32, "y"); + plotline (kapaGraph, &graphdata, graphdata.xmin, logHM, graphdata.xmax, logHM); + plotline (kapaGraph, &graphdata, 0.5*FWHM_MINOR, graphdata.ymin, 0.5*FWHM_MINOR, graphdata.ymax); + graphdata.ltype = 1; + plotline (kapaGraph, &graphdata, 0.5*MOMENTS_MINOR, graphdata.ymin, 0.5*MOMENTS_MINOR, graphdata.ymax); + graphdata.ltype = 0; - // 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; - } + graphdata.color = KapaColorByName ("green"); + graphdata.ptype = 0; + graphdata.size = 0.0; + graphdata.style = 0; + KapaPrepPlot (kapaGraph, rmod->n, &graphdata); + KapaPlotVector (kapaGraph, rmod->n, rmod->data.F32, "x"); + KapaPlotVector (kapaGraph, rmod->n, fmaj->data.F32, "y"); + plotline (kapaGraph, &graphdata, 0.5*FWHM_MAJOR, graphdata.ymin, 0.5*FWHM_MAJOR, graphdata.ymax); + graphdata.ltype = 1; + plotline (kapaGraph, &graphdata, 0.5*MOMENTS_MAJOR, graphdata.ymin, 0.5*MOMENTS_MAJOR, graphdata.ymax); + graphdata.ltype = 0; - // 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 peak - 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.deblend.sat")); + for (int i = 0; i < rmod->n; i++) { + rmod->data.F32[i] = log10(rmod->data.F32[i]); + } + + // ** loglog ** + KapaSelectSection (kapaGraph, "loglog"); + KapaGetGraphData (kapaGraph, &graphdata); + graphdata.color = KapaColorByName ("blue"); + graphdata.ptype = 0; + graphdata.size = 0.0; + graphdata.style = 0; + KapaPrepPlot (kapaGraph, rmod->n, &graphdata); + KapaPlotVector (kapaGraph, rmod->n, rmod->data.F32, "x"); + KapaPlotVector (kapaGraph, rmod->n, fmin->data.F32, "y"); + + graphdata.color = KapaColorByName ("green"); + graphdata.ptype = 0; + graphdata.size = 0.0; + graphdata.style = 0; + KapaPrepPlot (kapaGraph, rmod->n, &graphdata); + KapaPlotVector (kapaGraph, rmod->n, rmod->data.F32, "x"); + KapaPlotVector (kapaGraph, rmod->n, fmaj->data.F32, "y"); + + psFree (rmod); + psFree (fmin); + psFree (fmaj); + } + + if (!skipDisplay) { + KiiCenter (kapaImage, source->peak->xf, source->peak->yf, 1); + psphotVisualScaleImage (kapaImage, (psImage *) source->pixels->parent, NULL, "image", 1.0, 1); + KiiOverlay overlay; + overlay.x = source->peak->xf; + overlay.y = source->peak->yf; + overlay.dx = 5; + overlay.dy = 5; + overlay.angle = 0.0; + overlay.type = KiiOverlayTypeByName ("circle"); + KiiLoadOverlay (kapaImage, &overlay, 1, "red"); + overlay.x = source->moments->Mx; + overlay.y = source->moments->My; + overlay.dx = 8; + overlay.dy = 8; + overlay.angle = 0.0; + overlay.type = KiiOverlayTypeByName ("circle"); + KiiLoadOverlay (kapaImage, &overlay, 1, "blue"); + } + + psFree (rg); + psFree (Rg); + psFree (fg); + psFree (Fg); + psFree (rb); + psFree (Rb); + psFree (fb); + + if (skipDisplay) return true; + + // pause and wait for user input: + // continue, save (provide name), ?? + char key[10]; + fprintf (stdout, "[q]uit satstar? [e]rase and continue? [o]verplot and continue? [s]kip rest of stars? : "); + if (!fgets(key, 8, stdin)) { + psWarning("Unable to read option"); + } + if (key[0] == 'e') { + KapaClearPlots (kapaGraph); + } + if (key[0] == 'q') { + return false; + } + if (key[0] == 's') { + skipDisplay = true; + } return true; } + +bool psphotTestRadialModel (pmSource *source, psVector **logRmodelOut, psVector **logFmodelOut, psVector *logRad, psVector *logFlux, float Rmax) { + + // we have log(radius) & log(flux). find the median flux in log radial bins + + float logRmax = log10(Rmax); + float logRdel = 0.1; + + psVector *logRmodel = psVectorAlloc (1 + (int) (logRmax / logRdel), PS_TYPE_F32); + psVector *logFmodel = psVectorAlloc (1 + (int) (logRmax / logRdel), PS_TYPE_F32); + + pmSourceRadialProfileSortPair (logRad, logFlux); + + psStats *fluxStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); + psVector *flux = psVectorAllocEmpty (100, PS_TYPE_F32); + + int bin = 0; + for (int i = 0; i < logRmodel->n; i++) { + + // bin (i) has log radius range (i*logRdel : (i+1)*logRdel) + float lRmin = logRdel*(i + 0); + float lRmax = logRdel*(i + 1); + + // reset the flux vector + flux->n = 0; + psStatsInit (fluxStats); + + while (logRad->data.F32[bin] < lRmax) { + if (isfinite(logFlux->data.F32[bin])) { + psVectorAppend (flux, logFlux->data.F32[bin]); + } + bin ++; + } + + // we have the set of fluxes for this bin; find the median values + + float Rmin = pow(10.0, lRmin); + float Rmax = pow(10.0, lRmax); + + float Rmean = (2.0/3.0) * (pow(Rmax, 3.0) - pow(Rmin, 3.0)) / (PS_SQR(Rmax) - PS_SQR(Rmin)); + logRmodel->data.F32[i] = log10(Rmean); + + float Area = M_PI * (PS_SQR(Rmax) - PS_SQR(Rmin)); + if (flux->n < 0.25*Area) { + logFmodel->data.F32[i] = NAN; + continue; + } + + psVectorStats (fluxStats, flux, NULL, NULL, 0); + if (fluxStats->robustMedian > 0.0) { + logFmodel->data.F32[i] = log10(fluxStats->robustMedian); + } else { + logFmodel->data.F32[i] = -3.0; + } + // fprintf (stderr, "R: %f, F: %f +/- %f\n", Rmean, fluxStats->robustMedian, fluxStats->robustStdev); + } + + // now how do i use this to subtract a model?? + *logRmodelOut = logRmodel; + *logFmodelOut = logFmodel; + + // need to free stuff + psFree (fluxStats); + psFree (flux); + + return true; +} + +bool psphotTestRadialModelSub (pmSource *source, psVector *logRmodel, psVector *logFmodel, float Xo, float Yo, float Rmax, psImageMaskType maskVal) { + + int alt; + float logRdel = 0.1; + + for (int iy = 0; iy < source->pixels->numRows; iy++) { + for (int ix = 0; ix < source->pixels->numCols; ix++) { + if ((source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal)) continue; + + float radius = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ; + float logR = log10(radius); + + int bin = (int)(logR / logRdel); + + // we are going to interpolate if possible, or extrapolate if necessary + if (bin >= logRmodel->n - 1) bin = logRmodel->n - 1; + if (bin < 0) bin = 0; + + // interpolate to the current radial position + // XXX BIG HACK : skip nan bins for now + + float logF0 = logFmodel->data.F32[bin]; + float logR0 = logRmodel->data.F32[bin]; + if (!isfinite(logF0)) continue; + + // interpolate between closest two bins if possible, extrapolate on ends + if (logR < logR0) { + alt = (bin > 0) ? bin - 1 : bin + 1; + } else { + alt = (bin < logRmodel->n - 1) ? bin + 1 : bin - 1; + } + + float logF1 = logFmodel->data.F32[alt]; + float logR1 = logRmodel->data.F32[alt]; + if (!isfinite(logF1)) continue; + + // XXX use linear flux, not logFlux + float logF = InterpolateValues (logR0, logF0, logR1, logF1, logR); + float flux = pow (10.0, logF); + + // float flux = InterpolateValues (logR0, logF0, logR1, logF1, logR); + source->pixels->data.F32[iy][ix] -= flux; + } + } + return true; +} + +// only valid for F32 +bool psImageQuickStats (QuickStats *stats, psImage *image, psRegion region) { + + psAssert (image->type.type == PS_TYPE_F32, "unsupported image type"); + + int Npix = (region.y1 - region.y0) * (region.x1 - region.x0); + psVector *tmp = psVectorAllocEmpty (Npix, PS_TYPE_F32); + + int bin = 0; + for (int iy = region.y0; iy < region.y1; iy++) { + for (int ix = region.x0; ix < region.x1; ix++) { + if (!isfinite(image->data.F32[iy][ix])) continue; + tmp->data.F32[bin] = image->data.F32[iy][ix]; + bin ++; + } + } + tmp->n = bin; + + if (bin < 1) { + stats->min = NAN; + stats->lower20 = NAN; + stats->upper20 = NAN; + stats->max = NAN; + stats->Npts = 0; + psFree (tmp); + return true; + } + + + psVectorSortInPlace (tmp); + + int N20 = 0.2*tmp->n; + int N80 = 0.8*tmp->n; + int Nmax = tmp->n - 1; + + stats->min = tmp->data.F32[0]; + stats->lower20 = tmp->data.F32[N20]; + stats->upper20 = tmp->data.F32[N80]; + stats->max = tmp->data.F32[Nmax]; + stats->Npts = bin; + psFree (tmp); + + return true; +} Index: /branches/eam_branches/ipp-20120805/psphot/src/psphotReadout.c =================================================================== --- /branches/eam_branches/ipp-20120805/psphot/src/psphotReadout.c (revision 34367) +++ /branches/eam_branches/ipp-20120805/psphot/src/psphotReadout.c (revision 34368) @@ -128,10 +128,4 @@ } - // find blended neighbors of very saturated stars (detections->newSources) - if (!psphotDeblendSatstars (config, view, filerule)) { - psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis"); - return psphotReadoutCleanup (config, view, filerule); - } - // mark blended peaks PS_SOURCE_BLEND (detections->newSources) // XXX I've deactivated this because it was preventing galaxies close to stars from being @@ -147,4 +141,10 @@ psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications"); return psphotReadoutCleanup (config, view, filerule); + } + + // find blended neighbors of very saturated stars (detections->newSources) + if (!psphotDeblendSatstars (config, view, filerule)) { + psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis"); + return psphotReadoutCleanup (config, view, filerule); } Index: /branches/eam_branches/ipp-20120805/psphot/src/psphotVisual.c =================================================================== --- /branches/eam_branches/ipp-20120805/psphot/src/psphotVisual.c (revision 34367) +++ /branches/eam_branches/ipp-20120805/psphot/src/psphotVisual.c (revision 34368) @@ -1273,5 +1273,5 @@ } -static void plotline (int myKapa, Graphdata *graphdata, float x0, float y0, float x1, float y1) +void plotline (int myKapa, Graphdata *graphdata, float x0, float y0, float x1, float y1) { float x[2], y[2]; @@ -1285,7 +1285,9 @@ } -bool psphotVisualPlotRadialProfile (int myKapa, pmSource *source, psImageMaskType maskVal) { +bool psphotVisualPlotRadialProfile (int myKapa, pmSource *source, psImageMaskType maskVal, pmSourceMode showmode) { Graphdata graphdata; + + float Rmax = (showmode & PM_SOURCE_MODE_SATSTAR) ? 100.0 : 30.0; bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); @@ -1310,11 +1312,11 @@ Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS] - source->pixels->row0; } else { - Xo = source->moments->Mx; - Yo = source->moments->My; + Xo = source->moments->Mx - source->pixels->col0; + Yo = source->moments->My - source->pixels->row0; } for (int iy = 0; iy < source->pixels->numRows; iy++) { for (int ix = 0; ix < source->pixels->numCols; ix++) { - if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) { + if ((source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal)) { rb->data.F32[nb] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ; // rb->data.F32[nb] = hypot (ix - Xo, iy - Yo) ; @@ -1339,5 +1341,5 @@ // examine sources to set data range graphdata.xmin = -0.05; - graphdata.xmax = +30.05; + graphdata.xmax = Rmax + 0.05; graphdata.ymin = -0.05; graphdata.ymax = +8.05; @@ -1370,5 +1372,5 @@ // examine sources to set data range graphdata.xmin = -1.51; - graphdata.xmax = +1.51; + graphdata.xmax = log10(Rmax) + 0.02; graphdata.ymin = -0.05; graphdata.ymax = +8.05; @@ -1422,7 +1424,7 @@ params->data.F32[PM_PAR_YPOS] = 0.0; - psVector *rmod = psVectorAlloc(300, PS_TYPE_F32); - psVector *fmaj = psVectorAlloc(300, PS_TYPE_F32); - psVector *fmin = psVectorAlloc(300, PS_TYPE_F32); + psVector *rmod = psVectorAlloc(Rmax*10, PS_TYPE_F32); + psVector *fmaj = psVectorAlloc(Rmax*10, PS_TYPE_F32); + psVector *fmin = psVectorAlloc(Rmax*10, PS_TYPE_F32); psVector *coord = psVectorAlloc(2, PS_TYPE_F32); @@ -1569,5 +1571,25 @@ if (!(source->mode & showmode)) continue; - psphotVisualPlotRadialProfile (myKapa, source, maskVal); + psphotVisualPlotRadialProfile (myKapa, source, maskVal, showmode); + + if (pmVisualTestLevel("psphot.image", 1)) { + int display = psphotKapaChannel (1); + KiiCenter (display, source->peak->xf, source->peak->yf, 1); + KiiOverlay overlay; + overlay.x = source->peak->xf; + overlay.y = source->peak->yf; + overlay.dx = 5; + overlay.dy = 5; + overlay.angle = 0.0; + overlay.type = KiiOverlayTypeByName ("circle"); + KiiLoadOverlay (display, &overlay, 1, "red"); + overlay.x = source->moments->Mx; + overlay.y = source->moments->My; + overlay.dx = 8; + overlay.dy = 8; + overlay.angle = 0.0; + overlay.type = KiiOverlayTypeByName ("circle"); + KiiLoadOverlay (display, &overlay, 1, "blue"); + } // pause and wait for user input: