Changeset 34404 for trunk/psphot/src/psphotDeblendSatstars.c

Timestamp:

Sep 5, 2012, 4:19:30 PM (14 years ago)

Author:

eugene

Message:

adding SATSTAR profile code (but keeping it disabled for now)

Location:

trunk/psphot

Files:

• . (modified) (1 prop)
• src (modified) (1 prop)
• src/psphotDeblendSatstars.c (modified) (3 diffs)

trunk/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);
+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
@@ -16 +28 @@
 }
 
+/** this function does 3 things:
+
+    1) choose likely saturated stars
+       - mode | SATSTAR 
+       - if moments->Peak > 0.5*SATURATION, examine region of peak pixels (ensure we go down to 0.1*SAT) for SAT mask
+      
+    2) adjust the window for saturated objects
+
+    3) measure the radial profile
+
+    4) subtract the radial profile
+    
+    TBD:
+
+    * function to replace the radial profile for a source
+    * recenter the profile (based on cross-correlation / convolution with 1D profile)
+
+    * raise a bit somewhere for these super saturated stars (if they were so modeled)
+    * consider the subtraction bit : when to raise it?
+
+ **/
+
 bool psphotDeblendSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int fileIndex) {
+
+    int N;
+    pmSource *source;
+    bool status;
+
+    // XXX disable this function for now
+    return true;
+
+    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;
+
+        if (!psphotSatstarProfileModel (source, maskVal)) continue;
+
+        source->mode |= PM_SOURCE_MODE_SATSTAR; // yes, this source IS saturated
+        source->mode2 |= PM_SOURCE_MODE2_SATSTAR_PROFILE; // and we have in fact subtracted the profile
+
+        // 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;
@@ -219 +397 @@
     return true;
 }
+
+// create a model for the radial profile of a saturated star (is this actually more generic?)
+bool psphotSatstarProfileModel (pmSource *source, psImageMaskType maskVal) {
+
+    // XXX user define somewhere?
+    float Rmax = 320.0;
+
+    // XXX is this ever the case??
+    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 *logR = psVectorAllocEmpty (nPts, PS_TYPE_F32);
+    psVector *flux = psVectorAllocEmpty (nPts, PS_TYPE_F32);
+
+    int ng = 0;
+
+    // choose the best center for this profile
+    float Xo = NAN;
+    float Yo = NAN;
+
+    if (source->modelPSF) {
+        // XXX do we ever have a PSF model for a SATSTAR?
+        Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS];
+        Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS];
+    } else {
+        Xo = source->moments->Mx;
+        Yo = source->moments->My;
+    }
+
+    float Xc = Xo - source->pixels->col0 - 0.5;
+    float Yc = Yo - source->pixels->row0 - 0.5;
+
+    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;
+            // XXX do this faster by generating R^2 and returning 0.5*log10(R^2)?
+            logR->data.F32[ng] = log10(hypot (ix - Xc, iy - Yc));
+            flux->data.F32[ng] = source->pixels->data.F32[iy][ix];
+            ng++;
+        }
+    }
+    logR->n = ng;
+    flux->n = ng;
+
+    // XXX do something sensible here if there are no pixels
+
+    // measure the radial profile 
+    psVector *logFmodel = NULL;
+    psVector *logRmodel = NULL;
+    psphotSatstarProfileCreate (source, &logRmodel, &logFmodel, logR, flux, Rmax);
+    
+    // XXX do something sensible here if the profile is crap
+
+    source->satstar = pmSourceSatstarAlloc();
+    source->satstar->Xo = Xo;
+    source->satstar->Yo = Yo;
+    source->satstar->Rmax = Rmax;
+    source->satstar->logFmodel = logFmodel;
+    source->satstar->logRmodel = logRmodel;
+
+    // subtract the profile (false => subtract)
+    psphotSatstarProfileOp (source, maskVal, 1.0, 0, false);
+
+    psFree (logR);
+    psFree (flux);
+
+    return true;
+}
+
+// Take logR + flux and generate radial bins logRmodelOut, logFmodelOut
+bool psphotSatstarProfileCreate (pmSource *source, psVector **logRmodelOut, psVector **logFmodelOut, psVector *logR, psVector *flux, float Rmax) {
+
+  // we have log(radius) & log(flux).  find the median flux in log radial bins
+
+  float logRmax = log10(Rmax);
+  float logRdel = 0.1; // XXX user-define?
+
+  psVector *logRmodel = psVectorAlloc (1 + (int) (logRmax / logRdel), PS_TYPE_F32);
+  psVector *logFmodel = psVectorAlloc (1 + (int) (logRmax / logRdel), PS_TYPE_F32);
+
+  pmSourceRadialProfileSortPair (logR, flux);
+
+  psStats *fluxStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
+  psVector *fluxVals = 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
+    fluxVals->n = 0;
+    psStatsInit (fluxStats);
+
+    while (logR->data.F32[bin] < lRmax) {
+      if (isfinite(flux->data.F32[bin])) {
+          psVectorAppend (fluxVals, flux->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 (fluxVals->n < 0.25*Area) {
+      logFmodel->data.F32[i] = NAN;
+      continue;
+    }
+    
+    psVectorStats (fluxStats, fluxVals, 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 (fluxVals);
+
+  return true;
+}
+
+bool psphotSatstarProfileOp (pmSource *source, psImageMaskType maskVal, float FACTOR, pmModelOpMode mode, bool add) {
+
+    int alt;
+    float logRdel = 0.1;
+
+    float Xc = source->satstar->Xo - source->pixels->col0 - 0.5;
+    float Yc = source->satstar->Yo - source->pixels->row0 - 0.5;
+    psVector *logRmodel = source->satstar->logRmodel;
+    psVector *logFmodel = source->satstar->logFmodel;
+
+    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 - Xc, iy - Yc) ;
+            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);
+
+            if (mode & PM_MODEL_OP_NOISE) {
+                if (add) {
+                    source->variance->data.F32[iy][ix] += FACTOR * flux;
+                } else {
+                    source->variance->data.F32[iy][ix] -= FACTOR * flux;
+                }
+            } else {
+                if (add) {
+                    source->pixels->data.F32[iy][ix] += flux;
+                } else {
+                    source->pixels->data.F32[iy][ix] -= flux;
+                }
+            }
+        }
+    }
+    return true;
+}
+
+// Take logR + flux and generate radial bins logRmodelOut, logFmodelOut
+bool psphotSatstarPhotometry (pmSource *source) {
+
+    if (!source->satstar) return false;
+
+    psVector *logRmodel = source->satstar->logRmodel;
+    psVector *logFmodel = source->satstar->logFmodel;
+    
+    float fluxTotal = 0.0;
+    float logRdel = 0.1; // XXX user-define (or carry in satstar)
+
+    // integrate flux in radial bins 
+    for (int i = 0; i < logRmodel->n; i++) {
+        // just add up the mean flux in each bin and multiply by the area of the bin
+        
+        float logF = logFmodel->data.F32[i];
+        if (!isfinite(logF)) continue;
+        float flux = pow(10.0, logF); // this is the mean flux per pixel (ie, surface brightness)
+
+        // bin (i) has log radius range (i*logRdel : (i+1)*logRdel)
+        float lRmin = logRdel*(i + 0);
+        float lRmax = logRdel*(i + 1);
+
+        float Rmin = pow(10.0, lRmin);
+        float Rmax = pow(10.0, lRmax);
+
+        float Area = M_PI * (PS_SQR(Rmax) - PS_SQR(Rmin));
+
+        fluxTotal += flux * Area;
+    }
+
+    source->psfMag    = -2.5*log10(fluxTotal);
+    source->psfMagErr = 0.05; 
+
+    // XXX I have no idea of a realistic error on the photometry here.  the bottom line is that
+    // the error is totally dominated by the loss of charge due to saturation.  I am not
+    // modeling the profile in any real detail other than to follow the radial bins.  
+
+    source->extMag    = NAN;
+    source->apMag     = NAN;
+    source->apMagRaw  = NAN;
+    source->apFlux    = fluxTotal;
+    source->apFluxErr = 0.05*fluxTotal;
+
+    return true;
+}
+
+# if (0)
+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 {
+        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;
+        
+        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;
+        
+        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;
+}
+# endif
+
+// 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;
+}
+
+bool psphotAddOrSubSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int fileIndex, psMetadata *recipe, bool add) {
+
+    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->allSources;
+    psAssert (sources, "missing sources?");
+
+    if (!sources->n) {
+        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping satstar blend");
+        return true;
+    }
+
+    // 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);
+
+    // examine sources in decreasing SN order
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *source = sources->data[i];
+
+        if (!(source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE)) continue;
+
+        if (!psphotSatstarProfileOp (source, maskVal, 1.0, 0, add)) continue;
+    }
+
+    psLogMsg ("psphot", PS_LOG_DETAIL, "satstar op: %f sec\n", psTimerMark ("psphot.deblend.sat"));
+    return true;
+}
+