Changeset 34368

Timestamp:

Aug 30, 2012, 6:10:52 AM (14 years ago)

Author:

eugene

Message:

some generalization of visualization code; add radial profile fits to deblend satstars (still testing)

Location:

branches/eam_branches/ipp-20120805/psphot/src

Files:

• psphot.h (modified) (2 diffs)
• psphotDeblendSatstars.c (modified) (4 diffs)
• psphotReadout.c (modified) (2 diffs)
• psphotVisual.c (modified) (7 diffs)

branches/eam_branches/ipp-20120805/psphot/src/psphot.h

@@ -281 +281 @@
 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 +295 @@
 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);

branches/eam_branches/ipp-20120805/psphot/src/psphotDeblendSatstars.c

@@ -1 +1 @@
 # 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 +33 @@
 
 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 +230 @@
     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 +372 @@
 }
 
-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 (&section.name, "linlog");
+        KapaSetSection (kapaGraph, &section);
+        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 (&section.name, "loglog");
+        KapaSetSection (kapaGraph, &section);
+        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;
+}

branches/eam_branches/ipp-20120805/psphot/src/psphotReadout.c

@@ -128 +128 @@
     }
 
-    // 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 +141 @@
         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);
     }
 

branches/eam_branches/ipp-20120805/psphot/src/psphotVisual.c

@@ -1273 +1273 @@
 }
 
-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 +1285 @@
 }
 
-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 +1312 @@
         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 +1341 @@
     // 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 +1372 @@
     // 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 +1424 @@
         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 +1571 @@
         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: